Algorithmic Analysis of Finite-Source Multi-Server Heterogeneous Queueing Systems

: The paper deals with a ﬁnite-source queueing system serving one class of customers and consisting of heterogeneous servers with unequal service intensities and of one common queue. The main model has a non-preemptive service when the customer can not change the server during its service time. The optimal allocation problem is formulated as a Markov-decision one. We show numerically that the optimal policy which minimizes the long-run average number of customers in the system has a threshold structure. We derive the matrix expressions for performance measures of the system and compare the main model with alternative simpliﬁed queuing systems which are analysed for the arbitrary number of servers. We observe that the preemptive heterogeneous model operating under a threshold policy is a good approximation for the main model by calculating the mean number of customers in the system. Moreover, using the preemptive and non-preemptive queueing models with the faster server ﬁrst policy the lower and upper bounds are calculated for this mean value.


Introduction
The finite-source or finite-population queueing systems comparing to the ordinary markovian queues have no longer a Poisson arrival stream as in systems with an infinite source of customers, but rather have a finite source capacity N of possible customers. In such systems a customer can be inside the system, consisting in our case of one common queue with capacity N and K heterogeneous servers or outside the system in so-called arriving state. It is assumed that each customer outside arrives to the system in exponentially distributed time. After receiving the service a customer returns to the arriving area. Much attention by the study of finite-source queueing systems has been paid in terms of the machine repairman problem, see e.g., [1,2]. The customers outside the queueing system can be interpreted as unreliable machines with independent exponentially distributed life times. The queueing system represents then the repair facility where the failed machines must be recovered. Such systems are also used in various dispatching problems, they are appropriate queueing models for telephone registration systems, call centers, Ethernet systems, local-area networks, mobile communications, magnetic disk memory systems and so on.
The primary objective of this paper is to analyze optimal control and to evaluate performance metrics under fixed control policy for a non-preemptive finite-source queueing system with one class of customers and heterogeneous servers. In such a system, a customer that receives service on a slower server cannot change it if a faster server becomes available in the course of service. Unfortunately, performance analysis of this system in analytical form is limited firstly by the need to have a known allocation mechanism between the servers or control policy and secondly by the dimensionality of the corresponding Markov process, which is affected by the number of servers. To calculate the optimal allocation policy with the aim to minimize the long-run average number of customers in the system we formulate the Markov decision process (MDP) and apply the policy-iteration algorithm. This algorithm can be used not only for the optimal allocation policy calculation but also to obtain the mean number of customers in the system operating under that policy. Numerical experiments confirm our expectations that the optimal policy is of threshold type as in the models with an infinite source capacity [3]. According to this policy the fastest server must be activated whenever there is a customer in the system while the slower servers must be used only if the number of customers in the queue reaches some prespecified threshold level. The model of the non-preemptive queue operating under the optimal threshold policy (OTP) will referred to in the paper as the OTP-model.
The task of calculating other system performance characteristics for a given control policy remains unresolved. Furthermore, it should be taken into account that despite of advantages the policy-iteration algorithm has a limitation on the dimensionality of the random process for an arbitrary number of servers. In case of a threshold control policy for a particular states' ordering the corresponding Markov chain is a quasi-birth-and-death (QBD) process with a three-diagonal block infinitesimal matrix, where the blocks depend on the values of thresholds as it was shown in [4] for the infinite population system. In this case, matrix-analytic solution methods can be applied, but for a limited number of servers. This led us to discuss here in addition some simplified variants of the main model. The non-preemptive queueing system operating under a Fastest Server First (FSF) policy which prescribes for service the usage of the fastest idle server in each state and the preemptive queueing system (PS), where the service in a slower server can be interrupted if during the service time the faster server becomes idle. This system will operates according to a threshold policy, when the slower servers are activated or deactivated if the number of customers in the queue respectively exceeds or falls below a certain threshold level. Although these systems are simpler than the main model and have a low dimensions of the state-space, there are very few publications on such specific systems, especially those with analytical results.
Description of standard finite-source models with classical results, motivation examples and literature overview can be found in [5]. In [6] the author obtained the product form solution for the stationary state distribution for the finite population queueing model with a queue-dependent servers. A non-preemptive finite-population queueing system with heterogeneous classes of customers and a single server was studied in [7]. The problem of the throughput maximizing in a finite-source system with parallel queues was analyzed in [8], where some structural properties of the optimal control policy was proved. Heterogeneous multi-server finite-source queues with a FSF-policy and retrial phenomenon have been studied in [9], where numerical results were carried out by the help of the MOSEL tool. The model with machines having non-identical exponential service times was analysed in [10], where the repair policies which minimize the mean processing cost were considered. As for finite-source systems operating under a threshold policy, we know only research papers dedicated to single server systems with a unit threshold N-policy, see e.g., [11]. Here the threshold policy states that the server in a repair facility must be switched on only when the number of failed machines reaches to predefined threshold level N. In [12], the authors generalized the model to the case of two heterogeneous unreliable servers operating under policies N 1 and N 2 . Although there is a considerable amount of results on finite-source queues, the controlled systems with heterogeneous servers operating without preemption have not been investigated before. Therefore, the models presented in this paper and the corresponding results of the system analysis differ from those previously discussed.
The main contributions of paper are as follows. The structural or threshold properties of an optimal control policy are numerically established. We derive for the main model the matrix expressions used further by calculating different performance measures such as the mean number of waiting customers, the mean number of busy servers, the mean length of a busy period. The matrix-analytic solution for the stationary state distribution and performance measures are obtained for the FSF-model. Here we used the recursive definition of some blocks in the infinitesimal matrix. For the PS-models we obtain analytical product-form result for an arbitrary number of servers. Moreover, it is shown that performance characteristics of these systems in certain operation modes are the same or very close to those of the main system functioning under the optimal policy. Thus, these simplified models can be used under certain conditions to calculate upper and lower bounds for some performance characteristics and also as approximating models. We develop also the first step analysis to study the mean number of customers served in the system or by the kth server in a busy period and the probability of the maximum queue length observed during this period.
The rest of the paper is organized as follows. In Section 2, we describe the Markovdecision process of the main model and show that the system has a threshold-based optimal allocation policy. In this section we develop also the computational analysis for the mean values of performance measures including those characterizing the behaviour of the system in a busy period. The FSF-Model is presented and analysed in Section 3. Section 4 is devoted to the PS-model. Comparison analysis of the proposed models and illustrative examples are summarized in Section 5.
The following notations will be used throughout this paper: e(n), e j (n) and I n stand respectively for the unit vector of dimension n, for the basis vector of dimension n in R n with 0 ≤ j ≤ n − 1 or 1 ≤ j ≤ n depending on the context, and for the identity matrix of dimension n. If it is not necessary to specify a vector dimension, we will omit the corresponding argument. For example, e denotes a column unit vector of an appropriate dimension. The notation is used for the transpose. The notation ⊗ stands for the Kronecker product of two matrices, 1 {A} -for the indicator function, where 1 {A} = 1 if the condition A holds, and 0 otherwise. The notation |A| is used for the magnitude of a finite set A.

