Optimal Hysteresis Control via a Queuing System with Two Heterogeneous Energy-Consuming Servers

: A queuing system having two different servers is under study. Demands enter the system according to a Markov arrival process. Service times have phase-type distribution. Service of demands is possible only if the ﬁxed number of energy units, probably different for two servers, is available in the system at the potential service beginning moment. Energy units arrive in the system also according to a Markov arrival process and are stored in a stock (battery) of a ﬁnite capacity. Leakage of energy units from the stock can occur. Demands waiting in the inﬁnite buffer are impatient and can leave the buffer after an exponentially distributed waiting time. One server is the main one and permanently provides service when the buffer is not empty and the required number of energy units is available. The second server is the assistant server and is switched on or off depending on the availability of energy units and queue length according to the hysteresis strategy deﬁned by two thresholds. The assistant server is switched on when the queue length is not less than the greater threshold and is switched off when the queue length becomes smaller than the smaller threshold. The use of the assistant server has to be paid. Thus, the problem of the optimal selection of the thresholds deﬁning the control strategy naturally arises. To solve this problem, the study of the behavior of the system under any ﬁxed values of the parameters of the control strategy is necessary. Such a study is given in this paper. Numerical results are presented. They illustrate the feasibility of computer realization of the developed algorithms for computation of the stationary distribution of the system states and the main key performance indicators as well as the result of solving one of the possible optimization tasks.


Introduction
This paper contains the analysis of a queuing system that may be potentially useful for the optimal design of a variety of real-world systems.In brief, besides the consideration of arrival processes, which are essentially more general than the stationary Poisson, and service times having become more complicated than those popular in the literature's exponential distribution, the distinguishing features of the considered system are (i) the necessity of the use of energy units (e.u. for brevity) for the service of demands; (ii) the generation, consumption, and leakage of e.u.during the system operation with their temporal accumulation stored in the battery (stock) with limited capacity; (iii) the availability of two heterogeneous service devices having different parameters of service time distribution and different consumption of e.u.; (iv) the opportunity of the dynamic activation and deactivation of the server, which consumes more energy according to the hysteresis strategy, depending on the current queue length and amount of available energy; (v) the account of impatience pertaining to waiting demands.The simultaneous presence of all the listed features guarantees the novelty of the analyzed system.A very brief overview of the literature about the systems characterized by some of the above features is as follows.
1.1.Short Literature Overview 1.1.1.Single-Server Queues with Variable Service Rates The simplest kind of management of the service or arrival rate and the number of active service devices in queuing models is the threshold-type strategy.It is defined by a number called a threshold.Parameters of the system operation are changed when the number of demands in the system crosses up or down this threshold.It is considered quite a general single-server system with the batch Markov Arrival Process (MAP), the general service time distribution, and the multi-threshold control by the service time all being considered in [1].A disadvantage of the threshold strategy of control is the possibility of many changes in the service rate (oscillation) when the number of demands in the system admits values around the threshold.This may be undesirable in situations when the change in the rate is charged or is implemented non-instantaneously.The hysteresis-type control suggests that a threshold is replaced with two thresholds.The use of the higher service rate begins when the queue length increases up to the larger threshold.The use of the slower service rate begins when the queue length drops below the smaller threshold.When the length of the queue remains between two thresholds, no change in the current service rate is performed.Such a type of control gives a possibility to reduce the number of rate switching compared to the threshold-type control.The hysteresis-type control using the queuing system operation was introduced in 1960 and is often applied to different queuing systems.As early works in this direction, the following papers can be mentioned: [2][3][4].More references can be found, e.g., in [5][6][7][8] and the references therein.

Multi-Server Queues with a Variable Number of Homogeneous Active Service Devices
In the context of control via the number of active service devices, the following results can be mentioned.In [9], a two-service device queue with identical service devices and a threshold control of activation of an additional service device is analyzed.In [10,11], a more general system with many identical service devices and hysteresis control via the number of active service devices is considered.A similar model but with significantly less restrictive assumptions about the arrival and service processes was investigated in [12].A more detailed survey of the relevant research is given there.In [13], hysteresis-type control via service device reservation was applied to a multi-server priority queuing system describing the operation of a cognitive radio system.

Multi-Server Queues with a Variable Number of Heterogeneous Service Devices
The all above mentioned multi-server queues have identical service devices.The analysis of systems with distinct servers is essentially more complicated than the analysis of systems with identical servers.This is evidently explained as follows.In the case of identical servers, the service process in the system at an arbitrary moment is defined by the number of active servers.When service time distribution has exponential distribution, the number of active servers completely defines the total service rate at a given moment.In the case of non-identical servers, the service process at an arbitrary moment is not defined by the number of active servers.Keeping track of which servers namely are active and which are idle is required to build the Markov chain describing the behavior of the system.Consideration of a multi-dimensional process (number of active servers and indicators of which servers are active) instead of a one-dimensional process (just the number of active servers) is clearly a more difficult task.
The situation in the study of systems with distinct service devices is as follows.In [14], a system with two heterogeneous service devices and a threshold policy of control via the number of used service devices is considered.The practical motivation of model consideration is given and the invariant distribution of the system states is calculated.In [15], the optimality of the threshold-type strategy, which is used in [14], is proven.Nonoptimality of such a strategy in the case of more than two service devices is forecasted.In [16], the task of finding the optimal thresholds defining the hysteresis-type strategy for the two-server system providing service to the demands arriving in batches is solved.In [17], the optimality of the threshold-type control in a system with many heterogeneous service devices is proven and a numerical study of the model is implemented.The results from [18][19][20][21] are obtained as further development of the results from [17].In [11,22], the application of the hysteresis strategy in a system with many heterogeneous service devices is also discussed.

