An M/PH/1 queue with workload-dependent processing speed and vacations

Abstract

Motivated by the trade-off issue between delay performance and energy consumption in modern computer and communication systems, we consider a single-server queue with phase-type service requirements and with the following two special features: Firstly, the service speed is a piecewise constant function of the workload. Secondly, the server switches off when the system becomes empty, only to be activated again when the workload reaches a certain threshold. For this system, we obtain the steady-state workload distribution and its moments of any order. We use this result to choose the activation threshold such that a certain cost function, involving processing costs, activation costs and mean workload, is minimized.

1 Introduction

We consider an M/PH/1 queue with two special features: Firstly, the service speed is a piecewise constant function of the workload. Secondly, the server switches off when the system becomes empty, and it is reactivated when the workload reaches a certain threshold. In the remainder of this section we discuss the motivation for our study (Sect. 1.1), mention related literature (Sect. 1.2) and outline the rest of the paper (Sect. 1.3).

1.1 Motivation

Cloud services have become ubiquitous in our modern information society. Cloud services are more and more important in the current era where we are facing the Covid-19 pandemic for which remote work is strongly recommended to decrease the spread of the virus. Nowadays, most Internet users are familiar with some cloud-based storage services such as Dropbox and Google Drive, and with video-conferencing services such as Zoom and Microsoft Teams. These systems are supported by large-scale data centers where thousands of servers are available, consuming a huge amount of energy. Thus, it is crucial to have mechanisms that can balance energy consumption and performance.

While energy saving is indispensable to reduce CO${}_2$, most data centers are still designed for peak traffic. As a result, many servers are idle in off-peak periods – but they may still consume about 60% of their peak energy consumption [9, 15]. Autoscaling techniques are proposed to solve the trade-off between energy consumption and delay performance. In particular, an autoscaling algorithm adjusts the processing speed according to its workload. At the data center level, an autoscaling algorithm controls the number of active servers in the system in response to the workload [9, 15, 18, 21]. Furthermore, at the individual computer level, a CPU is also able to adjust the processing speed by either dynamic frequency scaling or dynamic voltage scaling techniques [28,29,30]. The processing speed is scaled up when the workload is high and scaled down under a low workload. As a result, less energy is consumed under off-peak periods while acceptable delay can be achieved in a heavy-load situation [7].

Queues with changeable service rate also fit the autoscaling mechanisms in 5G networks. In 5G networks, the key technology is network function virtualization (VNF) in which a physical resource can be virtualized to network functions. The operator can dynamically add or release these network functions to optimally construct cost-effective systems in response to their workload. As a result, these systems can be modeled using queueing models with changeable service speeds [20, 21].

Apart from the interest in power-saving computer and communication systems, queues with changeable service speed also naturally arise in service systems with human servers. In particular, in many real-world service systems such as call centers, staff are scheduled to meet the demands of customers. Also, a human server may serve at high speed when the workload is high and may spend more time on a job when the workload is low [5].

1.2 Related literature and our contribution

The topic of speed scaling in data centers, power-saving in CPUs and autoscaling of 5G networks provides motivation for our study. See [28] for an insightful discussion of speed scaling and [20, 21] for queueing analysis of an autoscaling algorithm in 5G networks. Recent papers which consider single server queues with speed scaling where the speed of the server is proportional to the number of jobs in the system are, for example [19, 29, 30]. Multiserver queues with ON-OFF control (turning idle servers off) have been extensively studied [9, 15, 18].

Our model is related to the literature on queues and dams with a level-dependent outflow rate. Influential early papers are [10, 11]; we refer to [2] for some more recent results and further references.

Our model is also related to the rich vacation literature: the server takes a vacation when the workload is zero, and returns to service when the workload reaches or exceeds a certain level. Such a D-policy has been extensively studied for the classical M/G/1 queue. See [8] for references and, in particular, for an optimality proof. For the case of switching costs and running costs, and with a holding cost per time unit which is a non-negative decreasing right-continuous function of the current workload, Feinberg and Kella [8] prove that D-policies are optimal for the average-cost-per-time-unit criterion: there is an optimal policy that either runs the server all the time or switches the server off when the system becomes empty and switches it on when the workload reaches or exceeds some threshold D.

Markov modulated queues with workload-dependent service rates have been extensively studied in the literature [4, 6, 13, 16, 26, 31]. All these studies first reduce the models to fluid models for which matrix analytic methods, spectral methods, and Schur decomposition methods are applied to derive numerical solutions for the distribution of the workload processes (see da Silva Soares and Latouche [4], Mandjes et al. [16] and Kankaya and Akar [13] for the methodologies). In addition to those studies, Bean and O’Reilly [1] consider a similar fluid model with a set of control levels, study the LST of various passage and sojourn times, and present algorithms for their evaluation. In contrast to the above, motivated by power-saving in modern computer and communication systems, we consider a model with Poisson input, phase-type service-time distribution, and vacation. For this model, using the level-crossing method and renewal theory, we obtain a direct solution for the workload distribution and its moments of any order. By virtue of the Poisson input and phase-type distribution, our solution is more direct in the sense that the workload distribution is expressed in terms of matrix exponentials whose components are explicitly written in terms of given parameters. These results are then used for an optimization problem balancing performance and energy-consumption tradeoff.

We finally would like to point out the relation to our recent paper [22]. There we propose and analyze an M/G/1-type queueing model that also features two power-saving mechanisms. The speed of the server is scaled according to the workload in the system. Moreover, the server is turned off when the system is empty and is activated again once the workload reaches a certain threshold. In the case of arbitrarily distributed service time and general service speed function, the stationary workload is expressed in terms of a series whose terms are recursively obtained by an integral formula. While the distribution of the workload can in principle be evaluated for this general case, the computation is highly complex. Simpler expressions are obtained for the case of exponential service requirement and the case where the service speed is a linear function of the workload.

In the present paper, our aim is to derive a computable solution for a relatively general model with phase-type service requirements and a piecewise constant service rate. Phase-type distributions lie dense in the class of distributions with positive support and thus can approximate any service requirement distribution with any accuracy. Furthermore, piecewise constant functions can also be used to approximate an arbitrary function. In addition, from a practical point of view, it is natural that the service rate is switched at discrete points [7].

1.3 Structure of the paper

The remainder of the paper is organized as follows: In Sect. 2 we derive the stationary workload distribution for the case of an M/PH/1 queue with vacations and with only two processing speeds. We first focus on this case, before tackling the case of an arbitrary number of different processing speeds in Sect. 3, because the analysis is quite technical; in this way, we improve the readability of the paper. We also obtain a computable form for the moments of any order for the stationary workload. In Sect. 4, we demonstrate the analysis of the active period for an M/M/1 with two processing speeds. In Sect. 5, we study a cost-optimization problem and show some numerical examples.

2 M/PH/1 queue with two processing speeds and vacations

In this section we consider the special case of a FCFS queue with a single server who, when active, works at one of two possible speeds: when the workload is below a threshold $d_1$ it works at speed $r_1$, and above it at speed $r_2$. Furthermore, when the system becomes empty, the server is switched off, only to be activated again when the workload exceeds some threshold level L. We assume that $L=d_1$ to keep the model for the moment as simple as possible while retaining its essential elements. In Sect. 3, we shall consider the more general case of piecewise constant service speed with $K_0$ different speed values, and in which L does not necessarily coincides with one of the thresholds at which the service speed in an active period changes. However, the analysis of that case is quite involved, and consideration of the simpler case of the present section will help the reader get acquainted with our approach.

The remainder of the section is organized as follows: Section 2.1 contains a more detailed model description, as well as a lemma about the computation of the convolution of matrix exponentials that will play a key role in the remainder of the paper. In Sect. 2.2, we derive the stationary workload density when the server is inactive. In Sects. 2.3 and 2.4, we successively determine the stationary workload density when the server is active while the workload is above, respectively below, $d_1$. The mean active period features in many of the formulas; we compute it explicitly in Sect. 2.5.

2.1 Model description

We consider a single-server FCFS queue, where the server has a single waiting line with infinite capacity. Customers arrive according to a Poisson process with rate $\lambda > 0$. The service requirements of the customers are independent and identically distributed (i.i.d.), generically indicated by B, with the following phase-type distribution (see, for example [14]):

$$\begin{aligned} B(x) = 1 - \varvec{\tau } \exp (\varvec{T} x) \varvec{1}, \quad x \ge 0, \end{aligned}$$

(2.1)

where $\varvec{\tau }$ is a $(1 \times N)$ probability row vector, and $\varvec{T}$ is an $(N \times N)$ defective transition rate matrix, where N is a positive integer, and $\varvec{1}$ is a column vector with ones whose size will be determined by the context. The tail of $B(\cdot )$ is denoted by $\overline{B}(x) := 1 - B(x)$.

Let $Z_{t}$, where $t \in \mathbb {R}_{+} := [0,\infty )$, denote the unfinished workload (workload, for short) in the system at time t. According to the value or history of the workload, the server is assumed to alternate between “inactive” and “active” states, which are referred to as modes 0 and 1, respectively. Let $S_{t} \in \{0,1\}$ denote the mode of the server at time t, which is defined as follows:

• (i) When the workload $Z_{t}$ hits 0 at a time t, the server enters mode 0, i.e., $S_{t} = 0$. After that, the server remains in mode 0 until the workload exceeds threshold $L=d_1 > 0$, i.e., $S_{u} = 0$ as long as $0 \le Z_{u} \le d_1$ for $u \ge t$.

• (ii) When the workload $Z_{t}$ exceeds $d_1$ at a time t while the mode was 0 at time $t-$, the server changes its mode from 0 to 1, i.e., $S_{t} = 1$. The server remains in that mode until the workload hits 0 again, i.e., $S_{u} = 1$ as long as $0< Z_{u} < \infty$ for $u \ge t$.

The processing speed of the server is assumed to depend on both the server’s mode and workload in the following way: Let $s_{i}(x)$, where $i = 0,1$ and $x \ge 0$, denote the processing speed of the server when $(Z_{t}, S_{t}) = (x,i)$; it is defined as follows:

$$\begin{aligned} s_{0}(x)&= 0, \quad 0 \le x \le d_1,\end{aligned}$$

(2.2)

