Computing transient probabilities in Markovian queues with balking conditional on a fixed number of joined customers

Abstract

We analyze the transient behavior of a Markovian queue with balking, conditional on exactly $K$ customers joining the system in a finite time interval $[0,T]$. A key quantity of interest is the cumulative number of balking customers, for which we consider the probability generating function (PGF) and derive an equation it satisfies. This approach circumvents the computational burden arising from dealing with high-dimensional Markovian system, which arises when attempting to directly compute the joint distribution of the cumulative number of balking customers together with the total number of joining customers and the number of customers in the system. To compute state transitions in $(t,T]$ efficiently, we examine the conditional state probability at time $T$ given the system state at time $T-u$, and show that it satisfies a linear differential equation in $u$. For piecewise-constant arrival rates, we develop a numerical procedure for evaluating the moment of the total number of balking customers, jointly with the cumulative number of joining customers and the number of customers in the system, under the condition that $K$ customers joined. We also present numerical examples that highlight counterintuitive behaviors arising from this conditioning, along with explanations of the underlying mechanisms.
Keyword balking, Markovian queue, fixed number of joined customers, transient analysis, computational procedure

1 Introduction

In many real-world queueing systems, such as waiting lines in cafeterias or at service counters in retail stores and offices, congestion can cause arriving customers to leave without receiving service. This general behavior can take two forms: when customers give up entering the system upon arrival, it is referred to as balking, while when they enter but leave before receiving service due to excessive waiting, it is referred to as reneging.

When customers abandon service opportunities through balking or reneging, the system loses potential revenue or throughput. Consequently, one of the primary performance metrics in such systems is the number of customers who arrive but do not receive service. Since these behaviors are typically influenced by factors that reflect the level of congestion, such as the number of customers already in the system or the expected waiting time, the congestion level at each time point also serves as an important measure.

Various queueing models incorporating balking and reneging have been studied extensively. Early works include [1, 2, 3, 4, 8, 9], which investigate customer abandonment in different settings. Reneging behavior under both random and first-come, first-served service disciplines is examined in [3, 4]. Balking behavior is modeled in [8], where stochastic thresholds are considered to represent each customer’s maximum acceptable number of customers in the system on arrival. This is extended in [9] to incorporate reneging based on individual waiting time tolerances. In [1, 2], probabilistic decisions to join the system based on the observed queue length are considered, where stationary distributions of several key system metrics are derived.

Subsequent studies have examined time-dependent arrival rates in queueing models. Fluid limits of queues with arrival rates depending on both time and congestion are analyzed in [20]. A discrete-time approximation to Poisson arrivals influenced by time and system state is adopted in [7]. For a comprehensive survey of queueing models with customer abandonment, see [17].

In this paper, we consider a Markovian queueing system with balking, under a nonstandard assumption that exactly $K$ customers join the system in a finite time interval $[0,T]$. This conditional framework departs from conventional models based on unconditioned arrival processes, and reflects realistic situations in which the total number of customers joining the system in a fixed period (such as one business day) is directly observable as a summary statistic. Although $K$ is not a controllable parameter of the system, it is often used as a reference quantity in performance evaluation and planning. Analyzing system behavior conditional on this observable quantity thus provides a complementary perspective to traditional time-averaged and transient analyses.

We consider a continuous-time Markovian queue with arrivals following a non-homogeneous Poisson process, where the system state at each time is given by the cumulative number of customers who have joined, the number of customers in the system, and a service phase, incorporating various service mechanisms. Among various performance metrics in this setting, a key quantity of interest is the cumulative number of balking customers. However, the computation of this quantity as a Markovian counting process requires analyzing a large underlying state-space Markov chain, leading to a substantial increase in computational complexity.

To obtain a computationally tractable characterization, we analyze the probability generating function (PGF) of the cumulative number of balking customers, and derive the corresponding factorial moments. We also consider the conditional state probability of the system at time $T$ given its state at time $T-u$, and show that it satisfies a linear differential equation in $u$. This result enables us to efficiently compute factorial moments of the cumulative number of balking customers conditional on the total number of joined customers.

Furthermore, for piecewise-constant arrival rates, we develop a numerical procedure for evaluating the factorial moments of the total number of balking customers under the condition that $K$ customers joined. We present numerical examples to illustrate how the system behavior under such conditioning can exhibit counterintuitive features, such as persistent time-dependent dynamics that arise even under time-homogeneous arrival rates, and we provide explanations for these dynamics.

We briefly review prior research on queueing models conditional on a fixed number of arrivals [5, 6, 10, 11, 12, 13, 15, 16]. In such models, the arrival times of $K$ customers are assumed to be independent and identically distributed (i.i.d.) in a finite interval $[0,T]$. This setup can be interpreted as conditioning a nonhomogeneous Poisson process to generate $K$ arrivals in $[0,T]$. Although balking is not considered in these models, several results have been established. In [16], a discrete-time model is analyzed, leading to equations for the PGFs of the number of customers in the system and the virtual waiting time. Continuous-time analogs are studied in [12, 13], where it is shown that the queue-length and workload processes converge weakly to Gaussian processes as $K$ increases. Approximation techniques have also been developed; for instance, fluid and diffusion limits are derived in [11], while reflected Brownian limits in critically loaded regimes are analyzed in [5, 6]. Furthermore, exact characterizations of the time-dependent workload and queue-length have been provided in [10] and [15].

These studies demonstrate the analytical richness of fixing the total number of arrivals, even without balking. In our setting, balking introduces additional complexity, particularly in computing the cumulative number of balked customers, as discussed later.

The remainder of this paper is structured as follows. In Section 2, we describe the model considered in this paper, explain the computational difficulties encountered with a naive analytical approach, and outline our approach. In Section 3, we present the analysis of the PGF of the cumulative number of balked customers conditional on the total number of arrivals, along with the corresponding results for the factorial moments. In Section 4, we provide a detailed construction of a computation algorithm for transient states probabilities in piecewise time-homogeneous systems, and in Section 5, we present numerical examples. Finally, we conclude the paper in Section 6.

2 Model and analytical method

2.1 Model

We consider a queueing system in which customers arrive in a finite time interval $[0,T]$ ($T\in(0,\infty)$). The arrival process of customers is assumed to follow a nonhomogeneous Poisson process with rate $\lambda(t)$, and for $t>T$, we set $\lambda(t)=0$. Additionally, arriving customers may balk based on the number of customers present just before their arrival. Specifically, customers who see $\ell$ customers ($\ell=0,1,\ldots$) present on arrival decides to join the system with probability $\beta_{\ell}$ ($0\leq\beta_{\ell}\leq 1$) and they leave immediately with probability $1-\beta_{\ell}$.

Let $A(t)$ ($t\geq 0$) denote the cumulative number of customers who have arrived at the system in the time interval $[0,t]$. Among these arriving customers, let $A^{\mathrm{join}}(t)$ denote the cumulative number of customers who have joined the system, and let $A^{\mathrm{balk}}(t)$ denote the cumulative number of customers who have balked. Let $D(t)$ denote the cumulative number of customers who have completed service and departed from the system in the time interval $[0,t]$. Also let $L(t)$ denote the number of customers in the system at time $t$. For simplicity, we assume that the system is initially empty, i.e., $L(0)=0$. By definition, we have the following relations:

$$A(t) =A^{\mathrm{join}}(t)+A^{\mathrm{balk}}(t),\quad t\geq 0,$$

$$A^{\mathrm{join}}(t) =L(t)+D(t),\quad t\geq 0.$$

We assume that the service mechanism of the system is Markovian, and we define $S(t)$ as the service phase at time $t$. The service phase $S(t)$ is assumed to contain sufficient information to describe the behavior of the system and it takes values in a finite set $\mathcal{S}$. Consequently, the triplet $(A^{\mathrm{join}}(t),L(t),S(t))_{t\geq 0}$ forms a continuous-time Markov chain defined on the state space $\Omega^{\mathrm{orig}}$, given by

$$\Omega^{\mathrm{orig}}\coloneqq\{(k,\ell,s);\ k=0,1,\ldots,\ \ell=0,1,\ldots,k,\ s\in\mathcal{S}\}.$$

The transition rate matrix of this continuous-time Markov chain is denoted by $\bm{Q}^{\mathrm{orig}}(t)$, where states are assumed to be ordered lexicographically. With this ordering, we regard $\bm{Q}^{\mathrm{orig}}(t)$ as the transition rate matrix of a bivariate Markov chain with level variable $A^{\mathrm{join}}(t)$ and phase variable $(L(t),S(t))$.

In this paper, we conduct a time-dependent analysis of the system conditional on the event that the cumulative number of joined customers by time $T$ equals a fixed constant $K$ ($K=1,2,\ldots$). Therefore, it suffices to consider only states where the cumulative number of joined customers does not exceed $K$. Accordingly, we define the restricted state space $\Omega$ as