Queuing Systems with Energy Consumption and Harvesting
In the above mentioned papers, no aspects related to the consumption of the e.u.during service of demands are touched.Meanwhile, this aspect is very important in many real-world systems because more e.u. can be required to provide a higher service rate and it is desirable to reduce its use due to economical and ecological reasons.In particular, this aspect is very important in the management of the operation of wireless sensor networks.Sensors designed for monitoring the state of some objects can be installed in hard-toreach areas and they can collect and transfer information only if they succeed in receiving (harvesting) enough e.u.The e.u. can be earned from sources, such as solar panels, and is accumulated in a small battery whose low charge makes it impossible for the transmission of information.Such situations are considered using queuing models with energy harvesting, see, e.g., [23][24][25][26][27][28][29].In such models, it is assumed that e.u.arrive in the system during its operation according to a stationary Poisson process.E.u. are stored in a finite stock and service of any demand can start only in the presence of a minimal required number of e.u. in the stock.Similar queues have been studied in the literature for about seventy years and they are called assembly-like systems; queues with paired demands; taxi/passenger queues; matching systems; double-ended queues; systems with passive service devices; systems with demands' arrival in sessions, etc. Queuing/inventory models (see, e.g., [30][31][32]) very close to such models are popular now.The distinction is that these models usually assume some standard replenishment policies for increasing the number of e.u.(inventory units), but not the arrival of these units in some random input flow.Aspects relating to the service rate control, depending on the amount of the harvested e.u., are considered very rare and, as such, they deserve attention.We can mention the paper [33] where the threshold strategy for service rate selection depending on the available amount of e.u. is analyzed.When the number of available e.u.increases (or decreases), higher (or lower) service rates can be used.

Queuing Systems with Demand Impatience
An important peculiarity of the model under study is that it takes into account the impatience of demands.Such impatience is typical in many real-world systems including communications systems, call centers, etc. Impatience means that a demand waiting for service in a buffer can renege from the system if the service is not started during some period.In the context of sensor networks, impatience is caused by the useless transmission of a signal to the coordinating node of the network due to the obsolescence of information contained in this signal.Lists of related research on queues with impatient demands can be found, e.g., in recent papers [34][35][36][37][38][39].Account of the impatience frequently causes essential complications in the analysis because the behavior of the random process describing the dynamics of the system becomes level-dependent and the corresponding process does not belong to the classes of QBD processes or M/G/1and GI/M/1-type structures that Markov chains were used to exhaustively investigate by M. Neuts, see [40,41].

Contributions of the Paper
As far as we know, this paper is the first one where the system with e.u.harvesting and heterogeneous service devices are considered.Two service devices differ by the service rate and the number of e.u.required for service providing.We suggest the hysteresis-type control via the activation and deactivation of the service device that consumes more energy.This service device is activated if the queue length is not below the larger threshold and enough e.u. is available.The service device is deactivated if the queue length becomes less than the smaller threshold or the number of available e.u. is not sufficient to begin service.Arrivals of demands and e.u. are described by the MAPs, which allows us to catch the possible correlation and high variability of inter-arrival times typical for many real-world systems.We assume the service times by both service devices have phase-type (PH) distribution that allows it to fit well with an arbitrary distribution, see [42], or, at least, several initial moments of service time distribution.Possible leakage of e.u. from the stock and the reneging of waiting demands from the buffer are taken into account.
The dynamic of the system is determined via a multidimensional continuous-time Markov chain that is a level-dependent QBD (Quasi-Birth-and-Death) process.The generator of this QBD is written down as the block matrix with sizes of blocks depending on the relation of the level of the chain with the thresholds of the control policy and the relation of the component of the chain describing both the number of e.u. and the number required for a demand service.Ergodicity conditions of the chain in the cases when demands are impatient and absolutely patient are presented and the invariant distribution of the chain is computed.Expressions for the computation of various performance measures are derived.The feasibility of the obtained algorithms for the computation of invariant distribution and performance measures is illustrated by several examples.

Possible Applications of the Model
Additional items required for demands' service, which are interpreted in the paper as e.u., may have a meaning of any other resource, e.g., computer memory, service assistants, any other kind of resources (workpiece, place in a warehouse, service device, money, etc.).Therefore, the analyzed model can be used for the performance evaluation and optimization of operation of any real-world object where service of demand requires the use, besides the service device, of some additional resources and available service devices require a different amount of resources, provide different service rates, and have a different cost.In particular, the model fits for the description of the work of the node of a sensor e.u.harvesting network, which has the main, low cost, and low e.u.consuming, transmitter, and reserve transmitter that should be used only in situations when there is enough energy, and there is also a risk of information loss due to its obsolescence.The results can be used for managing the system with the use of the best values of the thresholds of control policy by verifying the ability of the system to provide the required quality of service or the desired revenue in the case of the optimal choice of control strategy.

Structure of the Text
The mathematical model under study is fully described in Section 2. A multidimensional level-dependent Markov chain (MC) with floating cardinality of blocks that determines the dynamics of the system is introduced in Section 3. The infinitesimal generator of this chain is presented.The ergodicity condition for the system and the procedure for the computation of its stationary distribution in the case of absolute patient demands are presented in Section 4. Analogous results in the case of impatient demands are given in Section 5. Expressions for the computation of the performance indicators of the system based on the invariant distribution of MC are reported in Section 6. Section 7 deals with the results of some numerical experiments and a short illustration of the possibility of a solution to the optimization problem.Appendix A contains a brief explanation of the blocks of the generator.

