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

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.


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 socalled phase-type (PH) distribution.PH 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 Section 2, 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.

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 Figure 1.
Customers arrive at the system according to the Markovian arrival process.We will code this process as MAP  .Advantage of the MAP 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 MAP) is that the MAP allows taking into account correlation of the successive interarrival times and burstiness typical for modern telecommunication networks.Arrivals in the MAP  are directed by an irreducible continuous-time Markov chain   ,  ≥ 0, with the finite state space {0, 1, . . ., }.The sojourn time of the Markov chain   ,  ≥ 0, in the state  has an exponential distribution with the parameter   ,  = 0, .
Here, notation such as  = 0,  means that  assumes values from the set {0, 1, . . ., }.After this sojourn time expires, with probability   (,   ) the process   transits to the state   and  customers,  = 0,1, arrive at the system.The intensities of transitions from one state to another, which are accompanied by the arrival of  customers, are combined to the square matrices   ,  = 0, 1, of size  + 1.The matrix generating function of these matrices is () =  0 +  1 , || ≤ 1.The matrix ( 1) is an infinitesimal generator of the process   ,  ≥ 0. The stationary distribution vector  of this process satisfies the system of equations (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   (fundamental rate) of the MAP  is defined by   =  1 e.The variance V of customers interarrival time is calculated by the squared coefficient of the variation  var is calculated as and the coefficient of correlation  cor of intervals between successive arrivals is given as For more information about the MAPs 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 process MAP  .Arrivals in the MAP  are directed by an irreducible continuous-time Markov chain V  ,  ≥ 0, with the finite state space {0, 1, . . ., }.The MAP  is defined by the matrices  0 and  1 .Let us denote the average intensity of the MAP  as   .The buffer for energy has a finite capacity .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 a PH distribution with an irreducible representation (, ).This service time can be interpreted as the time until the underlying Markov process   ,  ≥ 0, with the finite state space {1, . . ., ,  + 1}, reaches the single absorbing state  + 1 conditioned on the fact that the initial state of this process is selected among the states {1, . . ., } according to the probabilistic row vector  = ( 1 , . . .,   ).The transition rates of the process   within the set {1, . . ., } are defined by the subgenerator , and the transition rates into the absorbing state (which leads to service completion) are given by the entries of the column vector S 0 = −e.The mean service time is calculated as  1 = (−) −1 e.
The class of PH distributions is dense (in the sense of weak convergence) in the set of all probability distributions of nonnegative variables.Thus, the PH 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 the MAP and PH 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 MAP and PH 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 threshold , 1 ≤  ≤ .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 , this customer occupies the server (see Figure 2).
Otherwise the customer joins the buffer (see Figure 3).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 threshold .
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 threshold  is a control parameter.Note that if the value of the threshold  is too small, many customers will be lost during service due to the lack of energy.If the value of the threshold  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  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  difficult.In some particular cases, this problem becomes a bit easier.For example, if PH distribution is the Erlangian of order , then the choice of  such that  >  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  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 threshold , 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.

Process of System States
It is easy to see that the behavior of the system under study is described by the following regular irreducible continuoustime Markov chain: where, during the epoch ,  ≥ 0, (i)   is the number of customers in the buffer,   ≥ 0; (ii)   is an indicator that indicates whether the server works or not:   = 0 corresponds to the case when the server does not work and   = 1 if the server works; (iii)   is the number of energy units in the buffer,   = 0, ; (iv)   is the state of the underlying process of the MAP  ,   = 0, ; (v) V  is the state of the underlying process of the MAP  , V  = 0, ; (vi)   is the state of the underlying process of the PH service process,   = 1, .
Analyzing all possible transitions of the Markov chain   ,  ≥ 0, during an interval of infinitesimal length and rewriting the intensities of these transitions in block matrix form we obtain the following result.
To illustrate the derivation of the form of the generator blocks, let us explain in brief the form of the blocks  ,−1 ,  ≥ 2. Decrease of the number of customers in the buffer from  to  − 1 without change of the status 0 of the server (the server is in the sleep mode) can only occur if one of  presenting in the buffer customers leaves the buffer due to impatience.The total intensity of such a transition is equal to .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    .As a result of this consideration we obtain formula  (0,0) ,−1 =    ,  ≥ 2. Decrease of the number of customers in the buffer from  to  − 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  customers were presenting in the buffer, the number of energy units in the system was equal to  − 1, and a new energy unit arrives in the system.The matrix  1 defines the transition intensities of the underlying process of energy arrival.The matrix Ẽ of size  × ( + 1) reflects the change of Mathematical Problems in Engineering the number of energy units in the system: this number was equal to  − 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  − 1. Size  × ( + 1) of the nonsquare matrix Ẽ is explained by the fact that, in the sleep mode, the possible number of energy units in the system varies from 0 to  − 1 while in the mode when the server works this number varies from 0 to .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   .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  (0,1) ,−1 = Ẽ ⊗   ⊗  1 ⊗ ,  ≥ 1.The order of multipliers in the Kronecker products is induced by the order of corresponding components of the Markov chain   .
Decrease of the number of customers in the buffer from  to  − 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  (+1)  ) 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  − .The transitions of the underlying processes of arrival of customers and energy are impossible, so their change is defined by the identity matrix    .As a result of this consideration, we obtain formula Decrease of the number of customers in the buffer from  to  − 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  ,−1 .Derivation of the presented above expressions for matrices  , and  ,+1 is implemented analogously.[16].