OTP-Model
Here we discuss the main model of the non-preemptive finite-source controlled queueing system of the type M/M/K/N//N illustrated in Figure 1. The system has K heterogeneous servers with different rates µ 1 ≥ µ 2 ≥ · · · ≥ µ K > 0 and N customers in a source. It operates under the optimal allocation policy which minimizes the mean number of customers in the system. It will be shown that this policy is defined through a sequence of threshold levels 1 = q 1 ≤ q 2 ≤ · · · ≤ q K < ∞ for the queue lengths which prescribe the activation of slower servers. The analysed system can be treated as a model for the machine-repairman problem, where N unreliable machines in a working area with exponential distributed life times and equal rates λ > 0 must be repaired by K heterogeneous repair stations. The machines fail independently of each other. The stream of failed machines can be treated as an arrival stream of customers to the queueing system. Hereafter, we will refer to the customer as a failed machine which enters the repair system and gets there a repair service. After the repair the machine becomes as good as a new one and it returns to the working area. The aim is to dynamically allocate the customers to the servers in order to minimize the long-run average number of customers in the system and to calculate the corresponding mean performance measures.

MDP Formulation
We formulate the optimal allocation problem in this machine-repairman system as a Markov Decision Process (MDP) in the following way. The behaviour of the system is described by a multi-dimensional continuous-time Markov-chain where Q(t) stands for the number of customers waiting in the queue at time t and D j (t) specifies the state of the jth server at time t, where  State space: The set E X consists of K + 1 dimensional row vectors, where q(x) denotes the number of customers in the queue and d j (x) -the status of the jth server in state x. The total number of states in the set E X is equal to |E X | = ∑ K j=0 ( K j )(N − j + 1). Decision epochs: The arrival and service completion epochs in the system with waiting customers. Action space: A = {0, 1, . . . , K}. To identify the group of idle and busy servers, the following sets are defined, With this notations the set of admissible control actions A(x) ⊆ A in state x ∈ E X can be defined as A(x) = J 0 (x) ∪ {0}. The action a ∈ J 0 (x) means that in state x a customer must be allocated to an idle server, while a = 0 means that the customer must be routed to the queue. At an arrival epoch, which occurs only if the number of customers in the system is less than N, the arrived customer joins the queue and simultaneously another one from the head of the queue must be routed to some idle server or returned back to the queue. At a service completion epoch the same happens, i.e. the customer from the head of the queue is routed either to one of idle servers or to the queue again. By service completion in a system without waiting customers no actions have to be performed. Immediate cost: The function l(x) specifies the number of customers in a state x ∈ E X , i.e., which is in fact independent of a control action a.
Transition rates: The policy-dependent infinitesimal matrix Λ f = [λ xy (a)] x,y∈E X of the Markov-chain (1) includes the rates to go from state x to state y given the control action is a defined as with λ x (a) = −λ xx (a) = − ∑ y =x λ xy (a), where S a and S −1 j stand for the shift operators applied to the vector state x in the following way, Due to the finiteness of the state space E X and boundedness of the immediate cost function l(x) ≤ N, a stationary average-cost optimal policy f : E → A exists with a finite constant mean value for the number of customers in the system which is independent of the initial state x. In this case the policy-iteration algorithm introduced in Algorithm 1 converges.
This algorithm consists of two main parts: Policy evaluation and Policy improvement. In the first part, a system of linear equations with immediate costs l( is solved for the unknown real-valued dynamic-programming value function v f : E X → R and mean value g f given a control policy is f . The second part of the algorithm is responsible for improving the previous policy, which for a given system consists in determining, for each system state, a control action a that minimizes the value function v(S a x). The improved control action in state x is defined then as Thus, the algorithm constructs a sequence of improved control policies until it finds one that minimizes the gain g f . In Algorithm 1 we perform a conversion of the K + 1-dimensional state space E X of the Markov chain (1) to one-dimensional equivalent state space using the function In one-dimensional state space the transitions due to arrivals and service completions can be defined then as For more details about derivation of the optimality equation for heterogeneous queueing systems the interested reader is referred to relevant publications, e.g., [3].
6: end for 7: Policy improvement else n ← n + 1, go to step 4 10: end if 11: Threshold evaluation Numerical analysis confirms our expectation that the optimal control policy in heterogeneous systems for a finite number of customers also belongs to a class of threshold policies, as in infinite population case. Theoretical justification of this statement is still difficult. For this purpose it is necessary to prove that the dynamic-programming operator B defined for our queueing model as where T 0 and T j are the events operators in case of a new arrival and a service completion at server j ∈ J 1 (x), preserves the monotonicity properties of the increments of the value function v: In proving the inequality (7) we encounter difficulty. This is due primarily to the form of the operator B in (5). There is a term describing arriving customers whose coefficient (N − l(x))λ depends on the system state x. Bringing the terms in inequality (7) to a common denominator by introducing fictitious transitions, we get terms which cannot be proved to be negative. We hope that we will be able to overcome these difficulties in our next paper, but to date we're basing our statement about a threshold structure of the optimal control policy f exclusively on the performed numerical experiments. The following example makes the case vividly. Example 1. Consider the system with K = 5, N = 60 and λ = 0.3. The service rates take the following values: (µ 1 , µ 2 , µ 3 , µ 4 , µ 5 ) = (20, 8, 4, 2, 1). The Table 1 of optimal control actions f (x) for selected system states x is of the form:  Threshold levels q k , 2 ≤ k ≤ K, are evaluated by comparing the optimal actions f (x) = 0 and In this example the optimal policy f is defined here through a sequence of threshold levels (q 2 , q 3 , q 4 , q 5 ) = (1,2,4,9) and g f = 4.91549.

Evaluation of System Performance Measures
We are concerned in calculation of the system performance measures for a given policy f . The state probabilities and performance characteristics defined here will refer to some particular fixed control policy f , so we will use in notations the corresponding upper index.
The states x of the set E X with q(x) = 0 are ordered according to the number of busy servers |J 1 (x)| while the states for q(x) > 0 are ordered with respect to the queue length, so that the infinitesimal matrix Λ f has a block three-diagonal structure for the fixed policy f . First we define the performance characteristics: • The probability that the kth server 1 ≤ k ≤ K is busy, The mean number of busy servers, The mean number of customers in the queue, The mean number of customers in the system,N f =C f +Q f . The following vectors of dimension |E X | − 1 comprise the policy-dependent values a f (x) and policy-independent values b(x), where the first elements of the vectors are respectively a f (e 1 (K + 1)) and b(e 1 (K + 1)). Denote byM where the vector a f is a solution of the system The matrixΛ f is obtained from Λ f by removing the first column and the first row, and Proof. We multiply the both sides of the equality (9) by the row-vector of the stationary x for the corresponding function b(x) is obviously equal to the performance measureM f 1 . The following sequence of relations validates the statement.
The following measures characterize the behaviour of the system in a busy period which we define as a duration starting when the arrived customer enters the empty system in state x 0 and finishes when the system visits x 0 again after a service completion.

•
The mean length of a busy period,L f = 1 • The mean number of customers served in a busy period by the kth server,N f L,k ; • The total mean number of customers served in a busy period, In the following proposition we describe a general way to calculate these characteristics for the fixed control policy f . Denote byM where the vector a f is a solution of the system The matrixΛ f is obtained from Λ f by removing the first column and the first row, and the Laplace-Stiltjes transform (LST) of the probability density function (PDF) ϕ f x (t) for the first passage time to state x 0 given that the initial state is x ∈ E X , the control policy is f and byL dt the corresponding first moment. According to the first step analysis we get for the LST the system , we can obtain from (13) the system for the conditional moments The system (14) for the transition rates (2) is of the form By expressing relations (15) in matrix form and taking into account thatL f := L f (e 1 (K + 1)) we obtain the expressions (11) for a f (x) =L f (x).
Denote now byψ x,k (i)z i , |z| ≤ 1, the probability generating function (PGF) of the PDF ψ f x,k (i) of the number of service completion at server k up to the end of busy period given that the initial state is x ∈ E X \ {x 0 }. With respect to the law of the total probability we get the following relations for the function ψ f x,k (i), The first term on the right hand side of (16) represents the transition to state u accompanied with an event we count, that is a service completion at server k. The second term stands for other possible transitions. The system (16) can be rewritten in terms of the PGF in the following form,ψ The expressions (17) can be modified using the propertyN in such a way that we get a system for the corresponding first moments, For the model under study the system (18) is of the form The last system can be also expressed in form (11) K + 1)). For the mean total number of customers servedN L the term d k (x)µ k on the right hand side of (19) must be replaced by ∑ K k=1 d k (x)µ k .
Finally, one more performance measure in this section is of our interest, namely, the distribution of the maximal queue length in a busy period for the given control policy f . Denote by Q f max the maximum number of customers waiting in the queue during a busy period. For each fixed value n ≥ 0 the event {Q f max ≤ n} is equivalent to the event that the process {X(t)} t≥0 starting in state e 1 (K + 1), where the first server is busy, hits the empty state x 0 before visiting the subset of states The probabilityQ f max,n = P[Q f max ≤ n] will be calculated by means of absorption probabilities for states in a set of absorbing states E max,n ∪ {x 0 } given that the initial state is x ∈ E X,n = E X \ E max,n ∪ {x 0 }. Denote by where n is fixed. Denote further byM where the vector a f is a solution of the system The matrixΛ f (n) is obtained fromΛ f by removing all columns and rows starting with the n + 1, andM Proof. Denote byQ f max,n (x) the probability of absorption into empty state x 0 starting in x ∈ E X,n , whereQ f max,n =Q f max,n (e 1 (K + 1)), where e 1 (K + 1) as before is the state after an arrival to an empty state x 0 . The following system can be obtained by conditioning on the next visited state Using again the first principles, For the queueing system operation under the control policy f the system (23) is of the form, Then after a routine of (block) identification the system (24) can be expressed in form (21), where a f (x) =Q f max,n (x), x ∈ E X,n .
As we can see, calculating the performance characteristics requires solving very similar systems of equations. Thus, the same algorithm can be used for this purpose by substituting appropriate values into vectors a f and b, This versatility of the proposed approach greatly simplifies the application of algorithmic types of analysis of complex controlled queueing systems. In principle, we assume that for a fixed control threshold policy, the structure of the infinitesimal matrix can be even fully defined for an arbitrary number of servers, as will be proposed in the next section for the special case of the control policy where all thresholds are equal to 1. Thus we believe that matrix expressions can be derived explicitly from the presented matrix systems for performance characteristics. We leave this problem for our research in the near future.