Mathematical Model
We consider a queuing model having two heterogeneous service devices and an infinite buffer for temporal storing of the demands.The scheme of the system operation is displayed in Figure 1.We assume that the MAP flow of demands arrives for service in the system.This flow is defined by the underlying process ν where 0 = (0, . . ., 0), e = (1, . . ., 1) T , and symbol T denotes the operation of transposition of the vector.
More information about the MAP, its characteristics, particular cases, and importance for an adequate description of flows in real-life systems can be found, e.g., in [43][44][45][46].The key benefit of the MAP over the widely used stationary Poisson process model is the possibility to fit real flow with any value of the coefficients of variation and correlation of successive inter-arrival times.The stationary Poisson process does not give any freedom in fitting, except on account of the mean arrival rate.The coefficients of variation and correlation of inter-arrival times in the stationary Poisson process are equal to 1 and 0, correspondingly.
If service devices are unavailable at the arrival moment, the demand joins the buffer of an unlimited capacity and will be processed later on.Otherwise, the demand may try to start the receiving service.
The start of service of any demand requires the presence of a fixed number of e.u.
E.u. arrive in the MAP flow defined by the underlying process ν (e) t , with the state space {1, . . ., W (e) } and matrices being D The e.u. are stored in the stock (battery) of a finite capacity N. If the stock is full upon the new e.u.arrival moment, this e.u. is permanently lost.Otherwise, it joins the stock.Leakage of e.u. from the stock occurs with rate γ n when n e.u.stay in the stock, n = 1, N. We assume that γ 0 = 0.
Service to demands can be provided via two service devices.The rth service device requires n (r) e.u. for a demand service.We assume that n (2) > n (1) , i.e., the second service device consumes more e.u.than the first one.The contrary case can be analyzed analogously.
Service time via the rth service device has the PH distribution defined by the underly- t with the state space being {1, . . ., M (r) }, r = 1, 2, of transient states.The state M (r) + 1 is the absorbing one.It is said that such service time has an irreducible represen- tation (β (r) , S (r) ).Transition rates to the absorbing state are defined by the column vector S (r) 0 = −S (r) e.The service rate is denoted by µ (r) = (β (r) (−S (r) ) −1 e) −1 , r = 1, 2. More useful information about the PH distribution and its properties can be found, e.g., in [40].
If the required service beginning number of e.u. is unavailable, the service device suspends the work until the required number of e.u. will be accumulated in the stock.After that, the service device may immediately begin service.At the moment of the service beginning, the used service e.u.immediately disappear from the stock.
We consider the first service device as the main one.It cannot stay idle if the buffer is not empty and at least n (1) e.u. are available in the stock.The second service device is the reserved, assistant service device.The mechanism of its activation or deactivation is as follows.The thresholds j and k are assumed to be fixed, 1 ≤ j ≤ k.An idle second service device is switched on if two conditions are simultaneously fulfilled: the number of demands in the queue is not less than k and at least n (2) e.u. are available.Interruption of service is not allowed.The second service device is switched off at the moment of service completion if (i) the queue length is below the value j or (and) (ii) n (2) e.u.required for service of the next demand are not available.
The switching-on of the idle second service device occurs when the queue length is not less than k and the number of available e.u. is not less than n (2) .In other words, there are two scenarios of switching on (activating) the idle second service device: (i) there were not less than n (2) e.u., where the number of demands in the system was equal to k − 1 and a new demand arrives; (ii) there were not less than k demands in the system, where there were n (2) − 1 e.u. and new e.u.arrives.We will suppose that n (1) + n (2) < N. The contrary case of a small buffer capacity N can be analyzed similarly to the analysis implemented here.
Demands residing in the buffer are, generally speaking, impatient.The patience time of a demand is exponentially distributed with rate ε.After this time finishes, the demand departs from the buffer, independently of other demands, and is considered lost.
We aim to analyze the behavior of this system and solve the optimization problem: to find thresholds j and k providing the minimum to the cost function.This function includes the charge for demands waiting in the buffer and the charges paid for both the use of the assistant service device during a unit of time and each activation of this device.

