**3. Time-Varying** *M***/***M***/2 System**

We start with the well-known time-varying *M*/*M*/2/∞ system with two servers and the infinite-capacity queue in which customers arrive one by one with the intensity *λ*(*t*). The service intensity of each server does not depend on the total number of customers in the queue and is equal to *μ*(*t*). The functions *λ*(*t*) and *μ*(*t*) are assumed to be nonrandom, nonnegative and locally integrable on [0, ∞) continuous functions. Let the integer-valued time-dependent random variable *X*(*t*) denote the total number of customers in the system at time *t* ≥ 0. Then *X*(*t*) is the CTMC with the state space {0, 1, 2 ... }. Its transposed time-dependent intensity matrix (generator) *A*(*t*)=(*aij*(*t*))<sup>∞</sup> *<sup>i</sup>*,*j*=<sup>0</sup> has the form

$$A(t) = \begin{pmatrix} -\lambda(t) & \mu(t) & 0 & 0 & \dots \\ \lambda(t) & -(\lambda(t) + \mu(t)) & 2\mu(t) & 0 & \dots \\ 0 & \lambda(t) & -(\lambda(t) + 2\mu(t)) & 2\mu(t) & \dots \\ 0 & 0 & \lambda(t) & -(\lambda(t) + 2\mu(t)) & \dots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{pmatrix}.$$

For all *t* ≥ 0 we represent the distribution of *X*(*t*) as a probability vector **p**(*t*), where **p**(*t*) = ∑<sup>∞</sup> *<sup>k</sup>*=<sup>0</sup> P(*X*(*t*) = *k*)**e***<sup>k</sup>* (as above, **e***<sup>k</sup>* denotes the *k*th unit basis vector). Given any proper initial condition **p**(0), the Kolmogorov forward equations for the distribution of *X*(*t*) can be written as

$$\frac{d}{dt}\mathbf{p}(t) = A(t)\mathbf{p}(t). \tag{5}$$

Assume that *X*(*t*) is null ergodic. The condition on the intensities *λ*(*t*) and *μ*(*t*), which guarantees null ergodicity will be derived shortly below, (clearly, if the intensities are constants, i.e., *λ*(*t*) = *λ* and *μ*(*t*) = *μ*, then the condition is simply *λ* > 2*μ*. If both are periodic and the smallest common multiple of the periods is *T*, then the condition is *T* <sup>0</sup> *λ*(*u*) *du* > 2 *T* <sup>0</sup> *μ*(*u*) *du*). Fix a positive number *d* > 1 and define the sequence {*δn*, *n* ≥ <sup>0</sup>} by *<sup>δ</sup><sup>n</sup>* <sup>=</sup> *<sup>d</sup>*−*n*. It is the decreasing sequence of positive numbers. By multiplying (5) from the right with Λ = *diag*(*δ*0, *δ*1,...), we get

$$\frac{d}{dt}\vec{\mathbf{p}}(t) = \vec{A}(t)\vec{\mathbf{p}}(t),\tag{6}$$

where **<sup>p</sup>**˜(*t*) = <sup>Λ</sup>**p**(*t*) and *<sup>A</sup>*˜(*t*) = <sup>Λ</sup>*A*(*t*)Λ−1. Denote by <sup>−</sup>*α*˜ *<sup>k</sup>*(*t*) the sum of all elements in the *k*th column of *A*˜(*t*). By direct inspection it can be checked that

$$\begin{aligned} \bar{\mathfrak{a}}\_0(t) &= \left(1 - d^{-1}\right) \lambda(t), \\ \bar{\mathfrak{a}}\_1(t) &= \left(1 - d^{-1}\right) (\lambda(t) - d\mu(t)), \\ \bar{\mathfrak{a}}\_k(t) &= \underbrace{\left(1 - d^{-1}\right) (\lambda(t) - 2d\mu(t))}\_{=\mathfrak{f}(t)}, k \ge 2. \end{aligned}$$

Since *α*˜ <sup>0</sup>(*t*) ≥ *β*(*t*) and *α*˜ <sup>1</sup>(*t*) ≥ *β*(*t*), the upper bound follows from (4) applied to (6):

$$\sum\_{k=0}^{\infty} d^{-k} p\_k(t) \le e^{-\int\_0^t \beta(u) \, du} \sum\_{k=0}^{\infty} d^{-k} p\_k(0). \tag{7}$$

