Compact Markov-modulated models for multiclass trace ﬁtting

Markov-modulated Poisson processes (MMPPs) are stochastic models for ﬁtting empirical traces for simulation, workload characterization and queueing analysis purposes. In this paper, we develop the ﬁrst counting process ﬁtting algorithm for the marked MMPP (M3PP), a generalization of the MMPP for modeling traces with events of multiple types. We initially explain how to ﬁt two-state M3PPs to empirical traces of counts. We then propose a novel form of composition, called interposition , which enables the approximate superposition of several two-state M3PPs without incurring into state space explosion. Compared to exact superposition, where the state space grows exponentially in the number of composed processes, in interposition the state space grows linearly in the number of composed M3PPs. Experimental results indicate that the proposed interposition methodology provides accurate results against artiﬁcial and real-world traces, with a signiﬁcantly smaller state space than superposed processes.

The MMPP is a special case of the Markovian Arrival Process (MAP) ( Neuts, 1979 ) in which the departure process is a modulation of Poisson processes and the active process is chosen according to the state of a continuous-time Markov chain.In this paper, we propose a generalization of the MMPP, which we call the marked MMPP (M3PP), and develop a scalable methodology for fit-ting M3PPs to empirical datasets.The M3PP can be regarded as a specialization of the MMAP, the marked extension of the MAP ( He & Neuts, 1998 ).A marked point process associates to each arrival a class label ( He & Neuts, 1998;He, 2001 ), thus it is useful to model traces with events of multiple classes (e.g., read and write requests in disk drives).Marked processes are also important in the analysis of priority queues and multiclass models ( Buchholz, Kemper, & Kriege, 2010;Horváth et al., 2009;Horváth, 2012;Houdt, 2012 ).However, few techniques exist for their fitting and they all focus on matching moments of the inter-arrival time process for marked MAPs (MMAPs) ( Buchholz et al., 2010;Horváth et al., 2009 ).Therefore, no techniques exist yet for fitting marked Markov-modulated processes to count data.Still, to limit monitoring overheads, only count data can be extracted from certain classes of computer and communications systems (e.g., network links).This motivates the investigation into methodologies to fit marked count data.
In this paper, we fill this gap by developing approximate fitting algorithms for the counting process of the M3PP.In particular, we first explain how to approximately fit a two-state MMPP counting process and then extend the idea to two-state M3PPs.The proposed approach is applicable to traces with aggregated or coarse-grained measurements, which cannot be analysed using approaches based on inter-arrival times, since these require moments from a trace recording all the arrivals.The main drawback is that counting processes are mathematically less tractable than interarrival processes, therefore one normally focuses on low-order moments of counts.
Two-state models are often insufficient to fit complex traces, therefore we also study the approximate fitting of large M3PPs.In the single class setting, a known limitation of MMPPs is the inability to simultaneously fit many statistical descriptors due to the non-linearity of their underlying equations ( Bodrog, Heindl, Horváth, & Telek, 2008;Heindl, Horváth, & Gross, 2006;Horváth & Telek, 2009 ).This has led to the definition of several approaches to fit complex traces by composing multiple small-sized MMPPs or MAPs using Kronecker operators ( Andersen & Nielsen, 1998;Casale, Zhang, & Smirni, 2010;Horváth & Telek, 2002 ).These methods employ composition operators for moment fitting, offering a different trade-off between computational cost and fitting accuracy compared to fitting methods based on the EM algorithm ( Breuer, 2002;Horváth & Okamura, 2013;Klemm, Lindemann, & Lohmann, 2003 ).In particular, the superposition operator allows one to describe a trace by the statistical multiplexing of several MMPPs, at the expense of an exponential growth of the number of states in the resulting process ( Sriram & Whitt, 1986 ).This state space explosion is an obstacle for the application of MMPPs and MAPs to modeling real systems; for example it considerably slows down, or even renders infeasible, the numerical evaluation of queueing models by matrix geometric methods ( Bini, Meini, Steffé, Pérez, & Houdt, 2012;Pérez, Velthoven, & Houdt, 2008 ).
In this paper, we tackle the state space explosion problem of superposition by showing that M3PPs admit a particular form of composition, which we call interposition , that enables several MMPPs to share the same state space without mutually affecting the marginal distributions of their counting processes.However, interposition introduces spurious covariances between class arrivals that may be seen as the cost of the state-space reduction.The interposition method defines an original approach to build large Markov-modulated processes, in which a set of J two-state M3PPs is merged into a single M3PP process with just J + 1 states and without affecting the marginal counting processes of the initial M3PPs.We identify general conditions for interposition to result in a valid MMAP and observe that these conditions can be satisfied by M3PPs, but not by general MMAPs.The ability to interpose processes makes a case for using M3PPs, instead of general MMAPs, for fitting count data.We then define a method to automatically identify the M3PPs to be interposed and a mixed-integer linear program (MILP) that can help in automatically identifying the order of interposition of the M3PPs.
We conclude the paper by reporting fitting experiments for a set of artificial and real-world traces.We show that the proposed M3PP fitting algorithms are widely applicable and run efficiently even in the case of approximate fitting.We also find that interposition is much more scalable than superposition, while retaining comparable accuracy.
Summarizing, our main contributions are as follows: • Section 3 defines fitting algorithms for the counting process of two-state M3PPs with arbitrary number of classes.Our formulas are in closed-form for exact fitting and use quadratic programming for approximate fitting.As a by-product, our approach also introduces the first infeasibility adjustment methodology for approximate fitting of unmarked MMPPs.• Section 4 introduces the new notion of interposition.This is an aid to compose multiple M3PPs, while preserving their statistical properties, as we rigorously establish in Theorem 1 .• Section 5 develops a methodology for fitting empirical traces using the interposition operator.
The paper reports in Section 6 an experimental study on random models and a real trace validating the effectiveness of the proposed models and algorithms.In addition to the above, a description of necessary background is given in Section 2 .Section 7 concludes the paper.