The Process Describing the Dynamics of the System and Its Generator
The dynamics of the described system are determined via the multidimensional MC with continuous time where at time t, t ≥ 0, • i t is the number of demands in the buffer, i t ≥ 0; • n t is the number of e.u. in the stock, n t = 0, N; t is the state of the underlying process of the MAP of demands, ν t is the state of the underlying process of the MAP of e.u., ν • vector process η t defines the states of the underlying processes of service in the devices.We will distinguish four macro-states of the process η t : the first macro-state means that both service devices are idle; the second macro-state contains the states (η t , 0) when the first service device is busy while the second service device is idle; the third macro-state contains the states (0, η (2) t ) when the second service device is busy while the first one is idle; and the fourth macro-state contains the states (η (1) t , η (2) t ) when both service devices are busy.Here, η (r) t denotes the current state of the underlying process of service in the rth service device, η (r) Here and below the notation n = 0, N means that the parameter n admits the values from the set {0, 1, . . ., N}.
Note that the state space of the vector process η t defined by Formula (1) depends on the values of the components {i t , n t } as follows.
If i t = 0 or i t > 0 and n t = 0, n (1) − 1, then all macro-states of the process η t are feasible and the cardinality of its state space is equal to M = (M (1) + 1)(M (2) + 1).
If i t > 0 and n t = n (1) , N, then the first and the third macro-states 0 of the process η t are not feasible because the first service device cannot stay idle when the queue is not empty and not less than n (1) e.u. are available.Correspondingly, the cardinality of the state space of the process η t is equal to M = M (1) (M (2) + 1).
Finally, if n t = n (2) , N, then for i t ≥ k, the second macro-state (η t , 0) of the process η t is also not feasible and the cardinality of the state space of the process η t is equal to M = M (1) M (2) .In the following, when we say that the process η t currently has four macro- states, we mean all four possible macro-states described above.When we say that the process η t has two macro-states, we mean the second and fourth macro-states.When we say that the process η t has one macro-state, we mean the fourth macro-state.
To simplify the formulas, we will use the following notation: is the elements' number of the state space of the components {ν 1 ⊗ I W 2 ; • ⊗ and ⊕ mean the operations of the Kronecker product and the sum for the matrices, see, e.g., [47]; • diag{. . .} is the diagonal matrix with the diagonal entries given in the brackets; describing the rates of the vector process η t transition when all of its states are feasible and which do not cause service completions; , describing the rates of the vector process η t transition when all of its states are feasible and which cause service completion in the first service device; describing the rates of the vector process η t transition when all of its states are feasible and which cause service completion in the second service device; describing the probabilities of the vector process η t transition when all of its states are feasible and service begins in the first service device; describing the probabilities of the vector process η t transition when all of its states are feasible and service begins in the second service device; • matrices A r,r defining the transitions of the vector process η t when it transits from the group of the states consisting of r macro-states to the group of the states consisting of r macro-states: • matrices B r,r defining the transitions of the vector process η t when service by the second service device starts and the process η t transits from the group of the states consisting of r macro-states to the group of the states consisting of r macro-states: r,r defining the transitions of the vector process η t when the first service device completes service and starts the new one and the process η t transits from the group of the states consisting of r macro-states to the group of the states consisting of r macro-states: Z (1) 0 β (1) 2,2 = S (1) 0 β (1) 0 β (1) 0 β (1) 0 β (1) r,r defining the transitions of the vector process η t when service by the second service device finishes service and starts the new one and the process η t transits from the group of the states consisting of r macro-states to the group of the states consisting of r macro-states: 2,2 defining the transitions of the vector process η t when service by the second service device finishes service and the new service does not start, so the process η t continues to stay in the group of the states consisting of two macro-states: Let also define the matrices F i,i , i, i ≥ 0, having blocks F , n, n = 0, N, defined by the following formulas: • square matrix F 0,0 of size (N + 1)W M has: the diagonal blocks given by the up-diagonal blocks F (n,n+1) 0,0 given by (n,n−n (1) ) 0,0 , n = n (1) , N given by F (n,n−n (1) ) 0,0 Remark 1.If n (1) = 1, we presented above two different expressions for the matrix F (n,n−1) 0,0 (formulas for the matrix F (n,n−1) 0,0 and F (n,n−n (1) ) 0,0 ).This happens because the number decrease in e.u. can occur either due to the leakage of e.u. or due to the service of a demand.In this and other similar cases, we merge two formulas into one formula via the summation of right-hand sides: • matrix F 0,1 has only diagonal blocks.
If k > 1, this matrix is defined by If k = 1, the formula is changed to are defined as follows.
If n = 0, n (1) − 1, then for all i ≥ 1 If n = n (1) , n (2) − 1, then for all i ≥ 1 If n = n (2) , N, then for 1 ≤ i < j, we have and for i ≥ k, we have The sub-diagonal blocks F For i ≥ k, these blocks have the form The up-diagonal blocks F (n,n+1) i,i for i = 1, k − 1 have the following form: The blocks F k−1,k−1 are already defined by the formulas above (when i = k − 1) and it has the complementary non-zero blocks • Matrices F i,i−1 , i ≥ 2 have the entries caused by the demands' departure from the queue due to impatience and also non-zero blocks which have different sizes and forms in the following three cases In the case of i = 2, k − 1, the matrix F i,i−1 is the sum of both the matrix iεI W N − and the matrix H − having the non-zero blocks
for n = n (2) , N are equal to zero.In the case of i > k, the matrix F i,i−1 is the sum of both the matrix iεI W N + and the matrix H + having the non-zero blocks Formulas for H (n,n−n (1) ) + have different forms depending on the relation between the numbers n (2) , and 2n (1) − 1.

•
Matrix F 1,0 has different forms for various values of the threshold k: k = 1 and k > 1.
Its form is not obtained from the above formula for F k,k−1 via formal setting k = 1 because level 0 of the chain has another cardinality than the other levels.
In the case of k = 1, the size of this matrix is W N + × W(N + 1)M and the matrix is equal to the sum of the matrix εR 1 where the matrix R 1 has the blocks R and the matrix has the following non-zero blocks: 2,4 , n = n (1) , n (2) − 1, If k > 1, the matrix F 1,0 of size W N − × W(N + 1)M is equal to the sum of the matrix εR 2 where (1) , N, and the matrix has the following non-zero blocks: F (n,n−n (1) ) 1,0 Blocks F (n,n−n (2) ) 1,0 are defined by the formula 1) .
• Matrix F i,i+1 , i ≥ 1 has different sizes and forms in the three cases is the block diagonal matrix of size W N − with the diagonal blocks defined by In the case of i = k − 1, it is the matrix of size W N − × W N + , with the blocks defined by In the case of i ≥ k, matrix F i,i+1 is a square matrix of size W N + , with the diagonal blocks defined by Let us arrange all states of the MC ζ t defined by Formula (1) in the direct lexicographic order of its components and call the group of the states of the MC having the value i of the first component as level i of the chain.The group of the states having the value (i, n) of the first two components is called the mini-level (i, n).
The cardinality of level 0 is equal to (N + 1)W M. The cardinality of the mini-levels Let F be the generator of MC ζ t , t ≥ 0.
Theorem 1.The infinitesimal generator F has the block-tridiagonal structure: where the square blocks, F i,i , i, i ≥ 0, |i − i | ≤ 1, have the following meaning.The matrix F i,i defines (except the diagonal elements of the matrix F i,i ) transition rates from level i to level i .The diagonal elements of the matrix F i,i are negative.The modules of these elements specify the departure rate of the MC from the respective state of MC belonging to level i.In turn, the matrix F i,i consists of the blocks (F i,i ) (n,n ) that define the rates of the chain ξ t transition from the mini-level (i, n) to the mini-level (i , n ).An explicit form of the blocks (F i,i ) (n,n ) is given by Formulas (2)- (41).
Proof of Theorem 1 is given in Appendix A. Having explicit formulas for the blocks of the generator of the MC ζ t , we have an opportunity to analyze its stationary behavior.

The Case of Absolutely Patient Demands-The Ergodicity Condition of the System and Its Stationary Distribution
Due to the control by the number of active service devices, the dynamics of the MC ζ t is level-dependent.However, if the demands are patient (i.e., ε = 0), this dependence disappears for levels higher than k.Therefore, the MC ζ t can be treated as level-independent QBD with many boundary states.It follows from [40] that the ergodicity criterion for this MC is formulated as follows. Let Theorem 2. The criterion for the ergodicity of MC ζ t is the fulfillment of the inequality where the vector y satisfies the system The vector y has a finite size W N + and the solution of the system ( 43) is easily implementable on the computer.After the computation of this vector, the validity of inequality can also easily be obtained.
Let us now suppose that the parameters of the system are such that inequality (42) holds good.
Then the limits ϕ{i, n, ν (c) , ν (e) , η} = lim which are called stationary probabilities exist and the minimal non-negative solution of the quadratic matrix equation for the matrix T exists as well.