If *d* is chosen such that *d* > 1 and <sup>∞</sup> <sup>0</sup> (*λ*(*t*) − 2*dμ*(*t*)) *dt* = +∞, then from (7) it follows that *pk*(*t*) → 0 as *t* → ∞ for each *k* ≥ 0 and thus *X*(*t*) is null ergodic. In such a case it is possible to extract more information from (7). Note that for any fixed *n* ≥ 0 it holds that

$$d^{-n}\sum\_{i=0}^{n}p\_i(t) \le \sum\_{k=0}^{n}d^{-k}p\_k(t) \le \sum\_{k=0}^{\infty}d^{-k}p\_k(t).$$

Thus, if *X*(0) = *N*, i.e., *pN*(0) = 1 then for any *n* ≥ 0 the following upper bound for the conditional probability P(*X*(*t*) ≤ *n*|*X*(0) = *N*), *N* ≥ 0, holds:

$$\mathbb{P}(X(t)\le n|X(0)=N)\le d^{n-N}e^{-\int\_0^t \mathbb{f}(u)\,du}.\tag{8}$$

Now assume that *X*(*t*) is weakly ergodic (the corresponding condition on the intensities *λ*(*t*) and *μ*(*t*) will be derived shortly below). Using the normalization condition *<sup>p</sup>*0(*t*) = <sup>1</sup> − <sup>∑</sup>*i*≥<sup>1</sup> *pi*(*t*) it can be checked that the system (5) can be rewritten as follows:

$$\frac{d}{dt}\mathbf{z}(t) = B(t)\mathbf{z}(t) + \mathbf{f}(t),\tag{9}$$

where the matrix *B*(*t*) with the elements *bij*(*t*) = *aij*(*t*) − *ai*0(*t*) has no probabilistic meaning and the vectors **f**(*t*) and **z**(*t*) are

$$\mathbf{f}(t) = (\lambda(t), 0, 0, \dots)^T, \ \mathbf{z}(t) = (p\_1(t), p\_2(t), \dots)^T.$$

Let **z**∗(*t*) and **z**∗∗(*t*) be the two solutions of (9) corresponding to two different initial conditions **z**∗(0) and **z**∗∗(0). Then for the vector **y**(*t*) = **z**∗(*t*) − **z**∗∗(*t*) = (*y*1(*t*), *y*2(*t*),...) *T*, with arbitrary elements we have the system

$$\frac{d}{dt}\mathbf{y}(t) = \mathcal{B}(t)\mathbf{y}(t). \tag{10}$$

The matrix *B*(*t*) in (10) may have negative off-diagonal elements. But it is straightforward to see, that the similarity transformation *TB*(*t*)*T*−<sup>1</sup> = *B*∗(*t*), where *T* is the upper triangular matrix of the form

