Transient Behavior of the MAP/M/1/N Queuing System

: This paper investigates the characteristics of the MAP/M/1/N queuing system in the transient mode. The matrix method for solving the Kolmogorov equations is proposed. This method makes it possible, in general, to obtain the main characteristics of the considered queuing system in a non-stationary mode: the probability of losses, the time of the transient mode, the throughput, and the number of customers in the system at time t . The developed method is illustrated by numerical calculations of the characteristics of the MAP / M /1/3 system in the transient mode.


Introduction
Queuing theory is used for the mathematical descriptions of a large number of problems in calculating the performance characteristics of telecommunications and distributed computing systems. As a rule, researchers are interested in the behavior of the system under consideration in a stationary mode, when the probabilities of the states of the system do not depend on time [1][2][3]. However, with the transition to high-speed optical information processing systems and 5G/6G communication systems, for a more accurate assessment of the characteristics of the system, it is necessary to take into account its behavior in a non-stationary mode. Such a mode of operation occurs, for example, when the switching equipment is rebooted, when the base station is turned on and off, when the equipment is switched from the main channel to the backup channel, in cases of equipment malfunction, and in other cases associated with changes in system states. In this case, one of the most important tasks is to determine the duration of the transient mode, that is, the time during which the system will go to a stationary state. This problem is also relevant in the simulation of queuing systems, when it is necessary to determine the moment of transition of the system to a stationary mode with a sufficiently high accuracy.
The transient operation mode of queuing systems was first studied in [4][5][6][7]. In [6], the author for the first time described the need to study the transient operation modes of telecommunication systems and introduced the concept of the time constant of the transient process as an important parameter for determining the rate of transition of the system to the stochastic equilibrium mode. He considered several examples that require taking into account the changes in the probabilities of system states from time to time. As an example, the problem of predicting the characteristics of the system when installing the system or disrupting its operation was presented.
In recent years, queuing systems with periodically varying parameters of input flows and service time for customers have been actively studied [8,9]. One paper [8] considered the problem of convergence for a nonstationary two-processor system with catastrophes, server failures, and repairs, when all parameters are harmonic functions of time. Another paper [9] proposed three different analytical methods for calculating upper bounds for the rate of convergence to a stationary mode of a process given by a continuous inhomogeneous Markov chain. The first method is based on the logarithmic norm of a linear operator function, the second uses the simplest Lyapunov functions, and the third uses differential inequalities. Note that the stability of processes described by Markov chains with periodic parameters was studied back in the 1980s [10,11].
Reference [12] presents an approximate approach to the analysis of the transient behavior of a M/M/1/N system. The time dependencies of the probabilities of the system states are presented for special cases only in numerical calculations. The papers [13,14] present analytical expressions for finding the dependencies of state probabilities on time for the M/M/1 queuing system. The solution to the Kolmogorov system of differential equations is sought using the Laplace transform. The papers [15,16] present an analysis of the transient mode of a single-line queuing system with catastrophes. Analytical expressions are given for finding the dependencies of the probabilities of the states of the systems under consideration on time, and their stationary values. The paper [17] investigates a transient operation mode of a queuing system with heterogeneous servers and impatient customers with a Poisson input stream. The most complete review concerning the study of the transient behavior of queuing systems with Poisson input flows is presented in [18]. In this paper, the queuing systems M/M/∞, M/M/1, M/E/1, and M/M/1/N with catastrophes are investigated.
Despite the described approaches to the analysis of the non-stationary modes of queuing systems with Poisson input flows, these results cannot be used in the design of various real telecommunication systems. This is due to the fact that the traffic of modern telecommunication systems is correlated, and for a more accurate description, one should use not the simplest flow, but Markov correlated MAP or BMAP streams. The monograph [1] presents a systematized presentation of research methods and estimations of stationary characteristics of queuing systems with correlated flows. However, nonstationary modes of queuing systems with MAP flows are poorly studied in the literature. Note, for example, the work [19], which considers a single-line queuing system with correlated input flows, for which a numerical calculation of the state probabilities versus time using the Runge-Kutta method is presented. This paper, for the first time, proposes an analytical method for studying the transient behavior of the MAP/M/1/N system. The paper is structured as follows. Section 2 gives the formulation of the considered problem of studying the unsteady modes of the MAP/M/1/N system. Section 3 presents an analytical method describing the behavior of the MAP/M/1/N system in a non-stationary mode. Section 4, for the first time, gives expressions for the main parameters of the system under consideration in the transient mode: the probability of losses, the probabilities of the system states, the time of the transient mode, the throughput, and the number of customers in the system at time t. Section 5 presents numerical calculations illustrating the proposed analytical approach for the MAP/M/1/3 system.