Model and notation
An m -state MMPP is specified by a continuous-time Markov chain (CTMC) with irreducible infinitesimal generator Q , having m states, and rate vector ( λ(1), . . ., λ( m )).When the CTMC is in state j , the MMPP generates arrivals according to a Poisson process with rate λ( j ).The effect of the modulating action of the underlying CTMC is to enable the modeling of non-Poisson, possibly nonrenewal, arrival processes.We assume the underlying CTMC to be initialized according to its steady-state distribution so that the process is time-stationary.For ease of comparison with the literature, we use throughout the paper the MAP ( D 0 , D 1 ) notation, where for We define a M3PP[ K ] as a MMPP in which arrivals are marked with one out of K available classes.This may be seen as a marking of the Poisson processes of the MMPP where one decomposes each rate λ( j ), 1 ≤ j ≤ m , into arrival rates q j , k λ( j ), 1 ≤ k ≤ K , subject to k q j,k = 1 , q j , k ≥ 0. When an arrival is generated by a Poisson process with rate q j , k λ( j ) it is said to be of class k .Equivalently, in matrix notation, augmenting a MMPP ( D 0 , D 1 ) with marks defines a set of matrices . Also, ( D 0 , D 1 ) is the embedded MMPP, which we refer to as the unmarked process .
In the rest of the paper we often deal with processes having m = 2 states.In this case, for readability, we use λ and λ in place of λ(1) and λ(2) and q k and q k in place of q 1, k and q 2, k .

Problem statement
Let us consider a M3PP[ K ] and let n k ( t ) ≥ 0 be the number of arrivals of class k at time t after initialization, subject to n k (0) = 0 for all classes.The counting process of the M3PP[ K ] is the CTMC with state n (t) = (n 1 (t) , n 2 (t ) , . . ., n K (t )) .We shall denote by n (t) = K k =1 n k (t) the aggregate arrival count.Given an initial state probability vector, the evolution over time of this process is characterized by a matrix P ( n , t ), with element p i , j ( n , t ) in row i and column j representing the probability that a M3PP initialized in state i is in state j at time t with an arrival count n ( t ).We refer to P ( n , t ) as the counting process matrix .The counting process evolves according to the Kolmogorov forward equations where e k is a column vector of zeros except for a one in position k .From this equation, it is simple to derive factorial moment functions, from which moments of the counting process can be obtained in closed form He and Neuts (1998) .For example, this method yields the mean arrival count of class k at time t as where π is the stationary distribution of the embedded CTMC, πQ = 0 , π1 = 1 , μ k is the mean arrival rate of class k in the timestationary process, 0 and 1 are respectively column vectors of zeros and ones.The variance of class k in counts (also called variancetime curve) is ( He & Neuts, 1998 ) Var where where Var [ n k (t) + n h (t)] can be obtained from (3) by replacing μ k with μ k + μ h and D 1, k with D 1 ,k + D 1 ,h .
The M3PP[ K ] fitting problem is to find a representation 4) match, to the best possible extent, the corresponding empirical moments at given timescales t .Such a representation needs to be valid, meaning that all rates need be non-negative real numbers and the generator Q of the embedded CTMC must be valid and irreducible.

Two-state MMPP fitting
In order to fit a two-state M3PP[ K ], we propose to separately fit the MMPP for the unmarked process and the class markings of the M3PP.Since MMPP fitting is well-understood, we just summarize the MMPP fitting method proposed in Heffes and Lucantoni (1986) .This algorithm receives in input the following descriptors of the counting process: where μ is the rate, I ( t ) is the index of dispersion for counts ( Sriram & Whitt, 1986 ) at a given timescale t , I is the asymptotic index of dispersion value for t → ∞ , and μ (3) ( t ) is the third central moment of counts at timescale t .The objective is to fit the rates used in the MMPP representation subject to all the parameters being non-negative and to λ + λ > 0 .Note in particular that we exclude the trivial cases λ = λ and r = r = 0 , where the MMPP degenerates into a Poisson process, which does not require a two-state model.
Fitting the above descriptors to a two-state MMPP requires to solve for the unknowns r , r , λ, and λ using a nonlinear system composed by the following four equations ( Heffes & Lucantoni, 1986 ) where x = r + r , x = r − r , t is an arbitrary finite timescale at which we want to fit the index of dispersion, and h = and g (3) (1, t ) is the third factorial moment of counts at timescale t .We point to Heffes and Lucantoni (1986 , Eq. (14d)) for explicit expressions of g (3) (1, t ).Heffes and Lucantoni (1986) propose the following fitting method.First, compute x = r + r , solving at an arbitrary timescale t = t 1 the fixed point equation obtained by rewriting the third equation appearing in (5) .Then, compute h at a second arbitrary timescale t = t 2 .If h = 0, the fitting formulas are then given by where y = (I − 1) μx 3 (2 h 2 ) −1 .If h = 0 , the explicit fitting formulas instead become The main drawback of this method is that the formulas can return negative rates for some combinations of the input parameters.In this case, approximate fitting techniques are required.Yet, to the best of our knowledge these are not available in the literature on MMPP counting process fitting.

