Numerically Representing A Stochastic Process Algebra

The syntactic nature and compositionality characteristic of stochastic process algebras make models to be easily understood by human beings, but not convenient for machines as well as people to directly carry out mathematical analysis and stochastic simulation. This paper presents a numerical representation schema for the stochastic process algebra PEPA, which can provide a platform to directly and conveniently employ a variety of computational approaches to both qualitatively and quantitatively analyse the models. Moreover, these approaches developed on the basis of the schema are demonstrated and discussed. In particular, algorithms for automatically deriving the schema from a general PEPA model and simulating the model based on the derived schema to derive performance measures are presented.


Introduction
Stochastic process algebras, such as PEPA [1], TIPP [2], EMPA [3], are powerful modelling formalisms for concurrent systems which have enjoyed considerable success over the last decade. A stochastic process algebra model is constructed to approximately and abstractly represent a system whilst hiding its implementation details. Based on the model, performance properties of the dynamic behaviour of the system can be assessed, through some techniques and computational methods. This process is referred to as the performance modelling of the system, which mainly involves three levels: model construction, technical computation and performance derivation. In order to derive performance measures from large scale stochastic process algebra models, many mathematical tools and approaches have been proposed to study the models. For instance, a fluid approximation method has been proposed in [4] to avoid the state space explosion problem encountered in the analysis of large scale PEPA models.
However, the syntactic nature of stochastic process algebras makes models easily understood by human beings, but not convenient for machines/computers (as well as for human beings) to directly employ these tools and approaches. In addition, the compositionality of the formalisms allow a model to be locally defined, but the analysis of the model or the underlying continuoustime Markov chain (CTMC) carries out in the global manner since it is the whole system rather than a part of it to be usually interested and considered. The syntactical and compositional qualities of the stochastic process algebras, which are advantages in model construction, turn to be disadvantages in model analysis.

Paper contributions
In order to overcome the obstacles in the direct and convenient application of the mathematical tools, we propose a new numerical representation schema for the formulism PEPA in this paper. In this schema, labelled activities are defined to cope with the difference between actions in PEPA and transitions in the underlying CTMC, so that the correspondence between them is one-to-one. Activity matrices based on the labelled activities are defined to capture structural information about PEPA models. Moreover, transition rate functions are proposed to capture the timing information. These concepts numerically describe and represent a PEPA model, and provide a platform for conveniently and easily simulating the underlying CTMC, deriving the fluid approximation, as well as leading to an underlying Place/Transition (P/T) structure. These definitions are consistent with the original semantics of PEPA, and a PEPA model can be recovered from its numerical representation. An algorithm for automatically deriving the schema from any given PEPA model has been provided. Some characteristics of this numerical representation are revealed. For example, using numerical vector forms the exponential increase of the size of the state space with the number of components can be reduced to at most a polynomial increase. The benefits of the schema embodies the three aspects of performance modelling, which are illustrated by Figure 1. At the first level, the proposed new schema numerically describes any given PEPA model and provides a platform to directly employ a variety of approaches to analyse the model. These approaches are shown at the second level. At this level, a fluid approximation method for the quantitative analysis of PEPA is established, as well as investigated, mainly with respect to its convergence and the consistency between this method and the underlying CTMC. In addition, a P/T structure-based approach is revealed, which can be utilised to qualitatively analyse the model. At the third level, both qualitative and quantitative performance measures can be derived from the model through those approaches. A stochastic simulation algorithm for the aggregated CTMC, which is based on the numerical representation schema, is proposed to obtain general performance metrics in this paper. As for the other two approaches, related investigation and analysis were given in [5] and [6] respectively, which will be briefly introduced in the next subsection.

Related work
Our work is motivated and stimulated by the pioneering work on the numerical vector form and activity matrix in [4], which was dedicated to the fluid approximation for PEPA. The P/T structure underlying each PEPA model, as stated in Theorem 1 in this paper, reveals tight connections between stochastic process algebras and stochastic Petri nets. Based on this structure and the theories developed for Petri nets, several powerful techniques for structural analysis of PEPA were presented in [5], including a structure-based deadlock-checking method which avoids the state space explosion problem. In [7], a new operational semantics was proposed to give a compact symbolic representation of PEPA models. This semantics extends the application scope of the fluid approximation of PEPA by incorporating all the operators of the language and removing earlier assumptions on the syntactical structure of the models amenable to this analysis. Moreover, the paper [6] shows how to derive the performance metrics such as action throughput and capacity utilisaition from the fluid approximation of a PEPA model.

