Perturbation analysis of Markov modulated fluid models

TL;DR

This paper investigates how small changes in the parameters of Markov modulated fluid models affect their behavior, especially focusing on the matrix of first return probabilities, by developing a new analytical approach.

Contribution

It introduces a novel method to analyze the impact of perturbations on the first return probability matrix in Markov modulated fluid models.

Findings

Developed a substitute for the first return probability matrix.

Analyzed the effect of perturbations on the matrix.

Provided insights into the stability of fluid models under perturbations.

Abstract

We consider perturbations of positive recurrent Markov modulated fluid models. In addition to the infinitesimal generator of the phases, we also perturb the rate matrix, and analyze the effect of those perturbations on the matrix of first return probabilities to the initial level. Our main contribution is the construction of a substitute for the matrix of first return probabilities, which enables us to analyze the effect of the perturbation under consideration.

Tables

Table 1. Table 1: Values of E + ​ ( ε ) subscript 𝐸 𝜀 E_{+}(\varepsilon) and E ⊕ ​ ( ε ) subscript 𝐸 direct-sum 𝜀 E_{\oplus}(\varepsilon) in Cases 1.a to 3.a, for ε 𝜀 \varepsilon equal to 10 − 4 superscript 10 4 10^{-4} and  10 − 2 superscript 10 2 10^{-2} .

	$\epsilon ={10}^{-4}$		$\epsilon ={10}^{-2}$
Case	$E_{+}$	$E_{\oplus }$	$E_{+}$	$E_{\oplus }$
1.a	5.37 ${10}^{-7}$	3.39 ${10}^{-6}$	4.60 ${10}^{-3}$	2.94 ${10}^{-3}$
2.a	1.92 ${10}^{-12}$	2.00 ${10}^{-12}$	2.08 ${10}^{-8}$	2.15 ${10}^{-8}$
3.a	3.77 ${10}^{-8}$	4.80 ${10}^{-8}$	3.66 ${10}^{-4}$	4.67 ${10}^{-4}$

Table 2. Table 2: Values of E ∞ ​ ( ε ) subscript 𝐸 𝜀 E_{\infty}(\varepsilon) in Cases 1.b to 3.b, for ε 𝜀 \varepsilon equal to 10 − 4 superscript 10 4 10^{-4} and 10 − 2 superscript 10 2 10^{-2} .

Case	$\epsilon ={10}^{-4}$	$\epsilon ={10}^{-2}$
1.b	1.08 ${10}^{-7}$	1.05 ${10}^{-3}$
2.b	5.11 ${10}^{-8}$	4.91 ${10}^{-4}$
3.b	1.33 ${10}^{-6}$	1.15 ${10}^{-2}$

Equations

$$X(t)$$

$$\text{where}Y(t)$$

$$A=\left[\begin{array}[]{ccc}A_{++}&A_{+0}&A_{+-}\\ A_{0+}&A_{00}&A_{0-}\\ A_{-+}&A_{-0}&A_{--}\end{array}\right],$$

$$C=\left[\begin{array}[]{ccc}C_{+}&&\\ &C_{0}&\\ &&C_{-}\end{array}\right]$$

$${\xi }_i=\underset{t\to \infty }{\lim }\mathbb{P}\left[\phi (t)=i|\phi \left(0\right)=j\right],$$

$${\Psi }_{ij}=\mathbb{P}\left[{\tau }_{-}<\infty ,\phi \left({\tau }_{-}\right)=j|Y\left(0\right)=0,\phi \left(0\right)=i\right],$$

$$C_{+}^{-1}{Q}_{+-}+C_{+}^{-1}{Q}_{++}X+X\left|C_{-}^{-1}\right|Q_{--}+X\left|C_{-}^{-1}\right|Q_{-+}X=0,$$

$$\displaystyle\left[\begin{array}[]{cc}Q_{++}&Q_{+-}\\ Q_{-+}&Q_{--}\end{array}\right]=$$

$${\hat{\Psi }}_{ij}=\mathbb{P}\left[{\tau }_{+}<\infty ,\phi \left({\tau }_{+}\right)=j|Y\left(0\right)=0,\phi \left(0\right)=i\right],$$

$$U$$

$$K$$

$$\Psi^{(1)}=\frac{\mathrm{d}\Psi(\varepsilon)}{\mathrm{d}\varepsilon}\bigg{|}_{\varepsilon=0}$$

$$A(\epsilon )=A+\epsilon \tilde{A},$$

$$\tilde{A}=\left[\begin{array}{ccc}{\tilde{A}}_{++}& {\tilde{A}}_{+0}& {\tilde{A}}_{+-}\\ {\tilde{A}}_{0+}& {\tilde{A}}_{00}& {\tilde{A}}_{0-}\\ {\tilde{A}}_{-+}& {\tilde{A}}_{-0}& {\tilde{A}}_{--}\end{array}\right],$$

$$C_{+}^{-1}{Q}_{+-}(\epsilon )+C_{+}^{-1}{Q}_{++}(\epsilon )X+X\left|C_{-}^{-1}\right|Q_{--}(\epsilon )+X\left|C_{-}^{-1}\right|Q_{-+}(\epsilon )X=0,$$

$$\displaystyle\left[\begin{array}[]{cc}Q_{++}(\varepsilon)&Q_{+-}(\varepsilon)\\ Q_{-+}(\varepsilon)&Q_{--}(\varepsilon)\end{array}\right]=\left[\begin{array}[]{cc}Q_{++}+\varepsilon\tilde{Q}_{++}&Q_{+-}+\varepsilon\tilde{Q}_{+-}\\ Q_{-+}+\varepsilon\tilde{Q}_{-+}&Q_{--}+\varepsilon\tilde{Q}_{--}\end{array}\right]+O(\varepsilon^{2}).$$

$$KX+XU=$$

$$\displaystyle F(\varepsilon,\mathcal{X})=\$$