FSF-Model
Here we discuss the FSF-Model which is a special case of the OTP-model, where q 1 = q 2 = · · · = q K = 1. The Markov-chain {X(t)} t≥0 operating under the FSF-policy has a state space The states in E X are divided in to levels y in the following way, Denote by s i,j = ( K−j+i i ) for K ≥ j, then |y| = s y,y for 1 ≤ y ≤ K and |y| = 1 for K + 1 ≤ y ≤ N. Within each level y, 1 ≤ y ≤ K, the states are ordered in the lexicographic order, where the rank of x in the level y with |J 1 (x)| = y and |J 0 (x)| = K − y can be evaluated by where n i (x) = |{j : d j (x) = 1, d i (x) = 0, j > i}| is the number of slower busy servers as the ith idle one. Note that this ordering of states differs from that defined in (4) and used in the policy iteration algorithm. In the lexicographic ordering within each level of states it is possible to obtain explicit matrix expressions for state probabilities in case of an arbitrary number of servers K. Denote further by L y , 1 ≤ y ≤ K, matrices whose rows consist of ordered elements of level y.

Proposition 4. The the system under FSF-policy is described by a QBD process with a block-three diagonal infinitesimal matrix of the form
The square blocks A 1,y of dimension s y,y for 0 ≤ y ≤ K − 1 and 1 for K ≤ y ≤ N consist of the rates to stay in the yth level, are defined as A 1,y = I s y,y (e (s y,y ) ⊗ [L y B 0,1 + (N − y)λe(s y,y )]), 0 ≤ y ≤ K − 1, The blocks A 0,y of dimension s y−1,y−1 × s y,y for 1 ≤ y ≤ K and of dimension 1 for K + 1 ≤ y ≤ N consist of the rates to move upwards from the level y − 1 to y due to arrivals and are defined as The blocks A 2,y of dimension s y,y × s y−1,y−1 for 1 ≤ y ≤ K and of dimension 1 for K + 1 ≤ y ≤ N consist of the rates to move downwards from the level y + 1 to y due to service completions and are defined as recursive matrices Proof. Analysing the transitions of the Markov-chain {X(t)} t≥0 we get a system of balance equations in form where π x = lim t→∞ P[X(t) = x], x ∈ E. Expressing Equation (30) for the sub-vectors π y , 1 ≤ y ≤ K − 1, and the scalars π 0 and π y , K ≤ y ≤ N, by means of defined blocks and taking into account the states' ordering (25) we get the system π 0 A 1,0 = π 1 A 0,1 , π y A 1,y = π y−1 A 0,y + π y+1 A 2,y , 1 ≤ y ≤ K − 2, Denote by π the macro-vector of the stationary state probabilities, i.e. π = (π 0 , π 1 , . . . , π K−1 , π K , . . . , π M ).
Compiling relations (31) to the the system πΛ = 0 we get then the infinitesimal matrix Λ is the form (26) with blocks defined by (27)-(29).