$$\Omega\coloneqq\{(k,\ell,s);\ k=0,1,\ldots,K,\ \ell=0,1,\ldots,k,\ s\in\mathcal{S}\}.$$

Under this condition, the transition rate matrix $\bm{Q}^{\mathrm{orig}}(t)$ of the continuous-time Markov chain $(A^{\mathrm{join}}(t),L(t),S(t))_{t\geq 0}$ can be expressed using the transition rate matrix $\bm{Q}(t)$ for transitions between states belonging to $\Omega$:

$$\bm{Q}^{\mathrm{orig}}(t)=\left[\begin{array}[]{cc}\bm{Q}(t)&\bm{Q}_{1,2}^{\mathrm{orig}}(t)\\ \bm{O}&\bm{Q}_{2,2}^{\mathrm{orig}}(t)\end{array}\right].$$

Let $p_{k,\ell,s}(t)$ denote the state probability at time $t$:

$$p_{k,\ell,s}(t)\coloneqq\Pr[A^{\mathrm{join}}(t)=k,L(t)=\ell,S(t)=s],\quad(k,\ell,s)\in\Omega.$$

Let $|\Omega|$ denote the size of the state space. We define $\bm{p}(t)$ as a $1\times|\Omega|$ vector consisting of the elements $p_{k,\ell,s}(t)$, arranged in lexicographic order. By definition, we have for $t=0$,

