# Analysis of Single-Server Queue with Phase-Type Service and Energy Harvesting

Authors: Sergey A. Dudin; Moon Ho Lee
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

We propose a queueing model suitable, for example, for modelling operation of nodes of sensor networks. The sensor node senses a random field and generates packets to be transmitted to the central node. The sensor node has a battery of a finite capacity and harvests energy during its operation from outside (using solar cells, wind turbines, piezoelectric cells, etc.). We assume that, generally speaking, service (transmission) of a packet consists of a random number of phases and implementation of each phase requires a unit of energy. If the battery becomes empty, transmission is failed. To reduce the probability of forced transmission termination, we suggest that the packet can be accepted for transmission only when the number of energy units is greater than or equal to some threshold. Under quite general assumptions about the pattern of the arrival processes of packets and energy, we compute the stationary distributions of the system states and the waiting time of a packet in the system and numerically analyze performance measures of the system as functions of the threshold. Validity of Little’s formula and its counterpart is verified.


---

## Body

## 1. Introduction

The wireless sensor networks technology has great advancement; small and smart wireless sensor networks now can be used for different complicated and challenging applications. Wireless sensor networks found applications for environmental monitoring, animal tracking and control, safety, security and military purposes, built environment, health, and so forth; for more details see, for example, Kausar et al. [1]. Nodes of sensor networks have small batteries with limited power and also have limited computational power and storage space. When the battery of a node is exhausted, it is not replaced and the node dies. When sufficient number of nodes die, the network may not be able to perform its designated task. Thus the life time of a network is an important characteristic of a sensor network and it is tied up with the life time of a node. Recent advances in energy harvesting technology have resulted in the design of new types of sensor nodes which are able to extract energy from the surrounding environment. The major sources of energy harvesting include solar, wind, sound, vibration, thermal, and electromagnetic power. The concept of extracting ambient energy is to convert harvested energy from existing environmental sources into electricity to power sensor nodes; an energy storage device is used to accumulate such energy. Energy harvesting sensor nodes have the potential to perpetuate the life time of a battery through continuous energy harvesting. Thus, a lot of attention in the literature is devoted to telecommunication systems with energy harvesting. Recent surveys in this subject are published in Cui et al. [2], Lu et al. [3], Zhang and Lau [4], and Ulukus et al. [5]. Because it is well known that many problems related to capacity planning, performance evaluation, and optimization of telecommunication networks can be effectively solved with help of queueing theory, there are a large number of papers devoted to application of this theory for analysis of telecommunication networks with energy harvesting. For background information see, for example, Sharma et al. [6], Tutuncuoglu and Yener [7], Yang and Ulukus [8], and Yang and Ulukus [9].In the existing literature, it is usually assumed that customers (packets) arrive according to the stationary Poisson process and are buffered if the server is not available at the customer arrival moment. It is also suggested usually that harvested energy is slotted to some discrete units where a unit is the amount of energy which is used to provide service to one customer. Arriving units of energy are accumulated in the buffer having a finite capacity. Similar to our model, a model was recently considered in Gelenbe [10]. The main contributions of our results in comparison with Gelenbe [10] are as follows.(i)
It is assumed in Gelenbe [10] and the overwhelming majority of other papers that the arrival flows of customers and units of energy are described by the stationary Poisson processes. This assumption greatly simplifies the analysis of the model, but it may be not true in many real life applications. The stationary Poisson process has a property that the intensity of arrival is constant. While in modelling operation of the sensor node of the network designed for monitoring the state of some object, arrivals of customers (signals from detectors of the node) may be quite rare in normal situation and may become quite frequent in security threat situation. The intensity of harvesting the units of energy may essentially vary depending on the state of environment of a sensor node. For example, it may essentially increase when the weather is sunny or windy and decrease in the opposite case if the energy is harvested from solar panels or wind turbines. In our paper, we suggest much more general processes of customers and energy arrivals, namely, Markovian arrival processes. This allows easily modelling situations when the intensities of arrivals essentially fluctuate.
(ii)
It is assumed in Gelenbe [10] that the service time is equal to zero. If a customer arrives and energy presents in the system, the customer and one unit of energy instantaneously leave the system. Here we assume more general situation. Service of a customer is not instantaneous. Service lasts during some random time. This time is assumed to be having so-called phase-type (P
H) distribution. P
H distribution is much more general than an exponential distribution, which is popular in the literature. Random time having such a distribution can be interpreted as the sum of known or random number of independent random times (phases), duration of which has an exponential distribution. Bearing in mind this interpretation, we assume that the unit of energy corresponds to amount of energy that is necessary for implementation of one phase of service. Thus, some number of energy units is required for service of a customer and if, at the beginning of some phase of service, the buffer of energy is empty, service is terminated and the customer in service is forced to permanently leave the system. Thus, the problem of control by service initiation arises. We assume that the control strategy is of the threshold type. Some threshold is fixed. Service of the customer does not start if the current number of energy units in the buffer is less than this threshold. Evident advantages of such a strategy are that, under the proper choice of the threshold, it allows decreasing the probability of forced termination of service and the waste of energy units spent during the elapsed service time. Presented in our paper, results are useful for solving the problem of the optimal choice of the threshold.It is worth noting that we assume in our model that the waiting time of a customer in the buffer is restricted. This allows taking into account the fact that the packet in a sensor node should be dropped if it cannot be successfully transmitted during a specific time interval.The problem of optimal (e.g., with respect to the loss probability of a customer) choice of the threshold is far from a trivial one. If this threshold is too small, many customers are lost due to the forced service termination caused by the lack of energy. If this threshold is too large, many customers are lost due to the restriction of allowed waiting time (impatience). Thus, careful mathematical analysis is necessary to solve an optimization problem.The rest of the paper is organized as follows. In Section2, the mathematical model is described. In Section 3, the process of the system states is defined and the problem of computation of the stationary probabilities of the system states is discussed. Formulas for computation of the performance measures of the system are presented in Section 4. The problem of computation of the stationary distribution of the waiting time is solved in Section 5. Numerical illustrations are given and briefly discussed in Section 6. Section 7 concludes the paper.


## 2. Mathematical Model