Paper organisation
The remainder of this paper is structured as follows: Section 2 gives a brief introduction to the PEPA formulism; In Section 3, 4 and 5, we respectively present the three combinators of the numerical schema, i.e. the numerical vector form, labelled activity and activity matrix, as well as the transition rate function. Computational approaches for performance derivation that are developed on the basis of the schema are demonstrated in Section 6. We finally conclude the paper in Section 7.
2 Introduction to PEPA PEPA (Performance Evaluation Process Algebra) [1], developed by Hillston in the 1990s, is a high-level model specification language for low-level stochastic models, and describes a system as an interaction of the components which engage in activities. In contrast to classical process algebras, activities are assumed to have a duration which is a random variable governed by an exponential distribution. Thus each activity in PEPA is a pair (α, r) where α is the action type and r is the activity rate. The language has a small number of combinators, for which we provide a brief introduction below; the structured operational semantics can be found in [1]. The grammar is as follows: where S denotes a sequential component and P denotes a model component which executes in parallel. C stands for a constant which denotes either a sequential component or a model component as introduced by a definition. C S stands for constants which denote sequential components. The effect of this syntactic separation between these types of constants is to constrain legal PEPA components to be cooperations of sequential processes.
Prefix: The prefix component (α, r).S has a designated first activity (α, r), which has action type α and a duration which satisfies exponential distribution with parameter r, and subsequently behaves as S.
Choice: The component S + T represents a system which may behave either as S or as T . The activities of both S and T are enabled. Since each has an associated rate there is a race condition between them and the first to complete is selected. This gives rise to an implicit probabilistic choice between actions dependent of the relative values of their rates.
Hiding: Hiding provides type abstraction, but note that the duration of the activity is unaffected. In P/L all activities whose action types are in L appear as the "private" type τ .
Cooperation:P L Q denotes cooperation between P and Q over action types in the cooperation set L. The cooperands are forced to synchronise on action types in L while they can proceed independently and concurrently with other enabled activities (individual activities). The rate of the synchronised or shared activity is determined by the slower cooperation (see [1] for details). We write P Q as an abbreviation for P L Q when L = ∅ and P [N ] is used to represent N copies of P in a parallel, i.e. P [3] = P P P .
Constant: The meaning of a constant is given by a defining equation such as A def = P . This allows infinite behaviour over finite states to be defined via mutually recursive definitions.
On the basis of the operational semantic rules (please refer to [1] for details), a PEPA model may be regarded as a labelled multi-transition system C, Act, where C is the set of components, Act is the set of activities and the multi-relation (α,r) −→ is given by the rules. If a component P behaves as Q after it completes activity (α, r), then denote the transition as P (α,r) −→ Q. The memoryless property of the exponential distribution, which is satisfied by the durations of all activities, means that the stochastic process underlying the labelled transition system has the Markov property. Hence the underlying stochastic process is a CTMC. Note that in this representation the states of the system are the syntactic terms derived by the operational semantics. Once constructed the CTMC can be used to find steady-state or transient probability distributions from which quantitative performance can be derived.