$${{\partial }_{\mathcal{X}}F(\epsilon ,\mathcal{X})|}_{\epsilon =0,\mathcal{X}=\Psi }(Y)=H,$$

$$KY+YU=H.$$

$$C(\epsilon )=C+\epsilon \tilde{C}$$

$$\tilde{C}=\left[\begin{array}[]{ccc}\tilde{C}_{+}&&\\ &\tilde{C}_{0}&\\ &&\tilde{C}_{-}\end{array}\right]$$

$$C_{+}^{-1}(\epsilon )Q_{+-}+C_{+}^{-1}(\epsilon )Q_{++}X+X\left|C_{-}^{-1}(\epsilon )\right|Q_{--}+X\left|C_{-}^{-1}(\epsilon )\right|Q_{-+}X=0.$$

$$\Psi (\epsilon )=\Psi +\epsilon {\Psi }^{(1)}+O({\epsilon }^2),$$

$$KX+XU=-\Psi |C_{-}^{-1}|{\tilde{C}}_{-}U-C_{+}^{-1}{\tilde{C}}_{+}\Psi U,$$

$$A=\left[\begin{array}[]{c|c|c}A_{++}&A_{+\oplus}&A_{+-}\\ \hline\cr A_{\oplus+}&A_{\oplus\oplus}&A_{\oplus-}\\ \hline\cr A_{-+}&A_{-\oplus}&A_{--}\end{array}\right].$$

$$A=\left[\begin{array}[]{cc|c}A_{++}&A_{+\oplus}&A_{+-}\\ A_{\oplus+}&A_{\oplus\oplus}&A_{\oplus-}\\ \hline\cr A_{-+}&A_{-\oplus}&A_{--}\end{array}\right]$$

$$\boldsymbol{\Psi}(\varepsilon)=\left[\begin{array}[]{c}{\Psi}_{+-}(\varepsilon)\\ {\Psi}_{\oplus-}(\varepsilon)\end{array}\right].$$

$$\displaystyle\left[\begin{array}[]{cc}C_{+}+\varepsilon\tilde{C}_{+}\\ &\varepsilon\tilde{C}_{\oplus}\end{array}\right]^{-1}\left(\left[\begin{array}[]{c}A_{+-}\\ A_{\oplus-}\end{array}\right]+\left[\begin{array}[]{cc}A_{++}&A_{+\oplus}\\ A_{\oplus+}&A_{\oplus\oplus}\end{array}\right]X\right)$$

$$\displaystyle+X\left|C_{-}+\varepsilon\tilde{C}_{-}\right|^{-1}\left(A_{--}+\left[\begin{array}[]{cc}A_{-+}&A_{-\oplus}\end{array}\right]X\right)=0.$$

$$\mathbf{\Psi }(\epsilon )=\overline{\Psi }+\epsilon {\Psi }^{(1)}+O({\epsilon }^2),$$

$$\overline{\Psi }$$

Full text

Keywords: Markov modulated fluid models; Perturbation analysis; First return probabilities.

1 Introduction

Most mathematical models have input parameters that are typically estimated from the real world data. Since the parameters in the modeled system represent quantities that can suffer from small errors, it is natural to analyze how the performance measures are affected by small changes in the parameters. Using the structural properties of the model, it becomes possible to assess the impact of perturbations on the key matrices of the underlying process by providing computationally feasible solutions along with probabilistic interpretation.

Markov modulated fluid models appeared in the 1960s to study the continuous-time behavior of queues and dams, an early paper being Loynes [11]. In the eighties, Markovian fluid models started to be more extensively investigated, in particular their stationary density, see for instance Rogers [14] and Asmussen [2]. The importance of the matrix of first return probabilities has been demonstrated in Ramaswami [13] and its computation has attracted much attention, see Bean et al. [3] and Bini et al. [4]. One may derive from $\Psi$, the matrix of first return probabilities form above, important performance measures of the model, such as the stationary density of the level of the fluid model.

The model $\{(X(t),\varphi(t)):t\in\mathbb{R}^{+}\}$ is described as follows: $\varphi(t)$ is a Markov chain, with finite state space $\mathcal{S}$, it is called the phase process; $X(t)$ is a continuous function, called the level. The evolution of the level is continuous and may be expressed as

[TABLE]

so that it varies linearly with rate $c_{i}$ when $\varphi(t)=i$, $i\in\mathcal{S}$. We partition $\mathcal{S}$ into $\mathcal{S}_{+}\cup\mathcal{S}_{0}\cup\mathcal{S}_{-}$ with $\mathcal{S}_{+}=\{i\in\mathcal{S}:c_{i}>0\}$, $\mathcal{S}_{0}=\{i\in\mathcal{S}:c_{i}=0\}$ and $\mathcal{S}_{-}=\{i\in\mathcal{S}:c_{i}<0\}$. The infinitesimal generator of the phase process is denoted by $A$ and is written, possibly after permutation of rows and columns, as

[TABLE]

and the rate matrix is denoted by

[TABLE]

with $C_{+}=\mbox{diag}(c_{i}:i\in\mathcal{S}_{+})$, $C_{-}=\mbox{diag}(c_{i}:i\in\mathcal{S}_{-})$ and $C_{0}$ is a null matrix. Throughout the paper, we make the following assumption.

Assumption 1.1

The Markov modulated fluid model is positive recurrent, that is, $\boldsymbol{{\boldsymbol{\xi}}}C{\boldsymbol{1}}<0$, where ${\boldsymbol{\xi}}$ is the stationary probability vector defined for $i,j\in\mathcal{S}$ by

[TABLE]

and is the unique solution of the equation $\boldsymbol{{\boldsymbol{\xi}}}A=0$ such that $\boldsymbol{{\boldsymbol{\xi}}1}=\boldsymbol{1}$, where $\boldsymbol{1}$ denotes the column vector of 1’s.

A key matrix for Markov modulated fluid models is the matrix $\Psi$ of first return probabilities to the initial level from above, with dimensions $|\mathcal{S}_{+}|\times|\mathcal{S}_{-}|$, and components