Denote
Theorem 3. The probability vectors ϕ i , i ≥ 0, are computed as follows.The vectors ϕ i , i > k, are computed as where the matrix T is the solution of Equation ( 44) having a spectral radius less than 1, while the vector (ϕ 0 , . . ., ϕ k ) is found as the solution to the finite system of equations and matrix F is given by Formula (45).

The Case of Impatient Demands-The Ergodicity Condition of the System and Its Stationary Distribution
Theorem 4. If the demands are impatient, i.e., ε = 0, then the MC ζ t is ergodic for any choice of the system parameters.
Proof.If ε = 0, then the dependence of the blocks of generator F on the level does not disappear as in the previous section.In this situation, it is feasible to prove that the MC ζ t belongs to the class of asymptotically quasi-Toeplitz Markov chains (AQTMC), which was defined in [48].
Indeed, the definition of AQTMC given in [48] requires the existence of limiting matrices where T i = −I • F i,i , i ≥ 0 and where • is the Hadamard product of matrices symbol, see, e.g., [49].
Using the above-obtained explicit expressions for the blocks of generator F , it is easy to demonstrate that if ε = 0, then for the MC ζ t , the matrices Y 0 , Y 1 , and Y 2 exist and are equal to I,O, and O, respectively.Thus, the MC ζ t belongs to the class of AQTMC.
It is proven in [48] that the condition sufficient for the ergodicity of the MC ξ t is the fulfillment of inequality wY 0 e > wY 2 e (46) where the row-vector w satisfies the system: Taking into account that Y 0 = I, Y 1 = O, and Y 2 = O, inequality (46) turns to the trivial inequality 1 > 0. Therefore, in the case of impatient demands, the MC ζ t is ergodic for any choice of the system parameters.Therefore, Theorem 4 is proven.
The statement of this theorem is intuitively clear because, due to impatience, the departure rate of demands from the system unboundedly increases when the length of a queue increases.Therefore, the queue length never becomes infinite.
In contrast to the case of the absolutely patient demands, in which generator F for large values has the block Toeplits-like structure, here the solution of this system does not have a matrix-geometric form.In general case, the solution of such systems is not possible.Some kinds of direct or soft truncation methods are applied to obtain approximate solutions.The most popular in the kind of soft truncation in the literature is described in [50].However, the truncation methods may not lead to success because a good quality of approximation can be achieved only for large values of the truncation levels and the existing computer memory is not enough to solve the corresponding finite truncated system.
As it was shown above, the Markov chain ζ t belongs to the class of AQTMC, of which the system of equilibrium equations can be solved via the corresponding algorithm proposed in [48].An effective and numerically stable algorithm for solving this system, which exploits the block tri-diagonal form of generator F , was elaborated in [51].We use this algorithm below for the numerical study of the system.

Performance Measures
With vectors ϕ i , i ≥ 0, computed, the values of several system performance metrics may be calculated.
The mean number L bu f f er of demands in the buffer is calculated as The mean number N energy of e.u. in the stock is calculated as

S
(1) The probability P 1,busy that the first service device is occupied at any time is calculated as The probability P 2,busy that the second service device is occupied at any time is calculated as The probability P busy that both service devices are occupied at any time is given by The probability P idle that both service devices are not occupied at any time is calculated as The probability P empty that the system is empty at any time is calculated as The mean number of demands in service at any time N service is calculated as The mean number of demands in the system at any time L system is calculated as The probability P energy−loss of any e.u.loss due to the battery overflow is calculated as The share P energy−leak of arriving energy, which is leaked from the battery, is calculated as It is worth noting here that the size of the vector e in this and similar formulas is not constant.It depends on the size of the vector ϕ(i, n), which is defined by the cardinality of the corresponding mini-level (i, n).
The share of the accumulated e.u., which is leaked from the battery Penergy−leak , is calculated as Penergy−leak = P energy−leak 1 − P energy−loss .
The probability of an arbitrary demand loss due to impatience P imp−loss is calculated as The departure rate λ (out) c of demands that received service is calculated as n=n (1)   ϕ(i, n) e W ⊗ S (1) 0

S
(1) Remark 2. Since the generator of the considered MC is quite cumbersome, it is necessary to manage some kind of control by using the accuracy of the theoretical derivations and computer implementation of the construction of the generator of the MC and the computation of its steady-state distribution.We present here a few of the existing possibilities used for such a control.
It is clear that the following relation must be fulfilled: The departure rate of demands that received service by the first service device λ n=n (1)   ϕ(i, n) e W ⊗ S (1) 0
The departure rate λ (2,out) c of demands that received service by the second service device is calculated as e M (1) ⊗ S (2) 0 The share Q 1 of demands that received service by the first service device is calculated as The share Q 2 of demands that received service by the second service device is calculated as .
The rate of switching off the second service device Q o f f is calculated as e M (1) ⊗ S (2) 0 ) .
The first summands here define the rate of switching off the second service device because the number of demands in the buffer at the service completion epoch is below the threshold j.The last summand is the rate of switching off the second service device because the number of demands in the buffer at the service completion epoch is not below the threshold j, but there is not enough e.u. to start a new service.
The rate Q on of switching on the second service device is calculated as Remark 3.According to the formula of the total probability, the mean number of e.u.spent for service of one demand is equal to n (1) It is clear that all which are not lost or leaked e.u. in the stationary operating system is spent for demands service.
Therefore, the following relation should be fulfilled and can be used for the accuracy of computation control: Also, it is obvious that in the stationary operating system, the relation