$$T = \begin{pmatrix} 1 & 1 & 1 & \cdots \\ 0 & 1 & 1 & \cdots \\ 0 & 0 & 1 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix}\_{\ell'} $$

gives the matrix *B*∗(*t*):

$$B^\*(t) = \begin{pmatrix} -\left(\lambda(t) + \mu(t)\right) & \mu(t) & 0 & 0 & \cdots\\ \lambda(t) & -\left(\lambda(t) + 2\mu(t)\right) & 2\mu(t) & 0 & \cdots\\ 0 & \lambda(t) & -\left(\lambda(t) + 2\mu(t)\right)2\mu(t) & \cdots\\ \vdots & \vdots & \vdots & \vdots & \ddots \end{pmatrix},\tag{11}$$

which off-diagonal elements are always nonnegative. Let **u**(*t*) = *T***y**(*t*) = (*u*1(*t*), *u*2(*t*),...) *T*. Then by multiplying both parts of (10) from the left by *T*, we get

$$\frac{d}{dt}\mathbf{u}(t) = B^\*(t)\mathbf{u}(t). \tag{12}$$

Fix a positive number *d* > 1 and define the increasing sequence of positive numbers {*δn*, *<sup>n</sup>* <sup>≥</sup> <sup>0</sup>} by *<sup>δ</sup><sup>n</sup>* <sup>=</sup> *<sup>d</sup>n*−1. Let *<sup>D</sup>* <sup>=</sup> *diag*(*δ*1, *<sup>δ</sup>*2, ...). By putting **<sup>w</sup>**(*t*) = *<sup>D</sup>***u**(*t*) in (12), we obtain the system of equations

$$\frac{d}{dt}\mathbf{w}(t) = B^{\*\*}(t)\mathbf{w}(t),\tag{13}$$

where the matrix *B*∗∗(*t*) = *DB*∗(*t*)*D*−<sup>1</sup> has nonnegative off-diagonal elements. Denote by −*αk*(*t*) the sum of all elements in the *k*th column of *B*∗∗(*t*) i.e.

$$\begin{aligned} \mathfrak{a}\_1(t) &= \mu(t) - (d-1)\lambda(t), \\ \mathfrak{a}\_2(t) &= \left(1 - d^{-1}\right)\mu(t) + \mu(t) - (d-1)\lambda(t), \\ \mathfrak{a}\_k(t) &= \underbrace{\left(1 - d^{-1}\right)(2\mu(t) - d\lambda(t))}\_{=\mathfrak{f}(t)}, \ k \ge 3. \end{aligned}$$

Note that if 1 < *d* ≤ 2 then *α*1(*t*) ≥ *β*(*t*) and *α*2(*t*) ≥ *β*(*t*). Now, remembering that **w**(*t*) = *D***u**(*t*) = *DT***y**(*t*), the upper bound for **y**(*t*) = **z**∗(*t*) − **z**∗∗(*t*) in the weighted norm due to (4) is (from (14) the purpose of the similarity transformation *DB*∗(*t*)*D*−<sup>1</sup> can be recognized: it is to make *β*(*t*) in the exponent as large as possible).

$$\|\|DT\mathbf{y}(t)\|\| \le e^{-\int\_0^t \beta(u) \, du} \|\|DT\mathbf{y}(0)\|\|. \tag{14}$$

The upper bound for **p**∗(*t*) − **p**∗∗(*t*) is obtained from (14). Firstly notice that **y**(*t*)- ≤ 2**p**∗(*t*) − **p**∗∗(*t*) since **y**(*t*) is the solution of (10)—the system with the excluded state (0). Secondly, it can be proved, (this is shown, for example, in [Equation (18)] of the [29]), that **x**- ≤ 2-*DT***x**for any vector **x**. Hence

$$\|\mathbf{p}^\*(t) - \mathbf{p}^{\*\*}(t)\| \le 4e^{-\int\_0^t \beta(u) \, du} \|DT\mathbf{y}(0)\|. \tag{15}$$

If *d* is chosen such that *d* > 1 and <sup>∞</sup> <sup>0</sup> (2*μ*(*t*) − *dλ*(*t*)) *dt* = +∞, then from (15) it follows that **p**∗(*t*) − **p**∗∗(*t*)- → 0 as *t* → ∞ for any initial conditions **p**∗(0) and **p**∗∗(0), i.e., *X*(*t*) is weakly ergodic. Note that it is sufficient to choose *d* ∈ (1, 2]: if the integral diverges for *d* > 2 it also diverges for *d* = 2 and this is sufficient for (14) to hold.

Sometimes it is also possible to obtain bounds similar to (15) for other characteristics of *X*(*t*). For example, denote by *E*(*t*, *k*) the conditional mean number of customers in the system at time *t*, given that initially there where *k* customers in the system, i.e., *<sup>E</sup>*(*t*, *<sup>k</sup>*) = <sup>∑</sup>*n*≥<sup>1</sup> *<sup>n</sup>*P(*X*(*t*) = *<sup>n</sup>*|*X*(0) = *<sup>k</sup>*). Then using [Equation (22)] of [29] it can be shown, that

$$|E(t,k) - E(t,0)| \le \frac{4(1-d^k)}{W(1-d)}e^{-\int\_0^t \beta(u) \, du}, \ k \ge 1, \ W = \inf\_n \frac{d^n}{n+1}.\tag{16}$$

The results obtained above for both, null and weak ergodic, cases can be put together in the single theorem.

**Theorem 1.** *Let there exist a positive <sup>d</sup>* <sup>=</sup> <sup>1</sup> *such that* <sup>∞</sup> 0 (<sup>1</sup> <sup>−</sup> *<sup>d</sup>*)*λ*(*t*) + <sup>2</sup>(<sup>1</sup> <sup>−</sup> *<sup>d</sup>*−1)*μ*(*t*) *dt* = +∞*. Then X*(*t*) *is null (weakly) ergodic if d* < 1 *(d* > 1*) and the ergodicity bounds (7) and (15) hold.*

Whenever the intensities *λ*(*t*) and *μ*(*t*) are constants or periodic functions stronger results can be obtained.

**Corollary 1.** *If in the Theorem 1 the intensities λ*(*t*) *and μ*(*t*) *are constants or* 1−*periodic, (i.e., λ*(*t*) *and μ*(*t*) *are periodic functions and the length of their periods is equal to one), then X*(*t*) *is exponentially null (weakly) ergodic if d* < 1 *(d* > 1*) and there exist R* > 0 *and a* > 0 *such that <sup>e</sup>*<sup>−</sup> *<sup>t</sup> <sup>s</sup> <sup>β</sup>*(*u*) *du* <sup>≤</sup> *Re*−*a*(*t*−*s*) *for* <sup>0</sup> <sup>≤</sup> *<sup>s</sup>* <sup>≤</sup> *t.*

We now consider the numerical example. Let *λ*(*t*) = 9(1 + sin 2*πt*) and *μ*(*t*) = 8(1 + cos 2*πt*). It is straightforward to check from the *Theorem 1* that if *d* = <sup>4</sup> 3 then *X*(*t*) is weakly ergodic. Then the ergodicity bounds follow from (15) and (16):

$$\|\mathbf{p}^\*(t) - \mathbf{p}^{\*\*}(t)\| \le 72e^{-t} \|DT(\mathbf{p}^\*(0) - \mathbf{p}^{\*\*}(0))\|,\tag{17}$$

$$|E(t,k) - E(t,0)| \le 162\left(\frac{4}{3}\right)^k e^{-t}, \ k > 0. \tag{18}$$

Figure 1 shows the graph of the probability *p*0(*t*) as *t* increases. It can be seen that for any initial condition **p**(0) there exists one periodic function of *t*, say *π*0(*t*) (i.e., *π*0(*t*) = *π*0(*t* + *T*), where *T* = 1 is the smallest common multiple of the periods of *λ*(*t*) and *μ*(*t*)), such that lim*t*→∞(*p*0(*t*) − *π*0(*t*)) = 0. Figure 2 shows the detailed behaviour of *π*0(*t*). Now consider (17). If *t* ≥ 37 then the right part of (17) does not exceed 10−<sup>3</sup> i.e., starting from the instant *t* = 37 = *t* ∗ the system "forgets" its initial state and the distribution of *X*(*t*) for *t* > *t* <sup>∗</sup> can be regarded as limiting. The error (in *l*1-norm), which is thus made, is not greater than 10−3. Moreover, since the limiting distribution of *X*(*t*) is periodic, it is sufficient to solve numerically the system of ODEs only in the interval [0, *t* ∗ + *T*]. The distribution of *X*(*t*) in the interval [*t* ∗, *t* ∗ + *T*] is the limiting probability distribution of *X*(*t*) (with error not greater than 10−<sup>3</sup> in *l*1-norm). Note that the system of ODEs contains infinite number of equations. Thus in order to solve it numerically one has to truncate it; this truncation was performed according to the method in [30]. The upper bound on the rate of convergence of the conditional mean *E*(*t*, *k*) is given in (18). If *t* ≥ *t* ∗ then the right part does not exceed 10−<sup>2</sup> i.e., starting from *t* = *t* ∗ the system "forgets" its initial state and the value of *E*(*t*, *k*) can be regarded as the limiting value of the conditional mean number of customers with the error not greater than 10−2. The rate of convergence of *E*(*t*, *k*) and the behaviour of its limiting value is shown in the Figures 3 and 4. Note that the obtained upper bounds are not tight: the system enters the periodic limiting regime before the instant *t* = *t* ∗.

**Figure 1.** Rate of convergence of the empty system probability *p*0(*t*) in the interval [0, 37] given two different initial conditions: *p*0(0) = 1 (**red line**), *p*189(0) = 1 (**blue line**).

**Figure 2.** Limiting probability *p*0(*t*) of the empty queue given two different initial conditions: *p*0(0) = 1 (**red line**), *p*189(0) = 1 (**blue line**).

**Figure 3.** Rate of convergence of the conditional mean *E*(*t*, *k*) number of customers in the system in the interval [0, 37]: *E*(*t*, 0) (**red line**), *E*(*t*, 189) (**blue line**).

**Figure 4.** Limiting conditional mean *E*(*t*, *k*) number of customers in the system: *E*(*t*, 0) (**red line**), *E*(*t*, 189) (**blue line**).