[TABLE]

where $\tau_{-}=\inf\{t>0:Y(t)<0\}$, $i\in\mathcal{S}_{+}$ and $j\in\mathcal{S}_{-}$. By Rogers [14, Theorem 1], $\Psi$ is the minimal nonnegative solution of the Riccati equation

[TABLE]

where $|C_{-}^{-1}|$ denotes the entrywise absolute value of $C_{-}^{-1}$ and

[TABLE]

Similarly, the matrix $\hat{\Psi}$ of first return probabilities to the initial level from below has components

[TABLE]

where $\tau_{+}=\inf\{t>0:Y(t)>0\}$, $i\in\mathcal{S}_{-}$ and $j\in\mathcal{S}_{+}$, it satisfies a Riccati equation similar to (6). The present article focuses on the perturbation analysis of $\Psi$ only, as the analysis for $\hat{\Psi}$ is similar.

Two other important matrices are

[TABLE]

The matrix $U$ is the infinitesimal generator of the process of downward record and is such that for $i,j\in\mathcal{S}_{-}$, $(e^{Ux})_{ij}$ is the probability that, starting from $(y,i)$, for any $y$, the process reaches level $y-x$ in finite time and that $(y-x,j)$ is the first state visited in level $y-x$. The matrix $K$ defined in (16) is also an important matrix for Markov modulated fluid models and appears in the sationary density of the fluid model, see Section 4.

For a long time there has been a recurrent interest in perturbation analysis, see for instance Cao and Chen [5], Heidergott, et al. [7], Antunes et al. [1]. In this paper, we analyze the perturbation of Markov modulated fluid models. When the infinitesimal generator (2) of the phases is perturbed into $A(\varepsilon)=A+\varepsilon\tilde{A}$, the analysis follows the usual path: the perturbed first return probability matrix $\Psi(\varepsilon)$ is shown to be analytic, and computable equations are readily obtained for the derivatives of $\Psi(\varepsilon)$. We focus on the first order derivative

[TABLE]

of a perturbed Markov modulated fluid model as it provides a good approximation of the effect of the perturbation on the system when compared to the unperturbed system. Furtermore, we are interested in the structures and going beyond the first derivative is rather computational and does not bring much more information.

We also analyze the effect on $\Psi$ of perturbations of the rate matrix (3). When $C$ is perturbed as $C(\varepsilon)=C+\varepsilon\tilde{C}$, phases of $\mathcal{S}_{0}$ may be transformed into phases of $\mathcal{S}_{+}$ or $\mathcal{S}_{-}$ in the perturbed model, with the consequence that a perturbation of the rates $c_{i}$ appearing in (1) may modify the structure of $\Psi(\varepsilon)$ as the dimensions are not the same as those of $\Psi$. Clearly, the comparison between the matrices $\Psi(\varepsilon)$ and $\Psi$ requires more care.

We do not consider cases where both the generator $A$ and the rate matrix $C$ are perturbed, as our results show that this may be done, at the cost of increased complexity in the expressions obtained.

In Section 2, we analyze perturbations of the infinitesimal generator of the phases. In Section 3, we analyze perturbations on the rate matrix $C$ in four different cases. In Section 3.1 we assume that the phases of $\mathcal{S}_{0}$ are unaffected by the perturbation. In Sections 3.2–3.4 we examine what happens when the phases of $\mathcal{S}_{0}$ are affected by the perturbation. We propose an adapted version of $\Psi$ which enables the analysis of the effect of the perturbation under consideration. We decompose the analysis in three subsections for the sake of clarity: firstly, we assume that all the phases in $\mathcal{S}_{0}$ become phases of $\mathcal{S}_{+}$ after perturbation, next, we assume that they all become phases of $\mathcal{S}_{-}$ after perturbation, finally, we assume that the phases in $\mathcal{S}_{0}$ are split between $\mathcal{S}_{+}$ and $\mathcal{S}_{-}$. The general approach is the same in the three cases but the details differ and become much more involved in the last. As an application, we derive in Section 4 the first order approximation of the stationary density of a perturbed fluid model. In Section 5, we provide a numerical illustration.

2 Perturbation of the infinitesimal generator

In this section, the infinitesimal generator $A$ is perturbed and becomes

[TABLE]

where

[TABLE]

$\widetilde{A}{\boldsymbol{1}}=0$, and we assume that $A(\varepsilon)$ is an irreducible infinitesimal generator for $\varepsilon$ sufficiently small in a neighborhood of 0.

The matrix $\Psi(\varepsilon)$ of first return probabilities for the perturbed model is the minimal nonnegative solution of the Riccati equation

[TABLE]

where $Q(\varepsilon)$ is defined by (14), with $A(\varepsilon)$ replacing $A$. We write

[TABLE]

Theorem 2.1

The matrix $\Psi(\varepsilon)$ of first return probabilities, minimal nonnegative solution to (19), for the perturbed model is analytic in a neighbourhood of zero. Furthermore, $\Psi^{(1)}$ is the unique solution of the Sylvester equation

[TABLE]

where $K$ and $U$ are defined in (16) and (15).

**Proof ** Define the continuous operator

[TABLE]

We have $F(0,\Psi)=0$ and ${\partial}_{\mathcal{X}}F(\varepsilon,\mathcal{X})$ exists in a neighborhood of $(0,\Psi)$ and is continuous at $(0,\Psi)$. For $Y,H\in\mathbb{R}^{|\mathcal{S}_{+}|\times|\mathcal{S}_{-}|}$, the equation

[TABLE]

is equivalent to the Sylvester equation

[TABLE]