Proposition 5. The elements of the stationary probability macro-vector π satisfy the relations
where the matrices M y satisfies the recursive relations (36) Proof. The probability π 0 and sub-vectors π y , 1 ≤ y ≤ K − 2, can be expressed from the balance equations (31) using a block forward elimination-backward substitution as We similarly obtain an expression for π K−1 , The relations (32) and (33) are obtained then through a successive substitution. The relation (34) is obtained by solving (31) for K ≤ y ≤ N recursively using π y = (N − y + 1)λ m K π y−1 starting from the last equation. The relation (35) for the probability π K is determined using the normalizing condition πe(N) = 1.
To calculate performance characteristics the expressions from the previous section applied to a control policy f (x) = argmax j∈J 0 (x) {µ j } can be used. As an alternative to the policy-iteration algorithm we can use the proposed matrix-analytic solution (32)-(35) to obtain the matrix expressions for the performance characteristics in an explicit form.

Corollary 1.
The probability that the kth server is busy and the mean number of busy servers, The mean number of customers in the queue, The mean number of customers in the system, Mean length of a busy period, Mean number of customers served in a busy period, The mean number of customers served by the kth server in a busy period and the distribution of the maximal queue length can be evaluated using the matrix systems (11) and (20) taking into account the structure (26) of the infinitesimal matrix Λ. Proposition 6. The mean numberN L,k of customers served in a busy period by the kth server satisfies the relationN L,k = e 1 (K) where