Corollary 2. The Markov chain 𝜉 𝑡 , 𝑡 ≥ 0, belongs to the class of continuous-time asymptotically quasi-Toeplitz Markov chains (𝐴𝑄𝑇𝑀𝐶); see Klimenok and Dudin
As follows from Klimenok and Dudin [16], a sufficient condition for the existence of a stationary distribution of the AQTMC   ,  ≥ 0, is expressed in terms of the matrices  0 ,  1 , and  2 defined as follows: where the matrix   is a diagonal matrix with the diagonal entries which are defined as the moduli of the corresponding diagonal entries of the matrix  , ,  ≥ 0.

Performance Measures
As soon as the vectors p  ,  ≥ 0, have been calculated, we are able to find various performance measures of the system.
The average number  customers of customers in the buffer is computed by The average number  energy of units of energy in the buffer is computed by The probability  busy where at an arbitrary moment the server is busy is computed by The probability  loss  that an arbitrary unit of energy will be lost is computed by The average intensity  out of flow of customers who receive service is computed by The probability  loss  that an arbitrary customer will be lost is computed by The probability  term of an arbitrary customer loss due to the termination of its service caused by the lack of energy is computed by The probability  imp that an arbitrary customer from the buffer will be lost due to impatience is computed by The probability  imm that an arbitrary customer succeeds to get immediate access to the server upon arrival is computed as

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 (LST).
To derive the expression for the LST () 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  as the intensity of some virtual stationary Poisson flow of catastrophes.Thus, () 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.
Let (, , , , V, ) 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 ,  ≥ 1, the state of the server is ,  = 0, 1, the number of units of energy in the buffer is equal to ,  = 0, , the states of the process V  are V, and the state of the process   is ,  ≥ 0.
Let us enumerate the probabilities (, , , , V, ) in the lexicographic order of components , V,  and form the column vectors z(, , , ) and z(, , ) from these probabilities.

Theorem 3. The 𝐿𝑆𝑇 𝑧(𝑠) of the distribution of an arbitrary customer's waiting time is computed by
The proof of the theorem evidently follows from the formula of total probability and probabilistic sense of the LST.
The statement of Theorem 3 defines the LST () up to the unknown vectors z(, , ).Let us solve the problem of computation of these vectors.
Using the probabilistic sense of the LST it can be shown that the vectors z(, , ) can be found from the following system of linear algebraic equations: Here,  , indicates the Kronecker delta.