Numerical Vector Form
The usual state representation in PEPA models is in terms of the syntactic forms of the model expression. When a large number of repeated components are involved in a system, the state space of the CTMC underling the model can be large. This is mainly because each copy of the same type of component is considered to be distinct, resulting in distinct Markovian states. The multiple states within the model that exhibit the same behaviour can be aggregated to reduce the size of the state space as shown by Gilmore et al. [8] using the technique based on a vector form. The CTMC is therefore constructed in terms of equivalence classes of syntactic terms. "At the heart of this technique is the use of a canonical state vector to capture the syntactic form of a model expression", as indicated in [4]. Rather than the canonical representation style, an alternative numerical vector form was proposed by Hillston in [4] for capturing the state information of models with repeated components. In the numerical vector form, there is one entry for each local derivative of each type of component in the model. The entries in the vector are the number of components currently exhibiting this local derivative, no longer syntactic terms representing the local derivative of the sequential component. Following [4], hereafter the term local derivative refers to the local state of a single sequential component. By adopting this model-aggregation technique, the number of the states of the system can be reduced to only increase (at most) polynomially with the number of instances of the components. According to Definition 1, m(C i j ) ≥ 0 for each C i j . At any time, each sequential component stays in one and only one local derivative. So the sum of m(C i ), i.e. d i j=1 m[C i j ], specifies the population of C i in the system. Notice that m(C i ) satisfies is an upper bound of the exact number of the states in terms of C i .
Notice that Therefore, it is easy to verify the following Proposition 1. Consider a system comprising n types of component, namely C 1 , C 2 , · · · , C n , with M i copies of the component of type C i in the system, where C i has d i local derivatives, for i = 1, 2, · · · , n. Then the size of the state space of the system is at most The upper bound given in Proposition 1 guarantees that the size of the state space increases at most polynomially with the number of instances of the components. Consider the following PEPA model. According to the semantics of PEPA originally defined in [1], the size of the state space of the CTMC underlying Model 1 is 2 M +N , which increases exponentially with the numbers of the users and severs in the system. According to Definition 1, the system vector m has four entries representing the instances of components in the total four local derivatives, that is Let M = N = 2, then the system equation of Model 1 determines the starting state (2, 0, 2, 0). By enabling activities or transitions, all reachable system states can be manifested, see Figure 2. The size of the state space is nine. The upper bound of the size given by Proposition 1, , is nine, coinciding with the size of the state space. The bound given in Proposition 1 is sharp and can be hit in some situations.