Fitting the counting process of a two-state M3PP[ K ]
In this section, we develop a method to fit two-state M3PP[ K ]s with representation As mentioned before, the availability of exact methods to fit unmarked MMPPs justifies an approach where we separately fit the embedded MMPP first, followed by the marking process that defines the M3PP.In practice, this means that we first fit D 0 and D 1 using a MMPP fitting algorithm applied to the unmarked trace.Then we determine the individual D 1, k matrices using the marked trace and subject to the condition D 1 = K k =1 D 1 ,k , which ensures that the embedded MMPP does not change.In Section 3.1 , we discuss how to perform approximate fitting of the unmarked process using an extension of the algorithm by Heffes and Lucantoni.In Section 3.2 , we discuss the fitting of the marked process.

Step 1: approximate MMPP fitting
The fitting algorithm of Heffes and Lucantoni described in Section 2.3 cannot be applied to cases where the input set of descriptors μ, I ( t ), I and μ (3) ( t ) is infeasible for the considered MMPP.
We have found this to happen frequently in real traces, and even though one may repeatedly attempt to fit different timescales t 1 and t 2 until finding a feasible set of descriptors, we believe that better solutions are needed and possible.We have not found previous works addressing this problem, at least for fitting counting processes.Our approximation consists in sacrificing the degree of freedom of the third moment of counts μ (3) ( t ) to restore feasibility of all second-order descriptors.As before, we focus on cases where the M3PP process is not Poisson, thus we assume λ = λ , r > 0, r > 0.
We begin by characterizing the feasibility region of the index of dispersion for a two-state MMPP.
Proposition 1.In a two-state MMPP with λ = λ , the index of dis- persion satisfies I > I ( t ) > 1 for any timescale t.
Proof.From the second equation in (5) , it readily follows that I > 1 since λ = λ .Moreover, since x > 0, it follows from (6) that I > by ( 5) and the constraint I > I ( t ) we conclude that I ( t ) > 1.
The previous proposition states the well-known fact that MMPPs can only fit traces that have greater variability than a Poisson process.We now show that asking for r > 0 and r > 0, to exclude the case where the MMPP is a Poisson process, is implied by the previous statement and thus always verified.
Proposition 2. For I > 1, the conditions r > 0 and r > 0 are always satisfied.
For the case h = 0 , the result follows immediately given that x > 0 since r = r = 0 would imply a Poisson process that would have I = 1 against the assumption.Therefore, provided that I > I ( t ) > 1, feasibility follows by ensuring that arrival rates are non-negative and λ + λ > 0 .Using (5) , we can reformulate this requirement as: Note that this expression implicitly characterizes the feasible region for the rate μ and, via the h term, for the third moment g (3) (1, t ).If one of these conditions is not met, we propose to relax the fitting method by ignoring the matching of the third moment g (3) (1, t ) as follows.Consider first the following lower bound on r .
Proposition 3. Without loss of generality, assume in which x = r + r is the solution of ( 6) .
Proof.We substitute the first equation in ( 5) into the second one and find after simple algebra Obtaining λ from the first equation in ( 5) and equating to (9) , we get the feasibility constraint The result follows using r = x − r and solving for r , which yields the lower bound u .
Any value u ≤ r < x leads to a feasible solution, thus we can for example choose the middle point of this interval r = u + ( u − x ) / 2 .
Our suggestion of the middle point is convenient to keep the formulas in closed-form, however other choices are possible.After choosing r , the other parameters are easily obtained by setting r = x − r and using ( 9) and ( 10) to obtain λ and λ.
Summarizing, provided that I > I ( t 1 ) > 1, the approximate method we have defined fits μ, I ( t 1 ) and I exactly at the expense of sacrificing an infeasible third moment μ (3) ( t 2 ), where t 1 and t 2 are arbitrary timescales.

Step 2: fitting the marking process
In the second step of the fitting algorithm we determine the D 1, k matrices.After Step 1, all the statistical properties of the unmarked process are constrained by the D 0 and D 1 matrices of the MMPP, thus the focus of this step is to fit the class-specific properties and the class covariances.We consider both exact and approximate fitting methods.
There are several ways of choosing which empirical descriptors to fit with a M3PP and each choice leads to equations that may differ in tractability compared to other choices.We have experimented with several possibilities, and we have found that fitting a set of central moments, such as mean class rates μ k , class variances, or covariances, typically leads to a difficult non-convex formulation.Conversely, we have found more efficient to fit a two-state M3PP[ K ] using the mean class rates and per-class contributions to the asymptotic index of dispersion I .Other efficient approaches are possible, such as fitting mean class rates and relative covariances or fitting mean class rates and differences between the class variances.However, we here focus on the first method, which we believe provides the best combination of efficiency (quadratic programming) and ease of interpretation.