Numerical Examples and Optimization Problem
Let us provide some findings from numerical experiments to demonstrate the influence of various system parameters on its performance metrics.
This MAP has the rate λ (c) = 1 and the coefficient of the correlation of successive inter-arrival times is equal to 0.0864406.The coefficient of the variation is 6.33639.
The MAP of e.u. is defined by the matrices 4.0297 0.0270042 0.0733051 0.0583402 .
This MAP has the rate λ (e) = 3 and the coefficient of the correlation of successive inter-arrival times is equal to 0.153007.The coefficient of the variation is 9.85627.To vary in the experiments of the arrival rates λ c and λ e ,, we multiply the matrices k , k = 0, 1 by the corresponding constants.
PH distribution of service time by the first service device has the irreducible representation (β (1) , S (1) ) where The mean service rate is µ (1) = 2.64706.The squared coefficient of variation is equal to 1.44291.
PH distribution of service time by the second service device has the irreducible representation (β (2) , S (2) ) where The mean service rate is µ (2) = 1.25.The squared coefficient of variation is equal to 0.5.
The number of e.u.n (1) required for a demand service by the first service device is equal to 1.The number of e.u.n (2) required for a demand service by the first service device is equal to 3. The rate ε of demands reneging due to impatience is equal to 0.01.The energy leakage rate γ n does not depend on the number of the e.u. in the stock and is equal to 0.05.
Let control parameters j and k be fixed as 1 and 3.This means that the second service device is activated (conditional that the required amount of e.u. is available) when the queue length reaches a value of 3.This service device is deactivated only when it finishes service when the queue is empty.
Figures 2-8 show the reliance of certain system performance indicators on the arrival rate λ c for different values of the e.u.arrival rate λ e .It is evident from Figures 2 and 3 that the output rate of serviced demands increases along with the increase in demands' arrival rate λ c .However, for a small rate λ e = 1 of e.u.arrival, this output becomes practically constant for λ c ≥ 1.5.It is explained by the lack of e.u. for the activation of the second service device.This is clearly seen in Figure 3 that the second service device is practically not used.This implies the sharp increase in the queue length observed in Figure 4.This increase, however, is smoothed by the mass exodus of demands from the system due to their impatience.
Figure 5, on the left, illustrates the behavior of the mean number of e.u. in the stock.This number decreases with the increase in demands' arrival rate (and higher consumption of energy).For a small value of the e.u.rate, this number is negligible, which implies the very rare use of the second service device and the stabilization of the output rate of serviced demands which we already commented on above.Figure 5, on the right, presents the mean number of demands in service, which also confirms this stabilization.Also, it illustrates the higher involvement of the second service device into demands' service when the e.u.flow has a higher rate.This involvement explains the increase in the rate of successfully serviced demands with the growth of λ c that we observed in Figure 2.
The essentially higher involvement of the second service device into demands' service with the increase in λ c when the e.u.flow has a higher rate is demonstrated in Figure 6. Figure 7 shows the decrease in the loss probability of e.u.due to the stock overflow when λ c increases and more e.u. is spent for service.For small λ e (equal to 1), this loss probability becomes practically equal to zero for λ c > 1.8 because all e.u. is quickly used for service.Figure 8, on the left, illustrates the higher involvement of the second service device to provide service when there is enough e.u. and the rate λ c of demands' arrival increases.
Figure 8, on the right, shows the behavior of the rate of the second service device deactivation.It is interesting to note here that this rate is higher for the middle value (3) of λ c .For λ e = 1, the rate of the second service device deactivation is small because it is rarely activated (due to the lack of energy).For λ e = 5, the rate of the second service device deactivation is relatively small because the queue length is not small and there is enough e.u. to use the second service device for a relatively long time.
Note that the form of the blocks of the generator given in Theorem 1 and the expressions for the system performance measures given in Section 6 are quite complicated.Thus, it is very important to control the accuracy of the analytical derivation of the obtained expressions of these blocks and the computer implementation of the computations of these blocks as well as the stationary probability vectors and the various performance measures of the system.During our experiments, we used many tests for the control, including the ones mentioned in Remarks 2 and 3 above, and all of them were satisfactory with an accuracy of order 10 −11 .This MAP has the rate λ (e) = 3 and the coefficient of the correlation of successive inter-arrival times is equal to 0.374795.The coefficient of the variation is 11.6831.

Illustration of the
To make clear the importance of correlation in the arrival processes' account, let us present a comparison of the shapes of dependence of three performance measures on the demand arrival rate λ c for the presented correlated MAPs and for the stationary Poisson arrival processes (having zero correlation) with the rates λ (c) = 1 and λ (e) = 3, correspondingly.
It is obvious from Figures 9 and 10 that the existence of correlation in arrival processes has a significant effect on the system performance indicators.As a result, if the flows in a real-world system described by the model are correlated or have a significant variance in inter-arrival time, using the stationary Poisson process leads to an overly optimistic assessment of the system's performance characteristics.This clearly motivates the necessity to analyze the systems with the MAP.

Illustration of the Effect of the First Service Device's Coefficient of Variation in Service Time
An account of the coefficient of variation in service time by the first service device is achieved via the assumption that this time has a PH distribution.It should be noted that if service time distribution would be assumed exponential (which is a very particular case of PH), the analysis implemented above would not become essentially simpler.This is because the main difficulties in the analysis consist of an account of the number of e.u.available at a given moment and the moment just after service begins and also on account of the states of two service devices and the change in cardinality of the levels of M depending on these states.Anyway, it is interesting to look at the influence of the variance of service time distribution.We illustrate such an influence by considering different variances of the service time distribution at the first service device.
We assume the same system's parameters as in the first experiment and the following three variants of PH distribution of service time at the first service device having the same mean service rate.
In the first variant, service time has an exponential distribution with the rate µ 1 .It is defined by the irreducible representation (β (1) , S (1) ) where β (1) = 1 and S (1) = −µ 1 .The squared coefficient of variation in this distribution is equal to 1.
In the second variant, service time has Erlangian distribution with the mean value µ −1 1 defined by the vector β (1) = (1, 0) and sub-generator The squared coefficient of variation in this distribution is equal to 0.5.
In the third variant, service time has the hyper-exponential distribution with the mean value µ −1 1 defined by the vector β (1) = (0.02, 0.98) and sub-generator The squared coefficient of variation in this distribution is equal to 24.5102.
In Figures 11 and 12, we show the dependencies of the output rate of the serviced demands, the mean number of the demands in the buffer, and the probability of an e.u.loss on the rate µ 1 for service times at the first service device with different coefficients of variation.Figures 11 and 12 show that the variance in the service time at the first service device has some impact, but, under the fixed sent of the system parameters, this impact is not very essential, but not negligible.For the mean number of the demands in the buffer, the difference is about 14 percent.A higher coefficient of variation worsens the quality of the system operation.