Labelled Activity and Activity Matrix
In the PEPA language, the transition is embodied in the syntactical definition of activities, in the context of sequential components. Since the consideration is in terms of the whole system rather than sequential components, the transition between these system states should be defined and represented. This section presents a numerical representation for the transitions between system states and demonstrates how to derive this representation from a general PEPA model. If a system vector changes into another vector after firing an activity, then the difference between these two vectors manifests the transition corresponding to this activity. Obviously, the difference is in numerical forms since all states are numerical vectors. Consider Model 1 and its transition diagram in Figure 2. Each activity in the model corresponds to a vector, called the transition vector. For example, task 1 corresponds to the transition vector l task 1 = (−1, 1, −1, 1) T . That is, the derived state vector by firing task 1 from a state, can be represented by the sum of l task 1 and the state enabling task 1 . For instance, (2, 0, 2, 0) T + l task 1 = (1, 1, 1, 1) T illustrates that s 1 = (2, 0, 2, 0) T transitions into s 2 = (1, 1, 1, 1) T after enabling task 1 .
Each activity in the model is represented by a transition vector -a column of the activity matrix, and each column expresses an activity. So the activity matrix is essentially indicating both an injection and a surjection from syntactic to numerical representation of the transition between system states. The concept of the activity matrix for PEPA was first proposed by Hillston in [4,10]. However, the original definition cannot fully reflect the representation mapping Table 1: Transition vectors form an activity matrix here. This is due to the fact of that the original definition is local-derivative-centric rather than transition centric. This results in some limitations for more general applications. For example, for some PEPA models (e.g. Model 2 in the following context), some columns of the originally defined matrix cannot be taken as transition vectors so that this definition cannot fully reflect the PEPA semantics in some circumstances. In the following, a modified definition of the activity matrix is given. The new definition is activity-or transition-centric, which brings the benefit that each transition is represented by a column of the matrix and vice versa. Figure 3: Transition diagram of Model 2 Let us consider a PEPA model, i.e. the following Model 2, in which there are multiple choices after firing some activities. In this model, firing α in the component P may lead to two possible local derivatives: P 2 and P 3 , while firing β may lead to P 1 and P 3 . In addition, firing γ may lead to P 1 , Q 1 . However, only one derivative can be chosen after each firing of an activity, according to the semantics of PEPA. But the original definition of activity matrix cannot clearly reflect this point. See the activity matrix of Model 2 given in Table 2. Moreover, the individual activity γ in this table, which can be enabled by both P 3 and Q 2 , may be confused as a shared activity.
In order to better reflect the semantics of PEPA, we modify the definition of the activity matrix in this way: if there are m possible outputs, namely {R 1 , R 2 , · · · , R m }, after firing either an individual or a shared activity l, then l is "split" into m labelled ls: l w 1 , l w 2 , · · · , l wm . Here respectively. Each l w i can only lead to a unique output R i . Here there are no new activities created, since we just attach labels to the activity to distinguish the outputs of firing this activity. The modified activity matrix clearly reflects that only one, not two or more, result can be obtained from firing l. And thus, each l w i can represent a transition vector.
For example, see the modified activity matrix of Model 2 in Table 3. In this activity matrix, the individual activity γ has different "names" for different component types, so that it is not confused with a shared activity. Another activity β, is labelled as β P 2 →P 1 and β P 2 →P 3 , to respectively reflect the corresponding two choices. In this table, the activity α is also split and attached with labels.
Before giving the modified definition of activity matrix for any general PEPA model, the pre and post sets for an activity are first defined. For convenience, throughout this paper any transition U (l,r) −→ V defined in the PEPA models may be rewritten as U 1. If a local derivative U can enable an activity l, that is U l −→ ·, then U is called a pre local derivative of l. The set of all pre local derivatives of l is denoted by pre(l), called the pre set of l.
2. If V is a local derivative obtained by firing an activity l, i.e. ·  Obviously, if l has only one pre local derivative, i.e. #pre(l) = 1, then l is an individual activity, like β in Model 2, whereafter the notation #A is defined as the cardinality of the set A, i.e. the number of elements of A. But l being individual does not imply #pre(l) = 1, see γ for instance. If l is shared, then #pre(l) > 1, for example, see #pre(α) = #{P 1 , Q 1 } = 2. For a shared activity l with pre(l) = k, there are k local derivatives that can enable this activity, each of them belonging to a distinct component type.
For an individual activity l, things are rather simple and easy: for U ∈ pre(l), l can be labelled as l U →post(U,l) [1] , l U →post(U,l) [2] , l U →post(U,l)[k U ] , where k U = #post(U, l). Varying U ∈ pre(l), there are U ∈pre(l) #post(U, l) labels needed to distinguish the possible transitions. See β P 2 →P 1 , β P 2 →P 3 , γ P 3 →P 1 , γ Q 2 →Q 1 in Model 2 for instance. Now we give the formal definition. 1. For any individual activity l, for each U ∈ pre(l), V ∈ post(U, l), label l as l U →V .
2. For a shared activity l, for each (V 1 , V 2 , · · · , V k ) ∈ post(pre(l) [1], l) × post(pre(l) [2], l) × · · · × post(pre(l)[k], l), label l as l w , where Each l U →V or l w is called a labelled activity. The set of all labelled activities is denoted by A label . For the above labelled activities l U →V and l w , their respective pre and post sets are defined as According to Definition 3, each l U →V or l w can only lead to a unique output. No new activities are created, since labels are only attached to the activity to distinguish the results after this activity is fired.
The impact of labelled activities on local derivatives can be recorded in a matrix form, as defined below. Definition 4. (Activity Matrix). For a model with N A label labelled activities and N D distinct local derivatives, the activity matrix C is an N D × N A label matrix, and the entries are defined as follows where l j is a labelled activity. The pre activity matrix C pre and post activity matrix C post are defined as otherwise. , The modified activity matrix captures all the structural information, including the information about choices and synchronisations, of a given PEPA model. From each row of the matrix, which corresponds to each local derivative, we can know which activities this local derivative can enable and after which activities are fired this local derivative can be derived. From the perspective of the columns, the number of "−1"s in a column tells whether the corresponding activity is synchronised or not. Only one "−1" means that this transition corresponds to an individual activity. The locations of "−1" and "1" indicate which local derivatives can enable the activity and what the derived local derivatives are, i.e. the pre and post local derivatives. In addition, the numbers of "−1"s and "1"s in each column are the same, because any transition in any component type corresponds to a unique pair of pre and post local derivatives. In fact, all this information is also stored in the labels of the activities. Therefore, with the transition rate functions defined in the next section to capture the timing information, a given PEPA model can be recovered from its activity matrix.
Moreover, the pre and post activity matrix indicate the local derivatives which can fire a labelled activity and the derived local derivative after firing a labelled activity respectively. The modified activity matrix equals the difference between the pre and post activity matrices, i.e. C = C Post − C Pre . Hereafter the terminology of activity matrix refers to the one in Definition 4. This definition embodies the transition or operation rule of a given PEPA model, with the exception of timing information. For a given PEPA model, each transition of the system results from the firing of an activity. Each optional result after enabling this activity corresponds to a relevant labelled activity, that is, corresponds to a column of the activity matrix. Conversely, each column of the activity matrix corresponding to a labelled activity, represents an activity and the chosen derived result after this activity is fired. So each column corresponds to a system transition. Therefore, we have the following proposition, which specifies the correspondence between system transitions and the columns of the activity matrix.
Proposition 2. Each column of the activity matrix corresponds to a system transition and each transition can be represented by a column of the activity matrix.