Fitting mean class rates and per-class contributions to the index of dispersion
We begin by considering the more general problem of fitting the mean class rates μ k and g k (t We then specialize the method to the index of dispersion, using the fact that Our goal is to compute the D 1, k matrices from the g k ( t ) values, given D 0 and The problem under consideration is to determine the values of the probabilities q k and q k that uniquely define the D 1, k matrices, given that λ and λ are known from the D 1 matrix fitted in Step 1.
Exact matching.We now give formulas to fit the probabilities q k and q k from μ k and g k ( t ).Note that the arbitrary timescale t used in the statement does not need to be equal to the timescales used for fitting the embedded MMPP.
Proposition 4. Given ( D 0 , D 1 ) and, for each class k , μ k and g k ( t ) at an arbitrary timescale t , the parameters of the M3PP can be computed as follows: and x = r + r is the solution of (6) .
Proof.The expressions of μ k and g k ( t ) for a two-state M3PP [ K ] process are found from the definitions to be where F 1 and F 2 are constant coefficients given D 0 , D 1 and the reference timescale t .Solving (13) for q k we obtain Substituting ( 15) into ( 14) , we obtain the fitting formulas ( 11) and (12) .
We are now ready to determine the contributions to the index of dispersion.Note that Proposition 4 may also be used to fit the contribution of class k to I ( t ) in the embedded MMPP, i.e., . This is because we can rewrite the fitting formulas as Combining these expressions with the requirements q k ≥ 0 and q k ≥ 0 , we find after simple passages this lower bound required to hold for the feasibility of the per-class contributions to the asymptotic index of dispersion For two-state M3PPs, the minimum value of the right-hand side is achieved for Poisson processes, where Approximate matching.Some values of the descriptors may lead to unfeasible values of the parameters (e.g., negative probabilities).In this case, we choose to fit the per-class arrival rates μ k exactly and find the feasible values ˜ G k that minimize the following L 2 -norm and subject to the constraints These constraints are derived using equations ( 16) and ( 17) for q k and q k .In particular, (20) stems directly from the constraints q k ≥ 0 and q k ≥ 0 , whereas (21) follows from the conditions K k =1 q k = 1 and K k =1 q k = 1 .The above optimization program may also be used for fitting mean class rates μ k and the contributions g k ( t ) to the index of dispersion I ( t ).In order to do so, it is sufficient to replace the coefficients c i , j with the coefficients f i , j , G k with g k ( t ), and ˜ G k with ˜ g k (t) .
In either case, the program has a convex quadratic objective function, thus its minimizer can be efficiently obtained using standard quadratic programming solvers.Once the problem is solved and the feasible values G k or ˜ g k (t) are obtained, the parameters q k and q k are readily fitted using ( 16) and (17) .

Compositional fitting of M3PP[ K ]s
In the previous section, we have defined a general purpose fitting method for two-state M3PPs.We now consider the problem of exploiting the additional degrees of freedom of a M3PP[ K ] to increase the flexibility of the fitting.In this case, exact fitting methods capable of exploiting all available degrees of freedom do not exist even for unmarked processes.Therefore we focus on compositional fitting, where one builds a large process by composition of smaller processes that are simpler to fit to count data.The drawback of compositional approaches of this kind is that they use just a few degrees of freedom of a general MMAP in return for ease of fitting.
We first review and generalize superposition methods for unmarked MMPPs.Afterwards, we introduce a novel form of composition, called interposition , which offers a different trade-off between accuracy and compactness of the representation.Note that we do not consider methods that are specific to inter-arrival processes, such as the Kronecker product composition ( Casale et al., 2010 ).
, in which denotes the Kronecker sum operator ( Brewer, 1978 ).Superposition naturally arises in networking to describe the traffic process obtained by merging separate traffic flows, each described by a MMPP.The superposed process defines the multiplexing of J channels, each with inter-arrival times described by the j -th MMPP.This process has mean arrival rate and variance-time curve equal to the sum of mean arrival rates and variance-time curves of the individual MMPPs.The index of dispersion is a weighted sum of the IDCs of the individual MMPPs, i.e., I(t ) = j ( μ j /μ) I j (t) , where μ = j μ j is the mean arrival rate of the superposition.Also, if MMPP j has m j states, the superposed process has J j=1 m j states, thus its size grows exponentially with J .
Marked case .Let K = { 1 , . . ., K} be a set of classes.We consider J M3PPs and assume that M3PP j generates arrival of classes In this case, some M3PPs contribute to arrivals of class k , thus where 0 is here a matrix of all zeros of order m j .The statistical properties of the embedded unmarked process are obtained as in an unmarked superposition.Lastly, note that if each M3PP has two states, the resulting process is a M3PP with 2 J states and K classes.Thus, the main drawback of the superposition method is the state-space explosion, which limits the composition to a small number of processes.

Interposition
We now propose a new form of composition for Markovmodulated processes that tackles the state-space explosion problem of superposition.Informally, our idea is to define an operator by which several M3PPs can by defined upon the same state space, without interfering with their respective marginal counting processes.This allows us to preserve in the interposition the counting process properties fitted in isolation for each composed M3PP.We characterize the main feature of this new composition operator in Theorem 1 , given later.Note that the results in this section are also applicable to unmarked MMPPs, and thus represent advances also for unmarked processes.
Consider a set of J two-state M3PPs where the i th process has K i classes.Assume without loss of generality that classes are labelled such that each class appears in one and only one M3PP.The M3PPs have representation where we define the probabilities k ≥ 0 , and the rates r i = J j= i α j and r i = i j=1 β j , for given constants α j ≥ 0 and β j ≥ 0. Equivalently, given the values of the r i and r i constants, we can define the rate differences α i = r i − r i +1 and β i = r i − r i −1 , with boundary values β 1 = r 1 and α J = r J .
We now define the interposition operator for M3PPs.Given the Please cite this article as: G.
, where I n is the identity matrix of order n , class k = i −1 j=1 K j + c is the class in the interposed process associated to class c of the ith composed M3PP, and the diagonal elements of D * 0 are such that Throughout this section, we denote by K i the set of class indexes in the interposed process associated to the i -th composed M3PP, and by The interposed process may be seen as a M3PP modulated by a CTMC with an infinitesimal generator ( Bolch, Greiner, de Meer, and Trivedi, 1998 , Chap 4).Specifically, for each of the initial M3PPs, we can define a partition of the states of the interposed process such that, by exact aggregation of the CTMC Q * one recovers the corresponding CTMC of the initial M3PP.For example, we may consider the aggregation where is the infinitesimal generator of the i th composed M3PP.Similarly, for each partition, the definition of D * 1 ,k ensures that the departure rates are identical to the ones in the composed M3PP.
In order to prove that the marginal counting process of each composed M3PP, for all classes k ∈ K i , is preserved by the interposition operator, we require additional notation.Let P ( n i , t ), n i = (n 1 , . . ., n K i ) , be the 2 × 2 counting process matrix for the i -th M3PP at time t .Similarly, let P * ( n , t ) be the (J + 1) × (J + 1) counting process matrix for the composed M3PP, where n = ( n 1 , . . ., n J ) .
Let 0 n and 1 n be row vectors of size n of all zeros and all ones, respectively, and define the aggregation matrix ( Buchholz, 1994 ) For an arbitrary (J + 1) × (J + 1) matrix X , the 2 × 2 matrix gives the aggregation of X into the two states associated to the i -th composed M3PP.Therefore is the counting process matrix aggregated onto the states of the ith composed M3PP.From this matrix, we can readily compute the marginal counting process matrix of the i -th M3PP as where the summation is over all classes not associated to the i -th M3PP.Using this notation, we prove the main characterization result for the interposed process, which states that the interposition operator does not affect the marginal counting process for the i -th M3PP.
Theorem 1. Assume all processes to be time-stationary, then Proof.For the composed M3PP it is simple to verify that V i has the following properties Note that conditions of this kind naturally arise in the study of minimal representations of Markovian and Rational Arrival Processes ( Buchholz & Telek, 2013 ).
Let us then consider the Kolmogorov forward equation for P * ( n , t ).Pre-multiplying (1) by W i and post-multiplying by V T i , we obtain where we have used ( 23) to simplify the expression and the notation n − e k indicates the removal of a job of class k from n .Now plugging the identity ,k and using again (23) we get We now consider the marginal probability P * i ( n i , t ) , for which the last expression implies Since the infinite summations in the last term of ( 24) are on the same class indexes ( K i sets), and P * ( n − e k , t ) = 0 when n k = 0 , by symmetry the double summation vanishes.Thus (24) reduces to which is identical to the Kolmogorov forward equation for the counting process of the i -th composed M3PP.In order to prove this statement, we just need to show that the initial conditions of the Kolmogorov forward equations are the same, i.e., P ( n i , 0) = Since we are considering time-stationary processes, the initial state of the interposed process is determined by the equilibrium distribution of the CTMC with generator Q * .Conversely, for the i -th M3PP this is given by the equilibrium distribution of the CTMC with generator Q i .By (23) we see that for an arbitrary vector π Please cite this article as: G. Therefore, if we choose the initial distribution to be the timestationary distribution πQ * = 0 , π1 = 1 , we readily find that its aggregation corresponds to the time-stationary distribution of the i -th M3PP, i.e., π i Q i = 0 and π i 1 = 1 , where the last condition holds since we use exact aggregation.This concludes the proof.
We remark that (23) provides a condition for a general MMAP to define a valid interposition.It is simple to see that the diagonal structure of the D 1, k matrices in a M3PP satisfies these assumptions.While our analysis does not exclude that other kinds of MMAPs may be used for interposition, the conditions required to satisfy (23) do not seem to readily suggest alternatives other than the M3PP.This motivates the use of M3PPs as a building block for interposition.

Feasibility of an interposition
We now give examples of valid and invalid interposed processes.Consider the following two-state M3PP[2]s: satisfies the assumptions on the non-negativity of the rates α j and β j and yields the composed 3-state M3PP[4] B and E are the same process, because their states are ordered in a different way and this affects the computation of the α j and β j coefficients.This last example provides intuition on the fact that, to obtain a feasible interposed process, one may need to re-order the states of the M3PPs.In the next section, we describe an algorithm to find a feasible interposition of a given set of M3PPs, if one exists.

Class covariances
While the interposed process preserves the marginal counting properties of the composed M3PPs, when compared to superposition this comes at the expense of introducing spurious covariances between arrivals of different M3PPs.This is because, by definition of the interposed process, a transition in the state space of a composed M3PP also changes the current state of the other composed M3PPs.In order to quantify the magnitude of these covariances, we look at the asymptotic covariances, which contribute to the index of dispersion.Observe that, by plugging (3) into (4) , we find after some simplifications and observe that this by construction is a stochastic matrix with equilibrium vector π.Using the fact that πD 1 ,k = μ k , for any class k , and after determining the structure of the P matrices for the interposed process, it is possible to compute their Jordan canonical form, which after algebraic simplifications yields the formula σ k,h = r i r j (q j,h λ j − q j,h λ j )(q i,k λ i − q i,k λ i )(x j + x i ) where it is assumed that k ∈ K i and h ∈ K j , and i < j .Some remarks on the formulas are as follows: • When any of the two M3PPs tends to a Poisson process, a pair of departure rates at the numerator of ( 25) annihilates and σ k , h goes to zero.• The order of the denominator suggests that a way to reduce the covariance introduced by interposition may consist in spending degrees of freedom to maximize x i and x j in the embedded MMPPs, for fixed arrival rates.When applying such a scheme, one should however take into account the bounds in Proposition 3 , since an increase of the value of x i also reduces fitting flexibility in the embedded MMPP.• Since x i > 0 for any i , we note that the sign of the covariance is determined only by the differences between the per-class arrival rates within each of the composed M3PPs.

Fitting the interposed process
In this section we consider two issues that arise in compositional fitting based on interposition.First, we consider the decision problem involving the mapping of a marked trace into a set of two-states M3PPs.Then we show that the problem of finding a feasible interposition of a set of J second-order M3PPs, where the i -th M3PP models any subset K i of the K classes, can be formulated as a MILP and hence solved using an integer programming solver.

Fitting a marked trace into a set of M3PPs
Given a trace with arrivals of K classes, if the class arrivals are independent it is possible to fit the trace using a superposition of K independent MMPPs, one for each class, and then reduce the size of the resulting process using interposition.However, if the classes have a significantly large covariance, this method may not produce good results.Therefore, we propose two heuristic fitting methods to address these two cases.We assess the effectiveness of these methods later in Section 6 .

Independent Method
In this method, we ignore class covariances and fit each class into a separate two-state MMPP.The method is similar to superposition, but returns a much smaller M3PP[K] with K + 1 states, instead of 2 K states.Interposition uses the order of composition obtained from the MILP method that we present in Section 5.2 .

Covariance-based method
We initially build the asymptotic co-variance matrix t) , n h (t)] between each pair of classes.We approximate the asymptotic timescale with the largest finite t for which we can average counts over at least 100 samples.The user is then requested to specify a covariance threshold δ.We then iterate on the classes, in decreasing order of variance σ k , k .For each class k , we first determine all the classes h = k with σ h , k ≥ δ and record the decision to fit these classes into the same M3PP [K].
The algorithm then continues by analyzing in a similar way the other classes not already planned for fitting in any M3PP.After this stage, each class k is mapped to a unique M3PP i .In order to obtain a representation for each M3PP i , we run the algorithm given in Section 5.2 , which returns the parameters of the embedded MMPPs and their order of composition in the interposition.
Please cite this article as: G.  [m5G;June 24, 2016;13:55 ] With these, we can conclude by fitting the M3PPs with the method in Section 3.2.1 and applying the interposition operator.