Statement of the Problem
This paper discusses an MAP/M/1/N queuing system with a Markov input flow of customers and a limited number of waiting places. The arrival of customers is controlled by an irreducible non-periodic Markov chain ν t , t ≥ 0 with continuous time and a state space 0, 1, . . . , M − 1. The residence time of the chain ν t , t ≥ 0 in a certain state ν has an exponential distribution with the parameter λ ν , ν = 0, M − 1. After the time spent by the process in this state has expired, with the probability p (1 (ν, ν ) the process ν t , t ≥ 0 goes to some state ν and a customer is generated, with the probability p (0) (ν, ν ) that the process makes a transition, but without generating a customer. It is assumed that a jump from the state ν to the same state is impossible without generating a customer, i.e., p (0) (ν, ν) = 0. It is also assumed that the probabilities p (k) (ν, ν ), k = 0, 1 satisfy the normalization condition The flow is described by two nonzero M × M-matrices D 0 and D 1 : The monograph [1] describes in detail the meanings of the elements of the matrices D 0 and D 1 . The service time of a customer in the system is exponentially distributed with the parameter µ.
The transition graph of the MAP/M/1/N system is shown in Figure 1. According to this graph, the queuing system MAP/M/1/N can be in one of K = M(N + 2) states. The vertices of the graph S 0 , S 1 , ..., S N+1 , designate the macrostates of the buffer and the server. Transitions between these macrostates are accompanied either by the generation of a new customer, or by servicing the old one. Thus, the system is in the S 0 macro state if there are no customers in it; that is, the server is idle and the buffer is empty. The system is in the S 1 macro state if the server processes one customer and the buffer is empty. The system is in the S N macrostate, if the server processes one customer, while there are N − 1 more customers in the buffer. The system is in the S N+1 macrostate if the server is busy, provided that there are N customers in the buffer, and the next incoming customer will be discarded. Each of the above macrostates S k , k = 0, N + 1 corresponds to M additional states S for the case M = 3. Note that transitions within macrostates S k , k = 0, N + 1 correspond to the transitions of the MAP with no customer arrivals, and these transitions are not shown in the generalized graph ( Figure 1). Additionally, the transitions between macrostates correspond to transitions of the MAP with customer arrival (M possible transitions) or the service completion with rate µ (also M possible transitions between states). It is easy to see that the resulting intensity of direct transitions between macrostates is defined as ij . In Figure 1, they are specified separately for each intensity λ i : ij . Reverse transitions between macrostates in the MAP/M/1/N system are determined both by the number of Markov flows M and by the service intensity µ : µ g = (M − 1)µ (see Figure 1).
Following [1], the Kolmogorov equation system for the MAP/M/1/N system under consideration can be written as in the matrix form where ) of the server and buffer, determined by the dimension of the matrices (2) and (3). Note that A T is an infinitesimal generator and A T − → e = − → 0 . Here − → e is the unit column; − → 0 is the zero column.

The Method of Studying the Non-Stationary Characteristics of the MAP/M/1/N System
First of all, we need to formulate the main theoretical provisions necessary for studying the non-stationary mode of the MAP/M/1/N queuing system.