From Rogers [14] and Govorun et al. [6], we have $\operatorname*{sp}(K)\in\{z\in\mathbb{C}:\mathrm{Re}({z})<0\}$ and $\operatorname*{sp}(-U)\in\{z\in\mathbb{C}:\mathrm{Re}({z})\geq 0\}$. Thus, $K$ and $-U$ have no common eigenvalue and, by Lancaser and Tismenetsky [10, page 414], (21) has a unique solution, so that $\left.{\partial}_{\mathcal{X}}F(\varepsilon,\mathcal{X})\right|_{\varepsilon=0,\mathcal{X}=\Psi(0)}$ is a nonsingular operator. We conclude that $\Psi(\varepsilon)$ is analytic at zero by the Implicit Function Theorem. $\square$

Remark 2.2

It immediately results from Xue et al. [15, Theorem 2.2] that small relative changes to the entries of $Q$ induce small relative differences between $\Psi$ and $\Psi(\varepsilon)$. The bounding coefficient matrix in [15, Eqn. (2.12)] is the solution of a Sylvester equation with the same coefficients $K$ and $U$ as in (20) and a different right-hand side.

3 Perturbation of the rate matrix

Define

[TABLE]

with $\tilde{C}=\operatorname*{diag}(\tilde{c}_{i}:i\in\mathcal{S})$, partitioned as

[TABLE]

where the orders of $\tilde{C}_{+}$, $\tilde{C}_{0}$ and $\tilde{C}_{-}$ are equal to those of $C_{+}$, $C_{0}$ and $C_{-}$, respectively. Assume that $\varepsilon$ is small enough so that the diagonal elements of $C_{+}(\varepsilon)$ are strictly positive and those of $C_{-}(\varepsilon)$ strictly negative.

We analyze separately the cases $\tilde{C_{0}}=0$ (in Section 3.1) and $\tilde{C_{0}}\neq 0$. If $\tilde{C_{0}}\neq 0$, the perturbation has the effect of changing null phases into non-null phases. To simplify the presentation, we suppose at first that all phases of $\mathcal{S}_{0}$ become phases of the same non-null subset $\mathcal{S}_{+}$ after perturbation. This is analyzed in Section 3.2. In Section 3.3, we treat the case where all the phases of $\mathcal{S}_{0}$ become phases of $\mathcal{S}_{-}$ after perturbation. Finally, we assume in Section 3.4 that the phases in $\mathcal{S}_{0}$ are split partially into $\mathcal{S}_{+}$ and into $\mathcal{S}_{-}$.

Clearly, Section 3.4 covers the cases analyzed in Sections 3.2 and 3.3. It is useful, nevertheless, to proceed through the special cases first, for which the results are easier to follow. In various remarks, we emphasize the unity of treatment.

The Implicit Function Theorem applies in all cases to prove the analyticity of $\Psi(\varepsilon)$, although details become more involved as we proceed from the simplest to the most general case. We show this in Theorem 3.2 and Theorem 3.5 and we omit the details for Theorem 3.7.

3.1 Phases in $\mathcal{S}_{0}$ unaffected

Assume that $\tilde{C}_{0}=0$ so that $C_{0}(\varepsilon)=0$ as well. The matrix $\Psi(\varepsilon)$ of first return probabilities for the perturbed model is the minimal nonnegative solution of the Riccati equation

[TABLE]

The next Theorem is proved by applying to (24) the same argument as in Theorem 2.1.

Theorem 3.1

Assume $C(\varepsilon)=C+\varepsilon\tilde{C}$, with $\tilde{C}_{0}=0$. The matrix $\Psi(\varepsilon)$ of first return probabilities for the perturbed model is analytic at zero and may be written as

[TABLE]

where $\Psi$ is the minimal non-negative solution to (6) and $\Psi^{(1)}$ is the unique solution of the Sylvester equation

[TABLE]

where $K$ and $U$ are defined in (16) and (15) respectively. $\square$

3.2 Migration of $\mathcal{S}_{0}$ to $\mathcal{S}_{+}$

Assume that $\tilde{c}_{i}>0$ for all $i$ in $\mathcal{S}_{0}$, this means that all phases of $\mathcal{S}_{0}$ become phases of fluid increase after perturbation. To make this explicit in our equations, we replace the subscript 0 by the subscript $\oplus$ and write $\mathcal{S}_{\oplus}$ instead of $\mathcal{S}_{0}$, etc. The infinitesimal generator of the phase process is written as

[TABLE]

After perturbation, it is partitioned as

[TABLE]

and the set of phases with positive rates is $\mathcal{S}_{+}\cup\mathcal{S}_{\oplus}$. The dimensions of the first return probability matrix become $(|\mathcal{S}_{+}|+|\mathcal{S}_{\oplus}|)\times|\mathcal{S}_{-}|$ after perturbation and $\Psi$ may not be directly compared to $\Psi(\varepsilon)$, the matrix of first return probabilities of the perturbed model, which is partitioned as

[TABLE]

The matrix ${\boldsymbol{\Psi}}(\varepsilon)$ is the minimal nonnegative solution of the Riccati equation

[TABLE]

As we show in the next theorem, comparisons are nevertheless possible, as $\Psi$ is immediately recognised in the limit $\overline{\Psi}=\lim_{\varepsilon\rightarrow 0}\boldsymbol{\Psi}(\varepsilon)$.

Theorem 3.2

The matrix (28) of first return probabilities for the perturbed model, minimal nonnegative solution of (37), is analytic near zero and may be written as

[TABLE]

where

[TABLE]

where $\Psi$ is given in (6), $\Psi_{\oplus-}=(-A_{\oplus\oplus}^{-1})(A_{\oplus-}+A_{\oplus+}\Psi)$, ${\Psi}_{+-}^{(1)}$ is the unique solution of the Sylvester equation

[TABLE]

and

[TABLE]

The matrices $K$ and $U$ are defined in (16) and (15), and

[TABLE]

with $K_{+\oplus}={C_{+}^{-1}A_{+\oplus}}+\Psi|C_{-}^{-1}|A_{-\oplus}$.