Solution of Optimization Problem
We may conceive and solve numerous optimization issues after we have solved the challenge of computing the stationary distribution of the system states and the main performance metrics of the system for an arbitrarily defined set (j, k) of the thresholds, e.g., let us consider the problem of selection of control parameters j, k that minimize the cost criterion (the economic function representing the system charges during a unit of time) of the form , where b is the charge paid per unit of time for an arbitrary demand waiting in the buffer, c 2 is the payment of the system due to the use (lease) of an assistant service device during a unit of time, and c o f f is the charge paid for each switching-off of the assistant service device.Evidently, the system obtains a profit from the service of demands and is charged for not enough quick access of demands to the service, for the use of an assistant service device, and for managing its switching on (or off).To maximize the revenue of the system, the cost criterion E(j, k) has to be minimized via the proper choice of the thresholds j and k.
Let us fix the following parameters of the system.Arrival flows of demands and e.u.
are defined by the same matrices D k , k = 0, 1 as in the first experiment scaled to provide mean arrival rates λ c = 5 and λ e = 7.The rate of service by the second service device is assumed to be five times higher (µ 2 = 6.25) and the impatience rate ε is assumed to be ten times higher than in the first experiment.The remaining system settings are the same as in the first experiment.
The costs are chosen as follows: b = 1, c 2 = 1, and c o f f = 15.Figure 13, on the left, depicts the cost criterion E(j, k) dependency on the variable thresholds j and k in the area 1 ≤ j ≤ k, 1 ≤ k ≤ 30. Figure 13, on the right, shows the same dependence when the arrival flows are replaced with the stationary Poisson flows having the same arrival rate.Some comments on these Figures are as follows.In the case of the MAP flows, the value E(1, 1) corresponding to the permanent use of the assistant service device when enough e.u. is available is equal to 21.9769.The value E(30, 30) corresponding to the switching on of the assistant service device when the queue reaches the value 30 and switching it off when the number of demands in the buffer drops below 30 is equal to 10.4861.The optimal value of the cost criterion is 9.04651.It is achieved when j = 1 and k = 30.This corresponds to the switching on of the assistant service device when the queue reaches the value of 30 and switching it off only when the number of demands in the buffer drops below 1.The absolute profit gained via control defined as ∆ = min{E(1, 1), E(30, 30)} − E(1, 30) is equal to 1.43959.The relative profit gained via control defined as ∆ E(1,30) is equal to 15.91%.
In the case of the stationary Poisson flows E(1, 1) = 23.3223,E(30, 30) = 7.03713, E(1, 30) = 5.72215, and ∆ = 1.31498, the relative profit is 22.98%, i.e., the optimal set of the thresholds is the same, but the relative profit is higher in the case of the stationary Poisson flows.Lower profit in the case of the correlated flows is explained by the fact that periods of the presence of a long queue occur due to the bursty arrivals.
Many other different conditional and unconditional optimization problems with various criteria of quality and restrictions on values of some performance measures can be solved with the use of the presented above analysis, e.g., the cost function may be assumed not user-centric and include the penalty of the system for waiting for the demands as suggested above, but it may be provider-centric and includes the penalty for the loss of demands leading to the decrease in both the throughput of the system and the profit of the service provider.
It is worth noting that we intentionally restricted ourselves from the early beginning to the consideration of a concrete control policy of the hysteresis type defined above.This allowed us to reduce, via consideration of a Markov chain with the generator having quite a complicated form, the problem of the optimization of the cost criterion to the problem of minimization of a function of two integer variables.If a more complicated policy would be considered, the solution of the optimization problem, e.g., via the use of Markov decision processes, looks to be extremely difficult, if feasible.

