Generalized Hypercube Queuing Models with Overlapping Service Regions

These states can be represented by a hyperlattice in dimension $I$. Each node of the hyperlattice corresponds to a state $B=(n_{i})_{i\in\mathcal{I}}$ represented by a tuple of numbers, where non-negative integer $n_{i}\in\mathbb{Z}_{*}$ indicates the status of server $i$. Server $i$ is idle if $n_{i}=0$ and busy if $n_{i}>0$. The value of $n_{i}-1$ represents the number of calls waiting in the queue of server $i$ when the server is busy. Figure 8 (a) gives an example of the state space of a hyperlattice queuing model for a system with two servers. It is worth noting that the state space in a hyperlattice queuing model can be divided into two parts, with one part consisting of states that possess at least one available server (referred to as non-saturated states) and the other part comprising states in which all servers are busy (referred to as saturated states), as shown in Figure 9.

To create the transition rate matrix, the states must be arranged in order. However, since the standard vector space lacks an inherent order, we employ a touring algorithm to traverse the entire hyperlattice. This algorithm generates a sequence of $I$-digit non-negative integers $B_{0},B_{1},\dots$, which contains an infinite number of members and fully describes the order of the states.

We simplify the tour problem by decomposing it into $K$ sub-problems. Each sub-problem $k\leq K$ first identifies the states in which exactly $k$ calls are staying in the system (either waiting in the queue or being served), and then proceeds to solve the enumerative combinatorics problem by searching for $\binom{k+I-1}{I-1}$ possible combinations in a pre-determined sequence. It is important to note that the sequence generated by each sub-problem is just one of numerous possible tours of the cutting plane in the hyperlattice. Figure 10 presents a hyperlattice with $I=3$, along with its corresponding tour of all the states. The combinatorics can be translated into the “stars and bars” [10], i.e., finding all possible positions of $I-1$ dashed circles that separate $k$ balls into $I$ groups, as illustrated in Figure 10 (c). The number of balls being separated in group $i$ can be regarded as $n_{i}$. Algorithm 1 summarizes the tour algorithm that generates the sequence $B_{0},B_{1},\dots$. For convenience, we index the states based on the order produced by the tour algorithm. The state index is represented by $u\in\mathcal{U}=\{0,1,2,\dots\}$.

Now we define the transition rate matrix for the hyperlattice queue $Q=(q_{uv})_{u,v\in\mathcal{U}}$, where $q_{uv}$ ($u\neq v$) denotes the transition rate departing from state $B_{u}$ to state $B_{v}$. Let $u$ represent the index of state $(n_{1},n_{2},\dots,n_{i},\dots,n_{I})$ and let $d_{uv}=\left\lVert B_{u}-B_{v}\right\rVert_{1}$ be the $\ell_{1}$ distance (or Manhattan distance) between two vertices $B_{u}$ and $B_{v}$ on the hyperlattice. We define $u^{i+}$ as the index of state $(n_{1},n_{2},\dots,n_{i}+1,\dots,n_{I})$; note that thus we have $d_{uu^{i+}}=1$ and $B_{u^{i+}}\succ B_{u}$ ($B_{u}$ is elementwise greater than $B_{v}$). Diagonal entries $q_{uu}$ are defined such that $q_{uu}=-\sum_{v:v\neq u}q_{uv}$ and therefore the rows of the matrix sum to zero. There are two classes of transitions on the hyperlattice: upward transitions that change a server’s status from idle to busy or add a new call to its queue; and downward transitions that do the reverse. The downward transition rate from state $B_{u^{i+}}$ to its adjacent state $B_{u}$ is always $q_{u^{i+}u}=\mu_{i}$. The upward transition rate, however, will depend on how regions overlap with each other and the dispatch rule when a call is received, which can be formally defined by the following proposition.