**Proof ** To remove the effect of $\varepsilon^{-1}$ in the left-most coefficient of (37), we pre-multiply both sides by $\operatorname*{diag}(I,\varepsilon I)$. For \mathcal{X}=\left[\begin{array}[]{c}\mathcal{X}_{+-}\\ \mathcal{X}_{\oplus-}\end{array}\right], we define the operator

[TABLE]

The equation

[TABLE]

is equivalent to the set of equations

[TABLE]

This is a non-singular system, so that $\boldsymbol{\Psi}(\varepsilon)$ is analytic, by the Implicit Function Theorem. From (37), we obtain the two equations:

[TABLE]

in which we take the limit for $\varepsilon\rightarrow 0$. The second equation gives

[TABLE]

and the first equation gives ${\Psi}_{+-}(0)$ as the solution of (6), so that ${\Psi}_{+-}(0)=\Psi.$ This proves (42).

Taking the coefficients of $\varepsilon$ in (46) and using (47) leads directly to (44). To prove (43), we note that $\lim_{\varepsilon\rightarrow 0}U(\varepsilon)=U$ so that, taking in (45) the limit for $\varepsilon\rightarrow 0$ and using (42), we obtain

[TABLE]

We take the coefficient of $\varepsilon$ in (45) and we use (48) to obtain

[TABLE]

with $K_{++}={C_{+}^{-1}A_{++}}+\Psi|C_{-}^{-1}|A_{-+}$. Using (44) and (16) gives then (43). $\square$

Remark 3.3

The components of the block $\Psi$ in $\overline{\Psi}$ are those defined in (5), for which one has a clear interpretation. The components of the second block have a probabilistic interpretation as well: the $ij$th entry, for $i\in\mathcal{S}_{\oplus}$ and $j\in\mathcal{S}_{-}$, is the sum of

• •

$[(-A_{\oplus\oplus}^{-1})A_{\oplus-}]_{ij}$, the probability that the phase process eventually goes from phase $i$ to phase $j$, after some time spent in $S_{\oplus}$ and

• •

$[(-A_{\oplus\oplus}^{-1})A_{\oplus+}\Psi]_{ij}$, the probability that the phase process leaves $\mathcal{S}_{\oplus}$ for a phase in $\mathcal{S}_{+}$ and later returns to the initial level in phase $j$.

Remark 3.4

The Sylvester equations (25) and (43) for $\Psi^{(1)}_{+-}$ are nearly identical. The only difference is the last term in the right-hand side of (43), reflecting the migration of all phases of $\mathcal{S}_{0}$ to phases of fluid increase.

3.3 Migration of $\mathcal{S}_{0}$ to $\mathcal{S}_{-}$

Assume that $\tilde{c}_{i}<0$ for all $i$ in $\mathcal{S}_{0}$, so that all the phases of $\mathcal{S}_{0}$ become phases of $\mathcal{S}_{-}$ after perturbation. The set of such phases is written $\mathcal{S}_{\ominus}$ and the infinitesimal generator of the phases is written as

[TABLE]

The matrix of first return probabilities of the perturbed model is partitioned as

[TABLE]

and it is the minimal nonnegative solution of a Riccati equation which we rewrite as the two equations

[TABLE]

Theorem 3.5

The matrix $\boldsymbol{\Psi}(\varepsilon)$ of first return probabilities, minimal nonnegative solution to (49) and (50) is near zero and may be written as

[TABLE]

where

[TABLE]

The matrix $\Psi$ is given in (6), ${\Psi}_{+-}^{(1)}$ is the unique solution of the Sylvester equation

[TABLE]

and

[TABLE]

the matrices $K$ and $U$ are defined in (16) and (15) and

[TABLE]

**Proof ** Here, to remove the effect of $\varepsilon^{-1}$ as a coefficient of $|\tilde{C}_{\ominus}^{-1}|$ in (49) and (50), we define $\Gamma(\varepsilon)=\varepsilon^{-1}\Psi_{+\ominus}(\varepsilon)$. We define the operator, for \mathcal{X}=\left[\begin{array}[]{cc}\mathcal{X}_{+\ominus}&\mathcal{X}_{+-}\end{array}\right],

[TABLE]

One shows that \left[\begin{array}[]{cc}\Psi_{+\ominus}^{(1)}&\Psi\end{array}\right] is a solution of $F(\varepsilon,\mathcal{X})=0$, where $\Psi_{+\ominus}^{(1)}$ is defined in (57). Next, we take the derivative of $F$ with respect to $\mathcal{X}$, evaluated at $\varepsilon=0$, \mathcal{X}=\left[\begin{array}[]{cc}\Psi_{+\ominus}^{(1)}&\Psi\end{array}\right]. The system is equivalent to the set of equations

[TABLE]

where, by (14), (15), (16),

[TABLE]

The system is non-singular so that \left[\begin{array}[]{cc}\Gamma(\varepsilon)&\Psi_{+-}(\varepsilon)\end{array}\right] is analytic.

The block components of $\overline{\Psi}$ are obtained as follows. As $\varepsilon\Gamma(\varepsilon)=\Psi_{+\ominus}(\varepsilon)$, we find that $\Psi_{+\ominus}(0)=0$. Next, define

[TABLE]

which is finite since $\Gamma(\varepsilon)$ is analytic. We rewrite (49) and find that

[TABLE]

Taking the limit as $\varepsilon\rightarrow 0$ in (50) and replacing $W$ by (61) leads to (6). Thus, $\lim_{\varepsilon\rightarrow 0}{\Psi}_{+-}(\varepsilon)=\Psi,$ and (53) is proved.

The block components of $\Psi^{(1)}$ are obtained as follows. Taking the coefficients of $\varepsilon^{0}$ in (49) gives directly (57). To show (56), we take the coefficients of $\varepsilon^{2}$ in (50) and get the equation

[TABLE]

We equate the coefficients of $\varepsilon$ in (49) and get

[TABLE]

By the Riccati equation (6) and the definition (58) of $U$, we have

