In this paper we consider a multi-server heterogeneous queueing system where servers di er in service rates and operating costs. The optimal allocation policy for this system is of threshold type. Gathering data for evaluated optimal threshold levels and system parameters is performed using a policy iteration algorithm. We study the possibility to use these data-sets to provide predictions for optimal thresholds through arti cial neural networks. The obtained results are accompanied by heuristic solution based on a uid approximation. Numerical examples illustrate the quality of provided predictions.
Many queueing systems are analyzed for their dynamic and optimal control
related to system access, resource allocation, changing service area characteristics
and so on. Sets of computerized tools and procedures provide large data-sets
which can be useful to expand potential of classical optimization methods. The
paper deals with a known model of a multi-server queueing system with
controllable allocation of customers between heterogeneous servers which are di
erentiated by their service and cost attributes. As it is known, see e.g. [1,9], the
optimal allocation policy which minimizes the long-run average cost per unit
of time for this queueing system belongs to a set of structural policies. For the
servers' enumeration (
Searching for the optimal values of q2; : : : ; qK by direct minimizing the average cost function can be expensive, especially when K is large. To calculate the optimal threshold levels we can use a policy iteration algorithm [4,6,10] which constructs a sequence of improved policies that converges to optimal one. Although this algorithm is a powerful tool for solving many optimization tasks, it has signi cant limitations on dimensionality of the model, number of states, convergence in a heavy tra c case. The contribution of the paper is two-fold. First, we provide a simple heuristic solution (HS) for a sub-optimal policy in order to avoid the search for the optimal one. Second, we investigate the possibility to use the data generated by a policy-iteration algorithm to provide a prediction for the optimal threshold levels with arti cial neural networks (NN) [3,7,8]. The trained network can be used then to calculate the optimal thresholds for those system parameters for which alternative numerical methods are di cult or impossible to use, for example, in heavy tra c case, or, in general, to reconstruct the areas of optimality without usage of time-expensive algorithms and procedures. We unsuccessfully tried to nd published works where machine learning methods would be used to solve a similar problem and therefore we consider this paper relevant.
This paper is organized as follows. In Section 2 we brie y discuss a mathematical model. Section 3 introduces some heuristic choices for threshold levels, that turn out to be nearly optimal. Section 4 presents results when the trained neural network was ran on veri cation data of the policy iteration algorithm. 2
We remind brie y the model under study. Consider an in nite-capacity M=M=K
queueing system with K heterogeneous servers and one common queue, see
Figure 1. The customers arrive to the system according to a homogeneous Poisson
process with a rate . The jth server has an exponentially distributed service
time with a rate j . The server j is called an available server if it is idle. The
service of customers is assumed to be without preemption, i.e. a customer being
served on a server can not change it. The inter-arrival and service times are
mutually independent. The system costs include an operating cost cj > 0 per
unit of time for busy server j and holding cost c0 > 0 per unit of time for any
customer waiting in the queue. Assume that the servers are enumerated in a way
(
The controller or decision maker, which has a full information about system
states, dispatchers or allocates customers according to a control policy f either
to one of available servers or to queue at a new arrival and service completion
epoch if it occurs with a nonempty queue. The system dynamics is common for
the systems with one queue and heterogeneous servers. At each arrival epoch the
customer joins the queue and the controller can allocate the customer staying at
the head of the queue to an available server j. At service completion epochs the
controller may decide to allocate the customer from the head of nonempty queue
to an available server or leave the customer in the queue. As it was mentioned
above, the optimal control policy, which minimizes the long-run average cost
per unit of time, belongs to a set of threshold policies de ned as a sequence of
threshold levels
1 = q1
q2
qK < 1:
(
We formulate the above optimization problem as a Markov decision problem associated with a multi-dimensional continuous-time Markov chain fX(t)gt 0 = fQ(t); D1(t); : : : ; DK (t)gt 0 with a set of admissible actions A = f0; 1; : : : ; Kg with elements a, where a = 0 means the allocation of the customer to the queue and a = j 6= 0 { to the jth server. The term Q(t) 2 N0 denotes the number of customers in the queue at time t, Dj (t) 2 f0; 1g { the number of customers at server j at time t. For any xed policy f , which is of threshold type with levels (q2; : : : ; qK ), we wish to guarantee that the process fX(t)gt 0 is an irreducible, positive recurrent Markov chain with a state space Ef = fx = (q(x); d1(x); : : : ; dK (x))g N0 f0; 1gK and in nitesimal generator f . The notations q(x) and dj (x) will be used further in the paper to specify the components of the vector state x 2 Ef , where q(x) denotes the queue length in state x and dj (x) { the state of the jth server in state x. The stability condition is obviously de ned through the inequality < PjK=1 j . For ergodic Markov chains with costs the long-run average cost per unit of time for the policy f coincides with the corresponding assemble average, i.e.
gf = lim sup 1 V f (x; t) = X c(y) yf ; t!1 t y2Ef where c(y) = c0q(y) + PK
j=1 cj dj (y) is an immediate cost in state y 2 Ef , V f (x; t) = E f h Z t 0
K
c0Q(t) + X cj Dj (t) dtjX(0) = xi
j=1
denotes the total average cost up to time t given initial state is x and yf =
Pf [X(t) = y] is a stationary state probability of the process under given policy
f . The policy f is said to be optimal when for gf de ned in (
min g(q2; : : : ; qK ):
q2;:::;qK
(
v(n)(e1) for x = (0; 1; 0; : : : ; 0) to (N; 1; 1; : : : ; 1) do
v(n)(x) end for + X
j2J1(x) + X j2J1(x)
1 + Pj2J1(x) j jv(n)(x hc(x)
g(n) + v(n)(x + ef(n)(x)) ej)1fq(x)=0g jv(n)(x ej e0 + ef(n)(x ej e0))1fq(x)>0g i f (n+1)(x)
argmina2A(x) v(n)(x + ea)
To evaluate optimal threshold levels and optimized value for the mean number of customers in the system the policy-iteration algorithm 1 is used. Here we use the notations
J0(x) = fj : dj(x) = 0g; J1(x) = fj : dj(x) = 1g to specify respectively a set of idle and busy servers in state x, A(x) = J0(x) [ f0g A the subset of admissible actions in state x and ej stands for a vector of dimension K + 1 with 1 in the jth position (j = 0; 1; : : : ; K) and 0 elsewhere. We convert the K + 1-dimensional state space Ef of the Markov decision process ordered in a certain way to a one-dimensional equivalent state space N0, : Ef ! N0, for state x = (q(x); d1(x); : : : ; dK (x)) 2 Ef ,
K
(x) = q(x)2K + X di(x)2i 1:
i=1
(
K 1)2K + X di(x)2i 1 = (x)
2K ; i=1
K ej) = q(x)2K + X di(x)2i 1 i=1 2j 1 = (x) 2j 1: For further details about derivation of the dynamic programming equation needed to evaluate the optimal policy the interested readers are referred to [1]. The in nite bu er queueing system is approximated by a nite bu er equivalent system in such a way that the loss probability does not exceed some speci ed small number " > 0.
Remark 1. For the bounded bu er size W the number of states is jEf j = 2K (W + 1): If the queue length q qK , all servers must be busy and the system behaves like a M=M=1 queueing system with a service rate PK j=1 j. The stationary state probabilities (q;1;:::;1), q qK , satisfy the di erence equation (q 1;1;:::;1)
K
+ X
j=1
j (q;1;:::;1) + X
j=1
K
j (q+1;1;:::;1) = 0;
which has a solution in a geometric form, (q;1;:::;1) = (qK;1;:::;1) q qK , q qK .
For details and theoretical substantiation see e.g. [2]. The threshold level qK can
be estimated using HS (
S =n( ; 1; : : : ; K ; c0; c1; : : : ; cK ) ! (q2; : : : ; qK ) : (
K < X Example 2. Some elements of the list S for the M=M=5 queueing system are (1; 20; 8; 4; 2; 1; 1; 1; 1; 1; 1; 1) ! (2; 5; 13; 30); (10; 20; 8; 4; 2; 1; 1; 1; 1; 1; 1; 1) ! (1; 4; 9; 21); (1; 20; 8; 4; 2; 1; 1; 5; 4; 3; 2; 1) ! (5; 12; 20; 20); (10; 20; 8; 4; 2; 1; 1; 5; 4; 3; 2; 1) ! (3; 8; 13; 13): 3
As it was mentioned above, the policy iteration algorithm has restrictions on dimensionality of the model, number of states, convergence in a heavy tra c case. In this section we derive a heuristic solution (HS) to estimate threshold levels qk, k = 2; : : : ; K, for the arbitrary K using a simple discrete uid approximation which is illustrated in Figure 2. Assume that qk is an optimal threshold to allocate the customer to server k in state (qk 1; 1; : : : ; 1; 0; : : : ; 0), where the | k{z1 } |K {kz+1} rst k 1 servers are busy. Now we compare the queues of the system given initial state is x0 = (qk; 1; : : : ; 1; 0; 0; : : : ; 0), where the kth server is not used for | k{z1 } | K{zk } a new customer, and y0 = (qk 1; 1; : : : ; 1; 1; 0; : : : ; 0), where the kth server is | k{z1 } | K{zk } occupied by a waiting customer. It is assumed that the stability condition holds. In Figure 2, the queue lengths are labeled by A = qk and B = qk 1. If the queue dynamics corresponded to the deterministic uid, it would decrease at the rate Pjk=11 j . When this rate is maintained until the queue is empty, it occurs respectively at points D = Pjk=1q1k j and C = Pjk=qk11 1j . The total holding times of customers in a queue with lengths qk and qk 1 are equal obviously to the areas FAOD = of triangles AOD and BOC. The mean operating cost of the rst k 1 servers until the queue is empty starting from state x0 is equal to qk where Pjk=i11 j is a probability to be served by the ith server, and starting from state y0 { is equal to (qk 1)PPjkjk==1111 cjj . According to a speci ed deterministic uid (a) (b) 1 2 3 4 5 6 7 8 9 10 4 5 6 7 8 9 10 11 12 13 1 9723 0 0 0 0 0 0 0 0 0 9723 4 4762 448 51 0 0 0 0 0 0 0 5261 2 17972431 0 0 0 0 0 0 0 0 4228 5 279 2929 192 6 0 0 0 0 0 0 3406 3 125 328 883 0 0 0 0 0 0 0 1336 6 68 179 1894 98 0 0 0 0 0 0 2239 4 0 41 136 337 0 0 0 0 0 0 514 7 0 19 106 1459 88 0 0 0 0 0 1672 ltcaua 56 00 00 306 3120 23123 1008 00 00 00 00 218510 ltcaua 89 00 00 101 8159 94364 78812 502 00 00 00 1819191 7 0 0 0 0 6 10 78 0 0 0 94 10 0 0 0 0 16 28 598 22 0 0 664 8 0 0 0 0 0 5 10 48 0 0 63 11 0 0 0 0 0 14 14 481 54 0 563 9 0 0 0 0 0 0 5 4 40 0 49 12 0 0 0 0 0 0 3 12 318 18 351 10 0 0 0 0 0 0 0 1 4 5 10 13 0 0 0 0 0 0 0 4 15 416 435 11645 2800 1055 379 251 231 93 53 44 5 predicted schema we formulate Proposition 1. The optimal thresholds qk, k = 2; : : : ; K, are de ned by Pk 1 Pk 1
j=1 j h ck j=1 cj i
qk q^k = min 1; :
c0 k Pk 1
j=1 j
(
V (x0)
V (y0) > 0: Pk 1
j=1 cj
V (x0) = c0FAOD + qk Pk 1
j=1 j
+ V (0; 1; : : : ; 1; 0; : : : ; 0);
| k{z1 } |K {kz+1}
(
| K{zk } = ck + c0FBOC + (qk k
Pk 1
1) Pkj=11 cj
j=1 j
+ V (0; 1; : : : ; 1; 0; : : : ; 0):
| k{z1 } |K {kz+1}
After substitution of (
K ; c0; c1; : : : ; cK ) and evaluate with HS the corresponding thresholds qk, k = 1; : : : ; K. Confusion matrices in Figure 3 visualize the performance of proposed heuristics respectively for the threshold levels (q2; q3; q4; q5). Each row of these matrices represents the instances in a predicted value while each column represents the instances in an actual value. The overall accuracies, i.e. the metric which describe the closeness of the measurements to a speci c value, as well as the accuracies for results with possible deviation of threshold values by 1 from the real value are summarized in Table 1.
HS q2 q3 q4 q5 Accuracy 0.8430 0.8778 0.7899 0.6282 Accuracy 1 0.9861 0.9884 0.9871 0.9769
Table 1. Accuracy for prediction with HS 4
Arti cial Neural Networks (NN) is a part of a supervised machine learning which
is most popular in di erent problems of data classi cation, pattern recognition,
regression, clustering, time series forecasting. Here we show that the NN can give
even more positive results comparing to the HS that indicates the possibility to
use it for predicting the structural control policies. The data-set S (
Example 4. The results of predictions in framework of some typical example are summarized in form of confusion matrices shown in Figure 4. The overall accuracy of classi cation and accuracies for the values with deviations are given in Table 2. We can see that the NN methodology exhibits more accurate predictions for the optimal thresholds comparing to the HS. Therefore we may conclude that classical queueing system analysis can and must be supplemented and extended by more active use of machine learning technologies. 5
We combine classic methodology of analyzing controllable queues with a heuristic solution and machine learning to study the possibility to forecast the values of optimal thresholds. When analyzing and comparing the results obtained by algorithms of the Markov decision theory and by supervised learning we may conclude that these methodologies can be seen as complementary rather than competitive.