Conclusions
In this paper, we analyzed a queuing system with two distinct service devices providing different service rates and consuming different numbers of e.u. for service of one demand.The main service device cannot stay idle if the buffer of the system is not full and the number of available e.u. is enough for service.The second, assistant service device can start service only when the queue length reaches the high threshold and the number of available e.u. is enough for service.The assistant service device can continue service after a given service completion if the queue length exceeds the low threshold and the number of available e.u. is enough for service.The parameters of the system are quite general: arrivals of demands and e.u.occur according to the MAPs, and the service time distribution is of the PH type.Demands and e.u. can depart the system after some random waiting time, having the exponential distribution.
Under any fixed pair of the low and high thresholds of the strategy of control via activation and deactivation of the assistant service device, the behavior of the system is described via a multidimensional Markov chain.In the case of absolutely patient demands, this chain belongs to the class of level-independent QBD processes with many boundary levels and its analysis is more or less standard.When the demands are impatient, this chain belongs to the class of Asymptotically Quasi-Toeplitz MC and its analysis was implemented based on the respective results from [48,51].The primary system performance measures are computed using non-trivial formulas.Their dependence on some parameters of the system is highlighted via numerical results.An example of solving the problem of optimal selection of the thresholds is presented.
The model under consideration can be utilized for analysis and optimization of realworld systems with service devices that require an online supply from some resource, e.g., it is suitable for modeling the operation of the node of the sensor network with e.u.harvesting.
The presented analysis can be extended to the cases of the correlated marked Markov arrival process of demands and energy, separate stocks for e.u.storing, finite buffer for demands and retrying of demands, disastrous disappearance of demands or (and) energy, service devices breakdowns, etc.An extension of the model to the controlled priority systems, see, e.g., [52], is planned.
block IεI W M appears due to impatience of demands and the block D(e) 1 ⊗ I M δ n,N naturally disappears because all values on n in the considered range n = 0, n (1) − 1 are not equal to N.
In the case n = n (1) , n (2) − 1, the state space of the vector process η t is reduced due to the impossibility of having an idle first service device in the presence of the queue and not less than n (1) e.u.The term Ẑ(2) appears in (10) instead of the term I W ⊗ ( S(1) 0 + S(2) 0 ) because (a) service completion by the first service device in case n ≥ n (1) implies the immediate start of the next service (and decreasing the value of i to i − 1); and (b) service completion by the second service device, in this case, is not accompanied by the immediate start of the next service due to the lack of the required e.u.(n < n (2) ).
Formula ( 11) is similar to Formula (10).Only here is the possibility of the e.u.loss in case a full stock is accounted for and a new service by the second service device cannot start, not due to the lack of e.u., but because in the considered case 1 ≤ i < j, the queue length is below the threshold value j.
In the case when j ≤ i < k, the second service device can start service, but this will again lead to the decrease in the value of i to i − 1.Thus, we obtain Formula (12).
In the case when i ≥ k and n = n (2) , N, both service devices are always busy.This implies the decrease in the cardinality of the state space of the process η t to M. As a result, we obtain Formula (13).
Formulas ( 14) and ( 15) describe the transitions of the considered MC within the level i caused by the e.u.leakage.Both formulas reflect that the other components of the MC ζ t have no transitions during the leakage moment.But, the cardinality of the state space of the process η t increases from M to M when the number of e.u.drops below the value n (1)  and this cardinality increases from M to M when the number of e.u.drops below the value n (2) when i ≥ k.The change in the cardinality from M to M is reflected by the matrix A 2,4 .The change in the cardinality from M to M is reflected by the matrix (O M×M (1) I M).This explains Formulas (14) and (15).Formulas ( 16) and ( 17) describe the transitions of the MC ζ t within the level i caused by an e.u.arrival.All other components of the MC ζ t have no transitions during the e.u.arrival moment.But, the cardinality of the state space of the process η t decreases from M to M when the number of e.u.reaches the value n (1) and this cardinality decreases from M to M when the number of e.u.reaches the value n (2) when i ≥ k.The change in the cardinality from M to M is reflected by the matrix A 4,2 .
The change in the cardinality from M to M is reflected by the matrix A 2,1 .
The form of the complementary non-zero block F (n,n−n (2) ) k−1,k−1 in the matrix F k−1,k−1 can be commented as follows.This block corresponds to the following scenario.The number of demands in the queue is equal to k − 1, the second service device is idle, and the number of available e.u. is enough to provide service by this service device.At this moment, a new demand arrives.According to the considered hysteresis strategy of control, this arrival triggers the switching on of the second service device and the service beginning.Simultaneously, n (2) e.u.disappear from the stock.Matrix B 2,4 defines the probabilities of the process η t transition in the case when n = n (2) , n (2) + n (1) − 1 and, consequently, the number of available e.u.drops below the level n (1) and the cardinality of the state space increases from M to M.
Matrix B 2,2 corresponds to the case when n = n (2) + n (1) , N, i.e., the number of available e.u.becomes not less than n (1) and cardinality M does not change.Summarizing, we obtain Formula (18).
Transitions of the MC ζ t from level i, i ≥ 2, to level i − 1 correspond to a demand departure from a queue.Such a departure can happen in the following five cases: (a) a demand abandons the queue due to impatience; (b) the first service device was idle, the number of available e.u. was n (1) − 1, new e.u.arrives, and the first service device starts work (the matrix F (n (1) −1,0) i,i−1 provides the related rates); (c) the first service device finishes service and starts the new one because the required e.u. is available (the matrix

Figure 1 .
Figure 1.Scheme of the system operation.

t 1 )
, t ≥ 0, with the state space {1, . . ., W (c) } and matrices being D , of size W(c) .The average arrival rate of demands λ c is defined by the formulaλ c = θ c D (c) 1 ewhere the invariant probability vector θ c of the process ν = 0, θ c e = 1 . The average arrival rate of e.u.λ e is defined by the formula λ e = θ e D (e) 1 e where the row vector θ e is defined as the unique solution to the system θ e (D (e) 0 + D (e)

7. 1 .
Dependence of Performance Measures on Demands' Arrival Rate λ c for Different Values of e.u.Arrival Rate λ e Let us fix the following set of system parameters.The MAP of demands is defined by the matrices D

Figure 2 .Figure 3 .Figure 4 .Figure 5 .Figure 6 .Figure 7 .
Figure 2. Dependence of the output rate λ (out) c of serviced demands on the arrival rate λ c for different values of e.u.arrival rate λ e .

Figure 8 .
Figure 8. Dependence of the share Q 1 of demands, which obtain service by the first service device, and the rate Q o f f of deactivation of the second service device on arrival rate λ c for different values of e.u.arrival rate λ e .

This
Impact of Correlation in Arrival Processes of Demands and e.u.In this experiment, we define the MAPs of demands and e.u. by the following matrices.The MAP of demands is defined by the matrices D MAP has the rate λ (c) = 1 and the coefficient of the correlation of successive inter-arrival times is equal to 0.312757.The coefficient of the variation is 10.0512.The MAP of e.u. is defined by the matrices D

Figure 9 .
Figure 9. Dependence of the output rate of serviced demands λ (out) c on the arrival rate λ c for different correlations in demand arrival processes.

Figure 10 .
Figure 10.Dependence of the mean number of e.u. in the stock N energy and the probability of an e.u.loss P energy−loss for different correlations in demand arrival processes on the mean arrival rate λ c .

Figure 11 .
Figure 11.Dependence of the output rate of serviced demands λ (out) c on the rate µ 1 for service times with different coefficients of variation.

Figure 12 .
Figure12.Dependence of the mean number of demands in the buffer L bu f f er and the probability of an e.u.loss P energy−loss on the rate µ 1 for service times with different coefficients of variation.

Figure 13 .
Figure 13.Dependence of the cost criterion E(j, k) on the thresholds for the correlated arrival processes and stationary Poisson flows of demands and e.u.having the same mean rates.