Finding a feasible interposition of M3PPs
The definition of the interposed process implies that the D 1, k matrices are always feasible by construction.Thus, infeasible M3PPs may arise only if the D 0 matrix has some negative elements.Since the entries of the D 0 matrix depend only on the MMPPs embedded in the two-state M3PPs, a method to determine a feasible interposition should be run prior to fitting the MMPPs.
To find a feasible interposition, we consider permutations of the mapping of the rates to the states in the embedded MMPP and the order in which the J M3PPs are composed.Furthermore, when deciding the structure of the MMPPs, we exploit the degree of freedom given by ignoring the fitting of the third moment of counts.Let x i and u i , be the values of x and u for the i -th embedded MMPP.Sacrificing the fitting of the third moment allows us to decide the value for r i , given x i , which then implies r i = x i − r i .Recall from Section 3 that when ignoring the third moment of counts, infinite feasible MMPPs exist as long as u i ≤ r i < x i and r i = x i − r i , where x i is found by ( 6) .These bounds have been obtained in Section 3 under the assumption that λ i > λ i .While this assumption does not have any impact on the feasibility of individual M3PPs, it has an impact on the feasibility of the interposition.
Making the opposite assumption λ i ≤ λ i , we obtain the specular constraints u i ≥ r i > x i and r i = x i − r i .Thus, the above bounds will need to be simultaneously considered and only one will be active depending on the relative value of λ i to λ i , which we decide using a binary variable b i that is responsible for state order.Summarizing, the above bounds allow us to obtain a feasible interposition by choosing from a continuous set of feasible MMPPs.
Before solving the MILP formulation that decides the order of states and the order of composition, we assume that the fixed point Eq. ( 6) has been solved for each M3PP i for the values of x i and u i .For the MILP formulation, we then consider the following decision variables: • the rates r i > 0 and r i > 0 of the i -th M3PP, ∀ i ∈ {1, ..., J }; • the integer variables p i ∈ N , ∀ i ∈ {1, ..., J }, deciding the position of the i -th M3PP in the interposed process, i.e., the order of composition.
Let M be a large constant.We consider the following MILP formulation: minimize 0 (feasibility problem) Strict bounds are obtained by adding small tolerance to the corresponding formulation with ≤ inequalities.Observe that the objective function is a constant, since we just need a feasible solution.Constraint ( 27) imposes the bounds on r i , with different bounds depending on whether the state order of the i -M3PP is inverted or not.Constraint (28) imposes specular bounds on r i .Antisymmetry constraints on the auxiliary variables z i , j are set in (29) .Constraint (30) ensures that z i, j = 1 if and only if p i ≤ p j .This constraint together with the uniqueness constraint on p i ensures the transitivity property for the binary variables z i , j .Ordering constraints on r i and r i are expressed by ( 31) and (32) , respectively.Two M3PP processes cannot be in the same position due to (33) ; note that inequality constraints can be handled for integer variables in MILP.
After solving the MILP, we have the parameters r i and r i of each M3PP.Applying the fitting formulas of Section 3 we obtain the remaining parameters of the embedded MMPP, λ i and λ i , taking care to swap the two values when states are in inverted order, i.e., b i = 1 .Then we compute the class parameters q i , k and q i,k for each class of the i -th M3PP and build the interposition of the J M3PPs as in Section 4 .