We consider a single-server queueing system with Markovian arrival processes of customers and energy. The structure of the system under study is presented in Figure1.Figure 1
Queueing system under study.Customers arrive at the system according to the Markovian arrival process. We will code this process asM
A
P
c. Advantage of the M
A
P pattern of the arrival process comparing to a popular in the queueing literature model of the stationary Poisson process (which is a very particular case of the M
A
P) is that the M
A
P allows taking into account correlation of the successive interarrival times and burstiness typical for modern telecommunication networks. Arrivals in the M
A
P
c are directed by an irreducible continuous-time Markov chain w
t, t
≥
0, with the finite state space {
0,1
,
…
,
W
}. The sojourn time of the Markov chain w
t, t
≥
0, in the state w has an exponential distribution with the parameter λ
w, w
=
0
,
W
¯. Here, notation such as w
=
0
,
W
¯ means that w assumes values from the set {
0,1
,
…
,
W
}. After this sojourn time expires, with probability p
k
(
w
,
w
′
) the process w
t transits to the state w
′ and k customers, k
=
0,1, arrive at the system. The intensities of transitions from one state to another, which are accompanied by the arrival of k customers, are combined to the square matrices D
k, k
=
0,1, of size W
+
1. The matrix generating function of these matrices is D
(
z
)
=
D
0
+
D
1
z, |
z
|
≤
1. The matrix D
(
1
) is an infinitesimal generator of the process w
t, t
≥
0. The stationary distribution vector θ of this process satisfies the system of equations θ
D
(
1
)
=
0, θ
e
=
1. Here and throughout this paper, 0 is a zero row vector, and e denotes unit column vector. In case the dimension of a vector is not clear from the context, it is indicated as a lower index.The average intensityλ
c (fundamental rate) of the M
A
P
c is defined by λ
c
=
θ
D
1
e. The variance v of customers interarrival time is calculated by (1)
v
=
2
λ
c
-
1
θ
-
D
0
-
1
e
-
λ
c
-
2
,the squared coefficient of the variation c
v
a
r is calculated as (2)
c
v
a
r
=
2
λ
c
θ
-
D
0
-
1
e
-
1
,and the coefficient of correlation c
c
o
r of intervals between successive arrivals is given as (3)
c
c
o
r
=
λ
c
-
1
θ
-
D
0
-
1
D
1
-
D
0
-
1
e
-
λ
c
-
2
v
.For more information about the M
A
P
s and their usefulness in modelling telecommunication networks see Chakravarthy [11]. Necessity of careful account of burstiness of the arrival process in the considered model is explained, as it was mentioned above, by the following essential feature of the nodes of the sensor networks: arrivals of customers (signals from detectors of the node) may be quite rare in normal situation (e.g., when there is no danger for the object under protection or the monitored animals do not move) and may become quite frequent in security threat situation or situation when the animals quickly relocate.Units of energy arrive at the system according to the Markovian arrival processM
A
P
e. Arrivals in the M
A
P
e are directed by an irreducible continuous-time Markov chain v
t, t
≥
0, with the finite state space {
0,1
,
…
,
V
}. The M
A
P
e is defined by the matrices H
0 and H
1. Let us denote the average intensity of the M
A
P
e as λ
e. The buffer for energy has a finite capacity N. An energy unit is lost if it arrives when the buffer is full. Necessity of account of burstiness of the energy arrival process in the considered model is quite clear and also was mentioned above.The service time of a customer for the server has aP
H distribution with an irreducible representation (
β
,
S
). This service time can be interpreted as the time until the underlying Markov process m
t, t
≥
0, with the finite state space {
1
,
…
,
M
,
M
+
1
}, reaches the single absorbing state M
+
1 conditioned on the fact that the initial state of this process is selected among the states {
1
,
…
,
M
} according to the probabilistic row vector β
=
(
β
1
,
…
,
β
M
). The transition rates of the process m
t within the set {
1
,
…
,
M
} are defined by the subgenerator S, and the transition rates into the absorbing state (which leads to service completion) are given by the entries of the column vector S
0
=
-
S
e. The mean service time is calculated as b
1
=
β
(
-
S
)
-
1
e.The class ofP
H distributions is dense (in the sense of weak convergence) in the set of all probability distributions of nonnegative variables. Thus, the P
H distribution is very general and can be used for approximation of an arbitrary distribution; see Asmussen [12].It is worth stressing that the assumptions made in our paper about the arrival processes of customers and energy units and the service process are quite realistic and are much weaker than the standard in the queueing literature assumptions that the flows are described by the stationary Poisson arrival processes and the service times are exponentially distributed. The single disadvantage of our assumptions is as follows. The stationary Poisson arrival process and the exponential distribution are defined by one parameter. Thus, estimation of this parameter based on real data is extremely simple. Estimation of parameters of theM
A
P and P
H is much more difficult problem. But this problem is successfully solved in the relevant literature. For references see, for example, Buchholz et al. [13], Buchholz and Panchenko [14], and Okamura and Dohi [15]. Importance of the use of the M
A
P and P
H for modelling of real sensor nodes operation where the arrival process of signals is correlated will be illustrated in the section devoted to numerical results. Correlation drastically changes the performance measures of the system and its ignorance, which occurs when the arrival process is assumed to be the stationary Poisson process, leads to poor prediction of the system performance measures.We assume that implementation of each phase of service requires a unit of energy. The unit, if any, leaves the buffer of energy at the phase beginning moment. If the buffer of energy is empty at such a moment, during service, service of a customer is interrupted and the customer is lost. The server switches to the sleep mode.We offer the following strategy of waking the server up. The staying of the server in the sleep mode finishes when the number of energy units becomes greater than or equal to the thresholdK, 1
≤
K
≤
N. Thus, if during an arbitrary customer arrival epoch the server is idle and the number of energy units is greater than or equal to the threshold K, this customer occupies the server (see Figure 2).Figure 2
Arriving customer occupies the server.Otherwise the customer joins the buffer (see Figure3).Figure 3
Arriving customer joins the buffer.The customers from the buffer are serviced in the order of arrival and the first customer from the buffer is accepted for service when the server becomes free and the number of energy units is greater than or equal to the thresholdK.The customers from the buffer are impatient; that is, the customer leaves the buffer and the system after an exponentially distributed time described by the parameterγ, 0
<
γ
<
∞. The impatience of a customer may be interpreted as the obsolescence of information collected by the sensor. After some amount of time, information becomes useless and may be dropped.We assume that the thresholdK is a control parameter. Note that if the value of the threshold K is too small, many customers will be lost during service due to the lack of energy. If the value of the threshold K is too large, which causes long sleep periods, then many customers will be lost due to impatience. Thus, the task of optimal choice of the threshold K has a major importance. In the general case, the number of units of energy that is required for service of one customer is random and unknown at the service beginning epoch. This makes the problem of suitable choice of the value of K difficult. In some particular cases, this problem becomes a bit easier. For example, if P
H distribution is the Erlangian of order m, then the choice of K such that K
>
m guarantees that service of any customer will not be interrupted because, at the customer’s service beginning epoch, there is enough energy to completely serve the customer. However, because the energy arrives during the service time, probably it does not make sense to wait until m energy units will be accumulated in the system to start service of a customer.Since the arrival of customers and energy units occurs at random moments and the service time is random as well, the number of customers and energy units in the system at an arbitrary moment is random. Thus, to have an opportunity to compute the performance measures of the system under arbitrary fixed value of the thresholdK, it is necessary to compute the stationary distribution of the system states. To this end, it is necessary to construct the Markovian process describing the dynamics of the system states and compute its stationary distribution. This will be done in the next section.


## 3. Process of System States