Numerical Example
To illustrate the effect of variation of the threshold  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   = 0.3 and energy   = 2, but different correlation and variation.
We assume that, in the case of System 1, the arrival process of customers is defined by the matrices  0 = (−0.3), 1 = (0.3).The arrival process of energy is defined by the matrices  0 = (−2),  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 ( 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  0 = ( −6.7954 0 0.002 −0.2204 ) ,  1 = ( 6.7246 0.0708 0.0243 0.1941 ) . (34) The coefficient of correlation of successive interarrival times is 0.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: PH service process of customers is characterized by the vector  = (1, 0, 0) and the matrix 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  var is equal to 0.653.Let us vary the threshold  in the range  ∈ [1; ].Figures 4-6 illustrate the dependence of the loss probability of an arbitrary customer  loss  , the loss probability  imp that an arbitrary customer from the buffer will be lost due to impatience, and the loss probability  term of an arbitrary customer due to the termination of its service in the case of lack of energy on parameter .
The first evident important conclusion following from Figures 4-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.
The second conclusion is that the loss probability  loss  decreases with growth of , reaches minimum at some point  * , and then starts the increase.This is easily explained by the figures for the probabilities  imp and  term , the sum of which is equal to  loss  .The probability  term is very high for small values of , and then it decreases and becomes very small for  > 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  = 1, 2 or 3 units of energy in the buffer, it is very likely that service will be terminated.When  becomes large enough, loss of a customer due to the service termination becomes almost impossible.Oppositely, the probability  imp increases when  grows because the increase of  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 probability  loss  but also the value  * of the threshold  where this probability is minimal.For System 1,  * = 8;  loss  = 0.0021.For System 2,  * = 5;  loss  = 0.3491.Thus, the problem of optimal choice of the threshold  (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 of  loss  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  = 1.For this strategy,  loss  = 0.007594 for System 1 and  loss  = 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  = .For this strategy,  loss  = 0.01859 for System 1 and  loss  = 0.44358 for System 2. It is evidently seen that the optimal choice of  provides essentially better value of  loss  comparing the trivial strategies.
Figures 7-9 illustrate the dependence of the loss probability of an arbitrary energy unit  loss  , the probability  imm that an arbitrary customer succeeds to start service immediately upon arrival, and the average waiting time  wait of an arbitrary customer on the parameter .
The probability of an arbitrary energy unit loss  loss  slightly depends on , but it increases when  becomes small or large.This agrees with Figure 4 because the probability of customer loss  loss  is also higher when  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  imm monotonically decreases when  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  wait of an arbitrary customer with growth of  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, as well as the formula 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) 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 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.
Table 1 contains the values of the key performance measures of the system for different values of the threshold .
It is easy to see that for this set of the system parameters the optimal value of the threshold is  * = 7, and the minimal value of the loss probability is  loss  = 0.06814.It is interesting to observe that the value of the average number  energy of energy units in the system is not monotone when the threshold  changes.The number  energy decreases when  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  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 value  loss  = 0.18410.The strategy "never finish sleep mode until the buffer of energy becomes full" provides value  loss  = 0.10801.It is evident that the optimal choice of ,  * = 7, provides essentially better value of  loss  = 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 threshold  * = 7 and present the values of the probability p(, 0)e that  customers are staying in the buffer and the server is in sleep mode and the probability p(, 1)e that  customers are staying in the buffer and the server provides service.We present here the values of the probabilities only for  = 0, 1, . . ., 20.Note that we show only the values of the sums of the components of the vectors p(, 0) and p(, 1), not the vectors themselves, because these vectors have quite high dimension.The size of the vectors p(, 0),  > 0, is equal to 104 while the size of the vectors p(0, 0) and p(, 1),  ≥ 0, is equal to 312 (see Table 2): (41) The distribution of the number of units of energy, which present in the system at an arbitrary moment, is characterized by Table 3.
It should be noted that two maximal values of probabilities are 0.39625689 for  =  = 25 and 0.0848896 for  = 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 point  = 6 is related with the fact that the optimal value of the threshold  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 correlation  cor = 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:  cor = 0,  cor = 0.2, and  cor = 0.4 also defined above.Figures 10-12 illustrate the dependence of the loss probability of an arbitrary customer  loss  , the probability  term of an arbitrary customer loss due to the termination of its service caused by lack of energy, and the average waiting time  wait of an arbitrary customer on the parameter  for arrival flows of energy with different coefficients of correlation.One can conclude from Figures 10-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  is equal to  * = 8 and the minimal value of a customer loss probability is equal to  loss  = 0.00903 for the process of energy arrival coded as MAP 0 , that is, having the coefficient of correlation equal to 0. In turn, for the MAP 0.2 energy arrival process,  * = 7 and  loss  = 0.06814, and for the MAP 0.4 energy arrival process  * = 7 and  loss  = 0.20245.Lack of monotonicity of some curves around the point  = 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  loss  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 correlation  cor = 0.2 and we compare the key performance measures  One can conclude from Figures 13-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  is equal to  * = 7 and the minimal value of a customer loss probability is equal to  loss  = 0.045256 for the process of customers arrival coded as MAP 0 .In turn, for the MAP 0.2 customers arrival process  * = 7   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.(43) Both these service processes have the same mean service time  1 = 1 and require exactly one energy unit for customer service but have different coefficients of variation.In the case of Exp, the coefficient of variation is equal to 1 and in the case of HExp, coefficient of variation is equal to 5.
Let us present the system performance measures for the considered cases and MAP 0.2 + MAP 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  * can be predicted without computation and  * = 1.
Table 4 presents the main performance measures of the system for Exp and HExp service processes in the case  = 1.
It is seen from Table 4 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  customers ,  loss  , and  wait and, thus, this variance has to be carefully taken into account.

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.

Figure 2 :
Figure 2: Arriving customer occupies the server.

Figure 9 :
Figure 9: Dependence of  wait on .

Figure 10 :
Figure 10: Dependence of  loss  on  for different arrival flows of energy.

Figure 11 :
Figure 11: Dependence of  term on  for different arrival flows of energy.

Figure 12 :
Figure 12: Dependence of  wait on  for different arrival flows of energy.

Figure 13 :
Figure 13: Dependence of  loss  on  for different arrival flows of customers.

Figure 14 :
Figure 14: Dependence of  term on  for different arrival flows of customers.

Figure 15 :
Figure 15: Dependence of  wait on  for different arrival flows of customers.
, is the matrix of size  ×  with all zero entries except the entries ( , ) , ,  = 0, min{ − 1,  − 1}, which are equal to 1;(viii)  − is the square matrix of size +1 with all zero entries except the entries ( − ) ,−1 ,  = 1, , which are equal to 1; i)  is the identity matrix, and  is a zero matrix of an appropriate dimension; (ii)   is the diagonal matrix with the diagonal entries equal to the corresponding diagonal entries of the matrix ; (iii) ⊗ and ⊕ indicate the symbols of Kronecker product and sum of matrices, respectively; (ix) Ẽ is the matrix of size  × ( + 1) with all zero entries except the entry ( Ẽ) −1,−1 , which is equal to 1; (x)  − is the square matrix of size +1 with all zero entries except the entries ( − ) ,−1 ,  = ,  which are equal to 1;

Table 2 :
Stationary distribution of the number of customers in the buffer and status of the server.

Table 3 :
Distribution of the number of units of energy.

Table 4 :
Main performance measures of the system for the service processes with different coefficients of variation in the case  = 1.