[TABLE]

We replace the first coefficient $\Psi_{+\ominus}^{(1)}$ in (63) by its expression (57), then we replace $\Psi_{+\ominus}^{(2)}|\tilde{C}_{\ominus}^{-1}|$ in (62) by the modified right-hand side of (63). We put together the coefficients of $\Psi_{+-}^{(1)}$, use (58), (59) and eventually obtain (56). $\square$

Remark 3.6

The physical justification of ${\Psi}_{+\ominus}(0)=0$ goes as follows: $({\Psi}_{+\ominus}(\varepsilon))_{ij}$ is the probability that the level moves to [math] in phase $j\in\mathcal{S}_{\ominus},$ given that the initial level is [math] and the phase is $i\in\mathcal{S}_{+}$, in the limit, when $\varepsilon$ approaches $0,$ this probability tends to [math] because the fluid can only return to level zero in a phase of $\mathcal{S}_{-}$.

3.4 General case

Assume $\tilde{c}_{i}\neq 0$ for $i$ in $\mathcal{S}_{0}$, so that all the phases of $\mathcal{S}_{0}$ disseminate in $\mathcal{S}_{+}$ and $\mathcal{S}_{-}$ after perturbation. The infinitesimal generator becomes

[TABLE]

We find here a superposition of the effects observed in the two special cases examined in Sections 3.2 and 3.3. The matrix of first return probabilities from above of the perturbed system takes the form

[TABLE]

it is the unique solution of the usual Riccati equation which may be rewritten as the following set of four equations:

[TABLE]

We show below that $\boldsymbol{\Psi}(\varepsilon)$ is analytic, thus we may write the matrices $U(\varepsilon)$ and $K(\varepsilon)$ as

[TABLE]

in particular, the blocks

[TABLE]

play an important role in what follows.

Theorem 3.7

The matrix $\boldsymbol{\Psi}(\varepsilon)$ of first return probabilities, minimal nonnegative solution to (66-69) for the perturbed model is near zero and may be written as

[TABLE]

where

[TABLE]

The block $\Psi$ is given in (6),

[TABLE]

$\Psi_{\oplus\ominus}$* is the minimal nonnegative solution to the Riccati equation*

[TABLE]

Furthermore,

[TABLE]

with

[TABLE]

${\Psi}_{\oplus\ominus}^{(1)}$* is the unique solution of the Sylvester equation*

[TABLE]

and with

[TABLE]

and $\Psi_{+-}^{(1)}$ is the unique solution of the Sylvester equation

[TABLE]

where

[TABLE]

**Proof **

To remove the effect of $\varepsilon^{-1}$ as $\varepsilon\rightarrow 0$, we need to combine the transformations of the previous two theorems. We pre-multiply the Riccati equation by $\operatorname*{diag}(I,\varepsilon I)$ and we use the matrix $\Gamma(\varepsilon)=\varepsilon^{-1}\Psi_{+\ominus}(\varepsilon)$. We obtain a new fixed-point equation, from which we eventually prove, by following the same steps as in Theorem 3.2 and Theorem 3.5, that the solutions are matrices of analytic functions.

Observe the terms in $\varepsilon^{-1}$ in the equations (66) to (69):

• •

we conclude from (66) that $\Psi_{+\ominus}=0$ by a similar argument to the proof of Theorem 3.5;

• •

multiply (68) by $\varepsilon$ and let $\varepsilon$ tend to zero to obtain the Riccati equation (80) satisfied by $\Psi_{\oplus\ominus}$;

• •

multiply (69) by $\varepsilon$ and let $\varepsilon$ tend to zero, gives (79), taking into account that $\lim_{\varepsilon\rightarrow 0}\Psi(\varepsilon)=\Psi$, an equality that is proved below.

To determine $\Psi_{+-}(0)$ is more involved. We proceed as follows. First, from (66), we take the terms in $\varepsilon^{0}$ and we find the expression (82) for $\Psi^{(1)}_{+\ominus}$ that we replace in (67). From (67), we take the terms in $\varepsilon^{0}$ and obtain $\Psi_{+-}$, after a reorganization of the terms, as the minimal nonnegative solution to

[TABLE]

with

[TABLE]

where

[TABLE]

To prove that the matrix $T$ is identical to the matrix $Q$ defined in (14), we only need to show that the matrix made up of the four blocks labeled with $D$s is equal to $(-A_{00}^{-1})$, partitionned according to (64), as

[TABLE]

where

[TABLE]

By (80), we have

[TABLE]

so that

[TABLE]

using (77). We write

[TABLE]

so that $B_{\oplus\oplus}=D_{\oplus\oplus}$.

Next, we have

[TABLE]

By (76), $-|\tilde{C}_{\ominus}|U_{\ominus\ominus}^{(-1)}+A_{\ominus\oplus}\Psi_{\oplus\ominus}$ simplifies to $-A_{\ominus\ominus}$ so that $B_{\ominus\oplus}=D_{\ominus\oplus}$.

Then, we have

[TABLE]

and we use (80) to replace $\tilde{C}_{\oplus}^{-1}A_{\oplus\ominus}$ in the first term to write

[TABLE]

We use (76), to write the second term as

[TABLE]

were we used (80) to replace $K_{\oplus\oplus}^{(-1)}\Psi_{\oplus\ominus}$. We find thus $B_{\oplus\ominus}=D_{\oplus\ominus}$.

Finally, we use the definition of $U_{\ominus\ominus}^{(-1)}$ to write

[TABLE]

We write

[TABLE]

by (77), so that $B_{\ominus\ominus}=D_{\ominus\ominus}$.

We find the block $\Psi_{\oplus-}^{(1)}$ of $\Psi^{(1)}$ given in (84) by observing the terms in $\varepsilon^{0}$ in (69). From (68), we obtain the Sylvester equation (84) for $\Psi_{\oplus\ominus}^{(1)}$. Taking the terms in $\varepsilon$ in (66) and (68) leads respectively to (3.7) and (3.7). $\square$