Transition Rate Function
The structural information of any general PEPA model is captured in the activity matrix, which is constituted by all transition vectors. However, the duration of each transition has not yet been specified. This section defines transition rate functions for transition vectors or labelled activities to capture the timing information of PEPA models.

Model 2 continued
Let us start from Model 2 again. As Table 3 shows, activity γ in Model 2 is labelled as γ P 3 →P 1 and γ Q 2 →Q 1 . For γ P 3 →P 1 , there are x[P 3 ] instances of the component type P in the local derivative P 3 in state x, each enabling the individual activity concurrently with the rate r γ . So the rate of γ P 3 →P 1 in state x is f (x, γ P 3 →P 1 ) = r γ x[P 3 ]. Similarly, the rate for γ Q 2 →Q 1 in state x is r γ x[Q 2 ]. This is consistent with the definition of apparent rate in PEPA, which states that if there are N replicated instances of a component enabling a transition (l, r), the apparent rate of the activity will be r × N .
In Model 2 activity β is labelled as β P 2 →P 1 and β P 2 →P 3 , to respectively reflect the corresponding two choices. According to the model definition, there is a flux of r β x(P 2 ) into P 1 from P 2 after firing β in state x. So the transition rate function is defined as f (x, β P 2 →P 1 ) = r β x[P 2 ]. Similarly, we can define f (x, β P 2 →P 3 ) = r β x[P 2 ]. These rate functions can be defined or interpreted in an alternative way. In state x, there are x[P 2 ] instances that can fire β. So the apparent rate of β is (r β + r β )x[P 2 ]. By the semantics of PEPA, the probabilities of choosing the outputs are r β r β +r β and r β r β +r β respectively. So the rate of the transition β P 2 →P 1 is while the rate of the transition β P 2 →P 3 is In Model 2, α is a shared activity with three local rates: r α , r α and r α . The apparent rate of α in P 1 is (r α +r α )x[P 1 ], while in Q 1 it is r α x[Q 1 ]. According to the PEPA semantics, the apparent rate of a synchronised activity is the minimum of the apparent rates of the cooperating components. So the apparent rate of α as a synchronisation activity is min{(r α + r α )x[P 1 ], r α x[Q 1 ]}. After firing α, P 1 becomes either P 2 or P 3 , with the probabilities r α r α +r α and r α r α +r α respectively. Simultaneously, Q 1 becomes Q 2 with the probability 1. So the rate function of transition (P 1 → P 2 , Q 1 → Q 2 ), represented by f (x, α (P 1 →P 2 ,Q 1 →Q 2 ) ), is The above discussion about the simple example should help the reader to understand the definition of transition rate function for general PEPA models, which is presented in the next subsection.

Definitions of transition rate function
In a PEPA model, as we have mentioned, we may rewrite any U  (5.6) The apparent rate of l in U in state x, denoted by r l (x, U ), is defined as The above definition is used to define the following transition rate function. −→ V , the transition rate function of labelled activity l U →V in state x is defined as Then the transition rate function of labelled activity l w in state x is defined as Remark 1. Definition 6 accommodates the passive or unspecified rate . If some r U →V l are , then the relevant calculation in the rate functions (5.8) and (5.9) can be made according to the inequalities and equations that define the comparison and manipulation of unspecified activity rates (see [1]). Moreover, we assume that 0 · = 0. So the terms such as "min{A , rB}" are interpreted as [11]: The definition of the transition rate function is consistent with the semantics of PEPA: The transition rate function in Definition 6 is consistent with the operational semantics of PEPA.
The proof is easy and omitted here. Since both the structural and timing information has been captured in the defined numerical representation schema, PEPA models can be therefore recovered from its representation schema. In addition, it is also easy to find that the transition rate function has the following homogenous property. This property will identify the CTMCs underlying a PEPA model to be density dependent (see Theorem 4 in the next section).