$$p_{k,\ell,s}(0)=\left\{\begin{aligned} &\Pr[S(0)=s],\quad&&k=0,\ \ell=0,s\in\mathcal{S},\\ &0,\quad&&\text{otherwise}.\end{aligned}\right.$$ (1)

Under this condition, $\bm{p}(t)$ is given by the solution of the following differential equation:

$$\frac{\mathrm{d}\bm{p}(t)}{\mathrm{d}t}=\bm{p}(t)\bm{Q}(t).$$ (2)

In this paper, we consider the time-dependent behavior of the system, conditional on exactly $K$ customers having joined before time $T$. Specifically, for each $t\in[0,T]$, we examine the joint probability of the cumulative number of joined customers, the number of customers in the system, the service phase, and the cumulative number of balked customers. We define this time-dependent conditional joint probability as $\pi_{k,\ell,s,b}(t\mid K)$:

$$\pi_{k,\ell,s,b}(t\mid K)\coloneqq\Pr[A^{\mathrm{join}}(t)=k,L(t)=\ell,S(t)=s,A^{\mathrm{balk}}(t)=b\mid A^{\mathrm{join}}(T)=K].$$ (3)

2.2 Analytical method

A straightforward approach to compute $\pi_{k,\ell,s,b}(t\mid K)$ is to analyze a continuous-time Markov chain defined on an expanded state space

$$\breve{\Omega}\coloneqq\{(k,\ell,s,b);\ k=0,1,\ldots,K,\ \ell=0,1,\ldots,k,\ s\in\mathcal{S},\ b=0,1,\ldots\},$$ (4)

where $k$, $\ell$, $s$, and $b$ represent the number of joined customers, the number of customers in the system, the service phase, and the number of balked customers, respectively. Let $p_{k,\ell,s,b}(t)$ denote the state probability of the expanded Markov chain defined on $\breve{\Omega}$:

$$p_{k,\ell,s,b}(t)\coloneqq\Pr[A^{\mathrm{join}}(t)=k,L(t)=\ell,S(t)=s,A^{\mathrm{balk}}(t)=b],\quad(k,\ell,s,b)\in\breve{\Omega}.$$

We also define $\breve{\bm{p}}(t)$ as the corresponding probability vector. Similarly to (2), the state probability vector $\breve{\bm{p}}(t)$ satisfies

$$\frac{\mathrm{d}\breve{\bm{p}}(t)}{\mathrm{d}t}=\breve{\bm{p}}(t)\breve{\bm{Q}}(t),$$ (5)

where $\breve{\bm{Q}}(t)$ denotes the transition rate matrix of the Markov chain defined on $\breve{\Omega}$.

We can then rewrite (3) as

$$\pi_{k,\ell,s,b}(t\mid K)$$

$$=\frac{\Pr[A^{\mathrm{join}}(t)=k,L(t)=\ell,S(t)=s,A^{\mathrm{balk}}(t)=b,A^{\mathrm{join}}(T)=K]}{\Pr[A^{\mathrm{join}}(T)=K]}$$

$$=\frac{p_{k,\ell,s,b}(t)\Pr[A^{\mathrm{join}}(T)=K\mid A^{\mathrm{join}}(t)=k,L(t)=\ell,S(t)=s,A^{\mathrm{balk}}(t)=b]}{\Pr[A^{\mathrm{join}}(T)=K]}$$

$$=\frac{p_{k,\ell,s,b}(t)\Pr[A^{\mathrm{join}}(T)=K\mid A^{\mathrm{join}}(t)=k,L(t)=\ell,S(t)=s]}{\Pr[A^{\mathrm{join}}(T)=K]},$$ (6)

where the last equation holds because $A^{\mathrm{join}}(T)$ and $A^{\mathrm{balk}}(t)$ are conditionally independent given the system state $(A^{\mathrm{join}}(t),L(t),S(t))$. Since the transition probability $\Pr[A^{\mathrm{join}}(T)=K\mid A^{\mathrm{join}}(t)=k,L(t)=\ell,S(t)=s]$ is also determined by the transition rate matrix $\breve{\bm{Q}}(t)$, the expression (6) fully characterizes $\pi_{k,\ell,s,b}(t\mid K)$.

As an example, consider a piecewise time-homogeneous system with

$$\breve{\bm{Q}}(t)=\breve{\bm{Q}}_{n},\quad T_{n-1}<t\leq T_{n},$$

for $T_{0}=0\leq T_{1}\leq T_{2}\leq\cdots\leq T_{N-1}\leq T_{N}=T$. The solution of (5) in this case is given by the following equation:

$$\breve{\bm{p}}(t)=\breve{\bm{p}}(0)\left(\prod_{i=1}^{n-1}\exp\bigl(\breve{\bm{Q}}_{i}\cdot(T_{i}-T_{i-1})\bigr)\right)\exp\bigl(\breve{\bm{Q}}_{n}\cdot(t-T_{n-1})\bigr),\;\;T_{n-1}<t\leq T_{n}.$$ (7)

Therefore, we obtain $\pi_{k,\ell,s,b}(t\mid K)$ from (6) with $p_{k,\ell,s,b}(t)$ (the elements of $\breve{\bm{p}}(t)$) given by (7).

However, there are two main challenges in computing $\pi_{k,\ell,s,b}(t\mid K)$ using this approach:

1. Accounting for the cumulative number of balked customers.

   • Since the Markov chain is defined on a countably infinite state space $\breve{\Omega}$, it is difficult to directly compute the transient state probabilities $\breve{\bm{p}}(t)$ using (7).

     One possible workaround is to treat $(A^{\mathrm{balk}}(t))_{t\geq 0}$ as a Markovian counting process driven by the underlying process $(A^{\mathrm{join}}(t),L(t),S(t))\in\Omega$, for which an efficient computational method is known in the literature [14, 18]. However, the state space $\Omega$ still grows as $|\Omega|=O(K^{2}|\mathcal{S}|)$, which remains computationally expensive. For instance, even with the modest values $K=100$ and $|\mathcal{S}|=2$, $|\Omega|$ already exceeds $10,\!000$, highlighting the difficulty of scaling to larger $K$.

2. Computing the transition probability $\Pr[A^{\mathrm{join}}(T)=K\mid A^{\mathrm{join}}(t)=k,L(t)=\ell,S(t)=s]$ in (6).

   • Although this transition probability is straightforward to evaluate in principle, it imposes a heavy computational burden when $\pi_{k,\ell,s,b}(t\mid K)$ must be computed across many time points $t$, due to repeated evaluations of this quantity.

To address challenge (i), we adopt an approach based on PGFs. While direct computation of the time-dependent probability $p_{k,\ell,s,b}(t)$ is computationally intensive, the use of PGFs enables efficient calculation of the moments of the cumulative number of balked customers. More specifically, we define $p_{k,\ell,s}^{*}(z,t)$ ($|z|\leq 1$, $t\in[0,T]$) as the (unconditional) PGF of the cumulative number of of balked customers, and $p_{k,\ell,s}^{(i)}(t)$ ($i=0,1,\ldots$ $t\in[0,T]$) as the $i$-th (unconditional) factorial moment of the cumulative number of balked customers:

$$p_{k,\ell,s}^{*}(z,t) \coloneqq\sum_{b=0}^{\infty}p_{k,\ell,s,b}(t)z^{b},\quad(k,\ell,s)\in\Omega,\ 0\leq t\leq T,$$

$$p_{k,\ell,s}^{(0)}(t) \coloneqq p_{k,\ell,s}^{*}(1,t)=p_{k,\ell,s}(t),$$ (8)

$$p_{k,\ell,s}^{(i)}(t) \coloneqq\frac{\partial^{i}p_{k,\ell,s}^{*}(z,t)}{\partial z^{i}}\bigg|_{z=1}$$

$$=\mathrm{E}\left[\prod_{h=0}^{i-1}(A^{\mathrm{balk}}(t)-h)\cdot\mathbbm{1}_{\{A^{\mathrm{join}}(t)=k,L(t)=\ell,S(t)=s\}}\right],\;\;i=1,2,\ldots,$$ (9)

where $\mathbbm{1}_{\{\cdot\}}$ denotes the indicator function. Similarly, we define the PGF and factorial moment of the number of balked customers conditional on exactly $K$ customers having joined before time $T$:

$$\pi_{k,\ell,s}^{*}(z,t\mid K) \coloneqq\sum_{b=0}^{\infty}\pi_{k,\ell,s,b}(t\mid K)z^{b},\quad(k,\ell,s)\in\Omega,\ 0\leq t\leq T,$$ (10)

$$\pi_{k,\ell,s}^{(0)}(t\mid K) \coloneqq\pi_{k,\ell,s}^{*}(1,t\mid K)$$

$$\>=\Pr[A^{\mathrm{join}}(t)=k,L(t)=\ell,S(t)=s\mid A^{\mathrm{join}}(t)=K],$$

$$\pi_{k,\ell,s}^{(i)}(t\mid K) \coloneqq\frac{\partial^{i}\pi_{k,\ell,s}^{*}(z,t\mid K)}{\partial z^{i}}\bigg|_{z=1}$$

$$=\mathrm{E}\left[\prod_{h=0}^{i-1}(A^{\mathrm{balk}}(t)-h)\cdot\mathbbm{1}_{\{A^{\mathrm{join}}(t)=k,L(t)=\ell,S(t)=s\}}\mid A^{\mathrm{join}}(T)=K\right],$$

$$\hskip 230.00035pti=1,2,\ldots.$$

We show that the unconditional PGF $p_{k,\ell,s}^{*}(z,t)$ (for each $|z|\leq 1$) satisfies a first-order differential equation in time $t$, and characterize $\pi_{k,\ell,s}^{*}(z,t\mid K)$ accordingly. Furthermore, using these results, we characterize both the unconditional factorial moments $p_{k,\ell,s}^{(i)}(t)$ and the conditional factorial moments $\pi_{k,\ell,s}^{(i)}(t\mid K)$. Based on these results, we develop a practical computational procedure applicable to piecewise time-homogeneous systems.

To address challenge (ii), we make use of the fact that the cumulative number of joined customers is fixed at time $T$, and consider the probability that $A^{\mathrm{join}}(T)=K$ given the system state at an earlier time $T-u$ for $0\leq u\leq T$. This backward-in-time formulation allows us to efficiently compute the transition probabilities $\Pr[A^{\mathrm{join}}(T)=K\mid A^{\mathrm{join}}(T-u)=k,L(T-u)=\ell,S(T-u)=s]$ for multiple values of $u$, as they satisfy a linear differential equation in $u$.

In the following sections, we analyze the conditional PGF $\pi_{k,\ell,s}^{*}(z,t\mid K)$ and develop an efficient method for computing it.

3 Analysis of the time-dependent PGF conditional on $K$ joined customers

We first consider the unconditional PGF. Let $\bm{p}^{*}(z,t)$ denote the row vector consisting of the elements $p_{k,\ell,s}^{*}(z,t)$, arranged in lexicographic order.

Lemma 1.

For each $|z|\leq 1$, the unconditional PGF $\bm{p}^{*}(z,t)$ is determined by the following linear differential equation:

$$\bm{p}^{*}(z,0) =\bm{p}(0),$$ (11)

$$\frac{\partial\bm{p}^{*}(z,t)}{\partial t} =\bm{p}^{*}(z,t)\bigl[\bm{Q}(t)-\lambda(t)(1-z)\bm{B}\bigr],$$ (12)

where $\bm{p}(0)$ is given by (1), and $\bm{B}$ denotes an $|\Omega|\times|\Omega|$ diagonal matrix representing state transitions due to customer balking:

$$[\bm{B}]_{(k,\ell,s),(k^{\prime},\ell^{\prime},s^{\prime})}=\left\{\begin{aligned} &1-\beta_{\ell},&&(k^{\prime},\ell^{\prime},s^{\prime})=(k,\ell,s),\\ &0,&&\text{otherwise}.\end{aligned}\right.$$

Proof.

Note that (11) is obvious from $A^{\mathrm{balk}}(0)=0$. By definition, as $\varDelta t\to 0+$,

$$\bm{p}^{*}(z,t+\varDelta t)=\bm{p}^{*}(z,t)\{\bm{I}+(\bm{Q}(t)-\lambda(t)\bm{B})\varDelta t\}+\bm{p}^{*}(z,t)\cdot z\lambda(t)\bm{B}\varDelta t+\bm{o}(\varDelta t),$$

which implies (12). ∎

We then characterize the conditional PGF $\pi_{k,\ell,s}^{*}(z,t\mid K)$ (defined as (10)) in terms of the unconditional PGF. Recall that the conditional and unconditional state probabilities $\pi_{k,\ell,s,b}(t\mid K)$ and $p_{k,\ell,s,b}(t)$ are related as (6). As explained in the previous section, the key quantity in developing an efficient computational method is the conditional probability $\Pr[A^{\mathrm{join}}(T)=K\mid A^{\mathrm{join}}(t)=k,L(t)=\ell,S(t)=s]$.

To facilitate efficient analysis of this quantity, we introduce a conditional probability $q_{k,\ell,s}(u,K)$, parametrized backward in time:

$$q_{k,\ell,s}(u,K)\coloneqq\Pr[A^{\mathrm{join}}(T)=K\mid A^{\mathrm{join}}(T-u)=k,L(T-u)=\ell,S(T-u)=s],\$$

$$(k,\ell,s)\in\Omega,\ 0\leq u\leq T.$$ (13)

By definition, $q_{k,\ell,s}(u,K)$ represents the conditional probability that exactly $K$ customers have joined before time $T$, given the system state $(k,\ell,s)$ at time $T-u$. Let $\bm{q}(u,K)$ denote a $|\Omega|\times 1$ vector consisting of the elements $q_{k,\ell,s}(u,K)$, arranged in lexicographic order.

Lemma 2.

$\bm{q}(u,K)$ is determined by the following linear equation

$$q_{k,\ell,s}(0,K) =\left\{\begin{aligned} &1,&&k=K,\ \ell=0,1,\ldots,K,\ s\in\mathcal{S},\\ &0,&&\text{otherwise},\end{aligned}\right.$$ (14)

$$\frac{\partial\bm{q}(u,K)}{\partial u} =\bm{Q}(T-u)\bm{q}(u,K).$$ (15)

Proof.

As the boundary condition (14) follows immediately from the definition of $\bm{q}(u,K)$, we prove (15) below. From the Markov property of the system state $(A^{\mathrm{join}}(t),L(t),S(t))_{t\geq 0}$, we have for $\varDelta t>0$,

$$q_{k,\ell,s}(u+\varDelta t)$$

$$=\Pr[A^{\mathrm{join}}(T)=K$$

$$\qquad\qquad\mid A^{\mathrm{join}}(T-u-\varDelta t)=k,L(T-u-\varDelta t)=\ell,S(T-u-\varDelta t)=s]$$

$$=\sum_{(k^{\prime},\ell^{\prime},s^{\prime})\in\Omega}\Bigl(\Pr[A^{\mathrm{join}}(T-u)=k^{\prime},L(T-u)=\ell^{\prime},S(T-u)=s^{\prime}$$

$$\hskip 60.00009pt\mid A^{\mathrm{join}}(T-u-\varDelta t)=k,L(T-u-\varDelta t)=\ell,S(T-u-\varDelta t)=s]$$

$$\hskip 40.00006pt\cdot\Pr[A^{\mathrm{join}}(T)=K\mid A^{\mathrm{join}}(T-u)=k^{\prime},L(T-u)=\ell^{\prime},S(T-u)=s^{\prime}]\Bigr).$$

Therefore, using the elements $[\bm{Q}(t)]_{(k,\ell,s),(k^{\prime},\ell^{\prime},s^{\prime})}$ of the transition matrix $\bm{Q}(t)$, we have

$$q_{k,\ell,s}(u+\varDelta t) =(1+[\bm{Q}(T-u-\varDelta t)]_{(k,\ell,s),(k,\ell,s)}\varDelta t)q_{k,\ell,s}(u,K)$$

$$\qquad+\sum_{(k^{\prime},\ell^{\prime},s^{\prime})\in\Omega\backslash(k,\ell,s)}[\bm{Q}(T-u-\varDelta t)]_{(k,\ell,s),(k^{\prime},\ell^{\prime},s^{\prime})}\varDelta tq_{k^{\prime},\ell^{\prime},s^{\prime}}(u,K)$$

$$\qquad+o(\varDelta t),$$

as $\varDelta t\to 0+$. This equation is further rewritten in matrix form as

$$\bm{q}(u+\varDelta t)=(\bm{I}+\bm{Q}(T-u-\varDelta t)\varDelta t)\bm{q}(u,K)+\bm{o}(\varDelta t),$$

which implies (15). ∎

The conditional PGF $\pi_{k,\ell,s}^{*}(z,t\mid K)$ is thus obtained in terms of $p_{k,\ell,s}^{*}(t)$ and $q_{k,\ell,s}(u,K)$:

Theorem 3.

The PGF $\pi_{k,\ell,s}^{*}(z,t\mid K)$ conditional on exactly $K$ customers joined in $[0,T]$ is given by

$$\pi_{k,\ell,s}^{*}(z,t\mid K)=\frac{p_{k,\ell,s}^{*}(z,t)q_{k,\ell,s}({T-t},K)}{\Pr[A^{\mathrm{join}}(T)=K]},\quad(k,\ell,s)\in\Omega,\ 0\leq t\leq T,$$ (16)

where $p_{k,\ell,s}^{*}({z},t)$ and $q_{k,\ell,s}(u,K)$ are given by Lemma 1 and Lemma 2.

Proof.

Theorem 3 follows immediately from (6). ∎

Theorem 3 provides the complete characterization of the state probabilities of the system, conditional on the event that exactly $K$ customers have joined in $[0,T]$. As $\pi_{k,\ell,s}^{*}(z,t\mid K)$ represents the number of balked customers in the form of PGF, we can derive its moments accordingly.

Letting $\bm{p}^{(i)}(t)$ denote a $1\times|\Omega|$ vector consisting of the elements $p_{k,\ell,s}^{(i)}(t)$ (cf. (8) and (9)) arranged in lexicographic order, we obtain the following results from Lemma 1 and Theorem 3:

Lemma 4.

$\bm{p}^{(i)}(t)$ ($i=0,1,\ldots$) is determined by the following differential equation:

$$\bm{p}^{(i)}(0) =\bm{0},\quad i=1,2,\ldots,$$ (17)

$$\frac{\mathrm{d}\bm{p}^{(i)}(t)}{\mathrm{d}t} =\bm{p}^{(i)}(t)\bm{Q}(t)+i\bm{p}^{(i-1)}(t)\lambda(t)\bm{B},\quad 0\leq t\leq T,\ i=1,2,\ldots.$$ (18)

Proof.

Differentiating both sides of (11) and (12) with respect to $z$ up to the $i$-th order and substituting $z=1$, we obtain (17) and (18). ∎

Theorem 5.

$\pi_{k,\ell,s}^{(i)}(t\mid K)$ are given by

$$\pi_{k,\ell,s}^{(i)}(t\mid K)=\frac{p_{k,\ell,s}^{(i)}(t)q_{k,\ell,s}(T-t,K)}{\Pr[A^{\mathrm{join}}(T)=K]},\quad 0\leq t\leq T,\ i=0,1,\ldots,$$ (19)

where $q_{k,\ell,s}(u,K)$ and $p_{k,\ell,s}^{(i)}(t)$ are given by Lemma 2 and Lemma 4.

Proof.

(19) follows immediately from (16) and the definitions of $\pi_{k,\ell,s}^{(i)}(t\mid K)$ and $p_{k,\ell,s}^{(i)}(t)$. ∎

The results above show that for given system parameters $\lambda(t)$, $\beta_{\ell}$, and $\bm{Q}(t)$, we can compute the conditional PGF $\pi_{k,\ell,s}^{*}(z,t\mid K)$ and the moments $\pi_{k,\ell,s}^{(i)}(t\mid K)$ by solving the linear differential equations in Lemma 1, Lemma 2, and Lemma 4. In the following sections, we develop a detailed computational algorithm focusing on the case that $\lambda(t)$ and $\bm{Q}(t)$ are piecewise constant in time $t$ and we present numerical examples, which demonstrate the feasibility of computation based on these general results.

4 Special case: Piecewise time-homogeneous systems

In this section, we develop a detailed computational algorithm for the case with piecewise constant $\lambda(t)$ and $\bm{Q}(t)$. More specifically, the time interval $(0,T]$ is divided into $N$ ($N\in\mathbb{N}$) subintervals $(T_{0},T_{1}],(T_{1},T_{2}],\ldots,(T_{N-1},T_{N}]$ with $T_{0}=0$ and $T_{N}=T$, and the arrival rate $\lambda(t)$ and transition rate matrix $\bm{Q}(t)$ are assumed constant during each interval:

$$\lambda(t)=\lambda_{n},\quad\bm{Q}(t)=\bm{Q}_{n},\qquad T_{n-1}<t\leq T_{n},\ n=1,2,\ldots,N.$$ (20)

In this case, the unconditional PGF $\bm{p}^{*}(z,t)$ characterized as Lemma 1 is given recursively by

$$\bm{p}^{*}(z,t)=\bm{p}^{*}(z,T_{n-1})\exp((\bm{Q}_{n}-\lambda_{n}(1-z)\bm{B})(t-T_{n-1})),\quad$$

$$T_{n-1}<t\leq T_{n},\ n=1,2,\ldots,N,$$ (21)

with $\bm{p}^{*}(z,0)$ given by (11). Similarly, the conditional probability $\bm{q}(u,K)$ characterized as Lemma 2 is given recursively by

$$\bm{q}(u,K)=\exp(\bm{Q}_{N-n+1}\cdot(u-T+T_{N-n+1}))\bm{q}(T-T_{N-n+1},K),\qquad$$

$$T-T_{N-n+1}<u\leq T-T_{N-n},\ n=1,2,\ldots,N,$$ (22)

with $\bm{q}(0,K)$ given by (14). Note here that the equality at right endpoints $t=T_{n}$ in (21) and $u=T-T_{N-n}$ in (22) follows from the continuity of $\bm{p}^{*}(z,t)$ and $\bm{q}(u,K)$. Therefore, we obtain the time-dependent PGF $\pi_{k,\ell,s}^{*}(z,t\mid K)$ from Theorem 3.

For the factorial moments $\pi_{k,\ell,s}^{(i)}(t\mid K)$, the following representation of unconditional factorial moments $\bm{p}^{(i)}(t)$ leads to efficient computation:

$$\widetilde{\bm{p}}^{(J)}(t)\coloneqq\left[\begin{array}[]{cccc}\displaystyle\frac{\bm{p}^{(0)}(t)}{0!}&\displaystyle\frac{\bm{p}^{(1)}(t)}{1!}&\cdots&\displaystyle\frac{\bm{p}^{(J)}(t)}{J!}\end{array}\right],\quad J=0,1,\ldots.$$ (23)

By definition, $\widetilde{\bm{p}}^{(J)}(t)$ denotes a vector of time-dependent joint probability $\bm{p}(t)=\bm{p}^{(0)}(t)$ and the factorial moments $\bm{p}^{(i)}(t)$ ($i=1,2,\ldots,J$), weighted according to their order. Under the assumption (20), $\widetilde{\bm{p}}^{(J)}(t)$ is obtained as follows:

Theorem 6.

For a fixed $J\in\{1,2,\ldots\}$, $\widetilde{\bm{p}}^{(J)}(t)$ is given by

$$\widetilde{\bm{p}}^{(J)}(0) =\underbrace{\left[\begin{array}[]{cccc}\bm{p}(0)&\bm{0}&\cdots&\bm{0}\end{array}\right]}_{\displaystyle J+1},$$

$$\widetilde{\bm{p}}^{(J)}(t) =\widetilde{\bm{p}}^{(J)}(T_{n-1})\exp(\widetilde{\bm{Q}}_{n}^{(J)}(t-T_{n-1})),\quad T_{n-1}<t\leq T_{n},\ n=1,2,\ldots,N,$$ (25)

where $\widetilde{\bm{Q}}_{n}^{(J)}$ is defined as

$$\widetilde{\bm{Q}}_{n}^{(J)} =\underbrace{\left[\begin{array}[]{ccccc}\bm{Q}_{n}&\lambda_{n}\bm{B}&\bm{O}&\cdots&\bm{O}\\ \bm{O}&\bm{Q}_{n}&\lambda_{n}\bm{B}&\cdots&\bm{O}\\ \bm{O}&\bm{O}&\bm{Q}_{n}&\cdots&\bm{O}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ \bm{O}&\bm{O}&\bm{O}&\cdots&\bm{Q}_{n}\end{array}\right]}_{\displaystyle(J+1)\times(J+1)}.$$

Proof.

As (6) follows immediately from (17), we prove (25) below. Using (20), we rewrite (2) and (18) as

$$\frac{\mathrm{d}\bm{p}^{(0)}(t)}{\mathrm{d}t} =\bm{p}^{(0)}(t)\bm{Q}_{n},\quad i=0,\ T_{n-1}<t\leq T_{n},\ n=1,2,\ldots,N,$$

$$\frac{\mathrm{d}\bm{p}^{(i)}(t)}{\mathrm{d}t} =\bm{p}^{(i)}(t)\bm{Q}_{n}+i\bm{p}^{(i-1)}(t)\lambda_{n}\bm{B},$$

$$\hskip 68.00012pti=1,2,\ldots,\ T_{n-1}<t\leq T_{n},\ n=1,2,\ldots,N,$$

which implies

$$\frac{\mathrm{d}}{\mathrm{d}t}\left[\frac{\bm{p}^{(0)}(t)}{0!}\right] =\frac{\bm{p}^{(0)}(t)}{0!}\cdot\bm{Q}_{n},\quad i=0,\ T_{n-1}<t\leq T_{n},\ n=1,2,\ldots,N,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\left[\frac{\bm{p}^{(i)}(t)}{i!}\right] =\frac{\bm{p}^{(i)}(t)}{i!}\cdot\bm{Q}_{n}+\frac{\bm{p}^{(i-1)}(t)}{(i-1)!}\cdot\lambda_{n}\bm{B},$$

$$\hskip 78.00014pti=1,2,\ldots,\ T_{n-1}<t\leq T_{n},\ n=1,2,\ldots,N.$$

It then follows from (23) and (6) that

$$\frac{\mathrm{d}\widetilde{\bm{p}}^{(J)}(t)}{\mathrm{d}t}=\widetilde{\bm{p}}^{(J)}(t)\widetilde{\bm{Q}}_{n}^{(J)},\quad T_{n-1}<t\leq T_{n},\ n=1,2,\ldots,N.$$

Therefore, we obtain (25) from this differential equation. ∎

Theorem 6 shows that the unconditional factorial moments $\bm{p}^{(i)}(t)$ ($i=0,1,\ldots,J$) up to the $J$-th order can be obtained simultaneously by computing $\widetilde{\bm{p}}^{(J)}(t)$. More specifically, we have

$$\bm{p}(t) =\widetilde{\bm{p}}^{(J)}(t)\widetilde{\bm{I}}_{0}^{(J)},$$ (31)

$$\bm{p}^{(i)}(t) =\widetilde{\bm{p}}^{(J)}(t)\widetilde{\bm{I}}_{i}^{(J)}i!, \quad i=1,2,\ldots,J,$$ (32)

where $\widetilde{\bm{I}}_{i}^{(J)}$ is given by

$$\widetilde{\bm{I}}_{i}^{(J)}\coloneqq\Bigl[\underbrace{\begin{array}[]{ccc}\bm{O}&\cdots&\bm{O}\end{array}}_{\displaystyle i}\begin{array}[]{c}\bm{I}\end{array}\underbrace{\begin{array}[]{ccc}\bm{O}&\cdots&\bm{O}\end{array}}_{\displaystyle J-i}\Bigr],\quad i=0,1,\ldots,J,\ J\in\mathbbm{N},$$

i.e., it has an identity matrix $\bm{I}$ at the $i$-th block, and zero matrices $\bm{O}$ elsewhere. Therefore, utilizing this representation, we can compute the conditional factorial moments $\pi_{k,\ell,s}^{(i)}(t\mid K)$ efficiently from Theorem 5.

To develop a stable computational algorithm, it is essential to use the uniformization technique [19, pp. 154–156] for computing the matrix exponentials. More specifically, we rewrite (22) as

$$\bm{q}(u,K)=\sum_{m=0}^{\infty}\mathrm{Poi}(\theta_{N-n+1}(u-T+T_{N-n+1}),m)(\bm{P}_{n})^{m}\bm{q}(T-T_{N-n+1},K),\qquad$$

$$T-T_{N-n+1}\leq u\leq T-T_{N-n},\ n=1,2,\ldots,N,$$ (33)

where $\theta_{n}$ and $\bm{P}_{n}$ are defined as

$$\theta_{n} \coloneqq\max_{i}|[\bm{Q}_{n}]_{i,i}|+\lambda_{n},$$ (34)

$$\bm{P}_{n} \coloneqq\bm{I}+\theta_{n}^{-1}\bm{Q}_{n},$$

and $\mathrm{Poi}(a,m)$ denotes the probability mass function of the Poisson distribution with mean $a$:

$$\mathrm{Poi}(a,m)\coloneqq\frac{a^{m}}{m!}e^{-a}.$$

Note that the right-hand side of (33) involves only additions and multiplications of non-negative numbers, which ensures numerical stability by avoiding the loss of significant digits. Although the choice of $\theta_{n}$ is arbitrary as long as it satisfies $\theta_{n}\geq\max_{i}|[\bm{Q}_{n}]_{i,i}|$, we adopt the specific form used here for convenience in computing $\widetilde{\bm{p}}^{(J)}(t)$, as discussed below.

For the computation of $\widetilde{\bm{p}}^{(J)}(t)$ using (25), it is important to note that $\widetilde{\bm{Q}}_{n}^{(J)}$ (defined as (6)) is not necessarily a proper transition rate matrix as its row sums may exceed zero. To improve numerical stability, we factor out a scalar exponential term:

$$\exp(\widetilde{\bm{Q}}_{n}^{(J)}\cdot(t-T_{n-1}))=e^{\lambda_{n}(t-T_{n-1})}\exp((\widetilde{\bm{Q}}_{n}^{(J)}-\lambda_{n}\bm{I})(t-T_{n-1})),$$

where $\widetilde{\bm{Q}}_{n}^{(J)}-\lambda_{n}\bm{I}$ represents a proper transition rate matrix. Using this relation, we rewrite (25) as

$$\widetilde{\bm{p}}^{(J)}(t)=e^{\lambda_{n}(t-T_{n-1})}\sum_{m=0}^{\infty}\mathrm{Poi}(\theta_{n}(t-T_{n-1}),m)\widetilde{\bm{p}}^{(J)}(T_{n-1})(\widetilde{\bm{P}}_{n}^{(J)}-\theta_{n}^{-1}\lambda_{n}\bm{I})^{m},\$$

$$T_{n-1}<t\leq T_{n},\ n=1,2,\ldots,N,$$ (35)

where

$$\widetilde{\bm{P}}_{n}^{(J)} \coloneqq\bm{I}+\theta_{n}^{-1}\widetilde{\bm{Q}}_{n}^{{(J)}}=\underbrace{\left[\begin{array}[]{ccccc}\bm{P}_{n}&\theta_{n}^{-1}\lambda_{n}\bm{B}&\bm{O}&\cdots&\bm{O}\\ \bm{O}&\bm{P}_{n}&\theta_{n}^{-1}\lambda_{n}\bm{B}&\cdots&\bm{O}\\ \bm{O}&\bm{O}&\bm{P}_{n}&\cdots&\bm{O}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ \bm{O}&\bm{O}&\bm{O}&\cdots&\bm{P}_{n}\end{array}\right]}_{\displaystyle(J+1)\times(J+1)}.$$

We can verify that the right-hand side of (35) involves only additions and multiplications of non-negative numbers.

Furthermore, to compute the term $\widetilde{\bm{p}}^{(J)}(T_{n-1})(\widetilde{\bm{P}}_{n}^{(J)}-\theta_{n}^{-1}\lambda_{n}\bm{I})^{m}$ in (35) efficiently, we leverage the block structure and sparsity of $\widetilde{\bm{P}}_{n}^{(J)}$. To this end, we define $\bm{x}_{n}(m)$ as

$$\bm{x}_{n}(m)\coloneqq\widetilde{\bm{p}}^{(J)}(T_{n-1})(\widetilde{\bm{P}}_{n}^{(J)}-\theta_{n}^{-1}\lambda_{n}\bm{I})^{m}=[\bm{x}_{n}^{(0)}(m),\bm{x}_{n}^{(1)}(m),\ldots,\bm{x}_{n}^{(J)}(m)],\quad$$

$$m=0,1,\ldots,$$

which satisfies the recurrence

$$\bm{x}_{n}(0)=\widetilde{\bm{p}}^{(J)}(T_{n-1}),\quad\bm{x}_{n}(m)=\bm{x}_{n}(m-1)(\widetilde{\bm{P}}_{n}^{(J)}-\theta_{n}^{-1}\lambda_{n}\bm{I}),\quad m=1,2,\ldots.$$

Each block $\bm{x}_{n}^{(i)}(m)$ of $\bm{x}_{n}(m)$ can then be computed recursively as

$$\bm{x}_{n}^{(0)}(m+1) =\bm{x}_{n}^{(0)}(m)(\bm{P}_{n}-\theta_{n}^{-1}\lambda_{n}\bm{I}), \quad i=0,$$ (41)

$$\bm{x}_{n}^{(i)}(m+1) =\bm{x}_{n}^{(i-1)}(m)\theta^{-1}\lambda_{n}\bm{B}+\bm{x}_{n}^{(i)}(m)(\bm{P}_{n}-\theta_{n}^{-1}\lambda_{n}\bm{I}), \quad i=1,2,\ldots.$$ (42)

Remark 7.

For the uniformization technique to work properly, both $\bm{P}_{n}$ in (33) and $\widetilde{\bm{P}}_{n}^{(J)}-\theta_{n}^{-1}\lambda_{n}\bm{I}$ in (35) must be substochastic matrices. The former condition is satisfied when $\theta_{n}\geq\max_{i}|[\bm{Q}_{n}]_{i,i}|$, while the latter requires $\theta_{n}\geq\max_{i}|[\bm{Q}_{n}]_{i,i}|+\lambda_{n}$. Therefore, we adopt the specific form of $\theta_{n}$ given in (34) to ensure that both conditions hold.

Figure 1 summarizes the computational procedure for evaluating the factorial moments $\pi_{k,\ell,s}^{(i)}(t\mid K)$. In this procedure, the truncation point $m=M$ for the infinite sums in (33) and (35) is treated as an input parameter.

    Input: $T$, $K$, $N$, $T_{n}$ ($n=1,2,\ldots,N$), $\lambda_{n}$ ($n=1,2,\ldots,N$), $\bm{B}$, $\bm{Q}_{n}$ ($n=1,2,\ldots,N$), $\bm{p}(0)$, $t$, $J$, $N^{*}(t)$, and $M$. Output: $\pi_{k,\ell,s}^{(i)}(t\mid K)$ ($(k,\ell,s)\in\Omega$, $i=0,1,\ldots,J$). Step 1: Computation of $\widetilde{\bm{p}}^{(J)}(T_{N^{*}(t)-1})$ and $\widetilde{\bm{p}}^{(J)}(T)$.  Set $\widetilde{\bm{p}}^{(J)}(0)$ as in (6).  for $n=1$ to $N$ do   Compute $\widetilde{\bm{p}}^{(J)}(T_{n})$ using (35) and the truncation point $m=M$.  endfor Step 2: Computation of $\bm{q}(T-T_{N^{*}(t)},K)$.  Set $\bm{q}(0,K)$ as in (14).  if $N^{*}(t)<N$ then   for $n=1$ to $N-N^{*}(t)$ do    Compute $\bm{q}(T-T_{N-n},K)$ using (33) and the truncation point $m=M$.   endfor  endif Step 3: Computation of output $\pi_{k,\ell,s}^{(i)}(t\mid K)$.  Compute $\widetilde{\bm{p}}^{(J)}(t)$ using (35) and $\bm{q}(T-t,K)$ using (33), with the truncation point $m=M$.  for $i=0$ to $J$ do   Compute $\bm{p}^{(i)}(t)$ by (31) or (32).   for $(k,\ell,s)\in\Omega$ do    Compute $\pi_{k,\ell,s}^{(i)}(t\mid K)$ by (19).   endfor  endfor    

Figure 1: Computational procedure for $\pi_{k,\ell,s}^{(i)}(t\mid K)$ ($t\in(0,T]$). $N^{*}(t)$ denotes a non-negative integer that satisfies $T_{N^{*}(t)-1}<t\leq T_{N^{*}(t)}$.

5 Numerical examples

In this section, we present numerical examples for an M/PH/1 queue and investigate how condition that the total number of joined customers is given affects the time-dependent behavior of the number of customers in the system. This setting corresponds to the case considered in Section 4 with $N=1$, and we fix the arrival rate as $\lambda_{1}=\lambda$ throughout this section. Service times are assumed to be i.i.d. according to a phase-type distribution. More specifically, we use a subclass of phase-type distributions [19, pp. 358–359] characterized by the mean and coefficient of variation $C_{V}$: (i) a mixture of Erlang distributions with $k$ and $k+1$ stages ($0<C_{V}<1$), (ii) an exponential distribution ($C_{V}=1$), and (iii) a balanced hyper-exponential distribution ($C_{V}>1$). Throughout the experiments, the mean service time is fixed to one. Unless otherwise mentioned, we set the observation horizon to $T=100$ and the total number of joined customers to $K=100$. We provide a detailed computational procedure for this model in Appendix A.

We begin by examining the time-dependent distribution of the number of customers in the system, conditional on the total number of joined customers. We set the entry probability $\beta_{\ell}$ as

$$\beta_{\ell}=\max\left(1-\frac{\ell}{50},0\right),\quad\ell=0,1,\ldots,100.$$ (43)

Figure 2 shows heatmaps of the distribution of the number of customers in the system for $C_{V}=1.0$ and different values of $\lambda$, conditional on exactly $K=100$ joined customers. Each panel also includes $10$ sample paths obtained from simulation. These sample paths tend to pass through regions of higher density in the heatmaps, highlighting the consistency between the simulated trajectories and the computed distributions.

These heatmaps and sample paths reveal behaviors that may initially seem counterintuitive. For example, in Figure 2 (d), the number of customers follows a unimodal trajectory (first increasing and then decreasing) despite the constant arrival rate $\lambda(t)=\lambda=2.0$. This illustrates how condition that the total number of joined customers ($K=100$) is given can induce dynamics that deviate fundamentally from those in conventional queueing models, where the system converges to the steady state over time.

(a) $\lambda=0.5$

(b) $\lambda=1.1$

(c) $\lambda=1.5$

(d) $\lambda=2.0$

Figure 2: The distribution of the number of customers in the system plotted with simulated sample paths ($C_{V}=1.0$ and $\beta_{\ell}$ is given by (43)).

To investigate the mechanism behind this behavior, Figure 3 presents the time evolution of several key system metrics, all conditional on exactly $K$ customers having joined by time $T$:

• the expected number of customers in the system $\mathrm{E}[L(t)\mid A^{\mathrm{join}}(T)=K]$,

• the increment in the expected cumulative number of joined customers $\mathrm{E}[\varDelta A^{\mathrm{join}}(t)\mid A^{\mathrm{join}}(T)=K]$,

• the increment in the expected cumulative number of balked customers $\mathrm{E}[\varDelta A^{\mathrm{balk}}(t)\mid A^{\mathrm{join}}(T)=K]$,

• the increment in the expected cumulative number of arrivals $\mathrm{E}[\varDelta A(t)\mid A^{\mathrm{join}}(T)=K]$,

• the increment in the expected cumulative number of departures $\mathrm{E}[\varDelta D(t)\mid A^{\mathrm{join}}(T)=K]$, and

• the system utilization $\Pr[L(t)\geq 1\mid A^{\mathrm{join}}(T)=K]$,

where for any time-dependent quantity $X(t)$, we define the increment $\varDelta X(t)\coloneqq X(t+1)-X(t)$. From Figure 3 (a), we observe that for $\lambda=1.3$, $1.5$, and $2.0$, the number of customers in the system first increases and then decreases over time. In contrast, for $\lambda=0.5$ and $0.8$, the number of customers remains nearly flat for most of the observation period, with slight increases near the beginning and end.

By definition, the change in the number of customers in the system satisfies $\varDelta L(t)=\varDelta A^{\mathrm{join}}(t)-\varDelta D(t)$, so the transient behavior of $L(t)$ can be explained in terms of the numbers of joined and departed customers. Figure 3 (b) shows that the mean joining rate $\mathrm{E}[\varDelta A^{\mathrm{join}}(t)\mid A^{\mathrm{join}}(T)=K]$ exhibits a bias under the conditioning $A^{\mathrm{join}}(T)=K$. In particular, larger values of the arrival rate $\lambda$ lead to more joining events in the early stage. This bias arises because, when $\lambda\gg K/T=1$, a larger number of balking events is necessary to limit the total number of joined customers to $K=100$. To satisfy this constraint, the system tends to admit more customers early on, causing $L(t)$ to increase rapidly and thereby creating more opportunities for balking.

Figure 3 (a), (c), and (d) further show that the bias in the conditional joining rate $\mathrm{E}[\varDelta A^{\mathrm{join}}(t)\mid A^{\mathrm{join}}(T)=K]$ originates mainly from a bias in the conditional arrival rate $\mathrm{E}[\varDelta A(t)\mid A^{\mathrm{join}}(T)=K]$, rather than from the balking behavior itself. In fact, the conditional balking rate $\mathrm{E}[\varDelta A^{\mathrm{balk}}(t)\mid A^{\mathrm{join}}(T)=K]$ closely follows the system congestion level $\mathrm{E}[L(t)\mid A^{\mathrm{join}}(T)=K]$. Therefore, we conclude that the counterintuitive time-dependent behavior of $\mathrm{E}[L(t)\mid A^{\mathrm{join}}(T)=K]$ mainly stems from a bias introduced in the arrival process due to the condition that the total number of joined customers is given.

(a) $\mathrm{E}[L(t)\mid A^{\mathrm{join}}(T)=K]$

(b) $\mathrm{E}[\varDelta A^{\mathrm{join}}(t)\mid A^{\mathrm{join}}(T)=K]$

(c) $\mathrm{E}[\varDelta A^{\mathrm{balk}}(t)\mid A^{\mathrm{join}}(T)=K]$

(d) $\mathrm{E}[\varDelta A(t)\mid A^{\mathrm{join}}(T)=K]$

(e) $\mathrm{E}[\varDelta D(t)\mid A^{\mathrm{join}}(T)=K]$

(f) $\Pr[L(t)\geq 1\mid A^{\mathrm{join}}(T)=K]$

Figure 3: Transient behavior of several performance metrics for different values of the arrival rate $\lambda$ ($C_{V}=1.0$ and $\beta_{\ell}$ is given by (43)).

Figure 4 complements this observation by showing the unconditional mean number of joined customers $\mathrm{E}[A^{\mathrm{join}}(T)]$ in $[0,T]$ as a function of $\lambda$. By examining Figure 3 (a) alongside Figure 4, we see how $\mathrm{E}[A^{\mathrm{join}}(T)]$ is related to the time-dependent behavior of $\mathrm{E}[L(t)\mid A^{\mathrm{join}}(T)=K]$: for $\mathrm{E}[A^{\mathrm{join}}(T)]>K$, the trajectory typically peaks in the middle of the interval, while for $\mathrm{E}[A^{\mathrm{join}}(T)]<K$, it tends to increase steadily over time.

Figure 4: $\mathrm{E}[A^{\mathrm{join}}(T)]$ with respect to $\lambda$, where $C_{V}=1.0$ and $\beta_{\ell}$ is given by (43).

(a) $\mathrm{E}[L(t)\mid A^{\mathrm{join}}(T)=K]$

(b) $\mathrm{E}[\varDelta A^{\mathrm{join}}(t)\mid A^{\mathrm{join}}(T)=K]$

(c) $\mathrm{E}[\varDelta A^{\mathrm{balk}}(t)\mid A^{\mathrm{join}}(T)=K]$

(d) $\mathrm{E}[\varDelta A(t)\mid A^{\mathrm{join}}(T)=K]$

(e) $\mathrm{E}[\varDelta D(t)\mid A^{\mathrm{join}}(T)=K]$

(f) $\Pr[L(t)\geq 1\mid A^{\mathrm{join}}(T)=K]$

Figure 5: Transient behavior of several performance metrics for different values of the coefficient of variation $C_{V}$ of service times ($\lambda=1.5$ and $\beta_{\ell}$ is given by (43)).

We next examine the effect of the coefficient of variation $C_{V}$ of service times on the transient behavior of the system. Figure 5 shows the dynamics of various performance metrics for several values of $C_{V}$, with the arrival rate fixed at $\lambda=1.5$. From Figure 5 (a), we observe that $\mathrm{E}[L(t)\mid A^{\mathrm{join}}(T)=K]$ increases initially and then decreases after reaching a peak, regardless of $C_{V}$. As $C_{V}$ increases, this peak becomes more pronounced and occurs later in time. Figure 5 (a) and (c) show that the timing of balking events remains closely aligned with the level of system congestion, as reflected in $\mathrm{E}[L(t)\mid A^{\mathrm{join}}(T)=K]$, across all cases.

Figure 5 (d) and (e) show that a higher $C_{V}$ leads to a stronger bias in the number of departures: for large values of $C_{V}$, the departure rate $\mathrm{E}[\varDelta D(t)\mid A^{\mathrm{join}}(T)=K]$ tends to be elevated during the early period, while the arrival rate $\mathrm{E}[\varDelta A(t)\mid A^{\mathrm{join}}(T)=K]$ remains relatively stable. In contrast, when $C_{V}$ is small, the bias shifts to the arrival process, with $\mathrm{E}[\varDelta A(t)\mid A^{\mathrm{join}}(T)=K]$ exhibiting a notable time dependence. In this case, the departure rate gradually approaches one, reflecting saturation of system capacity.

In this setting, the unconditional expected number of arrivals is $\lambda T=150$, which exceeds the fixed number of joined customers $K=100$ by a wide margin. When $C_{V}$ is large, frequent occurrences of long service times lead to heavy congestion, resulting in more balking events. In contrast, when $C_{V}$ is small, such long service times are rare, and the relative variability in the arrival process becomes more influential, shifting the bias toward fluctuations in arrivals.

Finally, we examine how different forms of the entry probability $\beta_{\ell}$ affect system behavior. In addition to the baseline defined in (43), we consider the following four alternative shapes:

$$\beta_{\ell} =\max\left((1-0.04\ell)\exp(0.04\ell),0\right), \quad\ell=0,1,\ldots,100,$$ (44)

$$\beta_{\ell} =\max\left((1-0.03\ell)\exp(0.017\ell),0\right), \quad\ell=0,1,\ldots,100,$$ (45)

$$\beta_{\ell} =0.4+0.6(1-0.02\ell)\exp(-0.025\ell), \quad\ell=0,1,\ldots,100,$$ (46)

$$\beta_{\ell} =0.6+0.4(1-0.015\ell)\exp(-0.09\ell), \quad\ell=0,1,\ldots,100.$$ (47)

These choices are calibrated so that, without condition that the total number of joined customers is given with $\lambda=1.5$ and $C_{V}=1.0$, the expected number of balking customers in $[0,T]$ satisfies $36<\mathrm{E}[A^{\mathrm{balk}}(T)]<37$. The corresponding entry probability functions are plotted in Figure 6.

Figure 6: Form of the entry probability $\beta_{\ell}$ given by (43)–(47).

(a) $\mathrm{E}[L(t)\mid A^{\mathrm{join}}(T)=K]$

(b) $\mathrm{E}[\varDelta A^{\mathrm{join}}(t)\mid A^{\mathrm{join}}(T)=K]$

(c) $\mathrm{E}[\varDelta A^{\mathrm{balk}}(t)\mid A^{\mathrm{join}}(T)=K]$

(d) $\mathrm{E}[\varDelta A(t)\mid A^{\mathrm{join}}(T)=K]$

(e) $\mathrm{E}[\varDelta D(t)\mid A^{\mathrm{join}}(T)=K]$

(f) $\Pr[L(t)\geq 1\mid A^{\mathrm{join}}(T)=K]$

Figure 7: Transient behavior of several performance metrics for different forms (43)–(47) of the entry probability $\beta_{\ell}$ ($\lambda=1.5$ and $C_{V}=1.0$).

Figure 7 presents the transient behavior of the system under the five different entry probability functions $\beta_{\ell}$, with the arrival rate $\lambda=1.5$ and the coefficient of variation of service times $C_{V}=1.0$. As shown in Figure 7 (a), the mean number of customers in the system exhibits a unimodal shape in all cases, and the peak tends to be lower when $\beta_{\ell}$ is convex. Figure 7 (b) and (c) show that when $\beta_{\ell}$ is concave, balking tends to occur more frequently in the later part of the interval, suggesting that customer entries are more concentrated early on. From Figure 7 (c), (d), and (e), we observe that the convex $\beta_{\ell}$ leads to fewer arrivals and more frequent departures, with fewer number of balking customers.

These observations can be explained by the structure of the entry probability $\beta_{\ell}$. When $\beta_{\ell}$ is concave, balking increases sharply with congestion, leading to a greater imbalance in the timing of entries. In contrast, the convex $\beta_{\ell}$ results in a more gradual response to congestion, making it less likely that the system compensates by increasing congestion to induce balking; instead, the arrival rate decreases and the departure rate becomes more similar to that in the unconditioned system.

6 Conclusion

In this paper, we analyzed a Markovian queueing system with balking, where arrivals follow a nonhomogeneous Poisson process over a finite time interval $[0,T]$, and customers may balk depending on the number of customers present upon arrival. Our analysis focused on the system behavior conditioned on exactly $K$ customers having joined the system during the interval.

We addressed two main computational challenges arising in the analysis of the cumulative number of balking events under the fixed-$K$ setting. First, directly analyzing the extended Markov chain leads to an intractably large state space. To overcome this, we introduced a PGF representation and derived a first-order differential equation satisfied by the unconditioned PGF. This formulation enables efficient computation without explicitly tracking the number of balking customers.

Second, to avoid computing the transition probabilities $\Pr[A^{\mathrm{join}}(T)=K\mid A^{\mathrm{join}}(t)=k,L(t)=\ell,S(t)=s]$ repeatedly for multiple time points, we considered a backward-in-time approach. We showed that the conditional probability $\Pr[A^{\mathrm{join}}(T)=K\mid A^{\mathrm{join}}(T-u)=k,L(T-u)=\ell,S(T-u)=s]$ satisfies a linear differential equation in $u$, allowing efficient evaluation across a range of values.

For the case where system parameters are piecewise constant, we developed a numerical procedure to compute the conditional factorial moments. Numerical examples for an M/PH/1 queue illustrated that conditioning on the total number of joined customers induces time-dependent dynamics that differ markedly from typical steady-state behavior. For instance, the expected number of customers in the system may increase and then decrease even under a constant arrival rate. We showed that such behaviors are driven by nontrivial biases in both arrivals and departures caused by the conditioning.

In this paper, we assumed that once customers join the system, they remain until service completion. Under this assumption, the total number of customers served equals the number of joined customers. However, in many real-world systems, customers may renege before completing service. In such settings, conditioning on the number of served customers no longer coincides with conditioning on the number of entries. Extending our approach to such cases, where early departures occur, remains an important direction for future work.

Appendix A Transition rate matrix and evaluation of (41) and (42) in the M/PH/1 queue

In this section, we present the specific form of the northwest corner block of the transition rate matrix in an M/PH/1 queue and describe the concrete calculation method for $\bm{x}_{n}^{(i)}(m)(\bm{P}_{n}-\theta_{n}^{-1}\lambda_{n}\bm{I})$ appearing on the right-hand side of equations (41) and (42). For the sake of notational simplicity, we omit subscripts and define:

$$\lambda_{1}\eqqcolon\lambda,\qquad\bm{Q}_{1}\eqqcolon\bm{Q}.$$

Additionally, we define:

$$\theta\coloneqq\max_{i}|[\bm{Q}]_{i,i}|+\lambda,\qquad\bm{P}\coloneqq\bm{I}+\theta^{-1}\bm{Q}.$$

The service time distribution follows a phase-type distribution characterized by the initial state $\bm{\gamma}$ and the transition rate matrix $\bm{\Gamma}$. The initial state $\bm{\gamma}$ and the transition rate matrix $\bm{\Gamma}$ for given mean and coefficient of variation can be obtained from [19, pp. 353, 358–359].

Let $\bm{L}(t)\coloneqq(L(t),S(t))$ denote the number of customers in the system and the service phase. The northwest corner block $\bm{Q}$ of the transition rate matrix of the continuous-time Markov chain $(A^{\mathrm{join}}(t),\bm{L}(t))_{t\geq 0}$ is then given by

$$\bm{Q}=\left[\begin{array}[]{cccccc}\bm{X}_{0}&\bm{Y}_{0}&\bm{O}&\cdots&\bm{O}&\bm{O}\\ \bm{O}&\bm{X}_{1}&\bm{Y}_{1}&\cdots&\bm{O}&\bm{O}\\ \bm{O}&\bm{O}&\bm{X}_{2}&\cdots&\bm{O}&\bm{O}\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ \bm{O}&\bm{O}&\bm{O}&\cdots&\bm{X}_{K-1}&\bm{Y}_{K-1}\\ \bm{O}&\bm{O}&\bm{O}&\cdots&\bm{O}&\bm{X}_{K}\end{array}\right].$$

where $\bm{X}_{k}$ and $\bm{Y}_{k}$ (for $k=0,1,\ldots$) are given by

$$\bm{X}_{k} ={\left[\begin{array}[]{cccccc}-\lambda\beta_{0}&\bm{0}&\bm{0}&\cdots&\bm{0}&\bm{0}\\ -\bm{\Gamma}\bm{e}&-\lambda\beta_{1}\bm{I}+\bm{\Gamma}&\bm{O}&\cdots&\bm{O}&\bm{O}\\ \bm{0}&-\bm{\Gamma}\bm{e}\bm{\gamma}&-\lambda\beta_{2}\bm{I}+\bm{\Gamma}&\cdots&\bm{O}&\bm{O}\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ \bm{0}&\bm{O}&\bm{O}&\cdots&-\lambda\beta_{k-1}\bm{I}+\bm{\Gamma}&\bm{O}\\ \bm{0}&\bm{O}&\bm{O}&\cdots&-\bm{\Gamma}\bm{e}&-\lambda\beta_{k}\bm{I}+\bm{\Gamma}\end{array}\right]},$$

$$\bm{Y}_{k} =\left[\begin{array}[]{cccccc}0&\lambda\beta_{0}\bm{\gamma}&\bm{0}&\bm{0}&\cdots&\bm{0}\\ \bm{0}&\bm{O}&\lambda\beta_{1}\bm{I}&\bm{O}&\cdots&\bm{O}\\ \bm{0}&\bm{O}&\bm{O}&\lambda\beta_{2}\bm{I}&\cdots&\bm{O}\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\ \bm{0}&\bm{O}&\bm{O}&\bm{O}&\cdots&\lambda\beta_{k}\bm{I}\end{array}\right].$$

Due to the sparsity of $\bm{Q}$, the computation of the matrix product $\bm{P}-\theta^{-1}\lambda\bm{I}$ from the right can be optimized as follows. Given a row vector $\bm{y}$, let $\bm{y}^{\prime}\coloneqq\bm{y}(\bm{P}-\theta^{-1}\lambda\bm{I})$. The elements $\bm{y}_{k,\ell}^{\prime}$ of $\bm{y}^{\prime}$, corresponding to the state where the cumulative number of joins is $k$ and the number of customers in the system is $\ell$, can be computed using the elements $\bm{y}_{k,\ell}$ of $\bm{y}$ as follows:

$$y_{0,0}^{\prime} =y_{0,0}(1-\theta^{-1}\lambda(\beta_{0}+1)),\hskip 108.20023ptk=0,\ \ell=0,$$

$$y_{k,0}^{\prime} =y_{k,0}(1-\theta^{-1}\lambda(\beta_{0}+1))+\bm{y}_{k,1}\theta^{-1}(-\bm{\Gamma}\bm{e}),\hskip 36.70003ptk=1,2,\ldots,K,\ \ell=0,$$

$$\bm{y}_{1,1}^{\prime} =y_{0,0}\theta^{-1}\lambda\beta_{0}\bm{\gamma}+\bm{y}_{1,1}(\bm{I}-\theta^{-1}(\lambda(\beta_{1}+1)\bm{I}-\bm{\Gamma})),\hskip 10.00002ptk=1,\ \ell=1,$$

$$\bm{y}_{k,1}^{\prime} =y_{k-1,0}\theta^{-1}\lambda\beta_{0}\bm{\gamma}+\bm{y}_{k,1}(\bm{I}-\theta^{-1}(\lambda(\beta_{1}+1)\bm{I}-\bm{\Gamma}))+\bm{y}_{k,2}\theta^{-1}(-\bm{\Gamma}\bm{e})\bm{\gamma},$$

$$\hskip 170.00026ptk=2,3,\ldots,K,\ \ell=1,$$

$$\bm{y}_{k,k}^{\prime} =\bm{y}_{k-1,k-1}\theta^{-1}\lambda\beta_{k-1}{\bm{I}}+\bm{y}_{k,k}(\bm{I}-\theta^{-1}(\lambda(\beta_{k}+1)\bm{I}-\bm{\Gamma})),$$

$$\hskip 170.00026ptk=2,3,\ldots,K,\ \ell=k,$$

$$\bm{y}_{k,\ell}^{\prime} =\bm{y}_{k-1,\ell-1}\theta^{-1}\lambda\beta_{\ell-1}{\bm{I}}+\bm{y}_{k,\ell}(\bm{I}-\theta^{-1}(\lambda(\beta_{\ell}+1)\bm{I}-\bm{\Gamma}))+\bm{y}_{k,\ell+1}\theta^{-1}(-\bm{\Gamma}\bm{e})\bm{\gamma},$$

$$\hskip 170.00026ptk=2,3,\ldots,K,\ \ell=2,3,\ldots,k-1.$$

$\bm{x}_{n}^{(i)}(m)(\bm{P}_{n}-\theta_{n}^{-1}\lambda_{n}\bm{I})$ appearing on the right-hand side of (41) and (42) can be computed using the above formulas. Moreover, for $\bm{P}_{n}\bm{q}(T-T_{N-n+1},K)$ in (33), a similar computation can be performed utilizing the sparsity of $\bm{Q}$.