Analysis of the Characteristic Polynomial Roots
This subsection analyzes the roots of the characteristic polynomial of the Kolmogorov equations (4) (the eigenvalues of the coefficient matrix (5)), which describe the MAP/M/1/N system. This analysis is of great importance for describing the behavior of this system in a transient mode, since these roots determine the duration of the transient mode and the kind of the dependence of the state probabilities on time. According to [20][21][22], the characteristic polynomial of the Kolmogorov system (4) of order K has the form Here b j are the coefficients expressed in terms of the elements of the matrix (5) [20][21][22][23][24], and γ are the roots of the characteristic polynomial of system (4).
Since the probabilities in (4) are linearly dependent, and accordingly, (D 0 + D 1 ) − → e = − → 0 [1], the determinant of the coefficient matrix (5) is equal to zero [24]. Taking into account the well-known equality b 0 = det A [24], it can be obtained from (6) that Thus, one of the roots γ 0 of the characteristic Equation (6) is necessarily equal to zero. Since matrix (5) is real, the roots of characteristic Equation (6) can be either real or complex conjugate [24].
In addition, based on Sylvester's theorem [25][26][27], the number of negative diagonal elements of matrix (5) is equal to the number of negative eigenvalues of this matrix. Since all diagonal elements (5), taking into account (2) and (3), are negative, all roots of the characteristic equation have negative real parts.
Thus, the roots of the characteristic polynomial of the Kolmogorov system of differential equations describing the behavior of the MAP/M/1/N system can be either real negative or pairwise complex conjugate with negative real parts.
Consequently, the transient mode in the MAP/M/1/N system always has a decaying exponential character (Figure 2a), if the characteristic numbers are real negative or of decaying wavy character [20][21][22] (Figure 2b), if at least two eigenvalues are complex conjugates with negative real parts.
(a) decaying exponential type (b) decaying wavy type In accordance with [20][21][22][23][24], the roots of the characteristic equation (Equation (4)) can be either simple or multiple. This study, for the sake of definiteness, considers only the case of simple roots, which is most often encountered in practice.

The Matrix Method for Solving the Kolmogorov Equations
The papers [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18] present various methods of studying the non-stationary mode of queuing systems. All of them are reduced either to a numerical solution of the corresponding Kolmogorov system, or to the use of Laplace transforms. All of them do not allow obtaining a general solution to the Kolmogorov system of equations.
This paper presents the analytical approach based on the use of the so-called probability translation matrix. Definition 1. The probability translation matrix L(t) of system states is the matrix that relates the state probabilities P(t) at an arbitrary moment of time t with the probabilities of the states P(t 0 ) of this system at the initial moment of time t 0 : Obviously, the dimension K of the matrix L in (8) is equal to the order of system (4). However, the rank of this matrix is less by one than the order of system (4), since the probability of one of the states is expressed in terms of the remaining probabilities. Note that in this case, it is not necessary to fulfill the condition t 0 = 0. In this case, the initial moment of time is understood as the moment of the beginning of observation over the behavior of a random process. Proof. First of all, it should be noted that system (4) is overdetermined under the condition that there exists a stationary state; therefore, the rank of the fundamental matrix is less by one than the dimensions of the coefficient matrix of the system (4), and one of the solutions of the Kolmogorov system is a linear combination of the others. Additionally, solutions to system (4) can be written in the form where ξ kl are the elements of the eigen-basis of the coefficient matrix of system (4) (the coefficients relating the probabilities of states); − → C is the column of integration constants determined by the initial conditions. Then for t = t 0 : From (11) we have Substituting (12) into (9), one can obtain where is the translation matrix. Thus, it follows from (14) that the matrix L(t) relates the states of the process at the time t with these states at the time t 0 . Moreover, if Y(t 0 ) = I, where I is the unit diagonal matrix, then L(t) = Y(t).
Obviously, the matrix (15) is neither a Cauchy matrix nor a fundamental matrix [22], since the system (4) is overdetermined, and the rank of these matrices is less by one than the dimensions of the matrix (15). Moreover, as t → ∞, determinant (15) tends to zero; and the determinant of the fundamental matrix is the Wronski determinant and is never equal to zero [24].
An important advantage of the presented approach is a decrease in the computational complexity associated with transformations that exclude the linear dependence of solutions when using fundamental matrices and calculating the Jordan forms of matrices. Moreover, this approach does not require the use of Laplace transforms and allows one to effectively solve synthesis problems (the so-called inverse problems). Now it is important to find the general analytical form of the solution for the queuing system under consideration. In the most general case, this solution is found as an exponent of the coefficient matrix of the initial differential equation [24]. However, the calculation of the matrix exponent requires finding the Jordan form of the original matrix and the matrix of eigenvectors, which for large sizes of the system is a very cumbersome operation and gives large errors. Moreover, this approach does not allow writing the final expression in an analytical form. This paper proposes the technique for finding this solution based on Theorem 1 without using the procedure for calculating the matrix exponent.  (7), has the form where the elements of the K × K matrix A j are written as ∆ lj is the determinant of the algebraic complement to the element ξ lj of the matrix ∆, ∆ is the proper basis of the matrix of coefficients of system (4) in K-dimensional space [22][23][24], and γ j is the j-th root of the characteristic polynomial (7) of system (4).
Proof. As noted above, the general solution of the system of linear homogeneous differential equations (4) with constant coefficients and simple roots of the characteristic equation has the form [20][21][22][23][24] P j (t) = One should look for the elements of the transformation matrix under the condition Y(t 0 ) = I (see (15)). Then the constants C j for the elements of the first column are found under the condition P i (t 0 ) = (1, 0, . . . , 0) T ; i.e., the system of equations to be solved has the form Hence, one can obtain Substitution of the C j1 found in (19) gives the elements of the first column of the transformation matrix. Similarly, the elements of the second column are found under the condition P i (t 0 ) = (0, 1, · · · , 0) T , as well as the elements of the K − 1 column under the initial conditions P i (t 0 ) = (0, 0, · · · , 1) T .