Fitting random MMAPs
We begin by examining the applicability of two-state M3PPs and the interposition process in fitting the characteristics of randomly generated processes.We consider random M3PP[K]s, with m ∈ {4, 8, 16} states and K ∈ {2, 3, 4} classes.For each choice of m and K , we generate 100 random M3PP[K] as follows.First, we generate a random infinitesimal generator Q = D 0 + D 1 using uniform random numbers.Then, for each class k , we first compute a vector of m uniform random numbers u and set The expression used to compute D 1, k provides a set of processes with index of dispersion I that is 3-6 times larger than the squared coefficient of variation (SCV), as shown in Table 1 , which gives statistics for the randomly generated M3PPs.The statistics in each row are averaged across the 100 random instances and we report mean and standard deviation.The third columns gives the average ratio which quantifies the relative magnitude of cross-covariances.
M3PPs are fitted as follows.We first compute the statistical descriptors of the random process using theoretical expressions.Two-state M3PPs are then obtained using the method described in Section 3.2.1 .The interposition process is fitted using the two methods described in Section 5.1 and the feasibility method in Section 5.2 .The latter is implemented in MATLAB using Please cite this article as: G.  [m5G;June 24, 2016;13:55 ] YALMIP ( Löfberg, 2004 ) and the CBC branch-and-cut solver ( Coinor branch & cut project ).Superposition is implemented as by the definition after fitting an independent MMPP for each class.Here and in the following sections, counts used to fit a class are obtained from inter-arrival times between successive arrivals of that same class.
Let n k ( t ) be the number of arrivals of class k in t time units for the randomly-generated M3PP and let ˜ n k (t) be the same descriptor in the fitted model.The metrics used to assess the quality of the fitting are the ability of the model to capture the asymptotic class variances and covariances, as quantified by the absolute relative errors We focus on second order class descriptors of the counting process since mean arrival rates can always be fitted exactly.We also quantify the absolute relative error on the variance of the underlying unmarked process For randomly-generated models, the values of v ar , cov and unmark are averaged across all the models.The corresponding standard deviations are indicated with σ v ar , σ cov and σ unmark , respectively.
Lastly, we assess the total number of states m fit of the fitted model and the time T required for the fitting algorithm to return a valid M3PP.In applying the M3PP fitting method of Section 3.2.1 , we choose t = 10 E[ X] as arbitrary timescale, being E [ X ] the mean of the random M3PP[K], and we approximate asymptotic values using the timescale t ∞ = 10 4 E[ X] .We have also repeated the experiments with t = E[ X] and t = 100 E[ X] , but the results closely resemble the ones reported in this section.