The following argument justifies the proposition. The first term of (1) suggests that, for an arbitrary server $i$, it must respond to a call received in its primary service region $i$. The second term of (1) means that the server $i$ needs to respond to a call received in the overlapping service region $e$ according to the dispatch policy $\eta_{e,u}$. It is important to emphasize that the dispatch policy $\eta_{e,u}$ is state-dependent and can be defined based on the user’s needs, offering a general framework that is apt for elucidating both random and deterministic dispatch strategies. For example, the one outlined in Section 2.1 represents a specific case of the dispatch policy, in which the assignment of server $i$ to the call received in overlapping service region $e$ can be discussed in two scenarios: (a) When all servers $i\in e$ are busy ($u:n_{i}>0,\forall~{}n_{i}\in B_{u}$), this call will join the queue of server $i$ based on the probability $\eta_{e,u}(i)$. (b) Conversely, when there is at least one server idle ($u:n_{i}=0,n_{i}\in B_{u}$), the server $i$ will be assigned to this call based on the probability $\eta_{e,u}(i)$. In Appendix, we also demonstrate how this proposition can be distilled to represent a rudimentary system where, at most, only two servers are permitted to patrol one overlapping service region.

To obtain the steady-state probabilities of the system, we can solve the balance equations given the transition rate matrix $Q$. Let $\mathbb{P}\{B\}$ denote the probability that the system is occupying state $B$ under steady-state conditions. The balance equations can be written as:

We also require that the probabilities sum to one, namely,

However, the above balance equations cannot be solved computationally as the number of states on the hyperlattice is infinite and grows exponentially with the number of servers even when queue capacity is limited. To overcome this difficulty, we introduce an efficient estimation strategy in the next section.

We first exploit the fact that the hyperlattice queuing model in the aggregate is simply a birth-death model whose states can be represented by $\{N_{k}\}_{k=0,1,\dots,\infty}$ and $k$ can be interpreted as the number of calls currently in the system. This equivalence is obtained by collecting together all states having an equal value of $k=\left\lVert B\right\rVert_{1}$ as shown in Figure 8 (b). Their summed steady-state probability is equal to the comparable probability of state occurring in the birth-death model, i.e.,

The birth and death rates of the aggregated model can be defined by the following proposition:

Now, rather than attempting to calculate the stationary distribution across an infinite number of states in the original state space, the problem can be simplified by focusing only on a finite subset of states. These states are selectively considered if their $k$ values are less than a user-defined hyper-parameter $K$, or $\left\lVert B\right\rVert_{1}\leq K$, and the sum of their steady-state probabilities remains constant. The choice of $K$ will discussed below and verified by an example in Table 2. Proposition 3 provides a closed-form formula for the steady-state probability of the aggregate state $N_{k}$ when $k\leq K$.

It is worth noting that the probabilities of the “tail” of the hyperlattice (states with large $k$ value) are usually low and can be negligible since their chance of happening decreases with the length of the queue. As a result, it is empirically benign to ignore the tail of the hyperlattice and focus on the states where the total number of calls in the system is less than $K$. We represent this set of states by $\mathcal{U}_{K}=\{0,\dots,U_{K}\}$, where $U_{K}$ denotes the size of the set. In the subsequent discussion, Lemma 1 presents the equation for $U_{K}$ in relation to $K$ and $I$. This demonstrates that the size of the set can expand exponentially relative to the number of servers $I$, and it can become quite substantial even for moderate values of $K$, particularly in systems with multiple servers.

It is evident that we need to strike a balance between estimating a larger number of states (or opting for a higher value of $K$) and the resulting computational efficiency. Table 2 presents a comparison of various metrics for the hyperlattice model under the different number of servers $I$ and choice of $K$. These metrics include the number of states required for estimation, the sum of their steady-state probabilities, and the corresponding running time in seconds. Our findings reveal a trade-off: while estimating a larger number of states increases the cumulative steady-state probability and provides a more in-depth view of the system dynamics, but at the same time, it negatively impacts computational efficiency.

Denote the steady-state probabilities and transition rate matrix for the states in set $\mathcal{U}_{K}$ as $\pi_{K}=\{\mathbb{P}\{B_{u}\}\}_{u\in\mathcal{U}_{K}}$ and $Q_{K}=(q_{uv})_{u,v\in\mathcal{U}_{K}}$, respectively. We can express the balance equations in (4) as a compact form:

where the right-hand side is a column vector of zeros. By using equation (5), the sum of probabilities in set $\mathcal{U}_{K}$ denoted by $\sigma_{K}$ is equal to the sum of probabilities in set $\{N_{k}\}$ up to $K$:

Therefore, we can estimate $\pi_{K}$ by solving (9) and (10). However, there is redundancy in these equations as the number of equations is $K+1$ and the number of unknown variables is $K$. To eliminate this redundancy, we replace one row of $Q_{K}$ with a row of ones, usually the last row, and denote the modified matrix as $Q^{\prime}_{K}$. Similarly, we modify the “solution” vector, which was all zeros, to be a column vector with $\sigma_{K}$ in the last row and zeros elsewhere, and denote it as $e_{K}$. In practice, we solve the resulting equation:

We note that the solution to the above balance equations requires a matrix inversion of $Q_{K}$ (or $Q^{\prime}_{K}$), which can become computationally challenging for large matrix sizes. By Lemma 1, the size of $Q_{K}$ can reach a staggering $U_{K}\approx(K+I)^{I}/I!$ (using Stirling’s approximation [19]) even for moderate values of $K$ when the system has $I>2$ servers, presenting a considerable computational hurdle. Fortunately, the sparsity of transition rate matrices (as illustrated in Figure 12) can be leveraged to iteratively find the steady-state distribution using the power method [16], by constructing a probability transition matrix of a discrete-time Markov chain through uniformization of the original continuous-time Markov chain. Let $\mathbb{I}$ represent the identity matrix and denote the maximum absolute value on the diagonal of the transition rate matrix $Q_{K}$ as $\gamma=\max|q_{uu}|$. By initializing with a vector $\pi_{K}^{(0)}$, we can compute $\pi_{K}^{(t)}=\pi_{K}^{(t-1)}(\mathbb{I}+Q_{K}/\gamma)$ at each iteration $t=1,2,\dots$, until the distance between $\pi_{K}^{(t)}$ and $\pi_{K}^{(t-1)}$ falls below a specified threshold. The convergence of this approach is guaranteed, and the steady-state distribution can be efficiently computed despite the potentially large size of the transition rate matrix.

The workload of server $i$, denoted by $\rho_{i}$, is simply equal to the fraction of time that server $i$ is busy serving calls. Thus $\rho_{i}$ is equal to the sum of the steady-state probabilities on part of the hyperlattice corresponding to $n_{i}>0$, i.e.,

where $\mathbbm{1}$ denotes indicator function and $n_{i}$ represents the status of server $i$ in state $B_{u}$. In practice, because the efficient model estimation derived in Section 2.3 only computes the steady-state probabilities for the states in $\mathcal{U}_{K}$, the workload of server $i$ is numerically approximated by

The individual workloads can be further used to calculate various types of system-wide workload imbalance defined in [13].

For the remainder of the system performance characteristics, it is necessary to compute the fraction of dispatches that send a server to each sub-region under its responsibility. We use $\rho_{i,e}$ to represent the fraction of dispatches from sending server $i$ to one of its primary overlapping service regions $\mathcal{O}_{e}$. We have

In addition, we use $\rho_{i,i}$ to represent the fraction of dispatches that send the server $i$ to its primary service region $i$, which is consistently equal to $\lambda_{i}/\lambda$ since this region can exclusively be handled by server $i$. We note that $\rho_{i,e}=0$ if $i\notin e$ and $\sum_{i,k}\rho_{i,e}+\sum_{i}\rho_{i,i}=1$. Additionally, in practice, we consider $u\in\mathcal{U}_{K}$ rather than $\mathcal{U}$ in (12) as we only calculate the steady-state probabilities for the first $U_{K}$ states according to Section 2.3.

We also test our models on the synthetic systems that allow overlapping patrolling. Specifically, we create a two-server system where the service region is a 100 by 100 square with a single overlapping service region in the middle ($I=2$ and $E=1$), as depicted in Figure 5 (a). We assume that calls are uniformly distributed with an arrival rate of $\lambda=1$ and that the service rate of each server is fixed at $\mu_{i}=1$. Changing the hyperparameters of the system, such as the size of the overlapping region and the dispatch policy, can significantly shift the queuing dynamics. Our numerical experiments suggest that the hyperlattice queuing models can accurately capture intricate queuing dynamics under different hyperparameter settings for such a system.