$$\begin{aligned} s_{1}(x)&= {\left\{ \begin{array}{ll} r_{1}, &{} 0< x \le d_{1},\\ r_{2}, &{} d_{1}< x < \infty , \end{array}\right. } \end{aligned}$$

(2.3)

where $r_{1}$ and $r_{2}$ are positive numbers. $\{(Z_{t}, S_{t}) ; t \ge 0\}$ is a Markov process. Figure 1 shows a sample path of this Markov process.

Fig. 1

The workload process for the case of two processing speeds

We assume that it is positive-recurrent, and denote its replica in steady state by (Z, S). The stability condition of this Markov process is given as follows:

$$\begin{aligned} \lambda \varvec{\tau }(-\varvec{T})^{-1}\varvec{1} < r_{2}, \end{aligned}$$

(2.4)

where the left-hand side is the mean amount of work arriving in a unit of time, and the right-hand side is the processing rate of the server when the workload is greater than $d_{1}$.

We introduce some notation. For $i=1,2$, let $n_{i}$ and $k_{i}$ be positive integers such that $k_{i} \le n_{i}$. For an $n_{1} \times n_{2}$ matrix $\varvec{A} := (a_{ij} ; 1 \le i \le n_{1}, 1 \le j \le n_{2})$, we denote the $k_{1} \times k_{2}$ northeast sub-matrix of $\varvec{A}$ by $[\varvec{A}]^{(k_{1}, k_{2})}$, i.e.,

$$\begin{aligned}{}[\varvec{A}]^{(k_{1}, k_{2})} = (a_{ij} ; 1 \le i \le k_{1}, n_{2}-k_{2} < j \le n_{2}). \end{aligned}$$

(2.5)

For a positive integer n, let $\varvec{I}_{n}$ denote the $n \times n$ identity matrix, and $\varvec{O}$ denote a zero matrix whose size will be determined by the context. Throughout this paper, the next lemma is useful when we compute the convolution of matrix exponentials. It is a direct consequence of Equation (10.40) of [12], and the proof is given in a similar way to that of Theorem 1 of [27] (see also [23]).

Lemma 2.1

For positive integers $n_{1}$ and $n_{2}$, let $\varvec{\alpha }$, $\varvec{X}$, $\varvec{b}$, and $\varvec{Y}$ be $1 \times n_{1}$, $n_{1} \times n_{1}$, $n_{1} \times 1$, and $n_{2} \times n_{2}$ matrices, respectively. For $x \ge 0$, the convolution of matrix exponentials

$$\begin{aligned} \int _{0}^{x} \varvec{\alpha } \exp (\varvec{X} u) \varvec{b} \exp (\varvec{Y}(x-u)) \mathrm{d}u \end{aligned}$$

(2.6)

is computed as follows: Let $\varvec{M}_{11} := \varvec{I}_{n_{2}} \otimes \varvec{X}$, $\varvec{M}_{12} := \varvec{I}_{n_{2}} \otimes \varvec{b}$, and $\varvec{M}_{22} := \varvec{Y}$, where $\otimes$ is the Kronecker product. Then (2.6) is given by

$$\begin{aligned} \int _{0}^{x} \varvec{\alpha } \exp (\varvec{X} u) \varvec{b} \exp (\varvec{Y}(x-u)) \mathrm{d}u = (\varvec{I}_{n_{2}} \otimes \varvec{\alpha }) [\exp (\varvec{M}x)]^{(n_{1}n_{2},n_{2})}, \end{aligned}$$

(2.7)

where

$$\begin{aligned} \varvec{M} := \left( \begin{array}{cc} \varvec{M}_{11} &{}\quad \varvec{M}_{12} \\ \varvec{O} &{}\quad \varvec{M}_{22} \end{array} \right) . \end{aligned}$$

(2.8)

By choosing $\varvec{\alpha } := 1$, $\varvec{X} := 0$, and $\varvec{b} := 1$ in Lemma 2.1, we obtain the next result.

Corollary 2.1

For an $n_{2} \times n_{2}$ matrix $\varvec{Y}$, we have

$$\begin{aligned} \int _{0}^{x} \exp (\varvec{Y}u) \mathrm{d}u = \left[ \exp \left( \left( \begin{array}{cc} \varvec{O} &{}\quad \varvec{I}_{n_{2}} \\ \varvec{O} &{}\quad \varvec{Y} \end{array} \right) x\right) \right] ^{(n_{2},n_{2})}. \end{aligned}$$

(2.9)

In the subsequent subsections, we present a computational procedure for the stationary density of the workload.

2.2 Stationary density in mode 0

In this subsection we determine the steady-state density of the workload during the times in which the server is in mode 0. Such an inactive period lasts from the instant the system becomes empty until the next instant in which the workload exceeds a certain threshold level $d_1$. Assume that $Z_{0} = 0$, i.e., there is no workload in the system and the server is in mode 0 at time 0. For $i \in \mathbb {N}_{0} := \{0,1,2,\ldots \}$, let $\theta _{i}$ denote the i-th arrival time after time 0, where $\theta _{0} := 0$. For $x \ge 0$, let

$$\begin{aligned} \nu (x) := \sup \{ i \in \mathbb {N}_{0} ; Z_{\theta _{i}} \le x \}, \end{aligned}$$

(2.10)

which is the number of customer arrivals until the workload becomes larger than x. We note that $m(x) := \mathsf {E}[\nu (x)]$ is the renewal function with the renewal interval distribution B(x) ($x \ge 0$) in (2.1). From Theorem 3.1.2 in [14], we have, for $0 \le x < d_1$,

$$\begin{aligned} m(x)&= \int _{0}^{x} \varvec{\tau } \exp (\varvec{D} u) \varvec{t} \mathrm{d}u = \varvec{\tau } \left[ \exp \left( \left( \begin{array}{cc} \varvec{O} &{}\quad \varvec{I}_{N} \\ \varvec{O} &{}\quad \varvec{D} \end{array} \right) x\right) \right] ^{(N,N)}\varvec{t},\end{aligned}$$

(2.11)

$$\begin{aligned} m'(x)&= \varvec{\tau } \exp (\varvec{D} x) \varvec{t}, \end{aligned}$$

(2.12)

where $\varvec{t} := -\varvec{T}\varvec{1}$ and $\varvec{D} := \varvec{T} + \varvec{t}\varvec{\tau }$, and the second equation in (2.11) follows from Corollary 2.1.

Remark 2.1

By Theorem 3.1.4 in [14], the renewal function and its derivative are given in more explicit formulas as follows: for $x \ge 0$,

$$\begin{aligned} m(x)&= \eta ^{-1} x + \varvec{\tau } (\exp (\varvec{D}x) - \varvec{I}) (\varvec{D} - \varvec{\Pi })^{-1} \varvec{t},\end{aligned}$$

(2.13)

$$\begin{aligned} m'(x)&= \eta ^{-1} + \varvec{\tau } \varvec{D} \exp (\varvec{D} x) (\varvec{D} - \varvec{\Pi })^{-1} \varvec{t}, \end{aligned}$$

(2.14)

where $\eta := \varvec{\tau } (-\varvec{T})^{-1} \varvec{1}$, $\varvec{\Pi } := \varvec{1}\varvec{\pi }$ and $\varvec{\pi }$ is the stationary distribution of $\varvec{D}$, i.e., $\varvec{\pi }$ is the nonnegative solution of $\varvec{\pi } \varvec{D} = \varvec{0}$ and $\varvec{\pi } \varvec{1} = 1$. In this paper, we need to compute some convolutions, including $m'(x)$. Therefore, we use (2.11) and (2.12) to obtain our formulas in a simpler form.

We are ready to consider the stationary density of the workload when the server is in mode 0. Let $m_{i}$ ($i=0,1$) be the mean length of a mode-i interval. From Subsection 2.1 of [22] we have

$$\begin{aligned} m_{0} = \lambda ^{-1} (1 + m(d_1)). \end{aligned}$$

(2.15)

Let $v(x | S = 0)$, where $x > 0$, denote the conditional stationary density of the workload when the server is in mode 0, which is given by (cf. Subsection 2.1 of [22])

$$\begin{aligned} v(x | S = 0) = {\left\{ \begin{array}{ll} \frac{ \lambda ^{-1} m'(x) }{ m_{0} }, &{} 0 < x \le d_1,\\ 0 , &{} x > d_1. \end{array}\right. } \end{aligned}$$

(2.16)

From (2.15) and (2.16), we have

$$\begin{aligned} v(x | S = 0) = {\left\{ \begin{array}{ll} \frac{m'(x)}{1 + m(d_1)}, &{} 0 < x \le d_1,\\ 0, &{} x > d_1. \end{array}\right. } \end{aligned}$$

(2.17)

For $i=0,1$, let $p_{i} := \Pr (S = i)$, which is the steady-state probability that the server is in mode i, and it is readily seen that $p_{i} = m_{i}/(m_{0} + m_{1})$, i.e.,

$$\begin{aligned} p_{0} = \frac{1+m(d_1)}{1+m(d_1) + \lambda m_{1}}, \quad p_{1} = \frac{\lambda m_{1}}{1+m(d_1) + \lambda m_{1}}, \end{aligned}$$

(2.18)

where $m_{1}$ later will be determined by using a normalizing condition. Let $v_{0}(x)$, $0 < x \le d_1$, denote the unconditional stationary density of the workload when the server is in mode 0, i.e., $v_{0}(x) = \frac{\mathrm{d}}{\mathrm{d}x} \Pr (Z \le x, S = 0)$. Then we have

$$\begin{aligned} v_{0}(x) = {\left\{ \begin{array}{ll} p_{0} \cdot v(x | S = 0), &{} 0 < x \le d_1,\\ 0, &{} x > d_1. \end{array}\right. } \end{aligned}$$

(2.19)

Furthermore, let $V(0) := \Pr (Z = 0, S = 0)$, which is the marginal probability that there are no customers in the system, and which is given by $V(0) = \lambda ^{-1} / (m_{0} + m_{1})$, i.e.,

$$\begin{aligned} V(0) = \frac{1}{1 + m(d_1) + \lambda m_{1}}. \end{aligned}$$

(2.20)

Combining (2.17), (2.18) and (2.19), the results in this subsection can be summarized as follows.

Lemma 2.2

The stationary density of the workload when the server is in mode 0 equals

$$\begin{aligned} v_{0}(x) ={\left\{ \begin{array}{ll} \frac{m'(x)}{1 + m(d_1) + \lambda m_{1}}, &{} 0 < x \le d_1,\\ 0, &{} x > d_1, \end{array}\right. } \end{aligned}$$

(2.21)

where m(x) is given by (2.11) and $m'(x)$ by (2.12); the probability that there is no customer in the system is given by (2.20).

In the next two subsections we determine the stationary density in mode 1, distinguishing between the cases where the workload is above $d_1$ (Sect. 2.3) and below $d_1$ (Sect. 2.4).

2.3 Stationary density in mode 1 when the workload is above $d_{1}$

Let $v_{1}(x)$, where $x > 0$, denote the stationary density of the workload when the server is in mode 1, i.e., $v_{1}(x) = \frac{\mathrm{d}}{\mathrm{d}x} \Pr (Z \le x, S = 1)$. We use the level crossing technique (see, for example [3]), which states that, in steady state, each workload level is crossed just as often from above and from below. Hence the stationary density of the workload, denoted by $v(x) := v_{0}(x) + v_{1}(x)$, where $x > 0$, satisfies the following relations:

$$\begin{aligned} r_{1} (v(x) - v_{0}(x))&= \lambda V(0) \overline{B}(x) + \lambda \int _{0}^{x} \overline{B}(x-y) v(y) \mathrm{d}y, \quad 0 < x \le d_{1}, \end{aligned}$$

(2.22)

$$\begin{aligned} r_{2} v(x)&= \lambda V(0) \overline{B}(x) + \lambda \int _{0}^{x} \overline{B}(x-y) v(y) \mathrm{d}y, \quad d_{1}< x < \infty , \end{aligned}$$

(2.23)

where $\overline{B}(x) = \varvec{\tau } \exp (\varvec{T} x) \varvec{1}$ for $x \ge 0$. We first solve (2.23), and show that the solution is given by a matrix exponential form. To this end, for $x > d_{1}$, we consider a $(1 \times N)$ row vector, which is denoted by $\varvec{v}(x)$, and satisfies the following equation:

$$\begin{aligned} r_{2} \varvec{v}(x) = \lambda V(0) \varvec{\tau } \exp (\varvec{T} x) + \lambda \int _{0}^{x} \varvec{v}(y) \varvec{1} \varvec{\tau } \exp (\varvec{T}(x-y)) \mathrm{d}y, \quad d_{1}< x < \infty . \end{aligned}$$

(2.24)

We note that $v_{1}(x) = v(x) = \varvec{v}(x) \varvec{1}$ for $x > d_{1}$ from (2.23) and (2.24). Similarly to [24], by taking the derivative of (2.24) with respect to x, we have

$$\begin{aligned} r_{2} \varvec{v}'(x)&= \lambda V(0) \varvec{\tau } \exp (\varvec{T} x) \varvec{T} + \lambda \varvec{v}(x) \varvec{1}\varvec{\tau } + \lambda \int _{0}^{x} \varvec{v}(y) \varvec{1} \varvec{\tau } \exp (\varvec{T}(x-y)) \mathrm{d}y \cdot \varvec{T} \nonumber \\&= r_{2} \varvec{v}(x) (\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}), \quad d_{1}< x < \infty , \end{aligned}$$

(2.25)

where the last equality follows by applying (2.24) to the last term of the first line of (2.25). The solution of (2.25) is given by the following matrix exponential form:

$$\begin{aligned} \varvec{v}(x) = \tilde{\varvec{u}} \exp ( (\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) (x - d_{1}) ), \quad d_{1}< x < \infty , \end{aligned}$$

(2.26)

where $\tilde{\varvec{u}}$ is a $(1 \times N)$ row vector, which will be determined later (see (2.40)).

The results in this subsection are summarized as follows.

Lemma 2.3

For $x > d_{1}$, the stationary density of the workload when the server is in mode 1 is given by $v_{1}(x) = \varvec{v}(x) \varvec{1}$, where

$$\begin{aligned} \varvec{v}(x) = \tilde{\varvec{u}} \exp ( (\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) (x - d_{1}) ), \quad x > d_{1}, \end{aligned}$$

(2.27)

where $\tilde{\varvec{u}}$ is given by (2.40).

2.4 Stationary density in mode 1 when the workload is below $d_{1}$

We next solve (2.22). Since $v(x) = v_{0}(x) + v_{1}(x)$, (2.22) is rewritten into

$$\begin{aligned} r_{1} v_{1}(x) = \lambda V(0) \overline{B}(x) + \lambda \int _{0}^{x} (v_{0}(y) + v_{1}(y)) \overline{B}(x-y) \mathrm{d}y, \quad 0 < x \le d_{1}. \end{aligned}$$

(2.28)

To solve (2.28), we consider $(1 \times N)$ row vectors $\varvec{v}_{0}(x)$ and $\varvec{v}_{1}(x)$ satisfying the following equation: for $0 < x \le d_{1}$,

$$\begin{aligned} r_{1} \varvec{v}_{1}(x) = \lambda V(0) \varvec{\tau } \exp (\varvec{T} x) + \lambda \int _{0}^{x} (\varvec{v}_{0}(y) + \varvec{v}_{1}(y)) \varvec{1} \varvec{\tau } \exp (\varvec{T}(x-y)) \mathrm{d}y. \end{aligned}$$

(2.29)

We note that $v_{0}(x) = \varvec{v}_{0}(x) \varvec{1}$ and $v_{1}(x) = \varvec{v}_{1}(x) \varvec{1}$ for $0 < x \le d_{1}$. Similar to (2.25), by taking the derivative of (2.29) with respect to x, and applying (2.29) to the derivative, we have

$$\begin{aligned} \varvec{v}_{1}'(x) = \varvec{v}_{1}(x) (\lambda r_{1}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) + \lambda r_{1}^{-1} v_{0}(x) \varvec{\tau }, \quad 0 < x \le d_{1}. \end{aligned}$$

(2.30)

Suppose that the solution of (2.30) is given by the following form:

$$\begin{aligned} \varvec{v}_{1}(x) = \varvec{w}(x) \exp ( (\lambda r_{1}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) x ), \quad 0 < x \le d_{1}, \end{aligned}$$

(2.31)

where $\varvec{w}(x)$ is a $(1 \times N)$ row vector and is determined as follows: By taking the derivative of (2.31) with respect to x, we have

$$\begin{aligned} \varvec{v}_{1}'(x) = \varvec{v}_{1}(x) (\lambda r_{1}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) + \varvec{w}'(x) \exp ( (\lambda r_{1}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) x ), \quad 0 < x \le d_{1}. \end{aligned}$$

(2.32)

By comparing (2.30) with (2.32), we have

$$\begin{aligned} \varvec{w}'(x) = \lambda r_{1}^{-1} v_{0}(x) \varvec{\tau } \exp (-(\lambda r_{1}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) x), \quad 0 < x \le d_{1}. \end{aligned}$$

(2.33)

We then have

$$\begin{aligned} \varvec{w}(x) = \int _{0}^{x}\lambda r_{1}^{-1} v_{0}(u) \varvec{\tau } \exp (-(\lambda r_{1}^{-1} \varvec{1} \varvec{\tau } + \varvec{T})u) \mathrm{d}u + \varvec{c}, \quad 0 < x \le d_{1}, \end{aligned}$$

(2.34)

where $\varvec{c}$ is a $(1 \times N)$ row vector and is determined as follows (see (2.37) below): By taking the limit $x \downarrow 0$ of (2.29), we have

$$\begin{aligned} r_{1} \varvec{v}_{1}(0+) = \lambda V(0) \varvec{\tau } = \frac{\lambda \varvec{\tau }}{1 + m(d_1) + \lambda m_{1}}, \end{aligned}$$

(2.35)

where the last equality follows from (2.20). From (2.31) and (2.34), we have

$$\begin{aligned} \varvec{v}_{1}(0+) = \varvec{c}. \end{aligned}$$

(2.36)

From (2.35) and (2.36), we obtain

$$\begin{aligned} \varvec{c}&= \frac{ \lambda r_{1}^{-1} \varvec{\tau } }{ 1 + m(d_1) + \lambda m_{1} }. \end{aligned}$$

(2.37)

We are now ready to formulate the following lemma.

Lemma 2.4

For $0 < x \le d_{1}$, the stationary density $v_{1}(x)$ when the server is in mode 1 is given by $v_{1}(x) = \varvec{v}_{1}(x) \varvec{1}$, where

$$\begin{aligned} \varvec{v}_{1}(x) = \frac{\lambda r_{1}^{-1} \varvec{\tau }}{1 + m(d_1) + \lambda m_{1}} \bigg \{ (\varvec{I}_{N} \otimes \varvec{\tau }) \left[ \exp \left( \overline{\varvec{M}} x \right) \right] ^{(N^{2},N)} + \exp ((\lambda r_{1}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) x) \bigg \}, \end{aligned}$$

(2.38)

with

$$\begin{aligned} \overline{\varvec{M}} := \left( \begin{array}{cc} \varvec{I}_{N} \otimes \varvec{D} &{}\quad \varvec{I}_{N} \otimes \varvec{t} \\ \varvec{O} &{}\quad \lambda r_{1}^{-1} \varvec{1} \varvec{\tau } + \varvec{T} \end{array} \right) . \end{aligned}$$

(2.39)

Furthermore,

$$\begin{aligned} \tilde{\varvec{u}} = \frac{ \lambda r_{2}^{-1} \varvec{\tau } }{ 1 + m(d_1) + \lambda m_{1}} \bigg \{ (\varvec{I}_{N} \otimes \varvec{\tau }) \left[ \exp \left( \overline{\varvec{M}} d_{1} \right) \right] ^{(N^{2},N)} + \exp ((\lambda r_{1}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) d_{1}) \bigg \}. \end{aligned}$$

(2.40)

Proof

From (2.21), (2.31), and (2.34) (see also (2.37)), $\varvec{v}_{1}(x)$ for $0 < x \le d_{1}$ is given by the sum of the following two terms:

$$\begin{aligned} \frac{\lambda r_{1}^{-1} \varvec{\tau }}{1 + m(d_1) + \lambda m_{1}} \int _{0}^{x} \varvec{\tau } \exp (\varvec{D} u) \varvec{t} \exp ((\lambda r_{1}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) (x - u)) \mathrm{d}u \end{aligned}$$

(2.41)

and

$$\begin{aligned} \frac{\lambda r_{1}^{-1} \varvec{\tau }}{1 + m(d_1) + \lambda m_{1}} \exp ((\lambda r_{1}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) x). \end{aligned}$$

(2.42)

From Lemma 2.1, the integral in (2.41) is given by

$$\begin{aligned} (\varvec{I}_{N} \otimes \varvec{\tau }) \left[ \exp \left( \overline{\varvec{M}} x \right) \right] ^{(N^{2},N)}. \end{aligned}$$

(2.43)

From (2.41), (2.42), and (2.43), we obtain (2.38). Using (2.26) and (2.38), and observing that $\varvec{v}_{1}(d_{1}+) = r_{1} r_{2}^{-1} \varvec{v}_{1}(d_{1}-)$ from (2.24) and (2.29), we obtain (2.40). $\square$

2.5 Computation of the mean active period

In this subsection, we compute the mean active period $m_{1}$ from the normalizing condition, i.e.,

$$\begin{aligned} V(0) + \int _{0}^{d_1} v_{0}(x) \mathrm{d}x + \int _{0}^{d_{1}} v_{1}(x) \mathrm{d}x + \int _{d_{1}}^{\infty } v_{1}(x) \mathrm{d}x = 1. \end{aligned}$$

(2.44)

The three integrals are successively obtained in Lemmas 2.5, 2.6 and 2.7. We first use Lemma 2.2 to obtain $\int _0^{d_1} v_0(x) \mathrm{d}x$.

Lemma 2.5

$$\begin{aligned} \int _{0}^{d_1} v_{0}(x) \mathrm{d}x&= \frac{m(d_1)}{1 + m(d_1) + \lambda m_{1}}, \end{aligned}$$

(2.45)

where

$$\begin{aligned} m(d_1) = \varvec{\tau } \left[ \exp \left( \left( \begin{array}{cc} \varvec{O} &{}\quad \varvec{I}_{N} \\ \varvec{O} &{}\quad \varvec{D} \end{array} \right) d_1\right) \right] ^{(N,N)}\varvec{t}. \end{aligned}$$

(2.46)

We next consider $\int _{0}^{d_{1}} v_{1}(x) \mathrm{d}x$. From Corollary 2.1, we have

$$\begin{aligned} \int _{0}^{d_{1}} \left[ \exp \left( \overline{\varvec{M}} x \right) \right] ^{(N^{2},N)} \mathrm{d}x&= \left[ \int _{0}^{d_{1}} \exp \left( \overline{\varvec{M}} x \right) \mathrm{d}x \right] ^{(N^{2},N)}\nonumber \\&= \left[ \left[ \exp \left( \left( \begin{array}{cc} \varvec{O} &{}\quad \varvec{I}_{N^{2}+N} \\ \varvec{O} &{}\quad \overline{\varvec{M}} \end{array} \right) d_{1} \right) \right] ^{(N^{2}+N,N^{2}+N)} \right] ^{(N^{2},N)}\nonumber \\&= \left[ \exp \left( \left( \begin{array}{cc} \varvec{O} &{}\quad \varvec{I}_{N^{2}+N} \\ \varvec{O} &{}\quad \overline{\varvec{M}} \end{array} \right) d_{1} \right) \right] ^{(N^{2},N)}, \end{aligned}$$

(2.47)

and

$$\begin{aligned} \int _{0}^{d_{1}} \exp ((\lambda r_{1}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) x) \mathrm{d}x = \left[ \exp \left( \left( \begin{array}{cc} \varvec{O} &{}\quad \varvec{I}_{N} \\ \varvec{O} &{}\quad \lambda r_{1}^{-1} \varvec{1} \varvec{\tau } + \varvec{T} \end{array} \right) d_{1} \right) \right] ^{(N,N)}. \end{aligned}$$

(2.48)

The next lemma immediately follows from Lemma 2.4, (2.47) and (2.48).

Lemma 2.6

The second integral in the left-hand side of (2.44) is given by

$$\begin{aligned} \int _{0}^{d_{1}} v_{1}(x) \mathrm{d}x&= \frac{\lambda r_{1}^{-1} \varvec{\tau }}{1 + m(d_1) + \lambda m_{1}} \bigg \{ (\varvec{I}_{N} \otimes \varvec{\tau }) \left[ \exp \left( \left( \begin{array}{cc} \varvec{O} &{}\quad \varvec{I}_{N^{2}+N} \\ \varvec{O} &{}\quad \overline{\varvec{M}} \end{array} \right) d_{1} \right) \right] ^{(N^{2},N)}\nonumber \\&\quad + \left[ \exp \left( \left( \begin{array}{cc} \varvec{O} &{}\quad \varvec{I}_{N} \\ \varvec{O} &{}\quad \lambda r_{1}^{-1} \varvec{1} \varvec{\tau } + \varvec{T} \end{array} \right) d_{1} \right) \right] ^{(N,N)} \bigg \} \varvec{1}, \end{aligned}$$

(2.49)

where $\overline{\varvec{M}}$ is given by (2.39).

We finally compute $\int _{d_{1}}^{\infty } v_{1}(x) \mathrm{d}x$ as follows.

Lemma 2.7

From Lemma 2.3 and the stability condition (2.4), we have

$$\begin{aligned} \int _{d_{1}}^{\infty } v_{1}(x) \mathrm{d}x = \tilde{\varvec{u}} \left( - (\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) \right) ^{-1} \varvec{1}, \end{aligned}$$

(2.50)

where $\tilde{\varvec{u}}$ is given by (2.40).

Proof

We first show that $\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}$ is invertible, i.e., $(\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T})^{-1}$ exists. To this end, we consider the following equation:

$$\begin{aligned} \varvec{x} (\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) = \varvec{0}, \end{aligned}$$

(2.51)

where $\varvec{x}$ is a $(1 \times N)$ row vector. By post-multiplying the left-hand side of (2.51) by $r_{2} (-\varvec{T})^{-1} \varvec{1}$, we have

$$\begin{aligned} \varvec{x} \varvec{1} \{\lambda \varvec{\tau } (-\varvec{T})^{-1} \varvec{1} - r_{2}\} = 0, \end{aligned}$$

(2.52)

and from the stability condition (2.4), we have

$$\begin{aligned} \varvec{x} \varvec{1} = 0. \end{aligned}$$

(2.53)

By applying (2.53) to (2.51) and noting the existence of $(-\varvec{T})^{-1}$, we obtain $\varvec{x} = \varvec{0}$, which implies that $(\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T})^{-1}$ exists.

We next show (2.50). From Lemma 2.3 and the existence of $(\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T})^{-1}$ we have

$$\begin{aligned} \int _{d_{1}}^{\infty } v_{1}(x) \mathrm{d}x&= \tilde{\varvec{u}} \left\{ \lim _{y \rightarrow \infty } \int _{0}^{y} \exp ( (\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) x ) \mathrm{d}x \right\} \varvec{1}\nonumber \\&= \tilde{\varvec{u}} \left( - (\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) \right) ^{-1} \left\{ \varvec{I}_{N} - \lim _{y \rightarrow \infty } \exp ( (\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) y ) \right\} \varvec{1}. \end{aligned}$$

(2.54)

It remains to show that

$$\begin{aligned} \lim _{y \rightarrow \infty } \exp ( (\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) y ) = \varvec{O}. \end{aligned}$$

(2.55)

Noting that $\varvec{\pi }$ is the stationary distribution of $\varvec{T} + \varvec{t}\varvec{\tau }$, i.e., $\varvec{\pi } \varvec{1} = 1$ and $\varvec{\pi } \varvec{T} = -\varvec{\pi } \varvec{t} \varvec{\tau }$, we have

$$\begin{aligned} \varvec{\pi } (\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T})&= (\lambda r_{2}^{-1} - \varvec{\pi } \varvec{t}) \varvec{\tau }\nonumber \\&= \left( \lambda r_{2}^{-1} - \frac{1}{\varvec{\tau }(-\varvec{T})^{-1}\varvec{1}} \right) \varvec{\tau }, \end{aligned}$$

(2.56)

where the last equation follows from $\varvec{\pi } \varvec{T} = -\varvec{\pi } \varvec{t} \varvec{\tau }$, i.e., $1 = \varvec{\pi } \varvec{1} = \varvec{\pi } \varvec{t} \varvec{\tau } (-\varvec{T})^{-1} \varvec{1}$. Since $\varvec{\tau }$ is a nonnegative and nonzero vector, (2.56) and the stability condition (2.4) imply that $\varvec{\pi } (\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) \le \varvec{0}$ and $\varvec{\pi } (\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) \not = \varvec{0}$. From Theorem 1.6(b) in [25], the Perron–Frobenius eigenvalue of $\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}$ is negative, which implies (2.55). $\square$

From Lemmas 2.5–2.7, the mean active period is given as follows.

Theorem 2.1

The mean active period for the case of two processing speeds is given by

$$\begin{aligned} m_{1}&= r_{1}^{-1} \varvec{\tau } \left\{ (\varvec{I}_{N} \otimes \varvec{\tau }) \left[ \exp \left( \left( \begin{array}{cc} \varvec{O} &{}\quad \varvec{I}_{N^{2}+N} \\ \varvec{O} &{}\quad \overline{\varvec{M}} \end{array} \right) d_{1} \right) \right] ^{(N^{2},N)} + \left[ \exp \left( \left( \begin{array}{cc} \varvec{O} &{}\quad \varvec{I}_{N} \\ \varvec{O} &{}\quad \lambda r_{1}^{-1} \varvec{1} \varvec{\tau } + \varvec{T} \end{array} \right) d_{1} \right) \right] ^{(N,N)} \right\} \varvec{1} \nonumber \\&\quad + r_{2}^{-1} \varvec{\tau } \left\{ (\varvec{I}_{N} \otimes \varvec{\tau }) \left[ \exp \left( \overline{\varvec{M}} d_{1} \right) \right] ^{(N^{2},N)} + \exp ((\lambda r_{1}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) d_{1}) \right\} \left( - (\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) \right) ^{-1} \varvec{1}. \end{aligned}$$

(2.57)

3 Extension to the case of multiple processing speeds

In this section, we consider an extension of the model of the previous section to the case of an arbitrary number of different (constant) processing speeds. Section 3.1 provides a detailed model description. In Sect. 3.2 we consider the stationarity workload density when the server is inactive (briefly, as this is very similar to the result in Sect. 2.2). Sections 3.3 and 3.4 are successively devoted to the stationary workload density when the server is active and the workload is above, respectively below, the threshold level $d_K$.

3.1 Model description

In this section we extend the assumptions (i) and (ii), and the service processing speed $s_{i}(x)$ ($i = 0,1$ and $x \ge 0$) in Sect. 2 as follows: Let $\{ d_{k} ; k = 0,1,2,\ldots , K_0\}$ be an increasing sequence such that $d_{0} := 0$.

• (i)’ When the workload $Z_{t}$ hits $d_{0} ( = 0)$ at time t, the server enters mode 0, i.e., $S_{t} = 0$. After that, the server remains in mode 0 until the workload exceeds a threshold $d_{K} > 0$, where $K \le K_0$ is a positive integer, i.e., $S_{u} = 0$ as long as $Z_{u} < d_{K}$ for $u > t$.

• (ii)’ When the workload $Z_{t}$ exceeds $d_{K}$ at time t, the server changes its mode from 0 to 1, i.e., $S_{t} = 1$. The server remains in that mode until the workload hits 0 again, i.e., $S_{u} = 1$ as long as $Z_{u} > 0$ for $u > t$.

The processing speed of the server depends on both the server’s mode and workload in the following way: Let $s_{i}(x)$, where $i = 0,1$ and $x \ge 0$, denote the processing speed of the server when $(Z_{t}, S_{t}) = (x,i)$; it is defined as follows:

$$\begin{aligned} s_{0}(x)&= 0, \quad 0 \le x \le d_{K},\end{aligned}$$

(3.1)

$$\begin{aligned} s_{1}(x)&= r_{k}, \quad x \in J_{k} := (d_{k-1},d_{k}], \ k = 1,2,3,\ldots ,K_0 , \end{aligned}$$

(3.2)

where $\{ r_{k} ; k = 1,2,\ldots , K_0 \}$ is a positive valued sequence and $d_{K_0}= \infty$. Note that $\{ s_{0}(x) ; x > d_{K} \}$ and $s_{1}(0)$ need not be specified in view of the definition of modes 0 and 1. Also note that $d_K$ coincides with one of the switching levels $d_k$ of mode 1. However, this is no restriction, as one could take $r_K=r_{K+1}$. A sample path of the workload is presented in Fig. 2.

Fig. 2

The workload process for the case of multiple processing speeds

We assume that the Markov process $(Z_{t}, S_{t})$ is positive-recurrent, and denote its replica in steady state by (Z, S). Similar to (2.4), the stability condition of this Markov process is given by

$$\begin{aligned} \frac{ \lambda \varvec{\tau }(-\varvec{T})^{-1}\varvec{1} }{ r_{K_0} } < 1. \end{aligned}$$

(3.3)

In the subsequent subsections, we present a computational procedure for the stationary density of the workload.

3.2 Stationary density in mode 0

Similar to Lemma 2.2, the stationary density in mode 0, and the marginal probability that there is no customer in the system, are given as follows.

Lemma 3.1

The stationary density of the workload when the server is in mode 0 is given by

$$\begin{aligned} v_{0}(x) = {\left\{ \begin{array}{ll} \frac{m'(x)}{1 + m(d_{K}) + \lambda m_{1}}, &{} 0 < x \le d_{K},\\ 0, &{} x > d_{K}, \end{array}\right. } \end{aligned}$$

(3.4)

where m(x) is given by (2.11) and $m'(x)$ by (2.12), and $m_{1}$ is the mean active period which will be determined later (see (3.48)). The marginal probability that there is no customer in the system is given by

$$\begin{aligned} V(0) = \frac{1}{1 + m(d_{K}) + \lambda m_{1}}. \end{aligned}$$

(3.5)

In the next two subsections we determine the stationary density in mode 1, distinguishing between the cases that the workload is above $d_K$ (Sect. 3.3) and below $d_K$ (Sect. 3.4).

3.3 Stationary density in mode 1 when the workload is above $d_{K}$

Similar to (2.22) and (2.23), the stationary density of the workload satisfies the following relations:

$$\begin{aligned} r_{k} (v(x) - v_{0}(x))&= \lambda V(0) \overline{B}(x) + \lambda \int _{0}^{x} \overline{B}(x-y) v(y) \mathrm{d}y, \quad x \in J_{k}, \ 1 \le k \le K, \end{aligned}$$

(3.6)

$$\begin{aligned} r_{k} v(x)&= \lambda V(0) \overline{B}(x) + \lambda \int _{0}^{x} \overline{B}(x-y) v(y) \mathrm{d}y, \quad x \in J_{k}, \ K+1 \le k \le K_0. \end{aligned}$$

(3.7)

To find a solution of (3.7), we consider the following equation:

$$\begin{aligned} r_{k} \varvec{v}(x) = \lambda V(0) \varvec{\tau } \exp (\varvec{T} x) + \lambda \int _{0}^{x} \varvec{v}(y) \varvec{1} \varvec{\tau } \exp (\varvec{T}(x-y)) \mathrm{d}y, \quad x \in J_{k}, \ K+1 \le k \le K_0. \end{aligned}$$

(3.8)

We note that $v_{1}(x) = v(x) = \varvec{v}(x) \varvec{1}$ for $x > d_{K}$ by (3.7) and (3.8). Similar to (2.26), by taking the derivative of (3.8) with respect to x, and applying (3.8) to the derivative, we have

$$\begin{aligned} \varvec{v}'(x) = \varvec{v}(x) (\lambda r_{k}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}), \quad x \in J_{k}, \ K+1 \le k \le K_0. \end{aligned}$$

(3.9)

The solution of (3.9) is given by the following matrix exponential form:

$$\begin{aligned} \varvec{v}(x) = \tilde{\varvec{u}}_{k} \varvec{U}_{k}(x), \quad x \in J_{k}, \ K+1 \le k \le K_0, \end{aligned}$$

(3.10)

where $\{ \tilde{\varvec{u}}_{k} ; K+1 \le k \le K_0 \}$ is a set of $(1 \times N)$ row vectors, which will be determined later (see (3.16)), and

$$\begin{aligned} \varvec{U}_{k}(x) := \exp ( (\lambda r_{k}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) (x - d_{k-1}) ), \quad x \in J_{k}, \ 1 \le k \le K_0. \end{aligned}$$

(3.11)

We note that

$$\begin{aligned} \lim _{x \downarrow d_{k-1}} \varvec{U}_{k}(x) = \varvec{I}_{N}, \quad 1 \le k \le K_0. \end{aligned}$$

(3.12)

The sequence $\{ \tilde{\varvec{u}}_{k} ; K+1 \le k \le K_0 \}$ is recursively determined as follows: By taking the two limits $x \uparrow d_{k}$ and $x \downarrow d_{k}$ of (3.8), we have

$$\begin{aligned} r_{k} \varvec{v}(d_{k}-) = r_{k+1} \varvec{v}(d_{k}+), \quad K+1 \le k \le K_0 -1. \end{aligned}$$

(3.13)

From (3.10) and (3.12), we have, for $K+1 \le k \le K_0 -1$,

$$\begin{aligned} \varvec{v}(d_{k}-) = \tilde{\varvec{u}}_{k} \varvec{U}_{k}(d_{k}), \quad \varvec{v}(d_{k}+) = \tilde{\varvec{u}}_{k+1}. \end{aligned}$$

(3.14)

From (3.13) and (3.14), we have

$$\begin{aligned} \tilde{\varvec{u}}_{k+1} = \frac{r_{k}}{r_{k+1}} \tilde{\varvec{u}}_{k} \varvec{U}_{k}(d_{k}), \quad K+1 \le k \le K_0 -1, \end{aligned}$$

(3.15)

which yields

$$\begin{aligned} \tilde{\varvec{u}}_{k} = \frac{r_{K+1}}{r_{k}} \tilde{\varvec{u}}_{K+1} \left\{ \Pi _{i=K+1}^{k-1} \varvec{U}_{i}(d_{i}) \right\} , \quad K+1 \le k \le K_0, \end{aligned}$$

(3.16)

where $\tilde{\varvec{u}}_{K+1}$ will be determined later (see (3.39)).

In what follows, we summarize the results in this subsection. To this end, we introduce the following notation:

$$\begin{aligned} \hat{\varvec{U}}_{k,l} := \Pi _{i=k}^{l} \varvec{U}_{i}(d_{i}), \quad 1 \le k, l \le K_0, \end{aligned}$$

(3.17)

where the empty product is an identity matrix, i.e., $\hat{\varvec{U}}_{k,l} = \varvec{I}_{N}$ for $k>l$.

Lemma 3.2

For $x > d_{K}$, the stationary density of the workload when the server is in mode 1 is given by $v_{1}(x) = \varvec{v}(x) \varvec{1}$, where

$$\begin{aligned} \varvec{v}(x)&= \frac{r_{K+1}}{r_{k}} \tilde{\varvec{u}}_{K+1} \hat{\varvec{U}}_{K+1,k-1} \varvec{U}_{k}(x),\end{aligned}$$

(3.18)

$$\begin{aligned} \varvec{U}_{k}(x)&= \exp ( (\lambda r_{k}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) (x - d_{k-1}) ) , \end{aligned}$$

(3.19)

for $x \in J_{k}$ and $K+1 \le k \le K_0$, and $\tilde{\varvec{u}}_{K+1}$ is given by (3.39). From the stability condition (3.3), we have $\varvec{U}_{K_{0}}(d_{K_{0}}) = \varvec{O}$.

3.4 Stationary density in mode 1 when the workload is below $d_{K}$

To solve (3.6), we consider $(1 \times N)$ row vectors $\varvec{v}_{0}(x)$ and $\varvec{v}_{1}(x)$ satisfying the following equation: for $x \in J_{k}, \ 1 \le k \le K$,

$$\begin{aligned} r_{k} \varvec{v}_{1}(x) = \lambda V(0) \varvec{\tau } \exp (\varvec{T} x) + \lambda \int _{0}^{x} (\varvec{v}_{0}(y) + \varvec{v}_{1}(y)) \varvec{1} \varvec{\tau } \exp (\varvec{T}(x-y)) \mathrm{d}y. \end{aligned}$$

(3.20)

We note that $v_{1}(x) = \varvec{v}_{1}(x) \varvec{1}$ for $0 < x \le d_{K}$. Similar to (2.30), by taking the derivative of (3.20) with respect to x, we have

$$\begin{aligned} \varvec{v}_{1}'(x) = \varvec{v}_{1}(x) (\lambda r_{k}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) + \lambda r_{k}^{-1} v_{0}(x) \varvec{\tau }, \quad x \in J_{k}, \ 1 \le k \le K. \end{aligned}$$

(3.21)

Suppose that the solution of (3.21) is given by the following form:

$$\begin{aligned} \varvec{v}_{1}(x) = \varvec{w}_{k}(x) \varvec{U}_{k}(x), \quad x \in J_{k}, \ 1 \le k \le K, \end{aligned}$$

(3.22)

where $\varvec{U}_{k}(x)$ is given by (3.11), and $\{ \varvec{w}_{k}(x) ; 1 \le k \le K \}$ is a set of $(1 \times N)$ row vectors. By taking the derivative in (3.22) and comparing with (3.21), we have

$$\begin{aligned} \varvec{w}_{k}'(x) = \lambda r_{k}^{-1} v_{0}(x) \varvec{\tau } \exp (-(\lambda r_{k}^{-1} \varvec{1} \varvec{\tau } + \varvec{T})(x - d_{k-1})), \quad x \in J_{k}, \ 1 \le k \le K. \end{aligned}$$

(3.23)

From (3.23), we have for $x \in J_{k}, \ 1 \le k \le K$,

$$\begin{aligned} \varvec{w}_{k}(x)&= \varvec{g}_{k}(x) + \varvec{c}_{k},\end{aligned}$$

(3.24)

$$\begin{aligned} \mathrm{with} ~ \varvec{g}_{k}(x)&:= \int _{d_{k-1}}^{x} \lambda r_{k}^{-1} v_{0}(u) \varvec{\tau } \exp (-(\lambda r_{k}^{-1} \varvec{1} \varvec{\tau } + \varvec{T})(u - d_{k-1})) \mathrm{d}u , \end{aligned}$$

(3.25)

where $\{ \varvec{c}_{k} ; 1 \le k \le K \}$ is a set of $(1 \times N)$ row vectors, which will be determined below (see (3.30) and (3.31)). We note that

$$\begin{aligned} \varvec{g}_{k}(d_{k-1}) = \varvec{0}, \quad \ 1 \le k \le K. \end{aligned}$$

(3.26)

In what follows, we summarize the results in this subsection. For simplicity of the exposition, let

$$\begin{aligned} \varvec{W}_{k}(x) := (\varvec{I}_{N} \otimes \varvec{\tau }_{k}) \left[ \exp \left( \overline{\varvec{M}}_{k} (x - d_{k-1}) \right) \right] ^{(N^{2},N)}, \quad x \in J_{k}, \ 1 \le k \le K, \end{aligned}$$

(3.27)

where

$$\begin{aligned} \varvec{\tau }_{k} := \varvec{\tau } \exp (\varvec{D} d_{k-1}), \quad \overline{\varvec{M}}_{k} := \left( \begin{array}{cc} \varvec{I}_{N} \otimes \varvec{D} &{}\quad \varvec{I}_{N} \otimes \varvec{t} \\ \varvec{O} &{}\quad \lambda r_{k}^{-1} \varvec{1} \varvec{\tau } + \varvec{T} \end{array} \right) . \end{aligned}$$

(3.28)

Lemma 3.3

For $0 < x \le d_{K}$, the stationary density $v_{1}(x)$ when the server is in mode 1 is given by $v_{1}(x) = \varvec{v}_{1}(x) \varvec{1}$, where, for $x \in J_{k}, \ 1 \le k \le K$,

$$\begin{aligned} \varvec{v}_{1}(x) = \frac{\lambda r_{k}^{-1} \varvec{\tau }}{1 + m(d_{K}) + \lambda m_{1}} \varvec{W}_{k}(x) + \varvec{c}_{k} \varvec{U}_{k}(x), \end{aligned}$$

(3.29)

where $\{ \varvec{c}_{k} ; 1 \le k \le K \}$ is given by

$$\begin{aligned} \varvec{c}_{k} = \frac{ \lambda r_{k}^{-1} \varvec{\tau } }{ 1 + m(d_{K}) + \lambda m_{1} } \left\{ \sum _{i=1}^{\min \{k-1,K\}} \varvec{W}_{i}(d_{i}) \hat{\varvec{U}}_{i+1,k-1} + \hat{\varvec{U}}_{1,k-1} \right\} . \end{aligned}$$

(3.30)

In particular, $\varvec{c}_{1}$ is given by

$$\begin{aligned} \varvec{c}_{1} = \frac{ \lambda r_{1}^{-1} \varvec{\tau } }{ 1 + m(d_{K}) + \lambda m_{1} }. \end{aligned}$$

(3.31)

Proof

From (3.22) and (3.24), we have

$$\begin{aligned} \varvec{v}_{1}(x) = \varvec{g}_{k}(x) \varvec{U}_{k}(x) + \varvec{c}_{k} \varvec{U}_{k}(x) , \quad x \in J_{k}, \ 1 \le k \le K. \end{aligned}$$

(3.32)

From (3.4), (3.11), and (3.25), the first term in the right-hand side of (3.32) is calculated as follows:

$$\begin{aligned} \varvec{g}_{k}(x) \varvec{U}_{k}(x)&= \frac{\lambda r_{k}^{-1} \varvec{\tau }}{1 + m(d_{K}) + \lambda m_{1}} \int _{0}^{x - d_{k-1}} \varvec{\tau } \exp (\varvec{D} d_{k-1}) \exp (\varvec{D}u) \varvec{t} \nonumber \\&\qquad \times \exp ((\lambda r_{k}^{-1} \varvec{1} \varvec{\tau } + \varvec{T})(x-d_{k-1}-u)) \mathrm{d}u\nonumber \\&= \frac{\lambda r_{k}^{-1} \varvec{\tau }}{1 + m(d_{K}) + \lambda m_{1}} \varvec{W}_{k}(x), \end{aligned}$$

(3.33)

where the last equation follows from Lemma 2.1, and then (3.29) is obtained.

It remains to show that $\{ \varvec{c}_{k} ; 1 \le k \le K \}$ is given by (3.30). By taking the limit $x \downarrow 0$ of (3.20) and applying (3.5) to it, we have

$$\begin{aligned} \varvec{v}_{1}(0+) = \frac{ \lambda r_{1}^{-1} \varvec{\tau } }{ 1 + m(d_{K}) + \lambda m_{1} }. \end{aligned}$$

(3.34)

From (3.11) and (3.27), we have $\varvec{U}_{1}(0+) = \varvec{I}_{N}$ and $\varvec{W}_{1}(0+) = \varvec{O}$, respectively, which imply $\varvec{v}_{1}(0+) = \varvec{c}_{1}$ from (3.34). Then we have (3.31).

Similar to (3.13), by taking the two limits $x \uparrow d_{k}$ and $x \downarrow d_{k}$ of (3.20) for $1 \le k \le K-1$, we have

$$\begin{aligned} r_{k} \varvec{v}_{1}(d_{k}-) = r_{k+1} \varvec{v}_{1}(d_{k}+), \quad \ 1 \le k \le K-1. \end{aligned}$$

(3.35)

For $1 \le k \le K-1$, by taking limits $x \uparrow d_{k}$ and $x \downarrow d_{k}$ of (3.29), we then have

$$\begin{aligned} \varvec{v}_{1}(d_{k}-)&= (\varvec{g}_{k}(d_{k}) + \varvec{c}_{k}) \varvec{U}_{k}(d_{k}),\end{aligned}$$

(3.36)

$$\begin{aligned} \varvec{v}_{1}(d_{k}+)&= (\varvec{g}_{k+1}(d_{k}) + \varvec{c}_{k+1}) \varvec{U}_{k+1}(d_{k})\nonumber \\&= \varvec{c}_{k+1}, \end{aligned}$$

(3.37)

where the last equality follows from (3.12) and (3.26). From (3.35), (3.36), and (3.37), we have

$$\begin{aligned} \varvec{c}_{k+1}&= r_{k+1}^{-1} \left\{ r_{k} \varvec{g}_{k}(d_{k}) \varvec{U}_{k}(d_{k}) + r_{k} \varvec{c}_{k} \varvec{U}_{k}(d_{k})\right\} \nonumber \\&= \frac{\lambda r_{k+1}^{-1} \varvec{\tau }}{1 + m(d_{K}) + \lambda m_{1}} \varvec{W}_{k}(d_{k}) + r_{k+1}^{-1} r_{k} \varvec{c}_{k} \hat{\varvec{U}}_{k,k} , \end{aligned}$$

(3.38)

for $\ 1 \le k \le K-1$, where the second equality sign follows from (3.33) and (3.17). By recursively applying (3.38) to itself, we obtain (3.30). $\square$

We now show that $\tilde{\varvec{u}}_{K+1}$ in (3.18) is equal to $\varvec{c}_{K+1}$ (here we extend the definition of $\varvec{c}_k$ in (3.30) to the case $k=K+1$).

Lemma 3.4

We have $\tilde{\varvec{u}}_{K+1} = \varvec{c}_{K+1}$, i.e.,

$$\begin{aligned} \tilde{\varvec{u}}_{K+1} = \frac{ \lambda r_{K+1}^{-1} \varvec{\tau } }{ 1 + m(d_{K}) + \lambda m_{1} } \left\{ \sum _{i=1}^{K} \varvec{W}_{i}(d_{i}) \hat{\varvec{U}}_{i+1,K} + \hat{\varvec{U}}_{1,K} \right\} . \end{aligned}$$

(3.39)

Proof

Similar to (3.13), from (3.8) and (3.20), we have

$$\begin{aligned} r_{K} \varvec{v}_{1}(d_{K}-) = r_{K+1} \varvec{v}_{1}(d_{K}+). \end{aligned}$$

(3.40)

From (3.29) and (3.18) (see also (3.12)), we have

$$\begin{aligned} \varvec{v}_{1}(d_{K}-)&= \frac{\lambda r_{K}^{-1} \varvec{\tau }}{1 + m(d_{K}) + \lambda m_{1}} \varvec{W}_{K}(d_{K}) + \varvec{c}_{K} \varvec{U}_{K}(d_{K}), \end{aligned}$$

(3.41)

$$\begin{aligned} \varvec{v}_{1}(d_{K}+)&= \tilde{\varvec{u}}_{K+1}. \end{aligned}$$

(3.42)

The preceding three equations and (3.30) imply that $\tilde{\varvec{u}}_{K+1} = \varvec{c}_{K+1}$. $\square$

3.5 Computation of the mean active period

We compute the mean active period $m_{1}$ by the normalizing condition, i.e.,

$$\begin{aligned} V(0) + \int _{0}^{d_{K}} (v_{0}(x) + \varvec{v}_{1}(x) \varvec{1}) \mathrm{d}x + \int _{d_{K}}^{\infty } \varvec{v}(x) \varvec{1} \mathrm{d}x = 1. \end{aligned}$$

(3.43)

From Lemma 3.1, we have

$$\begin{aligned} V(0) + \int _{0}^{d_{K}} v_{0}(x) \mathrm{d}x = \frac{1 + m(d_{K})}{1 + m(d_{K}) + \lambda m_{1}}. \end{aligned}$$

(3.44)

Note that from Lemma 3.2 and (3.39), for $x \in J_{k}$ and $K+1 \le k \le K_0$, we have

$$\begin{aligned} \varvec{v}(x)&= \frac{r_{K+1}}{r_{k}} \frac{ \lambda r_{K+1}^{-1} \varvec{\tau } }{ 1 + m(d_{K}) + \lambda m_{1} } \left\{ \sum _{i=1}^{K} \varvec{W}_{i}(d_{i}) \hat{\varvec{U}}_{i+1,K} + \hat{\varvec{U}}_{1,K} \right\} \hat{\varvec{U}}_{K+1,k-1} \varvec{U}_{k}(x)\nonumber \\&= \frac{ \lambda r_{k}^{-1} \varvec{\tau } }{ 1 + m(d_{K}) + \lambda m_{1} } \left\{ \sum _{i=1}^{K} \varvec{W}_{i}(d_{i}) \hat{\varvec{U}}_{i+1,k-1} + \hat{\varvec{U}}_{1,k-1} \right\} \varvec{U}_{k}(x)\nonumber \\&= \varvec{c}_{k} \varvec{U}_{k}(x), \end{aligned}$$

(3.45)

where the second and third equations follow from (3.17) and (3.30), respectively. From Lemmas 3.2 and 3.3, and (3.45), we have

$$\begin{aligned}&\int _{0}^{d_{K}} \varvec{v}_{1}(x) \varvec{1} \mathrm{d}x + \int _{d_{K}}^{\infty } \varvec{v}(x) \varvec{1} \mathrm{d}x \nonumber \\&\quad = \frac{\lambda \varvec{\tau }}{1 + m(d_{K}) + \lambda m_{1}} \left\{ \sum _{k=1}^{K} r_{k}^{-1} \int _{d_{k-1}}^{d_{k}} \varvec{W}_{k}(x) \mathrm{d}x + \sum _{k=1}^{K_0} r_{k}^{-1} \varvec{C}_{k} \int _{d_{k-1}}^{d_{k}} \varvec{U}_{k}(x) \mathrm{d}x \right\} \varvec{1}, \end{aligned}$$

(3.46)

where

$$\begin{aligned} \varvec{C}_{k} := \sum _{i=1}^{\min \{k-1,K\}} \varvec{W}_{i}(d_{i}) \hat{\varvec{U}}_{i+1,k-1} + \hat{\varvec{U}}_{1,k-1}, \quad 1 \le k \le K_0. \end{aligned}$$

(3.47)

Therefore, the mean active period is given as follows.

Theorem 3.1

From (3.43), (3.44), and (3.46), the mean active period is given by

$$\begin{aligned} m_{1} = \varvec{\tau } \left\{ \sum _{k=1}^{K} r_{k}^{-1} \int _{d_{k-1}}^{d_{k}} \varvec{W}_{k}(x) \mathrm{d}x + \sum _{k=1}^{K_0} r_{k}^{-1} \varvec{C}_{k} \int _{d_{k-1}}^{d_{k}} \varvec{U}_{k}(x) \mathrm{d}x \right\} \varvec{1}, \end{aligned}$$

(3.48)

where the integrals are calculated in the same way as (2.47) and (2.48), i.e.,

$$\begin{aligned} \int _{d_{k-1}}^{d_{k}} \varvec{W}_{k}(x) \mathrm{d}x&= (\varvec{I}_{N} \otimes \varvec{\tau }_{k}) \left[ \exp \left( \left( \begin{array}{cc} \varvec{O} &{}\quad \varvec{I}_{N^{2}+N} \\ \varvec{O} &{}\quad \overline{\varvec{M}}_{k} \end{array} \right) (d_{k} - d_{k-1}) \right) \right] ^{(N^{2},N)}, \quad 1 \le k \le K,\end{aligned}$$

(3.49)

$$\begin{aligned} \int _{d_{k-1}}^{d_{k}} \varvec{U}_{k}(x) \mathrm{d}x&= \left[ \exp \left( \left( \begin{array}{cc} \varvec{O} &{}\quad \varvec{I}_{N} \\ \varvec{O} &{}\quad \lambda r_{k}^{-1} \varvec{1} \varvec{\tau } + \varvec{T} \end{array} \right) (d_{k} - d_{k-1}) \right) \right] ^{(N,N)}, \quad 1 \le k \le K_0. \end{aligned}$$

(3.50)

For $k=K_0$, the last integral is given by $\left( - (\lambda r_{K_0}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}) \right) ^{-1}$.

3.6 Moments of the workload

In what follows, we show a computational procedure for the n-th moment of the workload, where n is a positive integer. To this end, we introduce some new notation. Let $\mathscr {M}$ be the set of all square matrices. For $\varvec{X} \in \mathscr {M}$, let $\sigma (\varvec{X})$ denote the order of $\varvec{X}$. For a positive integer n, let $\varvec{O}_{n}$ and $\varvec{I}_{n}$ denote the $n \times n$ zero and identity matrices, respectively. Let $\Psi : \mathscr {M}\rightarrow \mathscr {M}$ be a mapping defined as follows:

$$\begin{aligned} \Psi (\varvec{X}) := \left( \begin{array}{cc} \varvec{O}_{\sigma (\varvec{X})} &{}\quad \varvec{I}_{\sigma (\varvec{X})} \\ \varvec{O}_{\sigma (\varvec{X})} &{}\quad \varvec{X} \end{array} \right) , \quad \varvec{X} \in \mathscr {M}. \end{aligned}$$

(3.51)

From Corollary 2.1, we have, for $0 \le x \le y$,

$$\begin{aligned} \int _{x}^{y} \exp (\varvec{X} u) \mathrm{d}u = [ \exp (\Psi (\varvec{X}) y) - \exp (\Psi (\varvec{X}) x) ]^{(\sigma (\varvec{X}), \sigma (\varvec{X}))}. \end{aligned}$$

(3.52)

Lemma 3.5

For a positive integer n, $\varvec{X} \in \mathscr {M}$ and $d \ge 0$, we have

$$\begin{aligned} \int _{0}^{d} u^{n} \exp (\varvec{X}u) \mathrm{d}u = \sum _{i=0}^{n} (-1)^{i} \frac{n!}{(n-i)!} d^{n-i} [ \exp (\Psi ^{i+1}(\varvec{X}) d) ]^{(\sigma (\varvec{X}), \sigma (\varvec{X}))}, \end{aligned}$$

(3.53)

where $\Psi ^{i+1}$ is the $(i+1)$-th iterate of $\Psi$.

Proof

Noting that $u^{n} = \int _{0}^{u} n y^{n-1} \mathrm{d}y$, we have

$$\begin{aligned} \int _{0}^{d} u^{n} \exp (\varvec{X} u) \mathrm{d}u&= \int _{0}^{d} \left\{ \int _{y}^{d} \exp ( \varvec{X} u ) \mathrm{d}u \right\} n y^{n-1} \mathrm{d}y\nonumber \\&= \int _{0}^{d} n y^{n-1} \mathrm{d}y \times [ \exp (\Psi (\varvec{X}) d) ]^{(\sigma (\varvec{X}), \sigma (\varvec{X}))} \nonumber \\&\qquad - n \left[ \int _{0}^{d} y^{n-1} \exp ( \Psi ( \varvec{X} ) y ) \mathrm{d}y \right] ^{(\sigma (\varvec{X}), \sigma (\varvec{X}))}, \end{aligned}$$

(3.54)

where the last equality follows from (3.52). We then obtain the following recurrence formula:

$$\begin{aligned} \int _{0}^{d} u^{n} \exp (\varvec{X} u) \mathrm{d}u&= d^{n} [ \exp (\Psi (\varvec{X}) d) ]^{(\sigma (\varvec{X}), \sigma (\varvec{X}))}\nonumber \\&\quad - n \left[ \int _{0}^{d} u^{n-1} \exp ( \Psi ( \varvec{X} ) u ) \mathrm{d}u \right] ^{(\sigma (\varvec{X}), \sigma (\varvec{X}))}, \end{aligned}$$

(3.55)

which implies (3.53). $\square$

From Lemmas 3.1–3.3, we have

$$\begin{aligned} \mathsf {E}[Z^{n}]&= \frac{\varvec{\tau }}{1 + m(d_{K}) + \lambda m_{1}} \tilde{\varvec{D}}_{K}(n) \varvec{t} + \sum _{1 \le k \le K} \frac{\lambda r_{k}^{-1} \varvec{\tau }}{1 + m(d_{K}) + \lambda m_{1}} \tilde{\varvec{W}}_{k}(n) \varvec{1} \nonumber \\&\quad + \sum _{1 \le k \le K_{0}} \varvec{c}_{k} \tilde{\varvec{U}}_{k}(n) \varvec{1}, \end{aligned}$$

(3.56)

where $\varvec{t} := -\varvec{T}\varvec{1}$,

$$\begin{aligned} \tilde{\varvec{D}}_{K}(n)&:= \int _{0}^{d_{K}} x^{n} \exp (\varvec{D} x) \mathrm{d}x,\end{aligned}$$

(3.57)

$$\begin{aligned} \tilde{\varvec{W}}_{k}(n)&:= \int _{x \in J_{k}} x^{n} \varvec{W}_{k}(x) \mathrm{d}x, \quad 1 \le k \le K,\end{aligned}$$

(3.58)

$$\begin{aligned} \tilde{\varvec{U}}_{k}(n)&:= \int _{x \in J_{k}} x^{n} \varvec{U}_{k}(x) \mathrm{d}x,\quad 1 \le k \le K_{0}. \end{aligned}$$

(3.59)

By applying Lemma 3.5 to (3.56), a computable formula for the n-th moment of the workload is given as follows.

Theorem 3.2

For a positive integer n, the n-th moment of the workload is given by (3.56), where

$$\begin{aligned} \tilde{\varvec{D}}_{K}(n)&= \sum _{i=0}^{n} (-1)^{i} \frac{n!}{(n-i)!} (d_{K})^{n-i} [\exp ( \Psi ^{i+1}(\varvec{D}) d_{K} )]^{(N,N)}, \end{aligned}$$

(3.60)

$$\begin{aligned} \tilde{\varvec{W}}_{k}(n)&= (\varvec{I}_{N} \otimes \varvec{\tau }_{k}) \sum _{j=0}^{n} \left( {\begin{array}{c}n\\ j\end{array}}\right) (d_{k-1})^{n-j} \sum _{i=0}^{j} (-1)^{i} \frac{j!}{(j-i)!} (d_{k} - d_{k-1})^{j-i} \nonumber \\&\quad \times [ \exp ( \Psi ^{i+1}( \overline{\varvec{M}}_{k} ) (d_{k} - d_{k-1}) ) ]^{(N^{2}, N)}, \quad 1 \le k \le K, \end{aligned}$$

(3.61)

$$\begin{aligned} \tilde{\varvec{U}}_{k}(n)&= \sum _{j=0}^{n} \left( {\begin{array}{c}n\\ j\end{array}}\right) (d_{k-1})^{n-j} \sum _{i=0}^{j} (-1)^{i} \frac{j!}{(j-i)!} (d_{k} - d_{k-1})^{j-i} \nonumber \\&\quad \times [ \exp ( \Psi ^{i+1}( \lambda r_{k}^{-1} \varvec{1}\varvec{\tau } +\varvec{T} ) (d_{k} - d_{k-1}) ) ]^{(N, N)}, \quad 1 \le k \le K_{0}-1, \end{aligned}$$

(3.62)

$$\begin{aligned} \tilde{\varvec{U}}_{K_{0}}(n)&= \sum _{i=0}^{n} \frac{n!}{(n-i)!} (d_{K_{0}-1})^{n-i} (-\lambda r_{K_{0}}^{-1} \varvec{1}\varvec{\tau } - \varvec{T})^{-1-i}. \end{aligned}$$

(3.63)

Proof

Equation (3.60) is obtained by applying Lemma 3.5 to (3.57). Since $(-\lambda r_{K_{0}}^{-1} \varvec{1}\varvec{\tau } - \varvec{T})^{-1}$ exists, (3.63) is obtained by partially integrating (3.59) with $k := K_{0}$. Because the computation for (3.62) is similar to that of (3.61), we only show the derivation of the latter. From (3.58) and (3.27), we have

$$\begin{aligned} \tilde{\varvec{W}}_{k}(n)&= (\varvec{I}_{N} \otimes \varvec{\tau }_{k}) \int _{0}^{d_{k} - d_{k-1}} (x + d_{k-1})^{n} [ \exp (\varvec{\overline{M}}_{k} x) ]^{(N^{2}, N)} \mathrm{d}x \nonumber \\&= (\varvec{I}_{N} \otimes \varvec{\tau }_{k}) \sum _{j=0}^{n} \left( {\begin{array}{c}n\\ j\end{array}}\right) (d_{k-1})^{n-j} \left[ \int _{0}^{d_{k} - d_{k-1}} x^{j} \exp (\varvec{\overline{M}}_{k} x) \mathrm{d}x \right] ^{(N^{2}, N)}, \end{aligned}$$

(3.64)

where the last equality follows from the binomial theorem. We then obtain (3.61) by applying Lemma 3.5 to (3.64). $\square$

4 Example: the active period of an M/M/1 queue with two processing speeds and vacations

In previous sections, we have obtained a computable form for the mean active period; see (2.57) and (3.48). In this section, we further study the Laplace-Stieltjes transform (LST) of the active period for the special case with exponential service demand and multiple thresholds. From the LST, we derive the mean active period for an M/M/1 queue with two processing speeds and vacations. We also simplify the mean active period obtained in (2.57) and confirm that it is identical to the one obtained via the LST. Since this section is a supplement to support the previous formula for the mean active period, we only show key results and refer to our technical report [23] for their proofs.

4.1 LST of the active period for M/M/1 with multiple service speeds and vacations

We study the LST of the active period, denoted by A, of the queueing model in Sect. 3, where the service requirements are assumed to be exponentially distributed with mean $\mu ^{-1}$. For $s \ge 0$, let $\varphi (s) := \mathsf {E}[ e^{-s A} ]$ denote the LST of the active period, and $\varphi (s,x) := \mathsf {E}[e^{-s A} | \text {initial workload is } x]$ ($x > 0$) denote the LST of the active period when it starts with a workload at level x.

Lemma 4.1

For $x \in J_{k}$ and $\ 1 \le k \le K_{0}$, the LST of the active period starting at level x is given by

$$\begin{aligned} \varphi (s,x) = A_{k}(s) e^{\alpha _{k}(s) x} + B_{k}(s) e^{\beta _{k}(s) x}, \end{aligned}$$

(4.1)

where

$$\begin{aligned} \alpha _{k}(s)&= \frac{1}{2} \left\{ - \left( \frac{\lambda + s}{r_{k}} - \mu \right) + \sqrt{\left( \frac{\lambda + s}{r_{k}} - \mu \right) ^{2} + 4 \frac{\mu s}{r_{k}}} \right\} > 0, \end{aligned}$$

(4.2)

$$\begin{aligned} \beta _{k}(s)&= \frac{1}{2} \left\{ - \left( \frac{\lambda + s}{r_{k}} - \mu \right) - \sqrt{\left( \frac{\lambda + s}{r_{k}} - \mu \right) ^{2} + 4 \frac{\mu s}{r_{k}}} \right\} < 0, \end{aligned}$$

(4.3)

and $\{A_{k}(s), B_{k}(s) ; 1\le k \le K_{0}\}$ are constants which are determined by the following $2 K_{0}$ equations:

$$\begin{aligned} A_{K_{0}}(s) = 0, \quad A_{1}(s) + B_{1}(s) = 1, \end{aligned}$$

(4.4)

and, for $1 \le k \le K_{0}-1$,

$$\begin{aligned}&A_{k}(s) e^{\alpha _{k}(s) d_{k}} + B_{k}(s) e^{\beta _{k}(s) d_{k}} = A_{k+1}(s) e^{\alpha _{k+1}(s) d_{k}} + B_{k+1}(s) e^{\beta _{k+1}(s) d_{k}}, \end{aligned}$$

(4.5)

and

$$\begin{aligned}&r_{k} \left\{ A_{k}(s) \alpha _{k}(s) e^{\alpha _{k}(s) d_{k}} + B_{k}(s) \beta _{k}(s) e^{\beta _{k}(s) d_{k}} \right\} \nonumber \\&\qquad = r_{k+1} \left\{ A_{k+1}(s) \alpha _{k+1}(s) e^{\alpha _{k+1}(s) d_{k}} + B_{k+1}(s) \beta _{k+1}(s) e^{\beta _{k+1}(s) d_{k}} \right\} . \end{aligned}$$

(4.6)

We introduce some notation to express $\{A_{k}(s), B_{k}(s) ; 1\le k \le K_{0}\}$. For $1 \le k \le K_{0}-1$ and $i = 0, 1$, let

$$\begin{aligned}&a_{k+i,k} := e^{\alpha _{k+i}(s) d_{k}}, \quad b_{k+i,k} := e^{\beta _{k+i}(s) d_{k}},\end{aligned}$$

(4.7)

$$\begin{aligned}&a_{k+i,k}' := \alpha _{k+i}(s) e^{\alpha _{k+i}(s) d_{k}}, \quad b_{k+i,k}' := \beta _{k+i}(s) e^{\beta _{k+i}(s) d_{k}}, \end{aligned}$$

(4.8)

and then

$$\begin{aligned} \varvec{C}_{k} := \left( \begin{array}{cc} c_{00}^{(k)} &{}\quad c_{01}^{(k)} \\ c_{10}^{(k)} &{}\quad c_{11}^{(k)} \end{array} \right) , \qquad \overline{\varvec{C}}_{k} := \left( \begin{array}{cc} \overline{c}_{00}^{(k)} &{}\quad \overline{c}_{01}^{(k)} \\ \overline{c}_{10}^{(k)} &{}\quad \overline{c}_{11}^{(k)} \end{array} \right) = \prod _{l = k}^{K_{0}-1} \varvec{C}_{l}, \end{aligned}$$

(4.9)

where

$$\begin{aligned} c_{00}^{(k)}&:= \frac{a_{k+1,k}' b_{k,k} r_{k+1} - a_{k+1,k} b_{k,k}' r_{k}}{(a_{k,k}' b_{k,k} - a_{k,k} b_{k,k}' ) r_{k} }, \quad c_{01}^{(k)} := \frac{b_{k,k} b_{k+1,k}' r_{k+1} - b_{k+1,k} b_{k,k}' r_{k}}{(a_{k,k}' b_{k,k} - a_{k,k} b_{k,k}' ) r_{k} }.\end{aligned}$$

(4.10)

$$\begin{aligned} c_{10}^{(k)}&:= \frac{a_{k,k} a_{k+1,k}' r_{k+1} - a_{k+1,k} a_{k,k}' r_{k}}{(a_{k,k} b_{k,k}' - a_{k,k}' b_{k,k}) r_{k} }, \quad c_{11}^{(k)} := \frac{b_{k+1,k}' a_{k,k} r_{k+1} - b_{k+1,k} a_{k,k}' r_{k}}{(a_{k,k} b_{k,k}' - a_{k,k}' b_{k,k}) r_{k} }. \end{aligned}$$

(4.11)

We obtain the constants $\{A_{k}(s), B_{k}(s) ; 1\le k \le K_{0}\}$ as follows.

Lemma 4.2

For $1\le k \le K_{0}-1$,

$$\begin{aligned} A_{k}(s) = \frac{\overline{c}_{01}^{(k)}}{\overline{c}_{01}^{(1)} + \overline{c}_{11}^{(1)}}, \quad B_{k}(s) = \frac{\overline{c}_{11}^{(k)}}{\overline{c}_{01}^{(1)} + \overline{c}_{11}^{(1)}}, \end{aligned}$$

(4.12)

and

$$\begin{aligned} A_{K_{0}}(s) = 0, \quad B_{K_{0}}(s) = \frac{1}{\overline{c}_{01}^{(1)} + \overline{c}_{11}^{(1)}}. \end{aligned}$$

(4.13)

4.2 Mean active period in case of two service speeds, Approach 1: LST

We apply Lemma 4.1 to the case of two service speeds, i.e., $K_{0} := 2$, and derive the mean active period for the M/M/1 queue with two processing speeds and vacations. We denote the traffic intensities when the workload is in $[0,d_{1})$ and $[d_{1},\infty )$ by $\rho _{1} := \lambda /(r_{1}\mu )$ and $\rho _{2} := \lambda /(r_{2}\mu )$, respectively. The stability condition (see also (2.4)) is rewritten as

$$\begin{aligned} \rho _{2} < 1. \end{aligned}$$

(4.14)

By conditioning on the initial fluid level when the active period starts (see also Lemma 4.1), and from the memoryless property of the exponential distribution, the LST of the active period is given as follows: for $s \ge 0$, we have

$$\begin{aligned} \varphi (s) = \int _{0}^{\infty } \varphi (s, d_{1}+x) \mu e^{-\mu x} \mathrm{d}x = \frac{\mu e^{\beta _{2}(s) d_{1}}}{ \mu - \beta _{2}(s) } B_{2}(s), \end{aligned}$$

(4.15)

where the coefficients $\beta _{2}(s)$ and $B_{2}(s)$ are given as follows (see also (4.3) and (4.13)):

$$\begin{aligned} \beta _{2}(s)&= \frac{1}{2} \left\{ - \left( \frac{\lambda + s}{r_{2}} - \mu \right) - \sqrt{\left( \frac{\lambda + s}{r_{2}} - \mu \right) ^{2} + 4 \frac{\mu s}{r_{2}}} \right\} , \end{aligned}$$

(4.16)

$$\begin{aligned} B_{2}(s)&= \frac{1}{c_{01}^{(1)} + c_{11}^{(1)}}. \end{aligned}$$

(4.17)

Note that from (4.15) and (4.16), we have

$$\begin{aligned} \beta _{2}(0) = 0, \quad \varphi (0) = B_{2}(0) = 1, \end{aligned}$$

(4.18)

and then the mean active period $m_{1} = - \varphi '(0)$ is calculated as follows:

$$\begin{aligned} m_{1} = - B_{2}'(0) - \frac{1 + \mu d_{1}}{\mu } \beta _{2}'(0), \end{aligned}$$

(4.19)

where $\beta _{2}'(0) = - \{ r_{2} (1 - \rho _{2}) \}^{-1}$, and, for $\rho _{1} \not = 1$,

$$\begin{aligned} B_{2}'(0) = \frac{1}{1-\rho _{1}} \left\{ - \frac{d_{1}}{r_{1}} + \frac{1-e^{-(1-\rho _{1})\mu d_{1}}}{(1-\rho _{1}) \mu r_{1}} \right\} + \frac{1}{1 - \rho _{2}} \left\{ \frac{d_{1}}{r_{2}} - \frac{1 - e^{-(1-\rho _{1})\mu d_{1}}}{(1-\rho _{1}) \mu r_{1}} \right\} . \end{aligned}$$

(4.20)

From (4.19) and (4.20), the mean active period is given as follows.

Corollary 4.1

For $\rho _{1} \not = 1$, we have

$$\begin{aligned} m_{1} = \frac{\rho _{1}}{\lambda } \frac{(1 - \rho _{1}) \mu d_{1} - (1 - e^{-(1 - \rho _{1}) \mu d_{1}}) \rho _{1} }{(1 - \rho _{1})^{2}} + \frac{1 - \rho _{1} e^{-(1 - \rho _{1}) \mu d_{1}}}{\lambda (1 - \rho _{1})} \frac{\rho _{2}}{1- \rho _{2}}. \end{aligned}$$

(4.21)

Since the mean active period is a monotonic function of $\rho _{1}$, the mean active period for $\rho _{1} = 1$ is given by letting $\rho _{1} \rightarrow 1$ in (4.21), and then we have, for $\rho _{1} = 1$,

$$\begin{aligned} m_{1} = \frac{\mu d_{1}(2 + \mu d_{1})}{2 \lambda } + \frac{1 + \mu d_{1}}{\lambda } \frac{\rho _{2}}{1 - \rho _{2}}. \end{aligned}$$

(4.22)

4.3 Mean active period in case of two service speeds, Approach 2: application of Theorem 2.1

We apply Theorem 2.1 to the M/M/1 queue with two processing speeds and vacations. Since the service time is assumed to be exponentially distributed with mean $\mu ^{-1}$, we put

$$\begin{aligned} N:=1, \quad \varvec{\tau } := 1, \quad \varvec{T} := -\mu , \quad \varvec{t} := \mu , \quad \varvec{D} := 0, \end{aligned}$$

(4.23)

and apply

$$\begin{aligned}&\varvec{I}_{N} \otimes \varvec{\tau } := 1,\quad \varvec{I}_{N^{2}+N} := \varvec{I}_{2},\quad \lambda r_{1}^{-1} \varvec{1} \varvec{\tau } + \varvec{T} := -(1 - \rho _{1})\mu ,\nonumber \\&\quad \overline{\varvec{M}} := \left( \begin{array}{cc} 0 &{}\quad \mu \\ 0 &{}\quad -(1 - \rho _{1})\mu \end{array} \right) \end{aligned}$$

(4.24)

to (2.57), and then we have, for $\rho _{1} \not = 1$,

$$\begin{aligned} \left[ \exp \left( \left( \begin{array}{cc} \varvec{O} &{}\quad \varvec{I}_{N^{2}+N} \\ \varvec{O} &{}\quad \overline{\varvec{M}} \end{array} \right) d_{1} \right) \right] ^{(N^{2},N)}= & {} \frac{ (1 - \rho _{1}) \mu d_{1} - (1 - e^{-(1 - \rho _{1})\mu d_{1}}) }{ \mu (1 - \rho _{1})^{2} },\nonumber \\ \end{aligned}$$

(4.25)

$$\begin{aligned} \left[ \exp \left( \left( \begin{array}{cc} \varvec{O} &{}\quad \varvec{I}_{N} \\ \varvec{O} &{}\quad \lambda r_{1}^{-1} \varvec{1}\varvec{\tau } + \varvec{T} \end{array} \right) d_{1} \right) \right] ^{(N,N)}= & {} \frac{1 - e^{ -(1 - \rho _{1}) \mu d_{1} }}{ \mu (1 - \rho _{1}) }, \end{aligned}$$

(4.26)

$$\begin{aligned}{}[\exp (\overline{\varvec{M}} d_{1})]^{(N^{2},N)}= & {} \frac{1 - e^{ -(1 - \rho _{1}) \mu d_{1} } }{ 1 - \rho _{1} }, \end{aligned}$$

(4.27)

$$\begin{aligned} (-(\lambda r_{2}^{-1} \varvec{1} \varvec{\tau } + \varvec{T}))^{-1}= & {} \frac{1}{(1 - \rho _{2})\mu }. \end{aligned}$$

(4.28)

From (2.57) and (4.23)–(4.28), we obtain the same formulas given in Corollary 4.1 for the mean active period.

Corollary 4.2

For $\rho _{1} \not = 1$ and $\rho _{1} = 1$, the mean active period is given by (4.21) and (4.22), respectively.

5 Cost-optimization problem and numerical results

In this section, we study a cost-optimization problem and show some numerical examples. To this end, we define a cost function taking into account the trade-off issue between delay performance and energy consumption in Sect. 5.1, and then consider the optimization problem for minimizing the cost function in Sect. 5.2. In particular, we will find the activation threshold which minimizes the cost function. Furthermore, we also investigate the service curve which minimizes the cost function.

5.1 Formulation of optimization problem

We consider the following three types of costs per unit of time:

• $c_{p,k}$ ($1 \le k \le K_0$); A cost related to the power consumption when the server is in mode 1 and processes the workload with rate $r_{k}$.

• $c_{p,0};$ A cost related to the power consumption when the server is in mode 0.

• $c_{s};$ A cost when the server switches from mode 0 to 1.

For $1 \le k \le K_0$, let $m_{1,k}$ denote the mean active period when the workload is in $J_{k}$. The mean power consumption cost (energy consumption cost per time unit) is given by

$$\begin{aligned} \sum _{k=1}^{K_0} c_{p,k} \frac{m_{1,k}}{m_0+m_1} + c_{p,0} \frac{m_{0}}{m_0+m_1} + c_{s} \frac{1}{m_{0} + m_{1}}, \end{aligned}$$

(5.1)

where $m_{0} = \lambda ^{-1} (1 + m(d_{K}))$ (see Lemma 3.1) and $m_{1}$ is given by (2.57). From Theorem 3.1, $m_{1,k}$ ($1 \le k \le K_0$) is given as follows:

$$\begin{aligned} m_{1,k} = \varvec{\tau } \left\{ \left( r_{k}^{-1} \int _{d_{k-1}}^{d_{k}} \varvec{W}_{k}(x) \mathrm{d}x\right) \text{1 }\text{ l }(k \le K) + r_{k}^{-1} \varvec{C}_{k} \int _{d_{k-1}}^{d_{k}} \varvec{U}_{k}(x) \mathrm{d}x \right\} \varvec{1}. \end{aligned}$$

(5.2)

We note that each of $m_{0}$, $m_{1}$, and $m_{1,k}$ ($1 \le k \le K_0$) is given by a matrix exponential form (see Lemma 3.1 and Theorem 3.1), therefore it is easy to implement (5.1) in a numerical calculation.

In the cost function above, the cost $c_{p,k}$ will be set to $c_{p,k} = c_p r_k^2$. There is evidence that power consumption is a convex function of the processing speed and it is reasonable to set it as a quadratic function of the speed [17].

Also taking into account the performance, we consider the following cost function for our system:

$$\begin{aligned} \mathrm{Cost} = c_h \mathrm{E}[Z] + c_p \left( \sum _{k=1}^{K_0} r_k^2 \frac{m_{1,k}}{m_0+m_1} \right) + c_{p,0}\frac{ m_{0}}{m_0+m_1} + c_{s} \frac{1}{m_{0} + m_{1}}. \end{aligned}$$

(5.3)

In (5.3), the first term is related to the performance, i.e., the smaller the mean workload, the smaller the response time for jobs is. The second term (the summation) is related to the power consumption of the server in the active period. The smaller the processing speed $r_k$, the smaller the power consumption is; but it leads to a bigger $\mathrm{E}[Z]$. The third term is the holding cost (power consumption when the server is inactive). It should be noted that it is reasonable to set $c_{p,0} = 0.6 c_p$ as the server in inactive mode can consume about 60% power, compared to when it is busy processing a job. The last term is related to the switching cost. It should be noted that a CPU instantaneously consumes a large amount of energy once it is switched on.

5.2 Numerical results

We consider the effect of the service rate function r(x) and the activation threshold $d_K$ on the cost function. To this end, we fix the arrival rate and the job size distribution. In particular, $\lambda = 1$ and the job sizes follow a two-stage Erlang distribution with mean 2. The coefficients in the cost function are set as follows: $c_h = 0.1, c_s = 30, c_p = 1, c_{p,0} = 0.6 \times c_p$. We consider a service rate function of the form $r(x) = r_1 x^\alpha + r_0$, where we restrict the parameters as follows: $r_0 = 1$ and $r_1 = 0.1, 1, 10, 100$. The service rate function is approximated by the step function with step size 0.1 and $d_{K_{0} - 1} = 20$, and for $x \ge 20$ the service rate is approximated by r(20). We first consider the cost function against $d_K$ for some special service rate functions: $r(x) = r_1 x + r_0$, $r(x) = r_1 \sqrt{x} + r_0$ and $r(x) = r_1 x^2 + r_0$, where we fix $r_0 = 1$ and consider $r_1 = 0.1,1,10,100$.

Figure 3 shows the cost function against the threshold $d_K$ for $r(x) = r_1 \sqrt{x} + r_0$. In the case $r_1=0.1$, the stability condition is violated. We observe that the curves for $r_1 = 1$ and 10 are convex, implying the existence of a threshold that minimizes the cost function. In the third curve with $r_1 = 100$, the cost function monotonically increases, implying that $d_K = 0$ is the optimal threshold. These results suggest that when the service rate is large enough, it is optimal to switch the server on as soon as a job is available.

Fig. 3

Cost function against threshold ($r(x) = r_1 \sqrt{x} + r_0$)

Figure 4 shows the cost function against the threshold $d_K$ for $r(x) = r_1 x + r_0$. This figure shows that the cost function is minimized at positive values of $d_K$ for $r_1=0.1, 1$ and 10 while it is minimized at $d_K = 0$ for $r_1=100$.

Fig. 4

Cost function against threshold ($r(x) = r_1 x + r_0$)

Figure 5 shows the cost function against the threshold $d_K$ for $r(x) = r_1 x^2 + r_0$. We observe that the cost function is minimized at a positive value of $d_K$ for $r_1= 0.1$ and 1 while it is minimized at $d_K = 0$ for $r_1 = 10$ and 100.

Fig. 5

Cost function against threshold ($r(x) = r_1 x^2 + r_0$)

Observing all the graphs above, there exists a threshold $d_K$ which minimizes the cost function. Sometimes the cost function is minimized at $d_K = 0$, and in some other cases it is minimized at a non-trivial $d_K > 0$. The common trend is that for a fast service rate (larger $r_1$ and/or $\alpha$), the cost function is likely minimized at $d_K = 0$, whereas if the service rate is relatively small, the cost function is likely minimized at some positive value of the threshold $d_K$. For the case $r_1 = 100$, all the curves show that the cost function is minimized at $d_K = 0$.

Fig. 6

Case of $r_1 = 0.1$

Fig. 7

Case of $r_1 = 1.0$

This motivates us to have a closer look at the minimal cost for each fixed $\alpha \in [0.1,3]$, where the service curve is given by $r(x) = r_1 x^\alpha + r_0$. Figures 6, 7, and 8 show the optimal threshold $d_{K}$ (on the right y-axis) and the corresponding cost (on the left y-axis) against $\alpha$. We also find that the optimized cost function is minimized at some positive $\alpha$ for $r_1=0.1$, while it monotonically increases with $\alpha$ for $r_1 = 1, 10$. This implies that, for relatively large $r_1$, the value of $\alpha$ should be small so that the service rate is not too large to balance the power consumption.

From Figs. 6, 7, and 8, we observe a general rule that in most of the cases ($\alpha \ge 0.2$), the optimal threshold $d_K$ decreases with increasing $\alpha$. This suggests that for a fast service rate, it is better to set a low threshold.

Fig. 8

Case of $r_1 = 10.0$