Results
Results of the validation against random M3PPs are given in Table 2 .Remarks are as follows: • The results indicate that interposition has the same efficiency of superposition in capturing class variances, whereas two-state M3PPs are better at capturing class covariances.This aligns with expectations, given that interposition matches marginal counting processes, whereas the method presented in Section 3.2.1 for two-state M3PPs matches asymptotic covariances.
• It is fairly surprising to note the good performance of two-state M3PPs in fitting even large processes, with several states and classes.However, we conjecture this to be because the arrivals have the same inter-arrival time distribution of the embedded MMPP.In fact, we will show that when this assumption is removed in the fitting of real traces, two-state M3PPs perform worse.
• Interposition is on most cases superior to superposition.The number of states is significantly lower, without appreciable loss of accuracy in matching variances.It should be noted that, while in some cases interposition performs worse than superposition in matching covariances, superposition always return zero covariance, thus it is uninformative.Conversely, the error on covariance of interposition appears to decrease with growing number of states and classes.
• The covariance-based fitting method used for interposition typically fares better than the interposition of independent MMPPs.The approach has similar accuracy for variance matching, it is generally more accurate in describing covariances, and requires less states.• Computational times are small for all methods.However, there is an appreciable difference in the time to compute an interposition for large processes, due to the cost of solving the integer program introduced in Section 5.2 .Still, it should be noted that when the number of classes grows beyond K = 4 , superposition requires tens or even hundreds of states, thus it becomes far less tractable than interposition.