(38)
The column-vector b y = L y e k (K + 1)µ k for 1 ≤ y ≤ K and the scalar b y = µ k for K + 1 ≤ y ≤ N.

Proof. The system (11) can be rewritten for appropriate blocks in form
The elements of b y are equal to µ k if for some state x of the level y d k (x) = 1. This implies the relations for b y . Using a forward elimination -backward substitution we get the recursive relations where S y and T y are defined by (38). This statement follows through recurrence substitution taking into account thatN L,k = e 1 (K)a 1 , since the level 1 consists of K states.
The following statement for the matrix equation (20) can be proved in a similar way taking into account the structure (26) of the infinitesimal matrix Λ.

Proposition 7.
The probability of the maximum queue length in a busy period satisfies the relation where (40)

PS-Model
In this section we discuss a queueing system with a preemption operating under a general threshold policy f defined as a sequence of threshold levels (q 2 , . . . , q K ). The first server in this system is permanently available for service while the jth slower server must be used as soon as the number of customers in the system increases up to the value q j−1 + j − 2. This server must be removed from the system when the number of customers becomes again less as q j−1 + j − 2. Denote by {Y(t)} t≥0 the continuous-time Markov-chain with a state space E Y = {y : y ∈ N 0 }. All the rates are the same as in the model without preemption. The infinitesimal matrix Λ f Y = λ xy (q 2 , . . . , q K ) is then of the form: where m j = j ∑ i=1 µ i and q 1 = 1. The state transition diagram of the process {Y(t)} t≥0 is illustrated in Figure 2.