Remark 3.8

Not surprisingly, $\Psi_{+\ominus}=0$, as we found in (53).

As in Section 3.2, (79) is a function of $\Psi$ but also of the supplementary component $\Psi_{\oplus\ominus}$. This generalizes $\Psi_{\oplus-}$ given in (42). There is a probabilistic interpretation similar to the one given in (42), with, here, a correction term due to the introduction of $\mathcal{S}_{\ominus}$: $[\Psi_{\oplus-}]_{ij}$ is the sum of

• •

$[(-K_{\oplus\oplus}^{(-1)})^{-1}\tilde{C}^{-1}_{\oplus}A_{\oplus-}]_{ij}$, the probability that the phase process goes from $i$ to $j$, after some time spent in phases of $\mathcal{S}_{\oplus}$ or $\mathcal{S}_{\ominus}$,

• •

$[(-K_{\oplus\oplus}^{(-1)})^{-1}\tilde{C}^{-1}_{\oplus}A_{\oplus+}\Psi]_{ij}$, the probability that the process leaves $i$ for a phase in $\mathcal{S}_{+}$ and later returns to the initial level in $j$,

• •

$[(-K_{\oplus\oplus}^{(-1)})^{-1}\Psi_{\oplus\ominus}|\tilde{C}_{\ominus}^{-1}|A_{\ominus-}]_{ij}$ the probability that the process comes back to the initial level in a phase of $\mathcal{S}_{\ominus}$ and goes to $j$,

• •

$[(-K_{\oplus\oplus}^{(-1)})^{-1}\Psi_{\oplus\ominus}|\tilde{C}_{\ominus}^{-1}|A_{\ominus+}\Psi]_{ij}$ the process comes back to the initial level in a phase of $\mathcal{S}_{\ominus}$, goes to a phase of $\mathcal{S}_{+}$ and later returns to the initial level in $j$,

for $i\in\mathcal{S}_{\oplus}$, $j\in\mathcal{S}_{-}$.

Remark 3.9

Higher order terms (in particular, the coefficients of $\varepsilon^{2}$) may be of interest in some cases. It is clear that the principal difficulty lies in the necessity to deal with calculations that are steadily more cumbersome, but no more. We expect that coefficients of $\Psi_{+\ominus}$ or $\Psi_{\oplus-}$ will be given explicitly and that each successive coefficients of $\Psi_{+-}$ and $\Psi_{\oplus\ominus}$ will be solutions of Sylvester equations.

4 Impact on the stationary probability

For $j\in\mathcal{S}$ and $x\in\mathbb{R}^{+}$, we define the joint distribution function of the level and the phase at time $t$, $F_{j}(x,t)=\mathbb{P}\left[X\left(t\right)\leq x,\varphi\left(t\right)=j\right],$ and its density by

[TABLE]

The stationary density vector ${\boldsymbol{\pi}}(x)=(\pi_{j}(x):j\in\mathcal{S})$ of the fluid model, where, for $j\in\mathcal{S}$, $\pi_{j}(x)=\lim_{t\rightarrow\infty}f_{j}(x,t),$ exists if and only if the mean stationary drift is negative, that is, if and only if $\sum_{i\in\mathcal{S}}\xi_{i}c_{i}<0$, where $\xi_{i}$ is defined in (4) for all $i$. When the mean stationary drift of the fluid model is negative, from Govorun et al. [6], we have, for $x>0$,

[TABLE]

and the mass at zero is $[0\;;{\boldsymbol{p}}_{-}\;;{\boldsymbol{p}}_{0}]$ where

[TABLE]

and $[\boldsymbol{p}_{-}\;;{\boldsymbol{p}}_{0}]$ is the unique solution of the system

[TABLE]

Expression (88) is numerically stable and has a physical interpretation (da Silva Soares [38, Chapter 1, Section 1.3]). Furthermore, it appears clearly that all the quantities appearing in the expression of the stationary density are functions of $\Psi$.

The stationary density of (17) may be formulated as

[TABLE]

where $K(\varepsilon)$, $\Theta(\varepsilon)$ and ${\boldsymbol{q}}(\varepsilon)$ are defined similary to (89),(90) and (91) respectively. It is well known that the stationary density vector $\boldsymbol{\pi}(x,\varepsilon)$ is differentiable (see Kato [9, Section 2]) and such that $\boldsymbol{\pi}(x,\varepsilon)$ may be written as

[TABLE]

where

[TABLE]

for all $x\in\mathbb{R}^{+}$. We find

[TABLE]

where $\Psi^{(1)}$ is given in Theorem 2.1 and

[TABLE]

The vector ${\boldsymbol{q}}(\varepsilon)$ is differentiable by Kato [9, Section 2] and

[TABLE]

with

[TABLE]