It is easy to see that the behavior of the system under study is described by the following regular irreducible continuous-time Markov chain:(4)
ξ
t
=
i
t
,
r
t
,
n
t
,
w
t
,
v
t
,
m
t
,
t
≥
0
,where, during the epoch t, t
≥
0,(i)
i
t is the number of customers in the buffer, i
t
≥
0;
(ii)
r
t is an indicator that indicates whether the server works or not: r
t
=
0 corresponds to the case when the server does not work and r
t
=
1 if the server works;
(iii)
n
t is the number of energy units in the buffer, n
t
=
0
,
N
¯;
(iv)
w
t is the state of the underlying process of the M
A
P
c, w
t
=
0
,
W
¯;
(v)
v
t is the state of the underlying process of the M
A
P
e, v
t
=
0
,
V
¯;
(vi)
m
t is the state of the underlying process of the P
H service process, m
t
=
1
,
M
¯.The Markov chainξ
t, t
≥
0, has the following state space: (5)
0,0
,
n
,
w
,
v
,
n
=
0
,
N
¯
i
,
0
,
n
,
w
,
v
,
i
>
0
,
n
=
0
,
K
¯
∪
i
,
1
,
n
,
w
,
v
,
m
,
i
≥
0
,
n
=
0
,
N
¯
,
w
=
0
,
W
¯
,
v
=
0
,
V
¯
,
m
=
1
,
M
¯
.To simplify the analysis of the Markov chainξ
t, t
≥
0, let us enumerate the states of this process in the direct lexicographic order of the components r
,
n
,
w
,
v
,
m, and refer to the set of the states of the Markov chain having values (
i
,
r
) of the first two components of the Markov chain as the macrostate (
i
,
r
). Let G be the generator of the Markov chain ξ
t, t
≥
0. It consists of the blocks G
i
,
j, i
,
j
≥
0, each of which consists of blocks G
i
,
j
(
r
,
r
′
), r
,
r
′
=
0,1, defining the intensities of transitions from the macrostate (
i
,
r
) to the macrostate (
j
,
r
′
).Analyzing all possible transitions of the Markov chainξ
t, t
≥
0, during an interval of infinitesimal length and rewriting the intensities of these transitions in block matrix form we obtain the following result.Theorem 1.
The infinitesimal generatorG
=
(
G
i
,
j
)
i
,
j
≥
0 of the Markov chain ξ
t, t
≥
0, has a block-tridiagonal structure: (6)
G
=
G
0,0
G
0,1
O
O
…
G
1,0
G
1,1
G
1,2
O
…
O
G
2,1
G
2,2
G
2,3
…
⋮
⋮
⋮
⋮
⋱
.The blocks of the matrix G
i
,
i, i
≥
0, have the following form:(7)
G
i
,
i
=
G
i
,
i
0,0
G
i
,
i
0,1
G
i
,
i
1,0
G
i
,
i
1,1
,
i
≥
0
,
G
0,0
0,0
=
I
N
+
1
⊗
D
0
⊕
H
0
+
E
N
+
⊗
I
W
-
⊗
H
1
,
G
i
,
i
0,0
=
I
K
⊗
D
0
⊕
H
0
+
E
K
+
⊗
I
W
-
⊗
H
1
-
i
γ
I
,
i
>
0
,
G
0,0
0,1
=
E
-
-
⊗
D
1
⊗
I
V
-
⊗
β
,
G
i
,
i
0,1
=
O
,
i
>
0
,
G
0,0
1,0
=
I
^
N
+
1
,
N
+
1
⊗
I
W
-
V
-
⊗
S
-
S
d
e
+
I
N
+
1
W
-
V
-
⊗
S
0
,
G
i
,
i
1,0
=
I
^
N
+
1
,
K
⊗
I
W
-
V
-
⊗
S
-
S
d
e
+
I
N
+
1
,
K
⊗
I
W
-
V
-
⊗
S
0
,
i
>
0
,
G
i
,
i
1,1
=
I
N
+
1
⊗
D
0
⊕
H
0
⊕
S
d
+
E
N
+
⊗
I
W
-
⊗
H
1
⊗
I
M
-
i
γ
I
+
E
-
⊗
I
W
-
V
-
⊗
S
-
S
d
,
i
≥
0
.
The nonzero blocks of the matrixG
i
,
i
-
1, i
≥
1, have the following form:(8)
G
i
,
i
-
1
=
G
i
,
i
-
1
0,0
G
i
,
i
-
1
0,1
O
G
i
,
i
-
1
1,1
,
i
≥
1
,
G
1,0
0,0
=
γ
I
K
,
N
+
1
⊗
I
W
-
V
-
,
G
i
,
i
-
1
0,0
=
i
γ
I
K
W
-
V
-
,
i
≥
2
,
G
i
,
i
-
1
0,1
=
E
~
⊗
I
W
-
⊗
H
1
⊗
β
,
i
≥
1
,
G
i
,
i
-
1
1,1
=
i
γ
I
N
+
1
W
-
V
-
M
+
E
-
-
⊗
I
W
-
V
-
⊗
S
0
β
,
i
≥
1
.
The nonzero blocks of the matrixG
i
,
i
+
1, i
≥
0, have the following form: (9)
G
i
,
i
+
1
=
G
i
,
i
+
1
0,0
O
O
G
i
,
i
+
1
1,1
,
i
≥
0
,
G
0,1
0,0
=
I
-
⊗
D
1
⊗
I
V
-
,
G
i
,
i
+
1
0,0
=
I
K
⊗
D
1
⊗
I
V
-
,
i
≥
0
,
G
i
,
i
+
1
1,1
=
I
N
+
1
⊗
D
1
⊗
I
V
-
M
,
i
≥
0
.
Here,(i)
I is the identity matrix, and O is a zero matrix of an appropriate dimension;
(ii)
S
d is the diagonal matrix with the diagonal entries equal to the corresponding diagonal entries of the matrix S;
(iii)
⊗ and ⊕ indicate the symbols of Kronecker product and sum of matrices, respectively;
(iv)
W
-
=
W
+
1, V
-
=
V
+
1;
(v)
E
K
+ is the square matrix of size K with all zero entries except the entries (
E
K
+
)
l
,
l
+
1, l
=
0
,
K
-
2
¯, which are equal to 1;
(vi)
E
N
+ is the square matrix of size N
+
1 with all zero entries except the entries (
E
N
+
)
l
,
l
+
1, l
=
0
,
N
-
1
¯, and (
E
N
+
)
N
,
N which are equal to 1;
(vii)
I
k
,
j is the matrix of size k
×
j with all zero entries except the entries (
I
k
,
j
)
l
,
l, l
=
0
,
m
i
n
{
k
-
1
,
j
-
1
}
¯, which are equal to 1;
(viii)
E
- is the square matrix of size N
+
1 with all zero entries except the entries (
E
-
)
l
,
l
-
1, l
=
1
,
N
¯, which are equal to 1;
(ix)
E
~ is the matrix of size K
×
(
N
+
1
) with all zero entries except the entry (
E
~
)
K
-
1
,
K
-
1, which is equal to 1;
(x)
E
-
- is the square matrix of size N
+
1 with all zero entries except the entries (
E
-
-
)
l
,
l
-
1, l
=
K
,
N
¯ which are equal to 1;
(xi)
I
^
m
,
l is the matrix of size m
×
l with all zero entries except the entry (
I
^
m
,
l
)
0,0 which is equal to 1;
(xii)
I
- is the matrix of size (
N
+
1
)
×
K with all zero entries except the entries (
I
-
)
l
,
l, l
=
0
,
K
-
1
¯, which are equal to 1.To illustrate the derivation of the form of the generator blocks, let us explain in brief the form of the blocksG
i
,
i
-
1, i
≥
2. Decrease of the number of customers in the buffer from i to i
-
1 without change of the status 0 of the server (the server is in the sleep mode) can only occur if one of i presenting in the buffer customers leaves the buffer due to impatience. The total intensity of such a transition is equal to i
γ. During an interval of infinitesimal length, only one event from the following ones can occur: customers departure due to impatience, a transition of the underlying process of customers arrival, and a transition of the underlying process of energy arrival. Because we already assumed that departure of a customer takes place, a transition of the underlying process of customers arrival and a transition of the underlying process of energy arrival cannot occur. So, the number of energy units in the buffer and the states of the underlying processes of customers and energy arrival do not change. Thus, the matrix defining the transition probabilities of these processes is equal to I
K
W
-
V
-. As a result of this consideration we obtain formula G
i
,
i
-
1
(
0,0
)
=
i
γ
I
K
W
-
V
-, i
≥
2.Decrease of the number of customers in the buffer fromi to i
-
1 with the change of the status of the server from 0 (server is in the sleep mode) to 1 (server is working) can only occur if i customers were presenting in the buffer, the number of energy units in the system was equal to K
-
1, and a new energy unit arrives in the system. The matrix H
1 defines the transition intensities of the underlying process of energy arrival. The matrix E
~ of size K
×
(
N
+
1
) reflects the change of the number of energy units in the system: this number was equal to K
-
1, and the arrived energy unit was taken for service of a customer which started at this moment. Thus, the number of energy units in the system remains to be equal to K
-
1. Size K
×
(
N
+
1
) of the nonsquare matrix E
~ is explained by the fact that, in the sleep mode, the possible number of energy units in the system varies from 0 to K
-
1 while in the mode when the server works this number varies from 0 to N. Because in the considered case a transition of the underlying process of energy arrival occurs, the underlying process of customers arrival cannot make any transitions. Thus, its transitions are described by the identity matrix I
W
-. Start of customer service is accompanied by the choice of the initial state of the underlying process of service with probabilities defined by the entries of the vector β. Symbol of Kronecker product of matrices is very useful to describe simultaneous transitions of several independent Markov chains. Finally, we obtain formula G
i
,
i
-
1
(
0,1
)
=
E
~
⊗
I
W
-
⊗
H
1
⊗
β, i
≥
1. The order of multipliers in the Kronecker products is induced by the order of corresponding components of the Markov chain ξ
t.Decrease of the number of customers in the buffer fromi to i
-
1 without the change of the current status 1 of the server can occur if one of the customers in the queue leaves the buffer due to impatience (corresponding intensities of transitions are given by the matrix i
γ
I
(
N
+
1
)
W
-
V
-
M) or a customer finishes service and new service starts immediately. The intensities of finishing the service are given by the entries of the vector S
0 while the probabilities defining the choice of the initial state of the underlying process of service are given by the entries of the vector β. Start of the new service causes the decrease of the number of energy units in the system by one which is reflected by the matrix E
-
-. The transitions of the underlying processes of arrival of customers and energy are impossible, so their change is defined by the identity matrix I
W
-
V
-. As a result of this consideration, we obtain formula G
i
,
i
-
1
(
1,1
)
=
i
γ
I
(
N
+
1
)
W
-
V
-
M
+
E
-
-
⊗
I
W
-
V
-
⊗
S
0
β, i
≥
1.Decrease of the number of customers in the buffer fromi to i
-
1 does not imply that the server becomes idle. Thus, the change of status 1 of the server to status 0 is not possible. Eventually, we explained the form of all blocks of the matrix G
i
,
i
-
1. Derivation of the presented above expressions for matrices G
i
,
i and G
i
,
i
+
1 is implemented analogously.Corollary 2.
The Markov chainξ
t, t
≥
0, belongs to the class of continuous-time asymptotically quasi-Toeplitz Markov chains (A
Q
T
M
C); see Klimenok and Dudin [16].As follows from Klimenok and Dudin [16], a sufficient condition for the existence of a stationary distribution of the A
Q
T
M
C
ξ
t, t
≥
0, is expressed in terms of the matrices Y
0, Y
1, and Y
2 defined as follows: (10)
Y
0
=
lim
i
→
∞
⁡
R
i
-
1
G
i
,
i
-
1
,
Y
1
=
lim
i
→
∞
⁡
R
i
-
1
G
i
,
i
+
I
,
Y
2
=
lim
i
→
∞
⁡
R
i
-
1
G
i
,
i
+
1
,where the matrix R
i is a diagonal matrix with the diagonal entries which are defined as the moduli of the corresponding diagonal entries of the matrix G
i
,
i, i
≥
0.It is easy to verify that in the considered case the matricesY
0, Y
1, and Y
2 have the following form: (11)
Y
0
=
I
,
Y
1
=
O
,
Y
2
=
O
.As proven by Klimenok and Dudin [16], the steady-state distribution of the Markov chain ξ
t, t
≥
0, exists if the inequality (12)
y
Y
0
e
>
y
Y
2
eholds true, where the vector y is the unique solution to the system (13)
y
Y
0
+
Y
1
+
Y
2
=
y
,
y
e
=
1
.It is easy to see that ergodicity condition (12) is transformed to the inequality 1
>
0 which is true for all possible values of the system parameters.Thus, the stationary probabilitiesp
(
i
,
r
,
n
,
w
,
v
,
m
), i
≥
0, r
=
0,1, n
=
0
,
N
¯, w
=
0
,
W
¯, v
=
0
,
V
¯, m
=
1
,
M
¯, of the system states exist. Let us form the row vectors p
(
i
,
r
,
n
) of these probabilities enumerated in the lexicographic order of the components w
,
v
,
m. Then let us form the row vectors (14)
p
0,0
=
p
0,0
,
0
,
p
0,0
,
1
,
…
,
p
0,0
,
N
,
p
i
,
0
=
p
i
,
0,0
,
p
i
,
0,1
,
…
,
p
i
,
0
,
K
,
i
>
0
,
p
i
,
1
=
p
i
,
1,0
,
p
i
,
1,1
,
…
,
p
i
,
1
,
N
,
i
≥
0
,
p
i
=
p
i
,
0
,
p
i
,
1
,
i
≥
0
.It is well known that the probability vectorsp
i, i
≥
0, satisfy the following system of linear algebraic equations: (15)
p
0
,
p
1
,
…
,
p
i
,
…
G
=
0
,
p
0
,
p
1
,
…
,
p
i
,
…
e
=
1
,where G is the infinitesimal generator of the Markov chain ξ
t, t
≥
0. To solve this system, the numerically stable algorithm that takes into account that the matrix G has a block-tridiagonal structure, which is presented in Dudina et al. [17], can be proposed.