System Characteristics in the Transient Mode
The main indicators of the system in the transient mode are the probability of losses, the time of the transient mode, throughput, and the number of customers in the system [17,18]. As a rule, these indicators for various systems are studied in a stationary mode. The paper [6] presents relations for calculating the so-called time constant of a random process, but provides no substantiation for this parameter, even for the M/M/1/N system. Moreover, there are no ratios for calculating the transient time. The paper [18] considers the problem of the maximum queue size on the interval [0, t], the probability of losses, and the number of customers in the system at time t for the systems M/M/1 and M/M/1/N. The author used the uniformization procedure, passing from continuous to discrete time, which significantly complicated the calculations. At the same time, the calculation of the throughput of queuing systems in the transient mode has not yet been considered in the world scientific literature known to the authors.

Loss Probability
The expression discovered for the transformation matrix (16)  N+1 in this moment in time [1].
where P where L ij (t) are elements of matrix (16) at the time t.

Transient Time
Determining the transient time is one of the most important unsolved problems in the theory of MAP flows. It was indicated above that the values of the characteristic numbers fully described the duration and nature of the transient mode. Now let us consider the duration of the transient mode. In the general case, based on the analysis carried out in Section 3.1, the probability of the i-th state is where γ j < 0, j = 1, K − 1. That is, the overall process is a superposition of K-1 damped processes. If all γ j except γ 0 are real negative, then at t → ∞, P i (t) → ξ i1 C 1 . Each of the exponential terms in (26) decreases by a factor of e over some time τ j . This time will here be called the time constant of the j-th component of a random process with an input MAP flow, by analogy with the time constant borrowed by Harison [6] from the theory of electrical circuits for M/M/n queuing systems. However, it should be noted that the physical meaning of this quantity was not described in [6]. This constant can be calculated for real τ j based on the equality exp(−|γ j |τ j ) = 1/e = e −1 , taking into account that γ 0 = 0, ∀γ j < 0, j = 1, K − 1, whence, τ j = 1/|γ j |. In this case, a decrease in γ j leads to an increase in τ j . Thus, the term in (26) with the minimum γ j has the maximum time for decreasing by a factor of e. It also determines the transient mode duration.