Algorithm for deriving activity matrix and transition rate functions
This subsection presents an algorithm for automatically deriving the activity matrix and transition rate functions from any PEPA model, see Algorithm 1. The lines 3-12 of Algorithm 1 deal with individual activities while lines 13 − 32 deal with shared activities. The calculation methods in this algorithm are essentially the embodiment of the definitions of labelled activity and apparent rate as well as transition rate function. So we do not give more introduction to this algorithm.

Computational approaches for PEPA
As a model being represented numerically, efficient techniques such as stochastic simulation and fluid approximation can be can be directly utilised to analyse the model. This section briefly introduces these approaches as well as technical foundations for employing them in the context of PEPA.

Place/Transition structure in PEPA models
Whilst the focus of stochastic process algebras has understandably been primarily quantitative analysis, qualitative analysis can also provide valuable insight into the behaviour of a system. In contrast, in Petri net modelling there are well-established techniques of structural analysis [12][13][14]. This subsection shows how the new representation schema helps to manifest the P/T structure underlying PEPA models, and makes it possible to readily adapt structural analysis techniques for Petri nets to PEPA. First, the relevant definitions are given below.

2.
A marking is a vector m : P → N that assigns to each place of a P/T net a nonnegative integer (number of tokens).
3. A P/T system is a pair S = N , m 0 : a net N with an initial marking m 0 .
By Definition 7, it is easy and direct to verify Theorem 1. There is a P/T system underlying any PEPA model, that is N , m 0 , where m 0 is the starting state; N = D, A label , C Pre , C Post is P/T net: where D is the set of local derivatives, A label is the labelled activity set; C Pre and C Post are the pre and post activity matrices respectively. Moreover, each state m of the PEPA model is a marking.
Based on the P/T structure underlying PEPA models and the theories developed for P/T nets, several powerful techniques and approaches for structural analysis of PEPA were established in [5]. For instance, the authors gave a method of deriving and storing the state space which avoids the problems associated with populations of components, and an approach to find invariants which can be used to qualitatively reason about systems. Moreover, a structure-based deadlock-checking algorithm was proposed, which can avoid the state space explosion problem. // Form a corresponding column of the activity matrix and the rate function 8: for all (V 1 , V 2 , · · · , V k ) ∈ K(l) do 28:

29:
A label = A label ∪ {l w } // Label l as l w

30:
// Form a column of M a and the rate function corresponding to l w

