A queueing model with servers disguised as customers

We propose a new queueing model, motivated by the phenomenon of pseudoprogression in cancer, in which the length of a queue appears to increase initially, before reducing to a steady state. We assume that servers arrive to the queue alongside the customers, i.e. `disguised' as customers. We derive the general equations for this model using matrix analytic methods, and demonstrate its behaviour with numerical simulations.


Introduction
Numerous variations of queueing models have been proposed which take into account heterogeneous servers, multiple types of customers and priority queues.In this paper, we propose a new queueing model, in which servers arrive disguised as customers.The motivation behind this model arose from the documentary "Jim Allison: Breakthrough" (Haney, 2020), which spoke about the phenomenon of pseudoprogression in cancer immunotherapy treatments.Pseudoprogression is defined as a misleading increase in tumor size or appearance in new lesions, followed by a decrease in tumor burden.The misleading increase is due to treatment which increases the apparent tumor size.A thorough review of pseudoprogression and its comparison with true progression of a tumor is given in Jia et al.( 2019) and Ma et al.(2019) respectively.In a similar manner, our proposed queueing model appears to increase in size for a period of time due to the (disguised) arrival of servers among the customers.When these disguised servers begin serving, the queue size (rapidly) decreases.The addition of servers to the queue, disguised as customers, may be thought of as adding heterogeneous customers to the system.In section (2) we review the existing literature on queues with heterogeneous customers and servers, as well as matrix geometric solutions for queueing theory.In section (3), we introduce notation and the phase diagram for the model, and derive the equations for our model.Section (4) involves numerical simulations for the simplest case of our model using matrix analytic methods, as well as the variation of the queue properties for different values of the parameters.In section (5), we present concluding remarks.

Literature Review
There is a considerable amount of literature on queueing models with heterogeneous customers.Klimenok et al.(2020a) present a queueing model with two types of customers, where the priority changes after a random amount of time.A variation of this can be found in Klimenok et al.(2020b) which details a multi-server priority queue with marked Markovian heterogeneous arrivals.A detailed description of a queueing system with heteregenous customers is given by Collings (1974), in which he derives the probability flow through states by performing "cuts" at different points in the system.Customers belong to different groups, having different arrival and service rates, and it is shown that approximating these rates by a single exponential distribution underestimates the queue length.Gurtov and Mazalov (2012) presented a queueing system with servers that arrive when the number of customers exceeds a threshold.Yechiali and Naor (1971) put forward a queue with both heterogeneous arrivals and service times.In their model, arrival intensities can "jump" between two different levels with differing rates.The time spent in each level is exponentially distributed.Kim et al. (2011) examine the problem of assigning customers to different servers (having varied service rates) so as to minimize the waiting time of a customer in the case of priority queueing.Li and Yang (2000) give a model in which the number of servers in the queue is dynamic, a concept that we use in our system, and depends upon the number of customers present.Servers are added or removed on the basis of a search and release process.
We use matrix analytic methods for our model, a review of which is given by Neuts (1984).The procedure for solving such systems is given in Nelson (1991) and Stewart (2021).

Notation
Before defining the model, we introduce some preliminary notation: Customers arrive into the system according to a homogeneous Poisson process with rate λ c , where they're served on a first-come-first-serve basis.Service times are independent, exponential random variables with rate µ.In our model, servers arrive in the system along with customers, with rate λ sn , and leave the system with rate µ sn .In the general case, we assume that the arrival and departure of servers can depend on the number of customers (n) present in the system -for example, the larger the number of customers present in the system, the faster a new server would arrive.
We denote the states of the system by (n, k), where n ∈ {0, 1, 2, 3, ...} represents the number of customers and k ∈ {0, 1, 2, ...K}, the number of servers.We restrict the number of servers in the system to K, and for the sake of simplicity, we derive the equations for K = 3 servers.For each value of n, we have 0, 1, 2 or 3 servers in the system.

Phase Diagram
For each state (n, k), n is known as the level state, while k is known as the interlevel state.For each level, we have 4 interlevel states.Below is a phase diagram of the process: The dashes beside the numbers in the above diagram represent the number of servers.Here, we assume that a server leaves as soon as it is not needed, as can be seen from the transition from (1, 2) → (1, 1).Thus, we do not have the states (0, 2) as we reason that a transition cannot occur from (1, 2) → (0, 2).We will now derive the time dependant probability equations for the model.

Matrix Analytic Solution of the Model
We now proceed to derive the equations for the most general case of the model.The generator matrix for the system is given by: We can write Q in the form of a block matrix as follows: where: are known as the boundary conditions.We also have the matrices: where n ≥ 3. Notice that the matrix A 0 represents the movement to a lower level and A 1 represents the movement to a higher level.
Let π = (π 0 , π 1 , π 2 , ...) denote the limiting distribution of the system, where π i denotes the limiting distribution of the i th level.Then, we have the equations: where e denotes a column vector of 1's.Solving the above two equations gives us a set of boundary conditions: For i ≥ 3, we have the following recursive relation: To solve the above equations, we need to express A n in terms of the values we know.In the following section, we take the simplest case, that is, when λ sn = λ s and µ sn = µ s and derive the properties of the corresponding queue.
3.4 Simplest Case: λ sn = λ s and µ sn = µ s In this case, A n becomes a constant matrix given by: where λ = λ c + λ s .Our generator matrix becomes: The recurrence relation equation ( 7) reduces to the form: Note that in the matrix Q, once we reach the third level (n ≥ 3), the matrix starts repeating, precisely like a Quasi Birthand-Death process.We can now solve for the limiting distributions.

Condition for Stability
Let the generator matrix A be given by A = A 0 + A 1 + A 2 .This forms the generator for the repeating section of the system.To derive a condition for the stability of the queue, we first find the stationary distribution for A, given by π A = (π A0 , π A1 , π A2 , π A3 ).For this, we use the equations: Then, we have: From Equation ( 12), and the normalization equation, we get: For ergodicity, we require the drift to lower levels to be greater than that to higher levels.That is, the system will be more inclined to go to lower level.

Solving for Limiting Distribution
Once we check whether the queue is stable for given values of λ c , µ, λ s and µ s , we can find the limiting distribution.From Equation ( 10), since all three matrices A 0 , A 1 and A 2 are constant, we can assume that for i ≥ 3, π i = π i−1 R, where R is a constant matrix.Repeatedly substituting in the recurrence relation, we get: The simplest method to solve for R is to iterate numerically until the elements of the matrix differ by less than a specified value.Substituting the value of π i−1 , π i and π i+1 , the iterative formula is found as follows: We now solve the boundary equations ( 6) for π 0 , π 1 and π 2 .Substituting the value of π 3 = π 2 R, we can represent the equations in matrix form as: = (0, 0|0, 0, 0|0, 0, 0, 0) Note that π 0 = (π 00 , π 01 ), π 1 = (π 10 , π 11 , π 12 ) and π 2 = (π 20 , π 21 , π 22 , π 23 )We can rewrite the above equation as follows: Once we have values for π 0 , π 1 and π 2 , we normalize them to get probabilities: In the above equations, we take e merely to be a vector of ones of required dimension for each term.Once we have the boundary limiting distributions we can find any π n = π 2 R n−2 .In the next section, we find the properties of the queue.

Properties of the Queue
Expected queue length: For our model, the expected queue length includes the number of servers, since an observer of the queue cannot distinguish between a customer and a disguised server.Before we find an expression for the expected queue length, we first define: Then, the expression for expected queue length will be given by: Simplifying, we obtain our final equation: Expected Number of Customers Waiting in Queue: We denote this value by E(L n ).The derivation for the expected number of customers is similar to that for the expected queue length.
By rearranging the terms and changing the limits, we obtain: The expected waiting time can correspondingly be found by using Little's Law: ).
Delay Probability : The probability a customer has to wait, known as delay probability, is denoted by Π w , and it is derived as follows:

Numerical Simulations
We simulate the model for different values of λ c , λ s , µ and µ s and display the results of the model below.
1 1 4.474 0.756 0.378 0.557 We can make a few observations from the above values: • For a given µ s , as µ increases, the values of E(L n ) decrease.That is, as the service rate increases, fewer customers are found in the system.Naturally, the expected waiting time E(W) and the delay probability Π w will also decrease with an increase of µ.
• For a fixed µ, the increase of µ s resulted in an increase of customers waiting in the queue (E(L n )).This is consistent with the assumption that as servers leave the queue at faster rates, more customers will be left waiting.An interesting remark is that when λ c < λ s , the delay probability and waiting time decrease even if µ s increases.This could be attributed to the fact that though the rate of departure of servers increases, their rate of arrival makes up for it.Notice that when λ c > λ s , the expected waiting time and delay probability increase with µ s .
Example 1: Here, ρ c = 0.5, ρ s = 0.75 and (ρ s (1 149, hence the ergodicity condition is satisfiedthe queue will be stable. Using Equation ( 17), we get π A = (0.366, 0.274, 0.206, 0.154) with condition (19) indicating that the system is stable.We found R iteratively, till the difference between elements of successive iterations are less than 10 −12 .This gave R as: Now we proceed to solve the boundary conditions (21).We set π 00 = 1, that is, we start at the (0,0) state, and find the remaining 8 variables.Given the equation in the form Ax = 0, we can re-write this as A'x' = b, where b = −A[, 1], and A' is the matrix obtained by taking the submatrix A[2:9,2:9] of A (since we require 8 variables, we need only 8 equations, hence we take an 8x8 submatrix of A).Once we solve the matrix equation, we normalize the probabilities using Equation (18.Our final probabilities are: π 0 = (0.6557041, 0.1519753) π 1 = (0.052086275, 0.052086275, 0.005271742) π 2 = (0.030662973, 0.017641404, 0.006589678, 0.004503371) Once we have π 0 , π 1 and π 2 , we can find all other probabilities of the system as well as the properties: • The total expected queue length of the system E(L) = 0.576480174113887 • The expected waiting time of a customer W = 0.576480174113887 • The expected queue length (customers only) E(L n ) = 0.206882719667415 • The probability an arriving customer has to wait Π w = 0.182545493796219 A plot of the simulation of the instantaneous queue length is given below: The ergodicity condition is not satisfied here, and the queue length will blow up.This is illustrated in the simulation plot of instantaneous queue length below:

Figure 4 .
Figure 4. Instantaneous queue length for unstable queue