## 4. Performance Measures

As soon as the vectorsp
i, i
≥
0, have been calculated, we are able to find various performance measures of the system.The average numberN
c
u
s
t
o
m
e
r
s of customers in the buffer is computed by (16)
N
c
u
s
t
o
m
e
r
s
=
∑
i
=
1
∞
i
p
i
e
.The average numberN
e
n
e
r
g
y of units of energy in the buffer is computed by (17)
N
e
n
e
r
g
y
=
∑
i
=
1
∞
∑
n
=
1
K
-
1
n
p
i
,
0
,
n
e
+
∑
n
=
1
N
n
p
0,0
,
n
e
+
∑
i
=
0
∞
∑
n
=
1
N
n
p
i
,
1
,
n
e
.The probabilityP
b
u
s
y where at an arbitrary moment the server is busy is computed by (18)
P
b
u
s
y
=
∑
i
=
0
∞
p
i
,
1
e
.The probabilityP
e
l
o
s
s that an arbitrary unit of energy will be lost is computed by (19)
P
e
l
o
s
s
=
1
λ
e
p
0,0
,
N
I
W
-
⊗
H
1
e
+
∑
i
=
0
∞
p
i
,
1
,
N
I
W
-
⊗
H
1
⊗
I
M
e
.The average intensityλ
o
u
t of flow of customers who receive service is computed by (20)
λ
o
u
t
=
∑
i
=
0
∞
p
i
,
1
e
N
+
1
W
-
V
-
⊗
S
0
.The probabilityP
c
l
o
s
s that an arbitrary customer will be lost is computed by (21)
P
c
l
o
s
s
=
1
-
λ
o
u
t
λ
c
.The probabilityP
t
e
r
m of an arbitrary customer loss due to the termination of its service caused by the lack of energy is computed by (22)
P
t
e
r
m
=
1
λ
c
∑
i
=
0
∞
p
i
,
1,0
I
W
-
V
-
⊗
S
-
S
d
e
.The probabilityP
i
m
p that an arbitrary customer from the buffer will be lost due to impatience is computed by (23)
P
i
m
p
=
P
c
l
o
s
s
-
P
t
e
r
m
.The probabilityP
i
m
m that an arbitrary customer succeeds to get immediate access to the server upon arrival is computed as (24)
P
i
m
m
=
1
λ
c
∑
n
=
K
N
p
0,0
,
n
D
1
⊗
I
V
-
e
.