Corollary 2.
The kth server is busy with a probability, The mean number of busy servers,C The mean number of customers in the queue, The mean number of customers in the systemN f =C f +Q f .
The mean length of busy period, The mean number of customers served in a busy period, Further we use a similar methodology that has been employed in previous section to derive expressions forN f L,k andQ f max,n with the knowledge that all levels y consist now of only one state, and hence in the sequel we omit some details. where Proof. The proof follows directly from (11) by forward elimination -backward substitution taking into account the structure (41) of the infinitesimal matrix Λ f . where The last result can be rewritten in explicit form as well.
Proposition 11. The distribution function of the maximum queue length observed in a busy period is given byQ where the function F(n) has the following product form, where the function F(y) has a product form (48). Summing (51) for y = 1, . . . , n yields whereQ f max,n (n + 1) = 0 andQ f max,n (0) = 1. ExpressingQ f max,n (1) we obtain the explicit formula (47).

Comparison Analysis
In this section we discuss the results after having computed the performance metrics for the following finite-source heterogeneous queueing models: Non-preemptive queueing system operating under the optimal threshold policy (OTP-Model), non-preemptive queueing system with a fastest server first policy (FSF-Model), preemptive queueing system (PS1-Model), where the kth server is used when at least k customers present in the system, and preemptive queueing system (PS2-Model) operating according to a given threshold policy. This policy we calculate using a similar heuristic formula obtained in [13], which can be rewritten in form where q 1 = 1 and (N −N PS1 )λ is an average arrival rate in the PS1-Model which is derived in explicit form.
In our experiments we fix the number of servers K = 5, the source capacity N = 60 and service intensities (µ 1 , µ 2 , µ 3 , µ 4 , µ 5 ) = (20, 8,4,2,1). The rate λ will be varied in the interval [0.01, 0.7]. The choice of this interval is not random. At higher values of λ, the analysed functions become indistinguishable, since the corresponding queueing systems will have similar stochastic behaviour in a so-called heavy-traffic mode.
In Figure 3, we display the functionsN f (figure labeled by "a") andQ f (figure labeled by "b") for different models as λ varies. We observe that the functionsN FSF andQ PS1 models are the natural upper and low bounds forN OTP . It is clear that the FSF-model is a particular case of the OTP and the queue with a preemption is always superior in performance comparing to the non-preemptive case.
Differences between the functionsN OTP andN PS2 are almost not visible. This effect is also observed for other values of system parameters. It allows the preemptive system under a threshold policy to be used as an approximation of the original OTP-model. In contrary, the PS2-model exhibits the higher values of queue lengths while the PS1-modelthe shortest. The OTP-model also has in average more waiting customers as in FSF-model which is not surprising, since the optimal policy minimizes in our case the mean number of customers in the system but not in the queue. It should also be noted that the higher the degree of heterogeneity of the servers, the greater the differences in performance functions for different models become. 12  Figure 4 illustrates the influence of λ and model types on the functionsC f (figure labelled by "a") andL f (figure labelled by "b"). The functions of the mean number of busy servers for the OTP-and PS1-models are very close to each other. Thus, by subtracting the mean number of busy servers in PS1-model from the mean number of customers in PS2-model, an approximation can be obtained for the mean queue length of the OTP-model. The functionsC FSF andC PS2 represent the upper and low bound forC OTP . The longest busy period appears in FSF-model. In this case the slower servers can be occupied with higher probability and then these servers remain busy for a very long time. As expected, the shortest busy period exhibits the preemptive PS1-model.    Figure 5 shows the effect of the service speed of kth server, 1 ≤ k ≤ K, to the mean number of customersN L,k served in a busy period (figures are labeled respectively by "a"-"f"). We observe that the slow servers begin to contribute to the number of customers served as the intensity of λ increases. The functionsN f L,k are proportional to the rate λ, they are simply shifted to the right as λ is getting higher without changing their form. It can be observed also that the FSF-policy maximizes the number of customers served in a busy period at any server. This observation coincides with a statement in [14] that the fastest available server stochastically maximizes the number of service completions.     We now focus on the results obtained for the maximum queue length observed during a busy period. To study the influence of system parameters and model type we summarized the results in Table 2 for OTP-and FSF-models and in Table 3-for PS1-and PS2-models. In tables we vary the rate λ keeping as before other parameters constant. The results compiled and presented in tables correlate with the graphs for the mean length of the busy period. The longer the busy period is, the more likely there will be fewer waiting customers in the queue. In the FSF-model it is more likely that there is an empty waiting line. As λ increases, the queues grow and hence we observe that for all models that the 99th percentile increases.       We have also conducted various experiments where we analyzed the effect of the number of servers, the source capacity, the level of heterogeneity and so on to performance metrics of non-preemptive heterogeneous systems and possible approximations through their preemptive equivalents. Due to the space limitations of the paper, we omit these results. As a generalisation, we can state that the main observations we made in the presented examples remain valid also for other values of system parameters.

Conclusions
Finite-source multi-server heterogeneous systems without priority service interruption are described using a multivariate Markov-chains. For such a systems we have found the optimal threshold policy and calculated the corresponding performance measure. Both analytical and numerical studies of such a system face constraints on the dimensionality of the problem, i.e., on the number of servers. In this paper we have also tried to understand, whether there are simplified variations of the main model which are appropriate for boundary values calculation or even for approximation of the main model but without constraint on the number of servers. We have analyzed non-preemptive and preemptive queues and provided comparison analysis of the performance characteristics.

Acknowledgments:
The authors are very grateful to the reviewers of the paper for their comments which improved its quality.

Conflicts of Interest:
The authors declare no conflict of interest.