Transient analysis of piecewise homogeneous Markov fluid models

Piecewise homogeneous Markov fluid models are composed by homogeneous intervals where the model is governed by an interval dependent pair of generators and the model behaviour changes at the boundaries. The main difficulty of the transient analysis of piecewise homogeneous Markov fluid models is the appropriate description of the various boundary cases. The paper proposes an analytical approach to handle the wide variety of the possible boundary cases in a relatively simple to describe and implement manner.


Introduction
Markov fluid models (MFMs) gained significant popularity in modeling telecommunication systems in the 1980's (Anick et al. 1982). The first methodology to analyze the behaviour of such systems was based on spectral decomposition (Kulkarni 1997). In 1999, Ramaswami initiated a research line to analyze stochastic fluid models via matrix analytic methods (Ramaswami 1999), while Akar and Sohraby recommended the use of purely numerical matrix iterative methods (Akar and Sohraby 2004). Both methods provide numerically stable analysis, e.g., for the stationary distribution of the fluid level, but the approach based on matrix analytic methods gained more popularity due to the fact that it provides a stochastic interpretation of the considered performance measures. This work is partially supported by the OTKA K-123914 and the TUDFO/51757/2019-ITM grants.
In a series of consecutive papers the stationary da Silva Soares and Latouche 2006) and the transient (Ahn and Ramaswami 2004, 2005 analysis of homogeneous (finite and infinite) MFMs has been investigated. At the same time, the ingredients of the computational methods for various performance measures of fluid models has also been enhanced (Remiche 2005;Bean et al. 2009). Especially, the combination of two matrix exponential terms for describing the behaviour of finite buffer homogeneous fluid models got established.
Motivated by several practical examples, e.g. (Mandjes et al. 2003), the analysis of homogeneous MFMs has been extended to the analysis of piecewise homogeneous models, where the characterizing matrices of the model are constant in a region of the fluid level, but they might differ region by region. The terminology used to describe this set of models is rather diverse: "level dependent evolution" (da Silva Soares and Latouche 2009), "multi-layer" , "multi-regime" (Kankaya and Akar 2008), etc. We refer to such modes as piecewise homogeneous Markov fluid models (PHMFMs) and their homogeneous intervals as regions.
The main difficulty in the evaluation of PHMFMs comes from the potential sign change of the fluid rate at the region boundaries. As a result probability mass can develop at region boundaries (the fluid rate is positive below the boundary and negative above) and there might be cases when the evolution of the fluid process is not uniquely defined by the fluid rate (the fluid rate is positive above the boundary and negative below). These cases are also covered with a wide range of terminologies. Just to mention some, those subset of states are referred to as "absorbing" and "insulating" in (Gribaudo and Telek 2008), "absorbing" and "repulsive" in (Kankaya and Akar 2008), "sticky" and "repellent" in (da Silva Soares and Latouche 2009), "isolating" and "emitting" in (Chen et al. 2002).
There are several approaches to specify the boundary behaviour at internal boundaries. (da Silva Soares and Latouche 2009) assumed that the continuous region above a boundary determines the behaviour of the boundary in repulsive states, additional to the generator and fluid rate matrices describing the behaviour inside a region, (Kankaya and Akar 2008) introduces generator and fluid rate matrices associated with the boundaries. More general boundary behaviours are introduced and analyzed in  with the introduction of appropriately defined additional probability matrices characterizing the evolution of the fluid process at the environments of the boundaries. To keep the complexity of the analysis reasonably simple, we adopt the terminology and the boundary behaviour used in Kankaya and Akar (2008).
The stationary analysis of PHMFMs is based on performance analysis of individual regions and the solution of a linear system of equations, which can be described by a large coefficient matrix, e.g., (da Silva Soares and Latouche 2009), page 1048 and (Gribaudo and Telek 2008), Fig. 4 , and can be solved at once or by an iterative approach region by region.
Our Laplace transform domain transient analysis of PHMFMs follows a similar structure. The performance measures of the analysis of individual regions are available, e.g. from (Ahn et al. 2007;Bean et al. 2009). Based on that, an initial and final fluid level dependent counterpart of the stationary analysis is needed, since the stationary distribution is an initial condition independent measure, while the transient distribution depends also on the initial condition. The analysis approaches available for homogeneous (with infinite  and finite  buffer) MFMs describe the transient behaviour on the level of matrix blocks in Laplace transform domain using explicit expressions. The extension of this approach for PHMFMs gets prohibitively cumbersome because the proper description of the boundary behaviour at internal boundaries requires the consideration of all possible cases of sign changes separately. As an important contribution of this paper, to overcome the limitations of the explicit approach, we apply an equation systems based implicit description of the required performance measures, whose analytical description and implementation remain feasible due to easy to describe matrix block operations (see Theorems 1-3).
This paper focuses on the Laplace transform domain transient analysis of PHMFMs. Based on the Laplace transform domain transient description we calculate the time domain results using the concentrated matrix exponential (CME) based numerical inverse Laplace transform method (Horváth et al. 2019). We note that different transient analysis approaches are also considered in the literature. Chen et al. proposes a time domain numerical differential equation solution approach based on the finite difference method in (Chen et al. 2002), which discretizes the continuous fluid axis and its accuracy depends on discretization step. In Akar et al. (2020), Akar et al. adopts the approximate analysis approach of (Houdt and Blondia 2005) for computing the transient behaviour of PHMFMs based on the stationary analysis of an appropriately extended fluid model. While (Akar et al. 2020) builds on a stochastic interpretation based approximation approach, the current work is based on a exact Laplace transform domain analytical description.
The rest of the paper is organized as follows. Section 2 summarizes the basics of MFMs and provides the performance measures which are used later on. Section 3 introduces the considered class of PHMFMs with finite buffer and presents its transient analysis in Laplace transform domain. Section 4 discusses the model variant of PHMFMs with infinite buffer. Implementation details are provided in Sect. 5 and numerical examples in Sect. 6. The paper is concluded in Sect. 7.

Markov fluid models
MFMs are hybrid stochastic models composed of a continuous stochastic process X (t) ∈ R + , commonly referred to as fluid level, and a discrete stochastic process J (t) ∈ S • , commonly referred to as the state of the modulating Markov chain. Let us consider the MFM {X (t), J (t), t ≥ 0} defined by the generator of its background CTMC Q and the diagonal matrix of the fluid rates R. The subset of states with positive, negative and zero fluid rates are denoted by S + , S − and S 0 , while the set of states and the subset of states with non-zero fluid rates is denoted by S • and S • , respectively. That is S • = S + ∪ S − and S • = S • ∪ S 0 . The cardinality of the set S a , a ∈ {•, •, +, −, 0}, is denoted by |S a |. To order the states according to the sign of the fluid rates we introduce the permutation matrix Z with the following properties where subset specific matrix blocks of Q and R are with both R + and R − containing only positive diagonal elements. Matrix Z contains a single non-zero element in each row and each column which equals to one. That is, throughout the paper the underlined matrices refer to the original state ordering and the matrices without underline refer to the fluid rate specific ordering of the states. We are interested in the transient density and the transient boundary probability defined byṼ The corresponding matrices and Laplace transforms areṼ(t, e −stP (t, x, y)dt.

Characterizing matrices of infinite buffer MFMs
For i ∈ S + and j ∈ S − , the state dependent measure of returning to level zero is defined as where γ (x) is the first time when fluid process reaches level x, i.e., γ (x) = min{t : where I {event} is the indicator of event that equals to one when the event is true and otherwise it equals to zero. The matrix of size S + × S − , composed by the Ψ (s) i, j elements is Ψ (s).
In the sequel we define the measures of interest directly in transform domain, as in (8), and avoid the time domain definition as in (7). The return measure Ψ (s) satisfies the non-symmetric algebraic Riccati equation (NARE) with There are efficient numerical solution methods to compute Ψ (s) (Bean et al. 2018).
The matrices characterizing the fluid increase and fluid decrease process can be obtained from Ψ (s) as follows (Ahn and Ramaswami 2006) The spatial inverse of the fluid process is obtained by reverting the sign of the fluid rate in all states. The associated characterizing matrices are Q and −R. The characterizing matrices of the spatial inverse process, that is the matrices computed from the Q and −R matrices are denoted byΨ (s),K(s),Ĥ(s).

Characterizing matrices of the fluid process between boundaries
To characterize the evolution of the fluid process between two boundary fluid levels, 0 and b, we define the following state dependent hitting measures The special cases when the initial fluid level is one of the boundaries are According to da Silva Soares and Latouche (2006); Bean et al. (2009), the matrix composed by these hitting measures can be computed as and Similarly, the fluid density between two boundary levels, 0 and b, before reaching any of the boundaries are defined as and their matrix Laplace transforms are F The main advantage of (17) and (18) is that they makes the matrix exponential dependence on the fluid level explicit. As one of the consequences, the integral of F can also be computed explicitly as follow

Finite buffer piecewise homogeneous Markov fluid models
There are many variants of PHMFM considered in the literature differing mainly in the behaviour of the stochastic process at the boundaries. Here we consider the PHMFM variant from (Kankaya and Akar 2008), whose infinite buffer version is discussed in Sect. 4. The size of the buffer of the MFM is B and it is composed of K regions with region boundaries where the generator and the fluid rate matrices of the MFM are Q (k) and R (k) . At region borders, we assume that the generators and the fluid rate matrices areQ 0 denote the set of states where the fluid level increases, decreases and remains constant inside region k, that is, state i ∈ S (k) To separate the states where the fluid level is changing from the ones where it is constant, we introduce S (k) → denote the set of states where the fluid level increases, decreases and remains constant at boundary T k , that is, Condition 1 To avoid undefined stochastic behaviour we impose the following natural requirements on these sets These requirements mean that the fluid level cannot decrease below 0 and cannot increase above B and the fluid level can increase above (decrease below) boundary T k only if the fluid rate is positive above (negative below) the boundary.
We introduce the region and boundary specific permutation matrices with the following properties where underlined quantities refer to the original state ordering and quantities without underline refer to the region or boundary specific ordering of the states. Based on (20)-(21), we assume the availability of the required matrices in any convenient state ordering. The x, y) = 0 and is neglected in the sequel in any other cases. Due to potential rate change at boundary, V i j (s, x, y) might be discontinuous at T k for k ∈ {1, . . . , K − 1}. We apply the following density definition at the region borders Whenever the fluid level reaches boundary T k from below it migh be in a state of S (k) or in S (k) → . In the former case the fluid level keeps increasing at T k , while in the latter case it remains T k for a positive amount of time. Similar behaviours apply when the fluid level reaches boundary T k from above. To separate these two behaviours we refine the subset classification as follows We also introduce the subset specific filtering and reordering matrices a contains at most a single nonzero element in each row and column which equals to one. To keep the description simple we are going to present subset indexes for matrix blocks in the sequel. Behind these notations we assume the appropriate use of the related filtering and reordering matrices. We exemplify the use of the filtering and reordering matrices and their implementation in Sect. 5.
The following subsections describe the main steps of the proposed transient analysis approach in the order of their execution in the implemented algorithm. (14), (15), as well as their spatial inversesΨ (k)

Characteristic matrices of region k
The necessity of computing (s, y) and U (k) (s, x) depends on the initial and final fluid levels of the transient analysis. If needed, they are computed from (18), (19) and (17) (s) and so on. In some special cases (c.f. Sect. 3.6), we need to compute U(s, x), F(s, x) andF(s, x) for region k such that b is different from T k − T k−1 . In these cases, we explicitly indicate the interval size, i.e., U (k,z) (s, x) is obtained from (17) by assuming b = z, Ψ (s) = (k) (s) and so on.

Return measures of boundary k
we define the upward return measure of boundary T k as

External boundaries
For the following boundaries, the return measures are computed based on Bean and O'Reilly The expression to computeŶ (1) (s) contains matrices associated with the boundary T 0 (e.g. Q (0) →→ ) and also matrices associated with the (T 0 , T 1 ) region (e.g. U (1) −− (s)). Since the sign of the fluid rate might be different at T 0 and in the (T 0 , T 1 ) region, the dimensional validity of the expression raises notational problems. The notations we used in the expressions are sloppy. We applied a notation which would like to emphasize, on the one hand, that the matrices are associated with different boundary/region, on the other hand, that the expressions are dimensionally valid. We found this notational approach to be reasonable compact and expressive on the one hand and reasonable accurate on the other hand. The feasibility of the matrix operations are always ensured by Condition 1, that is + in this case. The same notational solution appears also in the sequel.

Internal boundaries
The return measures of the remaining boundaries are also discussed in  using an explicit approach which results in analytically complex expressions to evaluate. Here we adopt an implicit description which results in an easier to describe and implement analysis approach. This approach is based on the definition of appropriate auxiliary matrices. Using this approach the analytical complexity of the required measures is hidden by the matrix inverse operation. This implicit approach is applied also in the consecutive steps of the analysis.
For computing the return measures at the remaining boundaries we define the transition measure to move one boundary down and up. For i, j ∈ S • , are provided by the following theorems. where and, for a ∈ { , , →} and k ∈ {0, . . . , K } Proof Based on the intuitive stochastic meaning of introduced matrices we have The matrix form of (26)-(28) is, whose solution is (25).
If S (k) → = ∅ for an internal boundary, then the third block row and column vanishes in the expressions above. Since S (0) = ∅ and S (K ) = ∅ the second block row and column vanish for k = 0, and the third block row and column vanish for k = K . Based on the boundary transition measures, B (k) (s) andB (k) (s), the return measures of the internal boundaries can be computed as follows.

Starting and ending at boundaries
We start the Laplace transform domain description of the transient behaviour with the case when the initial and the final fluid levels are boundary values.

Starting and ending at the same boundary
First we consider the case when the initial and final levels are the same boundary level.

Theorem 3 Using the S (k) , S (k) , S (k)
→ division of the states, for boundary T k , k ∈ {0, . . . , K }, we have and with Proof Using the return measures of boundary T k and the Markov generator characterizing the evolution at boundary T k , for the blocks of P •→ (s, T k , T k ) we have The solution of (36) is (34). Similar block-wise equations apply for V(s, T k , T k ).
The first term in V (s, T k , T k ) represents the fact that starting from level T k results in a unit impulse in the fluid density at level T k at t = 0. A similar term appears in V (s, T k , T k ). The solution of (37) is (35).

Starting and ending at different boundaries
For k ∈ {1, . . . , K }, and k ≤ we have and for k ∈ {1, . . . , K }, and k > we have

Starting from boundary and ending between boundaries
The first term of (42) represents the cases when the last boundary visited before reaching y is T −1 and the second term represents the cases when it is T . To compute the CDF of the fluid level the integral of V •• (s, T k , y) is needed. It can be obtained using C ( ) (s, y) and C ( ) (s, y) as follows

Starting between boundaries and ending at boundary
and

Starting and ending between boundaries
For T k−1 < x < T k and T −1 < y < T and k = we condition on the first visited boundary and write where the integrals of the right hand side are provided in (43). For T k−1 < x, y < T k we have three cases: x = y, x < y and x > y. For x = y, we write for x < y, we condition on the last occasion when the process visits the border of the (x, T k ) interval before reaching level y and for x > y, we condition on the last occasion when the process visits the border of the (T k−1 , x) interval before reaching level y The integral of (44) form x to T k and the integral of (45) form T k−1 to x can also be computed as in (43).

Starting from and going to S 0
For T k−1 < x < T k and T −1 < y < T and for T k−1 < x < T k , T −1 < y < T and x = y The case of T k−1 < x = y < T k and J (0) ∈ S (k) 0 , results in a probability mass between boundaries which we avoid considering directly. Our analysis approach can handle this case by introducing an additional boundary at x and using Theorem 3.

Infinite buffer case
The case when the buffer is infinite can be defined and analysed as follows. The buffer is composed of K regions with region boundaries where the generator and the fluid rate matrices of the MFM are Q (k) and R (k) . Since we are interested in transient analysis, the last region, characterized by Q (K ) and R (K ) , might have positive drift as well. At region borders, we assume that the generators and the fluid rate matrices areQ With this model definition we only need to modify the analysis of the last region compared to the finite buffer case. The rest of the section collects the elements of the analysis, which needs to be modified when the buffer is infinite.

Condition 2
For having a consistent model the following requirements need to be satisfied.
Since we assume finite initial and final fluid levels, in the infinite buffer case, we neglect all quantities of Sect. 3.3 which starts or ends at the upper bound of the buffer, T K . The analysis of region K , in Sects. 3.4, 3.5 and 3.6 simplifies as follows. For k ∈ {0, . . . , K − 1} and T K −1 < y, and its integral is Finally, the cases when both, the initial and the final levels are above T K −1 can be handled are follows. For T K −1 < x = y,

Implementation notes
The above explained numerical procedure for the transient analysis of PHMFMs is implemented in Mathematica and available at webspn.hit.bme.hu/˜telek/tools.htm.
In this section, we summarize some interesting features of the implementation.

Model description
To simply the model definition it is also possible to generate a model description based only on the region descriptions, Q (k) and R (k) for k ∈ {1, . . . , K }, and defining if the upper or the lower regions dominate at boundaries. If the upper regions dominate at boundaries (as in da Silva Soares and Latouche (2009)) theñ The case when the lower regions dominate at boundaries can be defined as its spatial inverse. Our code checks the input data according to Conditions 1 in case of finite buffer or Condition 2 in case of finite buffer. The region based boundary definition characterized by (51) satisfies Conditions 1 and 2.

Index list based matrix operations
Our implementation makes use of the index list based matrix definition available in many programming environments (including Matlab and Mathematica). During the model descrip-tion phase the program calculates the S (k) a sets for k ∈ {1, . . . , K } and a ∈ {+, −, 0, •, •, , , →, + , → + , − , − →}. These sets contains the associated state numbers according to the original state numbering. With the help of these lists the required matrix blocks can be easily obtained. E.g., Q (k) where "[" and "]" denote the indexing operator.

Computational complexity
Each of the matrices introduces in Sect. 3.1 are at most of size |S • | × |S • |, and the computational complexity of obtaining one of such matrices is O(|S • | 3 ). In Sects. 3.1-3.2, O(K ) of such matrices are computed, while in Sect. 3.3, O(K 2 ) matrices are computed. The remaining computations, in Sects. 3.4-3.6, have lower order computational complexity, hence the overall computational complexity of the procedure is O(K 2 |S • | 3 ).
The memory complexity of the procedure is dominated by storing the O(K 2 ) P(s, T k , T ) matrices, which results in a memory complexity of O(K 2 |S • | 2 ).

Numerical experiment
As the homogeneous finite and infinite MFMs are special cases of PHMFMs, we can verify our numerical procedure against some of the transient results available in the literature. Additionally, in case of finite and stable infinite PHMFMs we can compare the t → ∞ limiting behaviour of the transient analysis with the stationary results available in the literature.
The set of transient results we evaluated with our implementation includes Example 1 of Ahn et al. (2007), which was also considered in Akar and Sohraby (2004). The results reported in (Ahn et al. 2007), Table 3 (based on an unspecified order numerical inverse Laplace transformation method) are identical to our results (which are based on the CME numerical inverse Laplace transformation method of Horváth et al. (2019) with N = 21) up to 4 digits.
A multi regime PHMFM example was introduced in Mandjes et al. (2003), which models the behaviour of a network access protocol with finite buffer and congestion control. The same example was considered in . The stationary measures of this example (e.g., the probability of full buffer) are evaluated in Mandjes et al. (2003) and some transient measures (e.g., the return time distribution of the internal boundary) are provided in . This example is investigated in the rest of the section.
The state space of the model is defined as S • = {0, 1, . . . , 10}, and the buffer is of size 5 with fluid boundaries T 0 = 0, T 1 = 2, T 2 = 5. The characterising matrices of region 1 and 2 are The boundary behaviour at T 0 , T 1 , T 2 is characterized bỹ i j , 0).
These characterizing matrices result in the following set definitions for region 1 and 2, and for boundaries T 0 , T 1 , T 2 : These set definitions satisfy Conditions 1. Figures 1 and 2 depict the CDF of the fluid level distribution at time t = 0.1, 1, 10 starting from state 0 (with negative fluid rate) and state 10 (with positive fluid rate) and from level 0 (empty buffer) and level 5 (full buffer). The figures indicate that at time t = 10 the fluid model reaches its steady state and the fluid level distribution starting from the different initial conditions converge to the same stationary distribution.
The results for t = 10 also indicate that the t → ∞ limit of the probability of full buffer computed by our tool, converge to the limit reported in Mandjes et al. (2003) (2.01 · 10 −4 ) in all evaluated cases if the initial states and fluid levels.
While the primary performance measure of the current paper is the transient behaviour according to (5) and (6), the performance measures introduced in course of the computations include the distribution of the return times of the boundaries, according to (23)

Conclusions
The paper presents an analysis approach for Laplace transform description of the transient behaviour of PHMFMs with arbitrary fluid rate changes at the barriers. Due to the large number of possible boundary behaviours the proposed approach replaces the previously applied direct descriptions with an equation system based indirect one. The obtained relatively simple analytical description keeps also the implementation of the method at a feasible level of complexity. The explicit matrix exponential description of the transient fluid density based on (18) allows a closed form computation of the cumulated density as well as the moments of the fluid distribution. The earlier is discussed in the paper, the later is neglected, but follows a similar pattern.
Our Mathematica implementation, containing also the model descriptions of all presented examples, is available at webspn.hit.bme.hu/˜telek/tools.htm. We also compared this implementation (in Mathematica) with the publicly available implementation of the method presented in Akar et al. (2020) (in Matlab). In all evaluated cases, we obtained practically identical results.
Funding Open access funding provided by Budapest University of Technology and Economics.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.