Definition 2.
The time constant of a random process with an input MAP flow with constant parameters λ i , p kl , µ is the time during which the component of the process with the smallest α j decays by e times. Its value is equal to the inverse smallest real part of the nonzero eigenvalue . Based on the mathematical formulation of the problem, the transient mode ends at → ∞, but from a physical point of view, the time of the transient mode can be considered finite at |P j (t) − π j | ≤ , where is some infinitesimal value, and π j is the stationary probability of the j-th state. The choice of is determined by the practical requirements for specific queuing systems. Now, it will be assumed that the mode can be considered stationary when the value of the j-th term decreases in (26) with γ ( j) = α jmin + iβ j in e 3 . Then, the transient time is defined as τ tr = 3/|α jmin |. Similarly, for the damping of the transient process by e 5 times, one can write τ tr = 5/|α jmin |. Definition 3. In the general case, the time of the transient process can be defined as τ tr = k/|α jmin |, where ∀α j ∈ Γ : (Γ = α j , α j ≥ α min ⇒ α j = α min . Here k(k > 0, k ∈ R) is defined as the coefficient of change in the state probability.
The numerical simulation presented below has shown that for practical calculations, k = 5 ÷ 10 is sufficient, depending on the formulation of the problem.
When considering the nature of the transient mode, the ratio of α j and β j is important. Indeed, if two roots of the characteristic equation are complex conjugate, then one of the terms in the elements of matrix (16) has the form exp(αt) cos(β j t). Thus, β j determines the oscillation frequency of state probabilities in the transient mode. If β j α j , then oscillations are practically not observed. If β j < α j , then one half-period of oscillations can be observed. If β j ≈ α j , then one oscillation period can be observed. If β j α j , then the process in the transient mode has a wavy character.

Throughput
The found probability of losses P loss (t) (25) makes it possible to calculate the throughput of the system at a given time in accordance with the definition [3] as Here λ = θD 1 − → e is the average intensity of the input flow [1]. Taking into account matrix (16): the average value of the throughput during the transient process is

The Number of Customers in the System at the Time t
Knowing the probability that the system is idle, that is, both the buffer is empty and the server is idle, one can determine the mathematical expectation of the average number of customers in the system at the time t [17] by the formula: Indeed, the probability of the system macrostate S i at the time t is equal to ∑ M−1 j=0 p (j) i (t) (Section 2). In this case, i customers in the system correspond to the macrostate S i . Then the mathematical expectation in accordance with the definition can be calculated by formula (30). Then the average statistical value of the number of received customers in the system during the transient process is Obviously, the integral in (31) is calculated analytically, since the integrand is the sum of exponentials. Integral (31) also has an analytical solution. However, the resulting analytical expressions are not given here due to their cumbersomeness.

Examples of a Numerical Study of the Characteristics of the MAP/M/1/N System
This section considers the process of servicing three types of traffic with the arrival rates λ 0 = 5 packets/s, λ 1 = 583 packets/s, and λ 2 = 833 packets/s. Assuming that the input flows are correlated, the values of the flow rates can be set by the matrix Λ = diag{λ 0 , λ 1 , λ 2 } = diag{5, 583, 833}. In this problem, transitions from the transmission state of one traffic with intensity λ n to the transmission state of another type of traffic with intensity λ m without generating customers are given by the transition probability matrix It is taken into account that a jump from one state to the same state is impossible without generating a customer-i.e., p ν,ν = 0; and the probabilities satisfy the normalization condition ∑ 2 j=0 p Thus, the information flow from a multimedia device can be specified as a MAP-flow described by two nonzero matrices D 0 = Λ(P (0) − I), D 1 = ΛP (1) [1]. It follows from the problem statement that the number of states of the MAP governing process is M = 3. In addition, it can be assumed that the server is a switch with a buffer size of N = 1. The transition graph for the case under consideration is shown in Figure 3. Here the states S 2 for video, and S (2) traffic in the buffer. In view of the above, transitions from one state to another are given by the matrices For the research, a software package in the Python language was developed. Let us investigate the dependence of the state probabilities on the ratio of the intensities of packet arrival and traffic service. To do this, consider the following three cases for all real γ k : (1) The average service rate µ = 5 packets/s; the average rate of packet arrival in accordance with [1]  (2) The average service rate µ = 11.3 packets/s; the average packet arrival rate λ = 139.04 packets/s. For this case, the modulus minimum characteristic exponent is γ min = −14.7. The results of calculating the probabilities of states are presented in Figure 5. The probability of loss, varying from zero to 0.41 during the transient, is equal to the probability that the buffer and the server are idle in the stationary mode. Solid lines indicate the probabilities of states. The time constant of the transient is τ max = 1/|γ min | = 1/14.7 = 0.068 s, and the transient time is τ tr = 0.34 s for k = 5.
(3) The average service rate µ = 33.33 packets/s; the average packet arrival rate λ = 139.04 packets/s. As a result of the calculation, it is found that γ min = −28.55. The results of calculating the state probabilities are presented in Figure 6. Solid lines indicate the probabilities of system states. The probability of losses during the transient mode increases from zero at t = 0 s to 0.26 at t > 0.18 s (see Figure 6). The probability that the system is free changes from unity at t = 0 s to 0.64 at t > 0.18 s. The time constant for this case is τ max = 1/|γ min | = 1/28.55 = 0.035 s, and the transient time is τ tr = 0.18 s for k = 5.
The results of studying the throughput of the system in the transient mode at µ = 5 packets/s and λ = 33.33 packets/s are presented in Figures 7 and 8, respectively. The throughput in the first case decreases from A(0) = 139.04 packets/s to A(t tr ) = 59.9 packets/s, and in the second case from A(0) = 139.04 packets/s to A(t tr ) = 103.4 packets/s. Thus, with an increase in the service rate µ, the transition time decreases. It should be noted that the transient time is long for high-speed data transmission systems, amounting to tenths of a second. In addition, the throughput of the system at the beginning of the system startup or reboot significantly exceeds the stationary value, which requires large resources for its processing. It follows from this that analysis of the characteristics of the system in the transient mode is necessary for the correct choice of the system configuration that provides the required quality of service indicators.    Consider the cases when any of γ k , k = 0, K − 1 is complex. This case corresponds, for example, to λ 0 = λ 1 = λ 2 = 833 packets/s, λ = 472.3 packets/s, µ = 11.3 packets/s ( Figure 9). In this case, there are four pairwise complex conjugate roots of the characteristic equation: γ 1,2 = −983.25 ± i81.86, γ 3 , 4 = −894.92 ± i74.1. The transient time is τ tr = 0.01 s. Note that for these values of τ k in the transient mode, there are bursts of probabilities that the server and the switch buffer are idle, along with the likelihood that the server is busy and the buffer is free for all types of traffic. At the same time, the values of some probabilities differ significantly from their values in the stationary mode. Thus, only the analysis of the behavior of the MAP/M/1/N system in a stationary mode can lead to an incorrect assessment of the system state.
Increasing the size of the switch buffer obviously leads to a decrease in the probability of losses; however, it also causes a significant increase in the transient time. Thus, Figure 10 shows the results for the previous values of the intensities and transition probabilities, and N = 8. The transient time is already τ tr = 0.025 s with a slight decrease in the probability of losses.

Conclusions
This work was devoted to an analytical study of the transient mode of the MAP/M/1 /N system. First of all, based on the analysis of the Kolmogorov system of equations, it was shown that in such a system, there will always be a stationary mode. The paper presented for the first time an analytical method for studying the MAP/M/1/N system in a transient mode, based on the concept of the so-called transformation matrix. Numerical methods were used only when calculating the eigenvalues of the coefficient matrix of the Kolmogorov system of equations for large-scale systems. Indeed, an algebraic equation higher than the fourth order has no analytical solution. The paper presents analytical expressions for calculating the time of the transient mode and the analysis of the possible behavior of the system in the transient mode, depending on the values of the roots of the characteristic equation of the system of Kolmogorov equations. Analytical expressions for calculating the probabilities of system states in a transient mode were presented. Analytical expressions were also presented for calculating the main indicators of the system operation, such as the probability of losses, throughput, the number of applications in the system at an arbitrary point in time, and the average values of throughput and the number of applications in the transient mode. Numerical modeling was carried out using the example of the MAP/M/1/1 system.

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