Fitting real traces
The analysis for random MMAPs is now repeated for fitting the BC-pAug89 and LBL-TCP-3 traces from the Internet Traffic Archive 1 .While these traces are commonly used in the literature of Markovmodulated processes, they are unmarked and thus cannot be readily used for validation of marked processes.In order to do so, we introduce a labelling that associates each arrival to a different bin of the histogram of packet sizes.We first consider the case where arrivals belong only to one of K = 2 classes, putting the separator at the p -th percentile of the inter-arrival time distribution, where p ∈ {25, 50, 75, 90}.In addition, we consider a separation of the histogram into five bins, defined by the same set of percentiles, which leads to a model with K = 5 classes.
Results are given in Tables 3 and 4 .The trends in the two tables are qualitatively similar and indicate that interposition can in most cases fit better the traces than a 2-state model or superposition.There are also several instances on which the composed M3PPs have better cov than the 2-state M3PP.This suggests that the good results on random instances for 2-state models may be biased by the fact that inter-arrival times of different classes are identically distributed in the random instances.This does not happen with the real traces, where the histogram bins do not overlap with each other.
It can also be noted that, as the number of classes grows to K = 5 , the performance of the different methods get closer to each other, presumably due to the increased difficulty in matching class covariances.

Queuing analysis
Lastly, we consider a queueing analysis application.We use again the LBL-TCP-3 trace, but with the class marking defined in Buchholz et al. (2010) .This marking defines 4 classes, each associated to a different bin of the histogram of packet sizes, and having service rates μ 1 = 300 , μ 2 = 250 , μ 3 = 200 , μ 4 = 100 .We simulate a queue with exponentially distributed service times, infinite buffer capacity, first-come first-served scheduling, and arrivals fed by the LBL-TCP-3 trace.In the simulations, we use 10 8 samples and record mean queue-lengths for each class.The results are compared with the predictions obtained from a MMAP/M/1 first-come first-served queue fed by a MMAP obtained by the following fitting methods: • 2-state: the fitting of a M3PP[4] (2 states).
• interpos-indep: interposition of the independent MMPPs used for the superposition, one per class (5 states).• interpos-cov: interposed process obtained by the covariancebased method, which interposes 3 M3PP[2], the first for classes 1 and 3, which have the largest covariances, the second for class 2, the third for class 4 (4 states).We have observed that different choices of the timescales used for fitting have a quite visible effect on the rate at which mean queue-lengths build up.Therefore, for complex traces such as LBL-TCP-3, we recommend to perform some calibration tests in order to find the best assignment of the arbitrary timescale t = t 1 = t 2 used in M3PP fitting.For this trace we perform calibration by varying t in {0.01, 0.05, 0.1, 0.5, 1, 5, 10, 50, 100} seconds and approximating the asymptotic timescale in each experiment as t ∞ = 10 t.Then we select the optimal timescale as the value for which superposition predicts the most accurate aggregated mean queue-length at 90 percent utilization.Using this calibration, we settle on t = 50 seconds and t ∞ = 500 seconds.
The resulting MMAP/M/1 queueing systems are analyzed using Q-MAM ( Bini et al., 2012 ).Results are illustrated in Fig. 1 .The simulation results indicate that mean queue-lengths grow exponentially, with classes 1-3 growing up to a few thousand jobs, whereas class 4 grows up to 86 jobs.The trends indicate that are all methods are accurate in heavy-load, above 60 percent uti-lization.Interposed processes are the most effective in capturing the job growths for the dominating classes 1-3.In heavy load interpos-cov is more accurate than interpos-indep, e.g., the mean class-2 queue length at 90 percent utilization predicted by interpos-cov is 3442 jobs, against a simulated value of 3881 jobs, whereas interpos-indep predicts 2196 jobs and superposition predicts 2024 jobs; the 2-state M3PP is also accurate in heavy load, with a prediction of 3416 jobs at 90 percent utilization, but this method is significantly worse than the other methods at lower loads.
Summarizing, by comparing against superposition as a baseline, the interposition method performs best in fitting real-world traces, delivering a more accurate estimate in loads that exceed 60 percent utilization.Combined with the previous validations, this outcome suggests that the interposition methods introduced in this paper can be helpful for real-world fitting and queueing analyses.
Please cite this article as: G.

Conclusion
In this paper, we have presented novel methods to fit multiclass arrival processes using marked Markov-modulated Poisson processes (M3PPs).We have defined exact and approximate algorithms to fit two-state M3PPs with arbitrary number of classes and introduced a new composition operator, called interposition, which enables composing several M3PPs while preserving their marginal counting processes.The state space of the interposed process grows linearly in the total number of composed M3PPs, instead than exponentially as in ordinary superposition.Experiments reveal that the interposed process can be effective in fitting real-Please cite this article as: G.
. The g k ( t ) terms have the simple interpretation of modeling the contribution of class k to the variance-time curve of the unmarked process, since σ 1) .The per-class arrival rates and variances of counts in the interposed process match the corresponding statistics of the two classes of A and the two classes of B .Conversely, the interposition of processes A and C =( C 0 , C 1 , 1 , C 1 , 2) is invalid because the non-diagonal rates of C 0 are both greater than the corresponding rates in A 0 .Similarly, the interposition of A and E

Table 1
Statistics for randomly generated M3PPs.Each row refers to 100 random instances.Mean μ and standard deviation σ are reported as ( μ, σ ).

Table 3
Fitting results for the BC-pAug89 trace for different percentile cut-points of the packet size histogram.