Stochastic simulation of PEPA models
By solving the global balance equations associated with the infinitesimal generator of the CTMC underlying a PEPA model, the steady-state probability distribution can be obtained, from which performance measures can be derived. According to the original definition of the PEPA language in which each instance of the same component type is considered distinctly, the size of the state space of this original CTMC may increase exponentially with the number of components. By adopting the numerical vector form to represent the system state which results in the aggregated CTMC, the size of the state space can thus be significantly reduced, as Proposition 1 shows, together with the computational complexity of deriving the performance by solving the corresponding global balance equations since, the dimension of the infinitesimal generator matrix is the square of the size of the state space.
Unless otherwise stated, hereafter the CTMC underlying a PEPA model refers to the aggregated CTMC, and the state of a model or a system is considered in the sense of aggregation. If the size of the state space is too large, it is not feasible to calculate the steady-state distribution and thus to get a performance measure R, which is usually expressed as R = s∈S ρ(s)π(s) where ρ and π defined on the state space S are the reward function and the steady-state probability distribution respectively. An alternative widely-used way to obtained performance is stochastic simulation.
As discussed previously, a transition between states, namely from x to x + l, is represented by a transition vector l corresponding to the labelled activity l (for convenience, hereafter each pair of transition vectors and corresponding labelled activities shares the same notation). The rate of the transition l in state x is specified by the transition rate function f (x, l). That is, x (l,f (x,l)) −→ x + l. Given a starting state x 0 , the transition chain corresponding to a firing sequence l 0 , l 1 , · · · , l, · · · is The above sequence can be considered as one path or realisation of a simulation of the aggregated CTMC, if the enabled activity at each state is chosen stochastically, i.e. is chosen through the approach of sampling. After a long time, the steady-state of the system is assumed to be achieved. Hence the average performance R = s∈S ρ(s)π(s) can be calculated. As one benefit of our numerical representation schema, it provides a good platform for directly and conveniently simulating the CTMC for PEPA, see Algorithm 2. In Algorithm 2, the states of a PEPA model are represented as numerical vector forms, and the rates between those states are specified by the transition rate functions which only depend on the transition type (i.e. labelled activity) and the current state. In this algorithm, the generated time τ in each iteration can be regarded as having been drawn from an exponential distribution with the mean 1 f (x) (see Example 2.3 in [15], page 38). That is, Line 9 in Algorithm 2 is in fact expressing: "generate τ from an exponential distribution with the mean 1 f (x) ". Line 10 determines which transition will be chosen, and consequently determines the next state that the system will transition into. Therefore, this algorithm is essentially to simulate the CTMC underlying a PEPA model. The choices for stopping the algorithm include a given large time, or the absolute or relative error of two continued iterations being small enough (since the output performance converges compute the transition rate function f (x, l j ), j = 1, 2, · · · , m 7: generate uniform random numbers r 1 , r 2 on [0, 1] 9: compute τ = −(1/f (x)) ln r 1 10: x = x + l µ // Update state vector of system 15: end while 16: Output performance: PerMeasure t as time goes to infinity as the following Theorem 3 states). Now we prove the convergence of the performance calculated using the algorithm. We need the following theorem. Moreover, for any bounded function ρ : S → R, we have where E[ρ(X)] = s∈S ρ(s)π(s).
Here is our conclusion (we assume the CTMCs underlying PEPA models to be irreducible and positive recurrent): Theorem 3. The performance measure calculated according to Algorithm 2 converges as time goes to infinity, that is, lim t→∞ Proof. Assume that n − 1 iterations have been finished and the time has accumulated to t n−1 . Suppose the current one is the n-th iteration and τ is the generated time in this iteration. After the n-th iteration is finished, the accumulated time will be updated to t n = t n−1 + τ . During the τ time interval, the simulated CTMC stays in the state x, that is, Therefore, after this n-th iteration, PerMeasure will be accumulated to = tn 0 ρ(X z )dz and PerMeasure t n = 1 t n tn 0 ρ(X z )dz.
According to Theorem 2, 1 t n tn 0 ρ(X z )dz tends to E[ρ(X)] as t n tends to infinity. So the performance obtained through Algorithm 2 converges to E[ρ(X)] as the simulation time goes to infinity.
Performance metrics, such as activity throughput of an activity and capacity utilisation of a local derivative that are discussed in [6], can be derived through this algorithm by choosing appropriate reward functions.

Fluid approximation of PEPA models
The weakness of the simulation method is its high computational cost, which makes it not suitable for real-time performance monitoring or prediction. Recently, a novel approach to get performance measures from PEPA models has been proposed in [4] and subsequently expanded in [7] and [17], making a continuous state space approximation as a set of ordinary differential equations (ODEs). In this section, we present a mapping semantics for this approaches, which is based on the numerical representation schema. In addition, a theoretical justification of this approach, mainly in terms of its consistency with the CTMCs, will be discussed.
In our representation schema, the transition rate f (x, l) reflects the intensity of the transition from state x to state x + l. The state space is inherently discrete with the entries within the numerical vector form always being non-negative integers and always being incremented or decremented in steps of one. As pointed out in [4], when the numbers of components are large these steps are relatively small and we can approximate the behaviour by considering the movement between states to be continuous, rather than occurring in discontinuous jumps. In fact, let us consider the evolution of the numerical state vector. Denote the state at time t by x(t). In a short time ∆t, the change to the vector x(t) will be x(·, t + ∆t) − x(·, t) = ∆t l∈A label lf (x(·, t), l).
Dividing by ∆t and taking the limit, ∆t → 0, we obtain a set of ODEs: lf (x, l). (6.12) Once the activity matrix and the transition rate functions are generated, the ODEs are immediately available. All of them can be obtained automatically by Algorithm 1.
For an arbitrary CTMC, the evolution of probabilities distributed on each state can be described using linear ODEs (see [18], page 52). For example, for the aggregated CTMC underlying a PEPA model, the corresponding differential equations are dπ dt = Q T π, (6.13) where each entry of π(t) represents the probability of the system being in each state at time t, and Q is an infinitesimal generator matrix corresponding to the CTMC. Clearly, the dimension of the coefficient matrix Q is the square of the size of the state space, which increases with the number of components.
The scale of (6.13), i.e. the number of the ODEs, depends on the size of the state space, so it suffers from the state sapce explosion problem. In contrast, the ODEs (6.12) reflect the evolution of the population of the components in each local derivative, so the scale is only determined by the number of local derivatives and is unaffected by the size of the state space. Therefore, it avoids the explosion problem. But the price paid is that the ODEs (6.12) are generally nonlinear due to synchronisations, whereas (6.13) is linear.
This paper emphasises the consistency between the fluid approximation and the aggregated CTMC. Obviously, the CTMC depends on the starting state of the given PEPA model. By altering the population of components presented in the model, which can be done by varying the initial states, we may get a sequence of aggregated CTMCs. Moreover, the homogenous property that the transition rate function satisfies, indicated in Proposition 4, identifies the aggregated CTMC to be density dependent.
Definition 8. ( [19]). A family of CTMCs {X n } n is called density dependent if and only if there exists a continuous function f (x, l), x ∈ R d , l ∈ Z d , such that the infinitesimal generators of X n are given by: q (n) x,x+l = nf (x/n, l), l = 0, where q (n) x,x+l denotes an entry of the infinitesimal generator of X n , x a numerical state vector and l a transition vector.
This allows us to immediately conclude the following conclusion: Theorem 4. ( [7,17]) Let {X n } be a sequence of aggregated CTMCs generated from a given PEPA model (by scaling the initial state), then {X n } is density dependent.
Since both ODEs and density dependent CTMCs can be derived from the same PEPA model through the same activity matrix and transition rate functions, it is natural to believe some kind of consistency between them. In fact, according to Kurtz's theorem [20], the complete solution of some ODEs can be the limit of a sequence of Markov chains. Such consistency in the context of PEPA has been previously illustrated for a particular PEPA model in [21], and subsequently generalised to general models in [7] and [17]. The result presented below is extracted from [17], in which the convergence is in the sense of almost surely rather than probabilistically as in [21] and [7].
Theorem 5. ( [17]) Let x(t) be the solution of the ODEs (6.12) derived from a given PEPA model with initial condition x 0 , and let {X n (t)} be the density dependent CTMCs with X n (0) = nx 0 . Then for any t > 0, lim n→∞ sup u≤t X n (u)/n − x(u) = 0 a.s. (6.14) This theorem justifies the fluid approximation by manifesting the consistency between this approach and the corresponding CTMCs for a general PEPA model. Furthermore, if there is no synchronisation contained in the model then the derived ODEs (6.12) becomes linear, and (6.13) and (6.12) coincide except for a constant factor. Moreover, the fundamental results on the fluid approximation of PEPA models such as the existence, uniqueness, boundedness and nonnegativeness of the ODEs' solution, as well as the solution's asymptotic behaviour, have been obtained. In particular, the convergence of the ODEs' solution as time tends to infinity, has been proved under a condition, which is revealed to relate to some famous constants of Markov chains such as the spectral gap and the Log-Sobolev constant. For more details about these stories, please refer to [22]. As for performance derivation via this approach, please see [6].

Conclusions
In this paper we have demonstrated a schema, which bridges the syntactic and numerical representation, as well as the local definition and global analysis for a PEPA model. Computational approaches and associated algorithms developed based on the schema have been presented, which can help to relieve the state space explosion problem for large scale models. For other stochastic process algebras, similar numerical representation schema can be established and expected to benefit relevant performance modelling.