To examine the impact of overlapping ratio on system dynamics, we first create 100 systems by varying the overlapping ratio $r$. This ratio represents the percentage of the area of the overlapping region with respect to the area of the entire service region of server 1, i.e., $r=|\mathcal{O}_{\{1,2\}}|/(|\mathcal{O}_{\{1,2\}}|+|\mathcal{P}_{1}|)$, as shown in Figure 13 (a-c). We then compare the estimated values of $\rho_{i}$ and $\rho_{i,e}$ to their true values (indicated by dashed lines) in Figure 13 (d) and (e), and find that the estimated values are a good approximation of their true counterparts. Clearly, the workloads of the two servers differ more with an increasing overlapping ratio, enlarging the area server 1 has to cover. Moreover, as the ratio increases, the probability of dispatching either server to the overlap region grows (indicated by green and orange lines), while the probability of sending server 2 to its primary area remains the same (blue line), and the chance of dispatching server 1 decreases (red line). This observation arises from our simplified synthetic setup, where we use a uniform dispatch policy and the dispatch probability to a region directly correlates with its size.

Furthermore, we also investigate the impact of dispatch policy by creating another 100 systems with different values of $\eta_{(1,2)}(1)$ and $\eta_{(1,2)}(2)$ (where $\eta_{(1,2)}(1)=1-\eta_{(1,2)}(2)$), while keeping the service regions the same, as shown in Figure 14 (a-c). For the sake of simplicity, we here omit the subscript $u$ in the notation of dispatch policy $\eta_{e}(i),~{}\forall i\in e$ as we assume it doesn’t depend on the state $u$ in this experiment. Figure 14 (d) and (e) show that the estimated $\rho_{i}$ and $\rho_{i,e}$ (indicated by solid lines) are able to closely approximate their true values suggested by the simulation. As we observe, when the dispatch probability for a call from the overlapping region to either server is 0.5, the workload distribution between the two servers is equal, i.e., $\eta_{1,2}(1)=\eta_{1,2}(2)=0.5$.

As discussed in Section 1, a main objective of the hyperlattice model is to extend the classic hypercube queue to capture the high-workload situations more accurately, thus allowing for a detailed analysis of over-saturated systems with queue lengths greater than one. In pursuit of this, we construct a synthetic system subjected to a traditional non-overlapping districting, defined by a $2\times 2$ grid encompassing four identical sub-regions served by their respective servers, as shown in Figure 15. For a meaningful comparison with the hypercube queuing model – where servers can attend to calls from any region if other servers are unavailable – we introduce an overlapping service zone comprised of all four sub-regions. This design ensures no primary service area is designated, granting every server the flexibility to address calls from any sub-region. Furthermore, we operate under the assumption that the call distributions are even across these sub-regions. This results in a call arrival rate of $\lambda/4$ = 1 for each. All servers in this model are consistent in their performance, exhibiting a uniform service rate denoted as $\mu$ = 1. Lastly, we adopt a standardized dispatch policy that is uniformly applied across all servers.

Figure 17 demonstrates the discrepancies in the estimation of individual workloads and dispatch fractions between the hypercube and hyperlattice models. In Figure 17 (a), it is evident that the hypercube model’s effectiveness diminishes substantially as the ratio escalates, whereas our method sustains commendable performance even at higher ratios. This advantage is primarily attributed to the hyperlattice’s capability to account for more detailed dynamics between queues, a feature not present in the hypercube model. Figure 17 (b) further reinforces the superior predictive accuracy of the hyperlattice model in estimating the fraction of dispatches, regardless of the $\lambda/\mu$ ratio. The success is attributed to the harmonized policy adherence between the simulation and that implemented by the hyperlattice model. Furthermore, this dispatch policy demonstrates robustness against fluctuations in the $\lambda/\mu$ ratio. Conversely, the hypercube model, which assumes a shared queue among all servers, fails to capture high-workload dynamics due to its simplified representation of the service process. Moreover, the hypercube model is constrained to a myopic policy (usually only dependent on travel distance) adopting a fixed preference for dispatch that becomes noticeably inadequate at lower $\lambda/\mu$ ratios, leading to a pronounced divergence from simulation results. However, an increase in the ratio sees a shift in the hypercube model towards a homogenous policy, akin to the one employed in our simulation, consequently reducing the gap in estimation accuracy. Such shift underscores the limitations in the hypercube model, demonstrating its restricted adaptability to different scenarios.