## 5. Distribution of the Waiting Time of an Arbitrary Customer in the System

We will derive the distribution of an arbitrary customer’s waiting time in terms of the Laplace-Stieltjes transform (L
S
T).To derive the expression for the LSTz
(
s
) of the distribution of an arbitrary customer’s waiting time we use the method of collective marks (method of additional events, method of catastrophes); for references, see, for example, Kesten and Runnenburg [18] and van Danzig [19]. To this end, we interpret the argument s as the intensity of some virtual stationary Poisson flow of catastrophes. Thus, z
(
s
) has the meaning of the probability that no catastrophe arrives during the waiting time of an arbitrary customer. Let us tag an arbitrary customer and keep track of its staying in the system.Letz
(
s
,
l
,
r
,
n
,
v
,
m
) be the probability that a catastrophe will not arrive during the rest of the tagged customer’s waiting time in the system conditioned on the fact that, at the given moment, the position of the tagged customer in the buffer is l, l
≥
1, the state of the server is r, r
=
0,1, the number of units of energy in the buffer is equal to n, n
=
0
,
N
¯, the states of the process v
t are v, and the state of the process m
t is m, t
≥
0.Let us enumerate the probabilitiesz
(
s
,
l
,
r
,
n
,
v
,
m
) in the lexicographic order of components n
,
v
,
m and form the column vectors z
(
s
,
l
,
r
,
n
) and z
(
s
,
l
,
r
) from these probabilities.Theorem 3.
TheL
S
T
z
(
s
) of the distribution of an arbitrary customer’s waiting time is computed by (25)
z
s
=
P
i
m
m
+
λ
c
-
1
∑
i
=
0
∞
∑
n
=
0
K
-
1
p
i
,
0
,
n
D
1
e
⊗
I
V
-
z
s
,
i
+
1,0
,
n
+
∑
i
=
0
∞
p
i
,
1
I
N
+
1
⊗
D
1
e
⊗
I
V
-
M
z
s
,
i
+
1,1
.The proof of the theorem evidently follows from the formula of total probability and probabilistic sense of theL
S
T.The statement of Theorem3 defines the L
S
T
z
(
s
) up to the unknown vectors z
(
s
,
l
,
r
). Let us solve the problem of computation of these vectors.Using the probabilistic sense of theL
S
T it can be shown that the vectors z
(
s
,
l
,
r
) can be found from the following system of linear algebraic equations: (26)
z
s
,
l
,
0
=
s
+
l
γ
I
-
I
K
⊗
H
0
-
1
1
-
δ
l
,
1
E
~
⊗
H
1
⊗
β
z
s
,
l
-
1,1
+
δ
l
,
1
E
~
⊗
H
1
e
+
E
K
+
⊗
H
1
z
s
,
l
,
0
+
l
-
1
γ
z
s
,
l
-
1,0
+
γ
e
,
z
s
,
l
,
1
=
s
+
l
γ
I
-
I
N
+
1
⊗
H
0
⊕
S
d
-
1
E
N
+
⊗
H
1
⊗
I
M
+
E
-
⊗
I
V
-
⊗
S
-
S
d
z
s
,
l
,
1
+
1
-
δ
l
,
1
E
-
-
⊗
I
V
-
⊗
S
0
β
z
s
,
l
-
1,1
+
δ
l
,
1
E
-
-
⊗
I
V
-
⊗
S
0
β
e
+
l
-
1
γ
z
s
,
l
-
1,1
+
γ
e
+
I
^
N
+
1
,
K
⊗
I
V
-
⊗
S
-
S
d
e
+
I
N
+
1
,
K
⊗
I
V
-
⊗
S
0
z
s
,
l
,
0
,
l
>
0
.Here,δ
i
,
j indicates the Kronecker delta.To find the solution to system (26), let us introduce the column vectors (27)
z
s
,
l
=
z
s
,
l
,
0
T
,
z
s
,
l
,
1
T
T
,and rewrite system (26) into the matrix form as (28)
-
s
+
l
γ
I
+
A
z
s
,
l
+
1
-
δ
l
,
1
B
+
l
-
1
γ
I
z
s
,
l
-
1
+
δ
l
,
1
a
+
γ
e
=
0
T
,
l
≥
1
,where (29)
A
=
I
K
⊗
H
0
+
E
K
+
⊗
H
1
O
A
1,0
A
1,1
,
A
1,1
=
I
N
+
1
⊗
H
0
⊕
S
d
+
E
N
+
⊗
H
1
⊗
I
M
+
E
-
⊗
I
V
-
⊗
S
-
S
d
,
A
1,0
=
I
^
N
+
1
,
K
⊗
I
V
-
⊗
S
-
S
d
e
+
I
N
+
1
,
K
⊗
I
V
-
⊗
S
0
,
B
=
O
K
V
-
E
~
⊗
H
1
⊗
β
O
E
-
-
⊗
I
V
-
⊗
S
0
β
,
a
=
E
~
⊗
H
1
e
T
,
E
-
-
⊗
I
V
-
⊗
S
0
e
T
T
.Thus, based on (28), it is easy to verify that the vectors z
(
s
,
l
), l
≥
1, can be recursively computed as (30)
z
s
,
1
=
s
+
γ
I
-
A
-
1
a
+
γ
e
,
z
s
,
l
=
s
+
l
γ
I
-
A
-
1
B
+
l
-
1
γ
I
z
s
,
l
-
1
+
γ
e
,
l
>
1
.Substituting expression (30) for the vectors z
(
s
,
l
), l
≥
1, in (25), it is easy to compute the L
S
T
z
(
s
) and the moments of the waiting time distribution.Corollary 4.
The average waiting timeV
w
a
i
t of an arbitrary customer is calculated as (31)
V
w
a
i
t
=
-
λ
c
-
1
∑
i
=
0
∞
∑
n
=
0
K
-
1
p
i
,
0
,
n
D
1
e
⊗
I
V
-
z
′
s
,
i
+
1,0
,
n
s
=
0
+
∑
i
=
0
∞
p
i
,
1
I
N
+
1
⊗
D
1
e
⊗
I
V
-
M
z
′
s
,
i
+
1,1
s
=
0
.
Here the column vectorsz
′
(
s
,
l
,
r
)
|
s
=
0 and z
′
(
s
,
l
,
r
,
n
)
|
s
=
0 are calculated as the blocks of the vector z
′
(
s
,
l
)
|
s
=
0 which can be calculated as follows: (32)
z
′
s
,
1
s
=
0
=
-
γ
I
-
A
-
2
a
+
γ
e
,
z
′
s
,
l
s
=
0
=
l
γ
I
-
A
-
1
-
e
+
B
+
l
-
1
γ
I
z
′
s
,
l
-
1
s
=
0
,
l
>
1
.


## 6. Numerical Example