where $M^{\#}$ denotes the group inverse of the matrix $M$. We find (4) by solving the Poisson equation (see Meyer [12]) satsified by $[{\boldsymbol{p}}_{-}^{(1)}\;;{\boldsymbol{p}}_{0}^{(1)}]$, deduced from (92), where $c$ is a normalisation found through (93). Finally,

[TABLE]

where $K^{(1)}=C_{+}^{-1}\tilde{Q}_{++}+\Psi^{(1)}|C_{-}^{-1}|Q_{-+}+\Psi|C_{-}^{-1}|\tilde{Q}_{-+}$. In order to actually compute $L^{(1)}(x)$, we refer the reader to Higham [8, Theorem 10.13, Equation (10.17a)].

5 Numerical Illustration

We evaluate and display the value of $E_{\infty}(\varepsilon)=\parallel\boldsymbol{\Psi}(\varepsilon)-\overline{\Psi}-\varepsilon\Psi^{(1)}\parallel_{\infty}$ for a few examples where only the phases in $\mathcal{S}_{0}$ are perturbed, either by a positive quantity as in Section 3.2, or by a negative quantity, as in Section 3.3. The narrative is as follows: assume that the rates $c_{i}$ in $\mathcal{S}_{0}$ are very small and positive (or very small and negative) and that they are set equal to 0. What is the effect on the matrix $\Psi$?

The controlling phase evolves as a birth-and-death process on the state space $\{1\ldots 3m\}$ where

• •

phases 1 to $m$ belong to $\mathcal{S}_{+}$ and all have the same positive rate $r_{+}$;

• •

phases $m+1$ to $2m$ belong to $\mathcal{S}_{0}$; when perturbed, they all have the same perturbation coefficient $\tilde{r}_{0}$, either positive or negative;

• •

phases $2m+1$ to $3m$ belong to $\mathcal{S}_{-}$ and all have the same negative rate $r_{-}$.

The parameters $r_{+}$ and $r_{-}$ are chosen in all cases such that the stationary drift of the non-perturbed fluid model is equal to $-0.1$.

Migration of $\mathcal{S}_{0}$ to $\mathcal{S}_{+}$.

Case 1.a

The infinitesimal generator is that of the M/M/1/N queue with $N=3m$, that is,

[TABLE]

We assume that $\lambda>\mu$, so that the process spends most of its time in $\mathcal{S}_{-}$ in this case. In our experimentation, we have noticed that the quantities

[TABLE]

and

[TABLE]

are significantly different sometimes, and for that reason we give separately their values in the figures of this section, the norm $E_{\infty}$ is easily found as the maximum of $E_{+}$ and $E_{\oplus}$.

The results on Figure 1 are obtained from $\lambda=2>\mu=1$, $m=5$, $r_{+}=0.4$, $r_{-}=-0.207$, $\tilde{r}_{\oplus}=r_{+}>0$. On the left, $E_{+}(\varepsilon)$ is displayed as a continuous line, marked with a ’+’ sign, $E_{\oplus}(\varepsilon)$ is displayed as a dashed line, marked with an ’o’ sign. On the right of Figure 1, we display the same functions on a logarithmic scale; in addition, we include for visual reference a function proportional to $\varepsilon^{2}$, as the lined marked with alternating dashes and dots.

The logarithm plot shows in a striking manner that the difference $\boldsymbol{\Psi}(\varepsilon)-\overline{\Psi}-\varepsilon\Psi^{(1)}$ is $O(\varepsilon^{2})$. This may also be seen in Table 1, where we give the values of $E_{+}$ and $E_{\oplus}$ for $\varepsilon=10^{-4}$ and $\varepsilon=10^{-2}$, for Cases 1.a, 2.a and 3.a.

Case 2.a

The infinitesimal generator for the phase is given by (97), the same as in Case 1.a, but here we take $\lambda<\mu$, so that the process spends most of its time in $\mathcal{S}_{+}$. The parameters in this case are $\lambda=1<\mu=2$, $m=5$, $r_{+}=0.4$, $r_{-}=-621$, $\tilde{r}_{\oplus}=r_{+}>0$. Notice that we must use a very small value for $r_{-}$, to compensate for the time spent in $\mathcal{S}_{+}$ and keep $-0.1$ as the stationary drift.

The results for this case are displayed on the left of Figure 2 and it appears that $E_{+}$ and $E_{\oplus}$ are very close to each other. Furthermore, they are much smaller that in Case 1.a. This is very clear from Table 1, where we observe that $E_{\infty}$ is several orders of magnitude smaller in Case 2.a than in Case 1.a.

Case 3.a

The infinitesimal generator in Case 3.a is that of a system of $N$ individuals independently alternating between two states. It is given by

[TABLE]

where $d_{i,j}=i\beta+(N-j)\alpha$ and $N=3m$. We take $\alpha=\beta$ so that the distribution is concentrated in the middle of the range $[1\ldots N]$, that is, in the region covered by $\mathcal{S}_{0}$.

The results are given on the right in Figure 2 and in the last row of Table 1, the parameters are $m=5$, $\alpha=\beta=1$, $r_{+}=0.4$, $r_{-}=-2.63$ and $\tilde{r}_{\oplus}=r_{+}$.

Migration of $\mathcal{S}_{0}$ to $\mathcal{S}_{-}$

In the second set of examples, labeled Cases 1.b to 3.b, we take the same parameters as in Cases 1.a to 3.a, except that the perturbation for the phases in $\mathcal{S}_{0}$ are negative, and we take in each case $\tilde{r}_{\ominus}=-r_{+}$. Here, there is only one set of functions $E(\cdot)$ to compute: the rows of $\Psi$ are all labeled by phases in $\mathcal{S}_{+}$ and so $E_{\infty}=E_{+}$, while $E_{\oplus}$ is not defined.

Graphically, we have observed results very similar to those in Figures 1 and 2 and we do not give the graphs here. Instead, we give in Table 2 the values of $E_{\infty}(\varepsilon)$ for $\varepsilon=10^{-4}$ and $10^{-2}$. The obtained values are similar to those obtained for Cases 1.a to 3.a and it appears clearly that $E_{\infty}(\varepsilon)$ is $O(\varepsilon^{2})$.

We might also compute the 1-norm of $\boldsymbol{\Psi}(\varepsilon)-\overline{\Psi}-\varepsilon\Psi^{(1)}$ instead of its $\infty$-norm, and compare the two partial norms

[TABLE]

and

[TABLE]

with $\parallel\boldsymbol{\Psi}(\varepsilon)-\overline{\Psi}-\varepsilon\Psi^{(1)}\parallel_{1}$ = $\max(E_{-}(\varepsilon),E_{\ominus}(\varepsilon))$. We would expect to observe differences similar to those between $E_{+}$ and $E_{\oplus}$ in Cases 1.a to 3.a.

Acknowledgement

This work was supported in part by the Ministère de la Communauté française de Belgique through the ARC grant AUWB-08/13-ULB 5 and in part by the Flemish Community of Belgium through the Methusalem program.