Discrete and Continuous Models and Applied Computational ScienceDiscrete and Continuous Models and Applied Computational Science2658-46702658-7149Peoples' Friendship University of Russia named after Patrice Lumumba (RUDN University)1837010.22363/2312-9735-2018-26-2-167-175Research ArticleAnalysis of Queueing Systems with an Infinite Number of Servers and a Small ParameterVasilyevS ACandidate of Physical and Mathematical Sciences, assistant professor of Department of Applied Probability and Informatics of Peoples’ Friendship University of Russia (RUDN University)vasilyev_sa@rudn.universityTsarevaG OPhD student of Department of Applied Probability and Informatics of Peoples’ Friendship University of Russia (RUDN University)gotsareva@gmail.comPeoples’ Friendship University of Russia (RUDN University)1512201826216717521042018Copyright © 2018, Vasilyev S.A., Tsareva G.O.2018In this paper we consider the dynamics of large-scale queueing systems with an infinite number of servers. We assume that there is a Poisson input flow of requests with intensity . We suppose that each incoming request selects two any servers randomly and at the next step of an algorithm is sending this request to the server with the shorter queue instantly. A share () of the servers that have the queues lengths with not less than can be described using a system of ordinary differential equations of infinite order. We investigate this system of ordinary differential equations of infinite order with a small real parameter. A small real parameter allows us to describe the processes of rapid changes in large-scale queueing systems. We use the simulation methods for this large-scale queueing systems analysis. The numerical simulation show that the solution of the singularly perturbed systems of differential equations have an area of rapid change of the solutions, which is usually located in the initial point of the problem. This area of rapid function change is called the area of the mathematical boundary layer. The thickness of the boundary layer depends on the value of a small parameter, and when the small parameter decreases, the thickness of the boundary layer decreases. The paper presents the numerical examples of the existence of steady state conditions for evolutions () and quasi-periodic conditions with boundary layers for evolutions ().countable Markov chainslarge-scale queueing systemssingular perturbed systems of differential equationsdifferential equations of infinite ordersmall parameterсчётные марковские цепикрупномасштабные системы массового обслуживаниясингулярные возмущённые системы дифференциальных уравненийдифференциальные уравнения бесконечного порядкамалый параметрIntroduction The current study of service networks with complex sending discipline in [1-3] transport networks [4-6] and the asymptotic behavior of Jackson networks [7] handled with the problematic of verifying the global convergence of the solutions of certain infinite systems of ordinary differential equations to a time-independent solution. In [8] the countable systems of differential equations with bounded Jacobi operators are studied and the necessary settings of global stability and global asymptotic stability was found. In [9] an infinite-server queuing system with a doubly stochastic Poisson input flow is considered. Assuming that the service time does not have expectation, limit theorems for the number of occupied servers is proven. As a consequence, limit theorems for systems is obtained in which the input flow intensity is a regenerative process. In [10] an infinite-server queueing system where customers come by groups of random size at random i.d. intervals of time is considered. Assuming that the number of requests in a group and intervals between their arrivals can be dependent, service times have a regularly varying distribution with infinite mean. Limit theorems for the number of customers in the system and prove limit theorems under appropriate normalizations are considered. In papers [11-13] the authors built various models of queueing systems and considered their dynamics. In paper [14] we examined the singular perturbed structures of ordinary differential equations of infinite order of Tikhonov-type ϵбє‹ = F x t,gx ,y t,gy ,t,бєЏ = f x t,gx ,y t,gy ,t with the initial conditions x(t0) = gx, y(t0) = gy, where x,gx О X, X М l1 and y, gy О Y , Y О Rn, t О t0,t1 (t0 t1), t0,t1 О T, T О R, gx and gy are given vectors, ϵ 0 is a small real parameter. In this paper we apply Dobrushin approaches from [1-3]. We consider the dynamics of large-scale queueing systems that consist of infinite number of servers with a Poisson input flow of requests of intensity Nl. We can use an algorithm that selects two any servers for each incoming request and sent it to one of the servers with the shorter queue instantly. We suppose that service time has mean 1в€•m with exponential distribution. In this case a share uk(t) of the servers that have the queues lengths with not less than k can be described using a system of ordinary differential equations of infinite order. We investigate this system of ordinary differential equations of infinite order with a small real parameter. A small real parameter allows us to describe the processes of rapid changes in large-scale queueing systems. Tikhonov type Cauchy problem for this system with small parameter ϵ and initial conditions is investigated. We investigate the truncation system of this ordinary differential equations of infinite order with a small real parameter order N. Tikhonov type Cauchy problem for this truncation system with small parameter ϵ and initial conditions is used for the simulation of behavior solutions and for analysis of large-scale queueing systems with taking into account parameters l, m, ϵ. 1. Queueing Systems with Infinite Number of Servers The basic model considered there is a queueing system SN, with N identical infinite-buffer FCFS (First-Come, First-Served) single-servers, with a Poisson arrival flow of rate Nl and with i.i.d. exponential service times of mean 1в€•m, where 0 l m. Upon its arrival each task chooses m servers at random (i.e., independently of the pre-history of the queueing system (QS) and with probability 1в€•(Nm)) and then selects, among the chosen ones, the server with the lowest queue-size, i.e., the lowest number of tasks in the buffer (including the task in service). If there happen to be more than one server with lowest queue-size, the task selects one of them randomly. One is interested in the “typical” behavior of a server in SN, as N ® Ґ. Formally, it means that "t і 0 and k = 0,1,†, we consider the fraction qk(t) = Mk(t)в€•N where Mk(t) is the (random) number of servers with the queue-size k at time t. Clearly, 0 Ј qk(t) Ј 1, е kqk(t) = 1; and Q(t) = (qk(t)), t і 0, forms a Markov process (MP). Technically, it is more convenient to pass to the tail probabilities rk(t) = е jіkQk(t); the state space of the corresponding MP UN(t) = (fk(t)),t і 0, is the set UN of non-increasing non-negative sequences u = (uk,k = 0,1,†) with u0 = 1, е k1uk Ґ and with the uk’s multiple of 1в€•N, which implies that uk = 0 for all k large enough. It is convenient to prolong the sequences u О UN to the negative k’s by the value 1. The generator of UN(t) is an operator A acting on functions f : UN ® C1 and given by ANf(u) = Nе k0 uk - uk+1 f u -ek N - f(u)+ + lNе k0 (uk-1)2 - (u k)2 f u + ek N - f(u).(1) Here, ek stands for the sequence with the k-th entry 1 and all others 0, the addition of the sequences is componentwise. Process UN(t) is positive-recurrent and thus possesses a unique invariant distribution, pN; given any initial distribution v, the distribution of UN(t) approaches pN as t ® Ґ. The main result of [1] is that, as N ® Ґ, the expected value EpNrk(t) converges to the value ak, where ak = l m(mk-1)в€•(m-1) ,k і 0. (2) Pictorially speaking, it means that, as N ® Ґ, an “average” server in the QS will have k or more tasks in the buffer with probability ak. It is interesting to compare SN with another queueing system L, where the arriving task chooses the server completely randomly (i.e., independently of the pre-history and with probability 1в€•N). Clearly, L is equivalent to an isolated Mв€•Mв€•Ґ queue with the arrival and service rates l and m, respectively, which justifies omitting subscript N in this notation. More precisely, the average server in L will have k or more tasks in the buffer with the geometrical probability ak0 = l mk,k і 1, (3) (independently of N), which is much larger than ak. In fact, as was shown in [1], the whole process UN(t) is asymptotically deterministic as N ® Ґ. More precisely, let U denote the set of the non-increasing non-negative sequences u = (uk,k О Z) with uk = 1 for k Ј 0 and е kЈ0uk Ґ.Then, if the distribution v of initial state UN(0) approaches a Dirac delta-measure concentrated at a point g = gk О U, the distribution of UN(t) is concentrated in the limit at the “trajectory” u(t) = uk(t), t і 0, giving the solution to the following system of differential equations uМ‡k(t) = m uk+1(t) - uk(t) + l (uk-1(t))2 - (uk(t))2 , u0(t) = 0,uk(0) = gk і 0,k = 1,2,†,t і 0, (4) where g = gk k=1Ґ is a numerical sequence (1 = g1, gk і gk+1, k = 1,2,†) [1]. Point a = (ak) (see (??)) is a (unique) fixed point for system (??) in U. These results illustrate the essence of the mean-field approximation for QS SN. Equations (??) describe a “self-compatible” evolution of vector u(t), or, equivalently, of the probability distribution q(t) = qk(t) defined by qk(t) = uk(t) - uk+1(t), t і 0, k = 0,1,†. As before, u(t) is simply the sequence of the tail probabilities for q(t). We can compare system (??) with the linear system бєЏk(t) = m yk+1(t) - yk(t) + l yk-1(t) - yk(t), (5) (where k і 1) describing the evolution of the probability distribution q0(t) = (qk0(t),qk0(t) = yk(t) - yk+1(t)) in a standard Mв€•Mв€•1в€•Ґ queue with the arrival and service rates l and m, respectively. The m-terms in (??) and (??) are the same; they correspond with the departure of the tasks and ’push’ the probability mass in q(t) and q(0)(t) towards k = 0. On the other hand, the l-terms (different in both SQ) correspond with the arrival of the tasks; these terms shift the probability mass to larger k’s. The l-term in (??) is smaller than the one in (??) when uk(t) is small; pictorially speaking, system (??) provides (for the same values of l and m) more “protection”, for large k, against the shift to the right, which may lead to an “explosion”, when the relation е k1uk(t) Ґ or е k1yk(t) Ґ may fail as t ® Ґ. Because of this, the entries ak of sequence a (see (??)) giving the fixed point of (??) decrease “super-exponentially”, in contrast with the exponential decay of the tail probabilities in the fixed point a0 = (ak0) of (??). 2. Queueing Systems with Infinite Number of Servers and a Small Parameter Let’s consider a system that consists of N servers with a Poisson input flow of requests of intensity Nl. Each request arriving to the system randomly selects two servers and is instantly sent to the one with the shorter queue. The service time is distributed exponentially with mean tМ„ = 1в€•m. Let uk(t) be a share servers that have the queues lengths with not less than k. It is possible to investigate the asymptotic distribution of the queue lengths as N ® Ґ and l 1 [1]. The considered system of the servers is described by ergodic Markov chain. There is a stationary probability distribution for the states of the system and if N ® Ґ the evolution of the values uk(t) becomes deterministic and the Markov chain asymptotically converges to a dynamic system the evolution of which is described by system of ordinary differential equations of infinite order uМ‡k(t) = m uk+1(t) - uk(t) + l (uk-1(t))2 - (u k(t))2 . (6) For this system of ordinary differential equations of infinite order we can formulate Cauchy problem in the form uМ‡k(t) = m uk+1(t) - uk(t) + l (uk-1(t))2 - (uk(t))2 , u0(t) = 0,uk(0) = gk і 0,k = 1,2,†,t і 0, (7) where g = gk k=1Ґ is a numerical sequence (1 = g1, gk і gk+1, k = 1,2,†) [1]. We can investigate Cauchy problem for system of ordinary differential equations of infinite order with small parameter such form uМ‡k(t) = m uk+1(t) - uk(t) + l (uk-1(t))2 - (uk(t))2 ,k = 0,1,†,n - 1, uМ‡n(t) = m Un+1(t) - un(t) + l (un-1(t))2 - (un(t))2 , ϵUМ‡k(t) = m Uk+1(t) - Uk(t) + l (Uk-1(t))2 - (Uk(t))2 ,k = n + 1,†, uk(0) = gk і 0,k = 0,1,2,†,n,Uk(0) = gk і 0,k = n + 1,†, (8) where ϵ 0 is a small parameter that bring a singular perturbation to the system (??), which allows us to describe the processes of rapid change of the systems. 3. Queueing Systems with Finite Number of Servers and a Small Parameter Using (??) we can rewrite system of differential equations of order N in the form uМ‡k(t) = m uk+1(t) - uk(t) + l (uk-1(t))2 - (uk(t))2 ,k = 0,1,†,n - 1, uМ‡n(t) = m Un+1(t) - un(t) + l (un-1(t))2 - (un(t))2 , ϵUМ‡k(t) = m Uk+1(t) - Uk(t) + l (Uk-1(t))2 - (Uk(t))2 ,k = n + 1,†,N. (9) For this truncation system of ordinary differential equations of order N we can formulate Cauchy problem in the form uМ‡k(t) = m uk+1(t) - uk(t) + l (uk-1(t))2 - (uk(t))2 ,k = 0,1,†,n - 1, uМ‡n(t) = m Un+1(t) - un(t) + l (un-1(t))2 - (un(t))2 , ϵUМ‡k(t) = m Uk+1(t) - Uk(t) + l (Uk-1(t))2 - (Uk(t))2 ,k = n + 1,†,N, uk(0) = gk і 0,k = 0,1,2,†,n,Uk(0) = gk і 0,k = n + 1,†,N. (10) The numerical analysis was carried out using the adaptive step Runge-Kutta integration method, which is one of the most commonly used methods for the numerical solution of the singularly perturbed system of differential equations. The numerical example is presented in the figure (see Fig. ??, ??) where n = 7, N = 10, l = 0.5, m = 1.0, g0 = 1, gk = 1 - 0.1k, k = 0,9ВЇ and a small parameter ϵ = 0.1 (Fig. ??), ϵ = 0.01 (Fig. ??), ϵ = 0.001 (Fig. ??). In these numerical examples we can see the existence of regularly perturbed solutions ui(t), i = 0,5ВЇ and singularly perturbed solutions ui(t), i = 6,10ВЇ with boundary layers. pict Figure 1. Evolution analysis of uk (ϵ = 0.1) pict Figure 2. Evolution analysis of uk (ϵ = 0.01) pict Figure 3. Evolution analysis of uk (ϵ = 0.001) The numerical simulation shows that the solution of the singularly perturbed systems of differential equations have an area of rapid change of the function, which is usually located in the initial point of the problem. This area of rapid function change is called the area of the mathematical boundary layer. The thickness of the boundary layer depends on the value of a small parameter, and when the small parameter decreases, the thickness of the boundary layer decreases. The integration area is divided into external (outside the boundary layer) and internal (inside the boundary layer). The solution of the singularly perturbed equation is sought in the form of a solution suitable for the outer domain, which is then refined in the vicinity of the boundary point where the boundary layer is located. The numerical examples have shown the existence of regularly perturbed solutions and singularly perturbed solutions with boundary layers for evolutions ui(t). Conclusions We investigate the dynamics of large-scale queueing systems that consists of infinite number of servers with a Poisson input flow of requests of intensity Nl. Each request arriving to the system randomly selects two servers and this request is instantly sent to the one with the shorter queue. We suppose that service time has mean 1в€•m with exponential distribution. In this case a share uk(t) of the servers that have the queues lengths with not less than k can be described using a system of differential equations of infinite order. Tikhonov type Cauchy problem for this system with small parameter ϵ. Tikhonov type Cauchy problem for this system with small parameter ϵ and initial conditions is investigated. We use the simulation methods for behavior solutions analysis with taking into account parameters l, m, ϵ. The numerical examples have shown the existence of steady state conditions for queueing systems with an infinite number of servers and singularly perturbed conditions with boundary layers for evolutions ui(t) with a small parameter.[N. D. Vvedenskaya, R. L. Dobrushin, F. I. Karpelevich, Queueing System with Selection of the Shortest of Two Queues: An Asymptotic Approach, Probl. Peredachi Inf. 32 (1) (1996) 20-34, in Russian.][N. D. Vvedenskaya, Yu. M. Suhov, Dobrushin’s Mean-Field Approximation for a Queue with Dynamic Routing, no. 3, 1997.][N. D. Vvedenskaya, Large Queueing System where Messages are Transmitted via Several Routes, Vol. 34, 1998, in Russian.][L. G. Afanassieva, G. Fayolle, S. Yu. Popov, Models for Transportation Networks, Journal of Mathematical Sciences 84 (3) (1997) 1092-1103.][D. V. Khmelev, V. I. Oseledets, Mean-Field Approximation for Stochastic Transportation Network and Stability of Dynamical System (Preprint No. 434), University of Bremen, Bremen, 1999.][D. V. Khmelev, Limit Theorems for Nonsymmetric Transportation Networks, Vol. 7, 2001.][V. V. Scherbakov, Time Scales Hierarchy in Large Closed Jackson Networks (Preprint No. 4), French-Russian A. M. Liapunov Institute of Moscow State University, Moscow, 1997.][V. I. Oseledets, D. V. Khmelev, Global Stability of Infinite Systems of Nonlinear Differential Equations, and Nonhomogeneous Countable Markov Chains, Vol. 36, 2000.][E. A. Chernavskaya, Limit Theorems for an Infinite-Server Queuing System, Vol. 99, 2015.][E. A. Chernavskaya, Limit Theorems for Queueing Systems with Infinite Number of Servers and Group Arrival of Requests, Vol. 71, 2016.][Yu. Gaidamaka, E. Sopin, M. Talanova, Approach to the Analysis of Probability Measures of Cloud Computing Systems with Dynamic Scaling, Vol. 601, 2016.][A. V. Korolkova, E. G. Eferina, E. B. Laneev, I. A. Gudkova, L. A. Sevastianov, D. S. Kulyabov, Stochastization of One-Step Processes in the Occupations Number Representation, 2016.][K. Samouylov, V. Naumov, E. Sopin, I. Gudkova, S. Shorgin, Sojourn Time Analysis for Processor Sharing Loss System with Unreliable Server, Vol. 9845, Springer Verlag, 2016.][G. O. Bolotova, S. A. Vasilyev, D. N. Udin, Systems of Differential Equations of Infinite Order with Small Parameter and Countable Markov Chains, Vol. 678, Springer Verlag, 2016.]