To illustrate the effect of variation of the thresholdK and importance of account of the coefficients of correlation and variation in the arrival processes of customers and energy units, let us consider the system with the stationary Poisson arrivals of customers and energy (System 1) and the system with the Markovian arrival processes of customers and energy (System 2) with the same average arrival rates of customers λ
c
=
0.3 and energy λ
e
=
2, but different correlation and variation.We assume that, in the case of System 1, the arrival process of customers is defined by the matricesD
0
=
(
-
0.3
), D
1
=
(
0.3
). The arrival process of energy is defined by the matrices H
0
=
(
-
2
), H
1
=
(
2
). These arrival flows are the stationary Poisson processes with the coefficient of correlation equal to 0 and the coefficient of variation equal to 1.In the case of System 2, we assume that the arrival process of customers is defined by the matrices(33)
D
0
=
-
1.0193
0
0.0003
-
0.033
,
D
1
=
1.0087
0.0106
0.0036
0.0291
.The coefficient of correlation of successive interarrival times is 0.4, and the coefficient of variation of interarrival times is 3.52.The arrival process of energy is defined by the matrices(34)
H
0
=
-
6.7954
0
0.002
-
0.2204
,
H
1
=
6.7246
0.0708
0.0243
0.1941
.The coefficient of correlation of successive interarrival times is0.4, and the coefficient of variation of interarrival times is 3.52.The rest of the system parameters are the same for each system under consideration and are assumed to be as follows:(35)
N
=
25
,
γ
=
0.01
.
P
H service process of customers is characterized by the vector β
=
(
1,0
,
0
) and the matrix (36)
S
=
-
6
6
0
1
-
6
5
0
1
-
6
.This means that service of a customer consists of at least three phases. Duration of each phase has an exponential distribution with parameter 6. After the first phase of service, immediately the second phase starts. After finishing this phase, with probability 5
/
6 the third phase starts, while with the complimentary probability service returns to the first phase, for example, due to not satisfactory implementation of the second phase. After finishing the third phase, with probability 1
/
6 service returns to the second phase while with the complimentary probability service is completed. The mean service time in this service process is equal to 0.646, and the coefficient of variation c
v
a
r is equal to 0.653.Let us vary the thresholdK in the range K
∈
[
1
;
N
].Figures4–6 illustrate the dependence of the loss probability of an arbitrary customer P
c
l
o
s
s, the loss probability P
i
m
p that an arbitrary customer from the buffer will be lost due to impatience, and the loss probability P
t
e
r
m of an arbitrary customer due to the termination of its service in the case of lack of energy on parameter K.Figure 4
Dependence ofP
c
l
o
s
s on K.Figure 5
Dependence ofP
i
m
p on K.Figure 6
Dependence ofP
t
e
r
m on K.The first evident important conclusion following from Figures4–6, as well as from Figures 7–9, is that correlation and variation in the arrival processes, under the same mean arrival rates, drastically change the values of the performance measures of the system. Higher correlation implies poorer performance. This effect is explained as follows. Positive correlation of two random variables implies that the smaller value of one variable causes the smaller value of another variable. The larger value of one variable causes the larger value of another variable. Thus, for the arrival flow with high positive correlation, time intervals, during which customers arrive rarely (and the server starves), alternate with time intervals, during which customers arrive very frequently and the system becomes congested. A lot of customers are lost or have a long waiting time. Note that the negative correlation in arrivals practically has no effect, comparing to the system with renewal arrival process.Figure 7
Dependence ofP
e
l
o
s
s on K.Figure 8
Dependence ofP
i
m
m on K.Figure 9
Dependence ofV
w
a
i
t on K.The second conclusion is that the loss probabilityP
c
l
o
s
s decreases with growth of K, reaches minimum at some point K
∗, and then starts the increase. This is easily explained by the figures for the probabilities P
i
m
p and P
t
e
r
m, the sum of which is equal to P
c
l
o
s
s.The probabilityP
t
e
r
m is very high for small values of K, and then it decreases and becomes very small for K
>
10. This fact stems from the fact that the minimal and the average number of phases, which should be implemented to provide service to a customer, are equal to 3 and 3.84, correspondingly, and the mean number of energy units arriving during one phase is equal to 0.333. Thus, if service to a customer starts at presence of K
=
1,2 or 3 units of energy in the buffer, it is very likely that service will be terminated. When K becomes large enough, loss of a customer due to the service termination becomes almost impossible.Oppositely, the probabilityP
i
m
p increases when K grows because the increase of K implies longer queue of customers and larger number of customers that leave the system due to impatience. Thus, behavior of curves in Figure 4 evidently matches the behavior of curves in Figures 5 and 6.It is worth noting that correlation significantly affects not only the value of the probabilityP
c
l
o
s
s but also the value K
∗ of the threshold K where this probability is minimal. For System 1, K
∗
=
8; P
c
l
o
s
s
=
0.0021. For System 2, K
∗
=
5; P
c
l
o
s
s
=
0.3491. Thus, the problem of optimal choice of the threshold K (as well as numerous other optimization problems related to this model) should be individually solved for any fixed set of the system parameters.It is interesting to compare the optimal value ofP
c
l
o
s
s with its values when more trivial strategies of control are applied. The first strategy is “never use sleep mode.” The server cannot be idle if there are customers in the system and energy in the system is not absent. This strategy corresponds to the choice K
=
1. For this strategy, P
c
l
o
s
s
=
0.007594 for System 1 and P
c
l
o
s
s
=
0.4808 for System 2. The second trivial strategy is “never finish sleep mode until the buffer of energy becomes full.” This strategy corresponds to the choice K
=
N. For this strategy, P
c
l
o
s
s
=
0.01859 for System 1 and P
c
l
o
s
s
=
0.44358 for System 2. It is evidently seen that the optimal choice of K provides essentially better value of P
c
l
o
s
s comparing the trivial strategies.Figures7–9 illustrate the dependence of the loss probability of an arbitrary energy unit P
e
l
o
s
s, the probability P
i
m
m that an arbitrary customer succeeds to start service immediately upon arrival, and the average waiting time V
w
a
i
t of an arbitrary customer on the parameter K.The probability of an arbitrary energy unit lossP
e
l
o
s
s slightly depends on K, but it increases when K becomes small or large. This agrees with Figure 4 because the probability of customer loss P
c
l
o
s
s is also higher when K becomes small or large. Loss of a customer implies that the customer does not go through all required phases of service, so less energy is consumed for its service, the number of energy units in the buffer becomes larger, and, thus, the probability of energy unit loss due to the buffer overflow becomes higher. The probability P
i
m
m monotonically decreases when K grows because an arriving customer has less chances to start service immediately upon arrival due to more strict requirement to the number of energy units in the buffer at the service initiation epoch and more long duration of sleep period. Increase of the average waiting time V
w
a
i
t of an arbitrary customer with growth of K is clear.It should be mentioned that the results of numerous numerical experiments show that the famous Little’s formula holds true for the system under study; that is,(37)
V
w
a
i
t
=
N
c
u
s
t
o
m
e
r
s
λas well as the formula (38)
V
w
a
i
t
=
P
i
m
p
γ
.These formulas as well as the results of computer simulations were used for validation of the presented analytical results and control of their computer implementation.Let us now assume that the arrival process of customers is defined by the matrices(39)
D
0
=
-
0.40549
0
0
-
0.01316
,
D
1
=
0.40279
0.0027
0.00733
0.00583
.The coefficient of correlation of successive interarrival times is 0.2, and the coefficient of variation of interarrival times is 3.51.The arrival process of energy is defined by the matrices(40)
H
0
=
-
2.70328
0
0
-
0.08774
,
H
1
=
2.6853
0.01798
0.04887
0.03887
.The coefficient of correlation of successive interarrival times is 0.2, and the coefficient of variation of interarrival times is 3.51. The rest of parameters is assumed to be the same as above.Table1 contains the values of the key performance measures of the system for different values of the threshold K.Table 1
Dependence of the main performance measures on the thresholdK in the case MAP0.2 + MAP0.2.
K
N
c
u
s
t
o
m
e
r
s
N
e
n
e
r
g
y
P
e
l
o
s
s
λ
o
u
t
P
c
l
o
s
s
P
t
e
r
m
V
w
a
i
t
P
i
m
m
K
=
1
0.86270
17.93
0.49658
0.24477
0.18410
0.1553438
2.88
0.61994
K
=
2
1.23368
17.28
0.47991
0.24992
0.16693
0.1258111
4.11
0.57533
K
=
3
1.55019
16.82
0.46626
0.26980
0.10068
0.0490042
5.17
0.53924
K
=
4
1.71485
16.77
0.46072
0.27378
0.08741
0.0302443
5.72
0.52086
K
=
5
1.83108
16.86
0.45780
0.27811
0.07297
0.0119331
6.10
0.50803
K
=
6
1.90329
17.06
0.45705
0.27897
0.07010
0.0066533
6.34
0.49947
K
=
7
1.96167
17.30
0.45702
0.27956
0.06814
0.0027496
6.54
0.49219
K
=
8
2.01069
17.57
0.45749
0.27945
0.06852
0.0014924
6.70
0.48564
K
=
9
2.05749
17.85
0.45814
0.27923
0.06922
0.0006383
6.86
0.47904
K
=
10
2.10305
18.14
0.45891
0.27886
0.07045
0.0003488
7.01
0.47240
K
=
11
2.14896
18.44
0.45974
0.27846
0.07178
0.0001529
7.16
0.46547
K
=
12
2.19611
18.75
0.46063
0.27801
0.07329
0.0000856
7.32
0.45799
K
=
13
2.24411
19.06
0.46155
0.27755
0.07484
0.0000381
7.48
0.45026
K
=
14
2.29365
19.39
0.46250
0.27706
0.07648
0.0000222
7.65
0.44182
K
=
15
2.34560
19.72
0.46351
0.27654
0.07820
0.0000090
7.82
0.43244
K
=
16
2.39824
20.05
0.46453
0.27602
0.07995
0.0000059
7.99
0.42277
K
=
17
2.45395
20.40
0.46561
0.27546
0.08180
0.0000027
8.18
0.41147
K
=
18
2.51290
20.77
0.46675
0.27487
0.08376
1.62982
E
-
06
8.38
0.39879
K
=
19
2.57221
21.12
0.46790
0.27428
0.08574
7.6309
E
-
07
8.57
0.38551
K
=
20
2.63865
21.51
0.46919
0.27361
0.08796
4.85818
E
-
07
8.80
0.36846
K
=
21
2.70996
21.91
0.47057
0.27290
0.09033
2.17321
E
-
07
9.03
0.34921
K
=
22
2.78556
22.30
0.47204
0.27214
0.09285
1.49155
E
-
07
9.29
0.32809
K
=
23
2.88689
22.75
0.47401
0.27113
0.09623
6.07256
E
-
08
9.62
0.29517
K
=
24
3.01828
23.18
0.47655
0.26982
0.10061
4.00805
E
-
08
10.06
0.25909
K
=
25
3.24017
23.59
0.48086
0.26760
0.10801
2.08683
E
-
08
10.80
0.21669It is easy to see that for this set of the system parameters the optimal value of the threshold isK
∗
=
7, and the minimal value of the loss probability is P
c
l
o
s
s
=
0.06814. It is interesting to observe that the value of the average number N
e
n
e
r
g
y of energy units in the system is not monotone when the threshold K changes. The number N
e
n
e
r
g
y decreases when K increases from 1 to 4 and then it starts the monotone increasing. This phenomenon is easily explained as follows. It was noted above that the average number of phases, which should be implemented to provide service of a customer, is equal to 3.84. Thus, if we start service when 4 energy units are staying in the buffer we have the high risk to lose all these units because service will be not successfully completed. When K is less than 4, the risk to interrupt service is even higher but the number of wasted energy units is less which results in larger number of energy units in the system.In this example the strategy “never use sleep mode” provides valueP
c
l
o
s
s
=
0.18410. The strategy “never finish sleep mode until the buffer of energy becomes full” provides value P
c
l
o
s
s
=
0.10801. It is evident that the optimal choice of K, K
∗
=
7, provides essentially better value of P
c
l
o
s
s
=
0.06814 comparing the trivial strategies.To give more information about the quantitative behavior of the system, let us consider the system under the optimal value in this example of the thresholdK
∗
=
7 and present the values of the probability p
(
i
,
0
)
e that i customers are staying in the buffer and the server is in sleep mode and the probability p
(
i
,
1
)
e that i customers are staying in the buffer and the server provides service. We present here the values of the probabilities only for i
=
0,1
,
…
,
20. Note that we show only the values of the sums of the components of the vectors p
(
i
,
0
) and p
(
i
,
1
), not the vectors themselves, because these vectors have quite high dimension. The size of the vectors p
(
i
,
0
), i
>
0, is equal to 104 while the size of the vectors p
(
0,0
) and p
(
i
,
1
), i
≥
0, is equal to 312 (see Table 2):(41)
∑
i
=
0
20
p
i
e
=
0.9881820310616707
.Table 2
Stationary distribution of the number of customers in the buffer and status of the server.
i
p
(
i
,
0
)
e
p
(
i
,
1
)
e
p
i
e
i
=
0
0.6239317218474825
0.1066112369136978
0.7305429587611788
i
=
1
0.01918369258400953
0.026417593940965396
0.04560128652497489
i
=
2
0.018435165982127483
0.009367941389389267
0.027803107371516724
i
=
3
0.01701956038133314
0.005685178227558685
0.022704738608891818
i
=
4
0.015589704990863436
0.004566368248038898
0.02015607323890231
i
=
5
0.014243214709919589
0.003955497221275474
0.01819871193119508
i
=
6
0.012980401351607156
0.0034751146180821662
0.01645551596968931
i
=
7
0.011793060615562114
0.003053511776530785
0.014846572392092897
i
=
8
0.010676341824961473
0.0026747366348585963
0.013351078459820058
i
=
9
0.009628050888483944
0.002333512545118889
0.011961563433602839
i
=
10
0.00864714408880838
0.002026836425554563
0.010673980514362942
i
=
11
0.007732854199632773
0.0017522708692883516
0.009485125068921121
i
=
12
0.006884290212401231
0.0015075471764451366
0.00839183738884637
i
=
13
0.006100279935315123
0.0012904635015116403
0.007390743436826769
i
=
14
0.005379323609073812
0.0010988627761308045
0.00647818638520462
i
=
15
0.004719595865365016
9.306354403460654
E
-
4
0.005650231305711074
i
=
16
0.004118968305954512
7.837299032657024
E
-
4
0.004902698209220216
i
=
17
0.003575040868797677
6.561646626477985
E
-
4
0.004231205531445474
i
=
18
0.0030851770217298546
5.460398016737795
E
-
4
0.003631216823403636
i
=
19
0.0026465407059894147
4.515470133235514
E
-
4
0.0030980877193129645
i
=
20
0.002256134121983131
3.709778645677343
E
-
4
0.0026271119865508663The distribution of the number of units of energy, which present in the system at an arbitrary moment, is characterized by Table3.Table 3
Distribution of the number of units of energy.
n
Probability thatn units of energy present in the system
n
=
0
0.005257905508569879
n
=
1
0.005509630208742114
n
=
2
0.014726179253861863
n
=
3
0.018781677169906343
n
=
4
0.056170706009369985
n
=
5
0.07240774731564238
n
=
6
0.08488960014753515
n
=
7
0.01537683795036771
n
=
8
0.011736083340386116
n
=
9
0.010285277312655848
n
=
10
0.009934953615597495
n
=
11
0.010251213367679358
n
=
12
0.010513483610080284
n
=
13
0.011259973932302178
n
=
14
0.01229433229048314
n
=
15
0.01286122668417683
n
=
16
0.014603533684450339
n
=
17
0.016238840375660396
n
=
18
0.017457901635544976
n
=
19
0.021154088085968227
n
=
20
0.024371579667587892
n
=
21
0.026382235045579013
n
=
22
0.038933398829995865
n
=
23
0.04079100939440461
n
=
24
0.041553695362992166
n
=
25
0.3962568902004601It should be noted that two maximal values of probabilities are 0.39625689 forn
=
N
=
25 and 0.0848896 for n
=
6.The first maximum may be explained by the following two facts: (i) the mean rate of customers arrival is equal to 0.3 and the average number of necessary energy units for an arbitrary customer service is equal to 3.84, so the mean number of energy units necessary per unit of time is 1.152; (ii) the mean rate of energy arrivals is equal to 2. Thus, the probability that the buffer of energy is full is high. The maximum at pointn
=
6 is related with the fact that the optimal value of the threshold K is 7. Thus, before and after sleep mode completion the number of energy units in the system is equal to 6.In all reported above experiments where we compared performance measures of the systems under different correlation in customers and energy streams, we compared the systems with the same coefficient of correlation for both streams.Let us assume that the arrival flow of customers is the one defined above having the coefficient of correlationc
c
o
r
=
0.2 and we compare the key performance measures of the system with this arrival process and the arrival flows of energy having different coefficients of correlation: c
c
o
r
=
0, c
c
o
r
=
0.2, and c
c
o
r
=
0.4 also defined above. Figures 10–12 illustrate the dependence of the loss probability of an arbitrary customer P
c
l
o
s
s, the probability P
t
e
r
m of an arbitrary customer loss due to the termination of its service caused by lack of energy, and the average waiting time V
w
a
i
t of an arbitrary customer on the parameter K for arrival flows of energy with different coefficients of correlation.Figure 10
Dependence ofP
c
l
o
s
s on K for different arrival flows of energy.Figure 11
Dependence ofP
t
e
r
m on K for different arrival flows of energy.Figure 12
Dependence ofV
w
a
i
t on K for different arrival flows of energy.One can conclude from Figures10–12 that correlation in the energy arrival process has a very essential impact on the system performance measures. In particular, the optimal value of the threshold K is equal to K
∗
=
8 and the minimal value of a customer loss probability is equal to P
c
l
o
s
s
=
0.00903 for the process of energy arrival coded as M
A
P
0, that is, having the coefficient of correlation equal to 0. In turn, for the M
A
P
0.2 energy arrival process, K
∗
=
7 and P
c
l
o
s
s
=
0.06814, and for the M
A
P
0.4 energy arrival process K
∗
=
7 and P
c
l
o
s
s
=
0.20245. Lack of monotonicity of some curves around the point K
=
3 was already mentioned above because the minimal number of energy units required for service of a customers is equal to 3 and the probability P
c
l
o
s
s has a maximum around this point. Thus, the other performance measures also have a bit irregular behavior around this point.Let us assume now that the arrival flow of energy is the one defined above having the coefficient of correlationc
c
o
r
=
0.2 and we compare the key performance measures of the system with this energy arrival process and the arrival flows of customers having different coefficients of correlation: c
c
o
r
=
0, c
c
o
r
=
0.2, and c
c
o
r
=
0.4 also defined above. Figures 13–15 illustrate the dependence of the loss probabilities P
c
l
o
s
s and P
t
e
r
m, and the average waiting time V
w
a
i
t on the parameter K for arrival flows of customers with different coefficients of correlation.Figure 13
Dependence ofP
c
l
o
s
s on K for different arrival flows of customers.Figure 14
Dependence ofP
t
e
r
m on K for different arrival flows of customers.Figure 15
Dependence ofV
w
a
i
t on K for different arrival flows of customers.One can conclude from Figures13–15 that correlation in the customers arrival process has a very essential impact on the system performance measures. In particular, the optimal value of the threshold K is equal to K
∗
=
7 and the minimal value of a customer loss probability is equal to P
c
l
o
s
s
=
0.045256 for the process of customers arrival coded as M
A
P
0
. In turn, for the M
A
P
0.2 customers arrival process K
∗
=
7 and P
c
l
o
s
s
=
0.06814; for the M
A
P
0.4 customers arrival process K
∗
=
8 and P
c
l
o
s
s
=
0.2806184.The purpose of the next experiment is to illustrate the impact of the coefficient of variation of the service time distribution on the system performance measures. Note that the service time distribution has a profound impact on the system performance because it determines not only the mean service time but also the distribution of the number of energy units required for one customer service. Thus, to investigate the impact of the coefficient of variation of the service process let us consider two service processes with the same mean service time and energy consumption for customer service. The first service process, coded asE
x
p, has an exponential distribution and is determined as follows: (42)
S
=
-
1
,
β
=
1
.The second service process, coded as H
E
x
p, has a hyperexponential distribution and is determined as follows: (43)
S
=
-
0.06208
0
0
-
4.882
,
β
=
0.05,0.95
.Both these service processes have the same mean service timeb
1
=
1 and require exactly one energy unit for customer service but have different coefficients of variation. In the case of E
x
p, the coefficient of variation is equal to 1 and in the case of H
E
x
p, coefficient of variation is equal to 5.Let us present the system performance measures for the considered cases andM
A
P
0.2
+
M
A
P
0.2 arrival processes. All the rest of parameters of the system are the same as above. Since the number of energy units required for one customer service in the considered cases is not random and is equal to 1, this is the unique situation when the optimal value of the threshold K
∗ can be predicted without computation and K
∗
=
1.Table4 presents the main performance measures of the system for E
x
p and H
E
x
p service processes in the case K
=
1.Table 4
Main performance measures of the system for the service processes with different coefficients of variation in the caseK
=
1.
N
c
u
s
t
o
m
e
r
s
N
e
n
e
r
g
y
P
e
l
o
s
s
λ
o
u
t
P
c
l
o
s
s
P
t
e
r
m
V
w
a
i
t
P
i
m
m
Exp
0.26745
19.16
0.8513
0.297
0.008915
0
0.8915
0.5939
HExp
2.27217
21.78
0.8585
0.283
0.056665
0
5.66
0.6235It is seen from Table4 that the variance of service times, under the same mean service time and energy consumption, has an essential impact on the system performance measures, in particular on N
c
u
s
t
o
m
e
r
s, P
c
l
o
s
s, and V
w
a
i
t and, thus, this variance has to be carefully taken into account.


## 7. Conclusion

A queueing model with energy harvesting and phase-type distribution of the service time under the threshold strategy of customers access is analyzed analytically. The presented numerical examples demonstrate the feasibility of analytical results and importance of the use of more general models of the arrival processes of customers and energy units than the stationary Poisson processes if the arrival processes in the real system under investigation exhibit correlation or high variability of interarrival times. Essential impact of variability of the service times is also demonstrated. It is numerically established that the famous Little’s formula holds true for the system under study.



---

*Source: 8142743-2016-03-10.xml*