Abstract

Translate

This paper analyzes a finite-capacity $G I / M / 2 / N$ queue with two heterogeneous servers operating under a multiple working-vacation policy, Bernoulli feedback, and customer impatience. Using the supplementary-variable technique in tandem with a tailored recursive scheme, we derive the stationary distributions of the system size as observed at pre-arrival instants and at arbitrary epochs. From these, we obtain explicit expressions for key performance metrics, including blocking probability, average reneging rate, mean queue length, mean sojourn time, throughput, and server utilizations. We then embed these metrics in an economic cost function and determine service-rate settings that minimize the total expected cost via the Bat Algorithm. Numerical experiments implemented in R validate the analysis and quantify the managerial impact of the vacation, feedback, and impatience parameters through sensitivity studies. The framework accommodates general renewal arrivals ( $G I$ ), thereby extending classical ( $M / M / 2 / N$ ) results to more realistic input processes while preserving computational tractability. Beyond methodological interest, the results yield actionable design guidance: (i) they separate Palm and time-stationary viewpoints cleanly under non-Poisson input, (ii) they retain heterogeneity throughout all formulas, and (iii) they provide a cost–optimization pipeline that can be deployed with routine numerical effort. Methodologically, we (i) characterize the generator of the augmented piecewise–deterministic Markov process and prove the existence/uniqueness of the stationary law on the finite state space, (ii) derive an explicit Palm–time conversion formula valid for non-Poisson input, (iii) show that the boundary-value recursion for the Laplace–Stieltjes transforms runs in linear time $O (N)$ and is numerically stable, and (iv) provide influence-function (IPA) sensitivities of performance metrics with respect to $(μ_{1}, μ_{2}, ν, α, ϕ, β)$ .

Full text

Turn on search term navigation

Translate

1. Introduction

Modern service systems routinely face finite buffers, rework loops, energy-aware operation, and behavioral customer responses. Treating these features jointly and without restricting arrivals to be Poisson is essential for credible design. Queuing theory (In 1909, Agner Krarup Erlang, a Danish engineer employed at the Copenhagen Telephone Exchange, published the pioneering paper that laid the foundation for what is now known as queuing theory [1]. Erlang modeled the arrival of telephone calls at an exchange as a Poisson process and subsequently provided solutions to the M/D/1 queue in 1917 and the M/D/k queuing model in 1920 [2])—aptly characterized as the mathematics of waiting lines—furnishes a rigorous probabilistic language for analyzing congestion, delay, and capacity allocation in service systems driven by stochastic demand. Its core analytic tools include Markovian and renewal models, embedded chains, transform and generating-function techniques, matrix-analytic methods for quasi-birth–death (QBD) structures, and diffusion/heavy-traffic limits; together, these enable the derivation of stability conditions, stationary distributions, and performance descriptors such as waiting times, loss probabilities, utilizations, and throughput (see the classics [3,4,5,6,7,8,9,10] and Little’s law [11]). For recent references, see [12,13,14,15,16,17,18,19,20,21]. This apparatus is not merely descriptive: it supports normative decisions in mission-critical environments—including emergency response, communication networks, cloud-computing platforms, healthcare delivery, transportation, and logistics—where service quality and economic efficiency must be optimized jointly under uncertainty (cf. [22,23,24]). By articulating principled trade-offs between limited resources and stochastic variability, queuing theory has become a cornerstone for designing resilient, cost-effective service systems.

A salient evolution in the literature is the explicit modeling of customer impatience, which manifests primarily as (i) balking, whereby an arrival observing unfavorable conditions (e.g., long visible queues or high predicted delay) elects not to join, and (ii) reneging, whereby a customer who has joined subsequently abandons after an excessive waiting time. From an operational perspective, impatience fundamentally reshapes effective load, stability, and welfare: it removes demand endogenously, introduces state- and history-dependent losses, and alters equilibrium congestion levels in ways that are highly sensitive to information structure and behavioral assumptions. Classical and modern studies have established its impact across telecommunications, web services, healthcare, and manufacturing [25,26,27,28,29]. From an analytical standpoint, impatience complicates standard embeddings (e.g., arrival-epoch Markov chains) and undermines the convenience of PASTA outside the Poisson regime; even under Poisson input, the state dependence of abandonment and balking may require careful Palm-calculus distinctions between time-average and arrival-average measures [4].

A complementary stream of work studies working-vacation (WV) policies, under which servers reduce—rather than suspend—their service rate during vacation periods. The $M / M / 1 / WV$ model introduced by [30] demonstrated that partial service during vacations can yield favorable sojourn-time and queue-length distributions relative to hard shutdowns, by tempering transient congestion while preserving energy or labor savings. Subsequent generalizations have spanned bulk arrivals [31,32], Erlang and general service distributions [33,34,35], multi-server configurations $M / M / c$ [36,37,38,39], and discrete-time analogues [40,41]. Methodologically, WV policies induce a modulated service mechanism that couples vacation phases with the service process; tractability often follows from phase-type representations or supplementary variables that capture the WV phase jointly with the queue state.

Real traffic is frequently over-dispersed, bursty, or correlated, rendering Poisson arrivals an inadequate approximation. Renewal-input models $G I / \cdot / \cdot$ capture such features parsimoniously. For $G I / M / 1$ and finite-capacity variants with WV, notable contributions include explicit sojourn-time characterizations, interruption mechanisms, and phase-type vacation structures [31,42,43,44,45,46]. The passage from Poisson to general renewal input removes PASTA and necessitates a careful separation of time-stationary quantities from arrival-stationary (Palm) quantities. In particular, performance evaluation at pre-arrival epochs typically requires either (i) Markovianization by tracking residual inter-arrival time as a supplementary variable, or (ii) matrix-analytic constructions on an augmented state space; passage between Palm and time averages demands size-biasing and residual-life arguments.

Combining impatience with WV dynamics raises subtle questions about (i) state-dependent stability and ergodicity, (ii) the structural properties of stationary laws, and (iii) optimal control under explicit economic objectives. For Markovian baselines, explicit probability-generating functions across vacation and busy periods have been obtained in [47], while finite buffers, policy comparisons, and cost optimization under impatience/WV interactions appear in [48,49,50,51,52,53,54,55]. A related and practically pervasive mechanism is feedback: after service completion, a customer returns to the queue with fixed probability, modeling rework loops or packet retransmissions. Originating with Takács [56], feedback has been analyzed in Markovian and renewal settings [57,58,59,60,61,62,63,64]. Analytically, feedback alters the effective arrival rate in a state-dependent manner and couples departure streams to subsequent congestion, generating nontrivial interactions with vacation phases and impatience.

Most multi-server studies posit homogeneous servers, an assumption natural for automated systems but frequently violated in human-service or mixed-technology environments. Heterogeneity ( $μ_{1} \neq μ_{2}$ ) destroys exchangeability, induces state-dependent departure rates tied to server occupancy configurations, and generally breaks product-form heuristics. Early work on heterogeneity includes [65,66,67,68,69], with focused developments for $M / M / 2$ (including WV) and transient analyses under impatience [70,71,72]. In contrast, renewal-input $G I / M / 2 / N$ queues with heterogeneous servers remain comparatively underexplored; to our knowledge, [73] addresses performance metrics but leaves open a broader structural and economic treatment that accommodates impatience, feedback, and WV control simultaneously.

This paper: model, analytic approach, and scope. We consider a finite-buffer $G I / M / 2 / N$ system with two heterogeneous exponential servers, a working-vacation policy, customer impatience in the form of balking and reneging, and Bernoulli feedback. The system captures operational features that co-occur in practice: (i) non-Poisson variability in demand (renewal inter-arrivals), (ii) endogenous demand attrition (balking, reneging), (iii) flexible, partially productive capacity modulation (working vacations), (iv) rework/retransmission loops (feedback), and (v) finite-capacity constraints (blocking). Formally, let ${A_{k}}_{k \geq 1}$ be i.i.d. inter-arrival times with cdf $A (u)$ , density $a (u)$ , and mean $1 / λ \in (0, \infty)$ ; service times are exponential with heterogeneous rates $μ_{1}, μ_{2} > 0$ . A WV policy modulates the active server’s rate according to a finite-state vacation phase (standard variants include single and multiple vacation schemes). Balking is captured by a state-dependent acceptance probability $\bar{θ} (i) \in [0, 1]$ when i customers are present upon arrival (with $\bar{θ} (N) = 1$ under finite capacity), while reneging is modeled via exponential patience with hazard $α (i)$ that may depend on the queue length. Feedback is Bernoulli with probability $\bar{β} \in [0, 1)$ , applied to each service completion. Customers are served under FCFS; ties in server selection are resolved by a fixed priority or randomized rule. The system state records the queue length (including any in service), server-occupancy configuration (which server(s) are busy), WV phase, and remaining inter-arrival time. Because $N < \infty$ , the augmented process is irreducible and aperiodic on a finite state space, hence admitting a unique stationary distribution for any fixed parameter vector; all performance measures below are computed in this stationary regime. We emphasize that finite capacity and exponential impatience together preclude null recurrence even under heavy traffic; positive recurrence is automatic in our augmented finite-state description. Please see Table 1 for details.

Contributions and organization. Our analysis proceeds via the supplementary-variable technique [77], which Markovianizes the renewal input by tracking the remaining inter-arrival time and couples it with the WV phase and occupancy configuration. The resulting piecewise-deterministic Markov process admits a system of linear balance equations whose structure allows a recursive solution for (i) pre-arrival (Palm) stationary probabilities and (ii) time-stationary probabilities at arbitrary epochs, together with explicit Palm–time conversions to reconcile performance metrics that require one or the other viewpoint. This framework facilitates the computation of loss probability, mean queue length, waiting- and sojourn-time descriptors, and server utilizations, while preserving the heterogeneity of $μ_{1}$ and $μ_{2}$ and the modulation induced by WV phases. We then embed the stationary descriptors into an economic objective that trades off (a) capacity costs (including heterogeneous service rates and WV control parameters), (b) congestion costs (delay and queue-length penalties), (c) abandonment costs (reneging-induced losses), and (d) rework costs arising from feedback. The resulting design problem is nonconvex due to the nonlinear dependence of performance measures on $(μ_{1}, μ_{2})$ , WV parameters, and impatience/feedback primitives; we therefore adopt the Bat Algorithm [74], a metaheuristic with robust exploration–exploitation behavior and favorable convergence in applied queuing design [75,76]. To our knowledge, this is the first study to combine renewal arrivals, heterogeneity, working vacations, impatience (balking and reneging), and Bernoulli feedback in the same finite-capacity two-server model with a full Palm–time reconciliation.

Unified state augmentation under renewal input. By including the remaining inter-arrival time in the state, we retain exactness for general $G I$ input without resorting to Poissonization or diffusion scaling. This preserves compatibility with balking and reneging rules that depend on instantaneous congestion and the WV phase, while enabling Palm-compatible pre-arrival descriptors. Heterogeneity-aware service dynamics. We treat $μ_{1} \neq μ_{2}$ explicitly and distinguish occupancy configurations—(0), (1 on server 1), (1 on server 2), or (2 in service)—which generate distinct departure hazards and interact with WV states; this non-exchangeability is retained in all performance expressions. Impatience and feedback coupling. Balking modifies the effective arrival stream in a state-dependent way at pre-arrival epochs; reneging introduces additional departures from the queue that depend on both the waiting population and the WV phase; and feedback couples completions to subsequent arrivals. Our formulation keeps these mechanisms explicit and compatible, rather than approximating by constant effective rates. Palm–time reconciliation. Many design metrics (e.g., loss on arrival, experienced waiting time) are naturally Palm quantities, whereas costs tied to time-average congestion require time-stationary laws. We make these distinctions explicit and provide the necessary conversions, avoiding implicit PASTA assumptions that fail under $G I$ input. Economic design under WV control. We formalize a cost functional that transparently prices (i) heterogeneity in service effort, (ii) responsiveness under WV modulation, and (iii) behavioral losses due to impatience; the resulting optimization illustrates how calibrated heterogeneity and WV tuning jointly improve responsiveness and cost.

Let $π^{(A)}$ denote the pre-arrival stationary law on the augmented state space and $π$ the time-stationary law. We first establish that, for any fixed parameter vector with finite buffer $N < \infty$ , the augmented process forms an irreducible and aperiodic Markov chain with a unique stationary distribution; hence, both $π^{(A)}$ and $π$ exist and are unique. For completeness, we further delineate Foster–Lyapunov conditions under which these conclusions extend to $N = \infty$ , while our numerical investigations focus on the finite-buffer setting. Next, we obtain closed recursive characterizations of $π^{(A)}$ and $π$ by balancing over renewal cycles and WV phases, and we develop Laplace–Stieltjes transform representations that render computations tractable for finite N. Building on these representations, we derive explicit performance formulas: (i) the loss probability arising from balking or blocking, (ii) the mean queue length together with the utilizations of heterogeneous servers, (iii) the waiting-time and sojourn-time distributions from both Palm and time perspectives, and (iv) sensitivities with respect to $μ_{1}, μ_{2}$ , the WV parameters, and the impatience/feedback primitives. Finally, within a convex–nonconvex composite cost framework, we demonstrate that a Bat Algorithm search identifies operating points that outperform homogeneous-server baselines and hard-vacation controls across broad regions of the parameter space [74,75,76].

Structure of the paper. Section 2 formalizes the model state space, WV mechanism, balking/reneging rules, feedback, and supplementary-variable augmentation for renewal input, with a precise definition of pre-arrival (Palm) versus time-average observables. Section 3 develops the stationary analysis: recursive balance equations, Laplace–Stieltjes transform solutions, and Palm–time conversions. Section 4 reports performance characteristics and sensitivity analyses, including heterogeneity-induced asymmetries in utilization and their interaction with WV phases. Section 5 formulates the economic objective and presents the optimization methodology. Section 6 provides numerical experiments and design maps that expose managerial trade-offs and robustness. Section 7 illustrates the implications on a flexible production facility benchmark. Section 8 concludes with research directions, including diffusion approximations for large N, control under partial state information, and learning-based calibration of impatience and feedback.

2. Model Mathematical Description

This study investigates a single-class $G I / M / 2 / N$ queue with two heterogeneous servers, a common FIFO queue, state-dependent balking, exponential reneging, Bernoulli feedback, and multiple working vacations taken by the slower server. A schematic of the system appears in Figure 1. The mathematical description relies on the following notation and assumptions.

Inter-arrival times. Inter-arrival times of successive arrivals form an i.i.d. sequence with cumulative distribution function $A (u)$ (density $a (u)$ for $u \geq 0$ ), Laplace–Stieltjes transform

$A^{*} (s) = \int_{[0, \infty)} e^{- s u} d A (u), s \geq 0,$

and mean inter-arrival time

$E [U] = - \frac{d}{d s} A^{*} (s) |_{s = 0} = \frac{1}{λ} \in (0, \infty) .$

Finite capacity. The total system capacity is $N \in N$ , counting customers both in service and waiting. Arrivals that would exceed capacity are blocked. Finite N ensures positive recurrence of the augmented Markov chain.

State-dependent balking. Upon arrival to a system containing $i \in {0, 1, \dots, N}$ customers, a customer is admitted with probability $θ_{i}$ and balks with probability ${\bar{θ}}_{i} : = 1 - θ_{i}$ . We assume

$θ_{0} = θ_{1} = 1, 0 \leq θ_{i + 1} \leq θ_{i} \leq 1 (2 \leq i \leq N - 1), θ_{N} = 0,$

so no one joins a full system and admission propensities are nonincreasing in congestion. (Feedback customers are subject to the same rule; see below) Alternative information structures (delayed or noisy observations) can be incorporated by redefining

(θ_{i})

; our analysis is unchanged.

Service structure and discipline. There is a single common FIFO queue feeding two heterogeneous servers. Server 1 and Server 2 provide exponential service with rates $μ_{1}$ and $μ_{2}$ , respectively, with $μ_{2} \leq μ_{1}$ . Server 1 is always available (no vacations).

Working vacations of Server 2. If Server 2 becomes idle when the queue is empty, it immediately enters a working vacation, whose duration is exponential with rate $ϕ$ . During a working vacation, Server 2 remains active but at a reduced exponential service rate $ν \leq μ_{2}$ . At the vacation’s end, Server 2 switches to regular service if at least one customer is waiting; otherwise, it initiates a new working vacation. This induces “multiple” working vacations separated by possible service epochs.

Reneging (impatience). Whenever both servers are busy and the system contains $i \geq 2$ customers, the queue length equals $i - 2$ . Each waiting customer runs an independent exponential impatience timer $T \sim Exp (α)$ . A customer whose service has not started before their timer expires abandons the system permanently. Impatience clocks are independent of the queue length process and of all other primitives.

Bernoulli feedback. Upon service completion (either during regular operation or while Server 2 is on a working vacation), a customer departs the system with probability $β$ and, with complementary probability $\bar{β} : = 1 - β$ , instantaneously feeds back to the input stream and behaves as a fresh arrival (i.e., is subjected to balking and capacity constraints in the current state). This mechanism models rework/retx loops; its coupling with balking is crucial in finite capacity.

All primitive sources of randomness—inter-arrival times, regular and working vacation service times, vacation durations, impatience timers, and feedback decisions—are mutually independent.

3. Steady-State Solution

This section analyzes the queuing system in steady-state using the supplementary variable technique to derive the system’s long-term behavior. The state of the system at time t is modeled as a continuous-time Markov process ${(L (t), J (t), W (t)), t \geq 0}$ , with the state space defined as

$Ω = {(0, 0, u) \cup (1, 0, u) \cup (i, 1, u) \cup (i, 2, u) : 1 \leq i \leq N, u \geq 0},$

where the system variables are

$L (t)$ : The number of customers in the system at time t, including those in service.

$J (t)$ : The state of the server at time t, defined as

$J (t) = \{\begin{matrix} 0, & if Server 2 is idle during the working vacation period, \\ 1, & if Server 2 is busy during the working vacation period, \\ 2, & if Server 2 is busy during the regular service period . \end{matrix}$

$W (t)$ : The remaining inter-arrival time for the next customer arrival at time t.

We define the joint probabilities for the system states as follows:

$\begin{matrix} P_{i, 0} (u, t) d u & = & P (L (t) = i, u \leq W (t) < u + d u, J (t) = 0), u \geq 0, i = 0, 1, \\ P_{i, j} (u, t) d u & = & P (L (t) = i, u \leq W (t) < u + d u, J (t) = j), u \geq 0, j = 1, 2, 1 \leq i \leq N . \end{matrix}$

In steady-state, as

t \to \infty

, we have

$\begin{matrix} P_{i, 0} (u) & = & lim_{t \to \infty} P_{i, 0} (u, t), i = 0, 1, \\ P_{i, j} (u) & = & lim_{t \to \infty} P_{i, j} (u, t), j = 1, 2, 1 \leq i \leq N . \end{matrix}$

By employing the probabilistic arguments and using the remaining inter-arrival time as the supplementary variable, we observe the state of the system at two consecutive time epochs, t and

t + d t .

Taking

lim_{t \to \infty}

and after some simplification, we have the following set of differential difference equations.

(1) $\begin{matrix} - P_{0, 0}^{(1)} (u) & = & β [μ_{1} P_{1, 0} (u) + ν P_{1, 1} (u) + μ_{2} P_{1, 2} (u)], \end{matrix}$

(2) $\begin{matrix} - P_{1, 0}^{(1)} (u) & = & - β μ_{1} P_{1, 0} (u) + β [ν P_{2, 1} (u) + μ_{2} P_{2, 2} (u)] + a (u) P_{0, 0} (0), \end{matrix}$

(3) $\begin{matrix} - P_{1, 1}^{(1)} (u) & = & - (ϕ + β ν) P_{1, 1} (u) + β μ_{1} P_{2, 1} (u), \end{matrix}$

(4) $\begin{matrix} - P_{2, 1}^{(1)} (u) & = & - [β (μ_{1} + ν) + ϕ] P_{2, 1} (u) + [β (μ_{1} + ν) + α] P_{3, 1} (u) \\ + a (u) [P_{1, 0} (0) + P_{1, 1} (0), \end{matrix}$

(5) $\begin{matrix} - P_{i, 1}^{(1)} (u) & = & - [β (μ_{1} + ν) + ϕ + (i - 2) α] P_{i, 1} (u) + [β (μ_{1} + ν) + (i - 1) α] P_{i + 1, 1} (u) \\ + a (u) [θ_{i - 1} P_{i - 1, 1} (0) + {\bar{θ}}_{i} P_{i, 1} (0)], 3 \leq i \leq N - 1, \end{matrix}$

(6) $\begin{matrix} - P_{N, 1}^{(1)} (u) & = & - [β (μ_{1} + ν) + ϕ + (N - 2) α] P_{N, 1} (u) + a (u) [θ_{N - 1} P_{N - 1, 1} (0) + P_{N, 1} (0)], \end{matrix}$

(7) $\begin{matrix} - P_{1, 2}^{(1)} (u) & = & - β μ_{2} P_{1, 2} (u) + ϕ P_{1, 1} (u) + β μ_{1} P_{2, 2} (u), \end{matrix}$

(8) $\begin{matrix} - P_{i, 2}^{(1)} (u) & = & - [β (μ_{1} + μ_{2}) + (i - 2) α] P_{i, 2} (u) + [β (μ_{1} + μ_{2}) + (i - 1) α] P_{i + 1, 2} (u) \\ + ϕ P_{i, 1} (u) + a (u) [θ_{i - 1} P_{i - 1, 2} (0) + {\bar{θ}}_{i} P_{i, 2} (0)], 2 \leq i \leq N - 1, \end{matrix}$

(9) $\begin{matrix} - P_{N, 2}^{(1)} (u) & = & - [β (μ_{1} + μ_{2}) + (N - 2) α) P_{N, 2} (u) + ϕ P_{N, 1} (u) \\ + a (u) [θ_{N - 1} P_{N - 1, 2} (0) + P_{N, 2} (0)], \end{matrix}$

where

P_{i, 0} (0), i = 0, 1

and

P_{i, j} (0), j = 1, 2, 1 \leq i \leq N

are the respective rate probabilities with the remaining inter-arrival time equal to zero denoting that an arrival is about to occur. Then, we introduce the following Laplace–Stieltjes transforms of the steady-state probabilities as follows:

$\begin{matrix} P_{i, 0}^{*} (s) & = & \int_{0}^{\infty} e^{- s u} P_{i, 0} (u) d u, i = 0, 1, \\ P_{i, j}^{*} (s) & = & \int_{0}^{\infty} e^{- s u} P_{i, j} (u) d u, j = 1, 2, 1 \leq i \leq N . \end{matrix}$

Let $P_{i, 0} = P_{i, 0}^{*} (0)$ for $i \in {0, 1}$ and $P_{i, j} = P_{i, j}^{*} (0)$ for $j \in {1, 2}$ and $1 \leq i \leq N$ denote the steady–state probabilities of having i customers in the system when the server is in state $j \in {0, 1, 2}$ , observed at an arbitrary epoch. For $2 \leq i \leq N$ , define

$δ_{i} : = β (μ_{1} + ν) + ϕ + (i - 2) α, ϑ_{i} : = β (μ_{1} + μ_{2}) + (i - 2) α .$

Multiplying Equations (1)–(9) by

e^{- s u}

and integrating with respect to u over

[0, \infty)

, we obtain the transformed relations

(10) $\begin{matrix} - s P_{0, 0}^{*} (s) & = & - P_{0, 0} (0) + β [μ_{1} P_{1, 0}^{*} (s) + ν P_{1, 1}^{*} (s) + μ_{2} P_{1, 2}^{*} (s)] . \end{matrix}$

(11) $\begin{matrix} (β μ_{1} - s) P_{1, 0}^{*} (s) & = & - P_{1, 0} (0) + β [ν P_{2, 1}^{*} (s) + μ_{2} P_{2, 2}^{*} (s)] + A^{*} (s) P_{0, 0} (0) . \end{matrix}$

(12) $\begin{matrix} (ϕ + β ν - s) P_{1, 1}^{*} (s) & = & P_{1, 1} (0) + β μ_{1} P_{2, 1}^{*} (s) . \end{matrix}$

(13) $\begin{matrix} (δ_{2} - s) P_{2, 1}^{*} (s) & = & - P_{2, 1} (0) + (δ_{3} - ϕ) P_{3, 1}^{*} (s) \\ + A^{*} (s) [P_{1, 0} (0) + P_{1, 1} (0) + {\bar{θ}}_{2} P_{2, 1} (0)] . \end{matrix}$

(14) $\begin{matrix} (δ_{i} - s) P_{i, 1}^{*} (s) & = & - P_{i, 1} (0) + (δ_{i + 1} - ϕ) P_{i + 1, 1}^{*} (s) \\ + A^{*} (s) [θ_{i - 1} P_{i - 1, 1} (0) + {\bar{θ}}_{i} P_{i, 1} (0)], 3 \leq i \leq N - 1 . \end{matrix}$

(15) $\begin{matrix} (δ_{N} - s) P_{N, 1}^{*} (s) & = & - P_{N, 1} (0) + A^{*} (s) [θ_{N - 1} P_{N - 1, 1} (0) + P_{N, 1} (0)] . \end{matrix}$

(16) $\begin{matrix} (β μ_{2} - s) P_{1, 2}^{*} (s) & = & - P_{1, 2} (0) + ϕ P_{1, 1}^{*} (s) + β μ_{1} P_{2, 2}^{*} (s) . \end{matrix}$

(17) $\begin{matrix} (ϑ_{i} - s) P_{i, 2}^{*} (s) & = & - P_{i, 2} (0) + ϕ P_{i, 1}^{*} (s) + ϑ_{i + 1} P_{i + 1, 2}^{*} (s) \\ + A^{*} (s) [θ_{i - 1} P_{i - 1, 2} (0) + {\bar{θ}}_{i} P_{i, 2} (0)], 2 \leq i \leq N - 1 . \end{matrix}$

(18) $\begin{matrix} (ϑ_{N} - s) P_{N, 2}^{*} (s) & = & - P_{N, 2} (0) + ϕ P_{N, 1}^{*} (s) + A^{*} (s) [θ_{N - 1} P_{N - 1, 2} (0) + P_{N, 2} (0)] . \end{matrix}$

Summing (10)–(18) and simplifying yields

(19) $\begin{matrix} - A^{*} (s) (\sum_{i = 0}^{1} P_{i, 0} (0) + \sum_{i = 1}^{N} (P_{i, 1} (0) + P_{i, 2} (0))) = s (\sum_{i = 0}^{1} P_{i, 0}^{*} (s) + \sum_{i = 1}^{N} (P_{i, 1}^{*} (s) + P_{i, 2}^{*} (s))) . \end{matrix}$

Differentiating (19) with respect to s, letting $s \to 0$ , and using the normalization condition

$\sum_{i = 0}^{1} P_{i, 0} + \sum_{i = 1}^{N} (P_{i, 1} + P_{i, 2}) = 1,$

we obtain

(20) $\sum_{i = 0}^{1} P_{i, 0} (0) + \sum_{i = 1}^{N} (P_{i, 1} (0) + P_{i, 2} (0)) = λ .$

Equation (20) expresses conservation of flow: the long–run average inflow to the system per unit time equals the mean arrival rate

λ

. For

i \in {0, 1}

, define

P_{i, 0}^{-}

as the steady–state probability (at a pre–arrival epoch) that an arriving customer sees i customers in the system and the server in state

j = 0

. For

j \in {1, 2}

and

1 \leq i \leq N

, define

P_{i, j}^{-}

analogously as the steady–state probability that an arrival sees i customers and server state j. Our first goal is to connect these pre-arrival probabilities to the stationary descriptors at an arbitrary time.

Let $L (t)$ denote the number of customers in the system at time t, let $S (t) \in {0, 1, 2}$ denote the server state, and let $U (t)$ denote the remaining inter-arrival time at t. By Bayes’ formula for conditional probabilities, for each admissible $(i, j)$ , we have

$\begin{matrix} P_{i, j}^{-} & = lim_{t \to \infty} P (L (t) = i, S (t) = j | U (t) = 0) \\ = lim_{t \to \infty} \frac{P (L (t) = i, S (t) = j, U (t) = 0)}{P (U (t) = 0)} . \end{matrix}$

Consequently,

(21) $P_{i, j}^{-} = \frac{1}{λ} P_{i, j} (0), \begin{matrix} j = 0, i \in {0, 1}, \\ j \in {1, 2}, 1 \leq i \leq N, \end{matrix}$

where

λ

is the arrival rate and

P_{i, j} (0)

is the corresponding boundary value at

s = 0

of the transform (defined below). To determine the pre-arrival probabilities, it suffices to compute

$\begin{matrix} P_{i, 0} (0) (i = 0, 1), P_{i, j} (0) (j = 1, 2, 1 \leq i \leq N), \\ P_{i, 0}^{*} (s) (i = 0, 1), P_{i, j}^{*} (s) (j = 1, 2, 1 \leq i \leq N), \end{matrix}$

recursively in terms of

P_{N, 0} (0)

(details below). Throughout,

A^{*} (\cdot)

denotes the relevant transform and

A^{* (k)}

its kth derivative.

Step 1: Eliminating $P_{N - 1, 1} (0)$ .

Setting $s = δ_{N}$ in (15) yields $P_{N - 1, 1} (0)$ in terms of $P_{N, 1} (0)$ :

(22) $P_{N - 1, 1} (0) = ψ_{N - 1} P_{N, 1} (0), ψ_{N - 1} : = \frac{1 - A^{*} (δ_{N})}{A^{*} (δ_{N}) θ_{N - 1}} .$

Step 2: An expression for $P_{N, 1}^{*} (s)$ .

Substituting (22) back into (15) gives

(23) $(δ_{N} - s) P_{N, 1}^{*} (s) = (A^{*} (s) (θ_{N - 1} ψ_{N - 1} + ψ_{N}) - ψ_{N}) P_{N, 1} (0),$

where we set

ψ_{N} \equiv 1

. Hence, for

s \neq δ_{N}

(24) $P_{N, 1}^{*} (s) = \frac{A^{*} (s) (θ_{N - 1} ψ_{N - 1} + ψ_{N}) - ψ_{N}}{δ_{N} - s} P_{N, 1} (0) .$

Step 3: The value at $s = δ_{N}$ and higher derivatives.

Differentiating (23) with respect to s and then setting $s = δ_{N}$ yields

(25) $P_{N, 1}^{*} (δ_{N}) = - A^{* (1)} (δ_{N}) (θ_{N - 1} ψ_{N - 1} + ψ_{N}) P_{N, 1} (0) .$

Differentiating (23) a total of l times with respect to s gives, for

l \geq 1

(26) $(δ_{N} - s) P_{N, 1}^{* (l)} (s) - l P_{N, 1}^{* (l - 1)} (s) = A^{* (l)} (s) (θ_{N - 1} ψ_{N - 1} + ψ_{N}) P_{N, 1} (0) .$

Step 4: Compact representation.

Combining (24)–(26), we obtain

(27) $P_{N, 1}^{*} (s) = ζ_{N, s} P_{N, 1} (0),$

where

$ζ_{N, s} = \{\begin{matrix} \frac{A^{*} (s) (θ_{N - 1} ψ_{N - 1} + ψ_{N}) - ψ_{N}}{δ_{N} - s}, & s \neq δ_{N}, \\ - A^{* (1)} (s) (θ_{N - 1} ψ_{N - 1} + ψ_{N}), & s = δ_{N}, \end{matrix}$

and

$ζ_{N, s}^{(l)} = \{\begin{matrix} \frac{A^{* (l)} (s) (θ_{N - 1} ψ_{N - 1} + ψ_{N}) + l ζ_{N, s}^{(l - 1)}}{δ_{N} - s}, & s \neq δ_{N}, \\ - \frac{A^{* (l + 1)} (s) (θ_{N - 1} ψ_{N - 1} + ψ_{N})}{l + 1}, & s = δ_{N}, \end{matrix}$

with

ζ_{N, s}^{(l)}

denoting the lth derivative of

ζ_{N, s}

with respect to s. Equations (21)–(27) provide the desired bridge from arbitrary-epoch quantities to pre-arrival probabilities and furnish the recursive ingredients needed to express all required terms in

P_{N, 0} (0)

. For

i = N - 1

, taking

s = δ_{N - 1}

in Equation (14) and using Equation (22), we obtain

P_{N - 2, 1} (0)

in terms of

P_{N, 1} (0)

as follows:

(28) $P_{N - 2, 1} (0) = ψ_{N - 2} P_{N, 1} (0),$

with

$ψ_{N - 2} = \frac{ψ_{N - 1} - (δ_{N} - ϕ) ζ_{N, δ_{N - 1}} - A^{*} (δ_{N - 1}) {\bar{θ}}_{N - 1} ψ_{N - 1}}{A^{*} (δ_{N - 1}) θ_{N - 2}} .$

By applying Equation (28) to Equation (14) for

i = N - 1

, we can express

P_{N - 1, 1}^{*} (s)

in terms of

P_{N, 1} (0)

as follows:

(29) $P_{N - 1, 1}^{*} (s) = ζ_{N - 1, s} P_{N, 1} (0),$

where

$ζ_{N - 1, s} = \{\begin{matrix} \frac{A^{*} (s) (θ_{N - 2} ψ_{N - 2} + {\bar{θ}}_{N - 1} ψ_{N - 1}) + (δ_{N} - ϕ) ζ_{N, s} - ψ_{N - 1}}{(δ_{N - 1} - s)}, & if s \neq δ_{N - 1}, \\ - (A^{* (1)} (s) (θ_{N - 2} ψ_{N - 2} + {\bar{θ}}_{N - 1} ψ_{N - 1}) + (δ_{N} - ϕ) ζ_{N, s}^{(1)}), & if s = δ_{N - 1}, \end{matrix}$

where

$ζ_{N - 1, s}^{(l)} = \{\begin{matrix} \frac{A^{* (l)} (s) (θ_{N - 2} ψ_{N - 2} + {\bar{θ}}_{N - 1} ψ_{N - 1}) + (δ_{N} - ϕ) ζ_{N, s}^{(l)} + l ζ_{N - 1, s}^{(l - 1)}}{(δ_{N - 1} - s)}, & if s \neq δ_{N - 1}, \\ - \frac{A^{* (l + 1)} (s) (θ_{N - 2} ψ_{N - 2} + {\bar{θ}}_{N - 1} ψ_{N - 1}) + (δ_{N} - ϕ) ζ_{N, s}^{(l + 1)}}{l + 1}, & if s = δ_{N - 1} . \end{matrix}$

For

i = N - 2, N - 3, \dots, 3

in Equation (14),

P_{i, 1} (0),

for

2 \leq i \leq N - 3,

and

P_{i, 1}^{*} (s),

for

3 \leq i \leq N - 2,

are expressed in terms of

P_{N, 1} (0)

as follows:

(30) $P_{i - 1, 1} (0) = ψ_{i - 1} P_{N, 1} (0), i = N - 2, N - 3, \dots, 3,$

where

$ψ_{i - 1} = \frac{ψ_{i} - (δ_{i + 1} - ϕ) ζ_{i + 1, δ_{i}} - A^{*} (δ_{i}) {\bar{θ}}_{i} ψ_{i}}{A^{*} (δ_{i}) θ_{i - 1}} i = N - 2, N - 3, \dots, 3,$

and

(31) $P_{i, 1}^{*} (s) = ζ_{i, s} P_{N, 1} (0), i = N - 2, N - 3, \dots, 3,$

where

$ζ_{i, s} = \{\begin{matrix} \frac{A^{*} (s) (θ_{i - 1} ψ_{i - 1} + {\bar{θ}}_{i} ψ_{i}) + (δ_{i + 1} - ϕ) ζ_{i + 1, s} - ψ_{i}}{(δ_{i} - s)}, & if s \neq δ_{i}, \\ - (A^{* (1)} (s) (θ_{i - 1} ψ_{i - 1} + {\bar{θ}}_{i} ψ_{i}) + (δ_{i + 1} - ϕ) ζ_{i + 1, s}^{(1)}), & if s = δ_{i}, \end{matrix}$

where

$ζ_{i, s}^{(l)} = \{\begin{matrix} \frac{A^{* (l)} (s) (θ_{i - 1} ψ_{i - 1} + {\bar{θ}}_{i} ψ_{i}) + (δ_{i + 1} - ϕ) ζ_{i + 1, s}^{(l)} - l ζ_{i, s}^{(l - 1)}}{(δ_{i} - s)}, & if s \neq δ_{i}, \\ - \frac{(δ_{i + 1} - ϕ) ζ_{i + 1, δ_{i}}^{(l + 1)} + A^{* (l + 1)} (s) (θ_{i - 1} ψ_{i - 1} + {\bar{θ}}_{i} ψ_{i})}{l + 1}, & if s = δ_{i} . \end{matrix}$

Substituting

s = δ_{2}

into Equation (13) yields

P_{1, 1} (0)

in terms of

P_{N, 1} (0)

and

P_{1, 0} (0)

as follows:

(32) $P_{1, 1} (0) = ψ_{1} P_{N, 1} (0) + ω P_{1, 0} (0),$

where

$ψ_{1} = \frac{ψ_{2} - (δ_{3} - ϕ) ζ_{3, δ_{2}} - A^{*} (δ_{2}) {\bar{θ}}_{2} ψ_{2}}{A^{*} (δ_{2})}, and ω = - \frac{A^{*} (δ_{2})}{A^{*} (δ_{2})} = - 1 .$

Using Equation (32) in Equation (13), we obtain

P_{2, 1}^{*} (s)

in terms of

P_{N, 1} (0)

as follows:

(33) $P_{2, 1}^{*} (s) = ζ_{2, s} P_{N, 1} (0),$

where

$ζ_{2, s} = \{\begin{matrix} \frac{- ψ_{2} + (δ_{3} - ϕ) ζ_{3, s} + A^{*} (s) (ψ_{1} + {\bar{θ}}_{2} ψ_{2})}{(δ_{2} - s)}, & if s \neq δ_{2}, \\ - ((δ_{3} - ϕ) ζ_{3, s}^{(1)} + A^{* (1)} (s) (ψ_{1} + {\bar{b}}_{2} ψ_{2})), & if s = δ_{2}, \end{matrix}$

and

$ζ_{2, s}^{(l)} = \{\begin{matrix} \frac{(δ_{3} - ϕ) ζ_{3, s}^{(l)} + A^{* (l)} (s) (ψ_{1} + {\bar{θ}}_{2} ψ_{2}) - l ζ_{2, s}^{(l - 1)}}{(δ_{2} - s)}, & if s \neq δ_{2}, \\ - \frac{(δ_{3} - ϕ) ζ_{3, s}^{(l + 1)} + A^{* (l + 1)} (s) (ψ_{1} + {\bar{b}}_{2} ψ_{2})}{l + 1}, & if s = δ_{2} . \end{matrix}$

From Equation (12), we have

(34) $P_{1, 1}^{*} (s) = ζ_{1, s} P_{N, 1} (0) + τ_{1, s} P_{1, 0} (0),$

where

$\begin{matrix} ζ_{1, s} & = & \{\begin{matrix} \frac{β μ_{1} ζ_{2, s} - ψ_{1}}{(ϕ + β ν - s)}, & if s \neq ϕ + β ν, \\ - β μ_{1} ζ_{2, s}^{(1)}, & if s = ϕ + β ν . \end{matrix} \\ τ_{1, s} & = & \{\begin{matrix} - \frac{ω}{(ϕ + β ν - s)}, & if s \neq ϕ + β ν, \\ 0, & if s = ϕ + β ν, \end{matrix} \end{matrix}$

where

$\begin{matrix} ζ_{1, s}^{(l)} & = & \{\begin{matrix} \frac{β μ_{1} ζ_{2, s}^{(l)} - l ζ_{1, s}^{(l - 1)}}{(ϕ + β ν - s)}, & if s \neq ϕ + β ν, \\ - \frac{β μ_{1} ζ_{2, s}^{(l + 1)}}{l + 1}, & if s = ϕ + β ν . \end{matrix} \\ τ_{1, s}^{(l)} & = & \{\begin{matrix} \frac{l τ_{1, s}^{(l - 1)}}{(ϕ + β ν - s)}, & if s \neq ϕ + β ν, \\ 0, & if s = ϕ + β ν . \end{matrix} \end{matrix}$

Putting

s = ϑ_{N}

in Equation (18) and using the relation of

P_{N, 1}^{*} (s)

, we obtain

P_{N - 1, 2} (0)

in terms of

P_{N, 1} (0)

and

P_{N, 2} (0)

as follows:

(35) $P_{N - 1, 2} (0) = η_{N - 1} P_{N, 2} (0) + γ_{N - 1} P_{N, 1} (0),$

where

$η_{N - 1} = \frac{1 - A^{*} (ϑ_{N})}{A^{*} (ϑ_{N}) θ_{N - 1}} and γ_{N - 1} = - \frac{ϕ ζ_{N, ϑ_{N}}}{A^{*} (ϑ_{N}) θ_{N - 1}} .$

Using Equation (35) in Equation (18), we get

(36) $P_{N, 2}^{*} (s) = ρ_{N, s} P_{N, 2} (0) + χ_{N, s} P_{N, 1} (0),$

where

$ρ_{N, s} = \{\begin{matrix} \frac{A^{*} (s) (θ_{N - 1} η_{N - 1} + η_{N}) - η_{N}}{ϑ_{N} - s}, & if ϑ_{N} \neq s, \\ - A^{* (1)} (s) (θ_{N - 1} η_{N - 1} + η_{N}), & if ϑ_{N} = s . \end{matrix}$

$χ_{N, s} = \{\begin{matrix} \frac{ϕ ζ_{N, s} + A^{*} (s) θ_{N - 1} γ_{N - 1}}{ϑ_{N} - s}, & if ϑ_{N} \neq s, \\ - (ϕ ζ_{N, s}^{(1)} + A^{* (1)} (s) θ_{N - 1} γ_{N - 1}), & if ϑ_{N} = s, \end{matrix}$

where

$ρ_{N, s}^{(l)} = \{\begin{matrix} \frac{A^{* (l)} (s) (θ_{N - 1} η_{N - 1} + η_{N}) + l ρ_{N, s}^{(l - 1)}}{(ϑ_{N} - s)}, & if ϑ_{N} \neq s, \\ - \frac{A^{* (l + 1)} (ϑ_{N}) (θ_{N - 1} η_{N - 1} + η_{N})}{l + 1}, & if ϑ_{N} = s . \end{matrix}$

$χ_{N, s}^{(l)} = \{\begin{matrix} \frac{ϕ ζ_{N, s}^{(l)} + A^{* (l)} (s) θ_{N - 1} γ_{N - 1} + l χ_{N, s}^{(l - 1)}}{ϑ_{N} - s}, & if ϑ_{N} \neq s, \\ - \frac{ϕ ζ_{N, ϑ_{N}}^{(l + 1)} + A^{* (l + 1)} (ϑ_{N}) θ_{N - 1} γ_{N - 1}}{l + 1}, & if ϑ_{N} = s, \end{matrix}$

where

η_{N} = 1

and

γ_{N} = 0

. Following the procedure carried out for obtaining

P_{i, 1} (0)

and

P_{i, 1}^{*} (s)

, the probabilities

P_{i, 2} (0),

1 \leq i \leq N - 2,

and

P_{i, 2}^{*} (s),

2 \leq i \leq N - 1,

can be obtained from Equation (17) in terms of

P_{N, 1} (0)

and

P_{N, 2} (0)

(37) $P_{i - 1, 2} (0) = η_{i - 1} P_{N, 2} (0) + γ_{i - 1} P_{N, 1} (0), 2 \leq i \leq N - 1,$

where

$\begin{matrix} η_{i - 1} & = & \frac{η_{i} - ϑ_{i + 1} ρ_{i + 1, ϑ_{i}} - A^{*} (ϑ_{i}) {\bar{θ}}_{i} η_{i}}{A^{*} (ϑ_{i}) θ_{i - 1}}, \\ γ_{i - 1} & = & \frac{γ_{i} - ϑ_{i + 1} χ_{i + 1, ϑ_{i}} - A^{*} (ϑ_{i}) {\bar{θ}}_{i} γ_{i} - ϕ ζ_{i, ϑ_{i}}}{A^{*} (ϑ_{i}) θ_{i - 1}} . \end{matrix}$

By using Equation (37) in Equation (17), we get

(38) $P_{i, 2}^{*} (s) = ρ_{i, s} P_{N, 2} (0) + χ_{i, s} P_{N, 1} (0),$

where

$\begin{matrix} ρ_{i, s} & = & \{\begin{matrix} \frac{- η_{i} + ϑ_{i + 1} ρ_{i + 1, s} + A^{*} (s) (θ_{i - 1} η_{i - 1} + {\bar{θ}}_{i} η_{i})}{(ϑ_{i} - s)}, & if ϑ_{i} \neq s, \\ - (ϑ_{i + 1} ρ_{i + 1, s}^{(1)} + A^{* (1)} (s) (θ_{i - 1} η_{i - 1} + {\bar{θ}}_{i} η_{i})), & if ϑ_{i} = s . \end{matrix} \\ χ_{i, s} & = & \{\begin{matrix} \frac{- γ_{i} + ϕ ζ_{i, s} + ϑ_{i + 1} χ_{i + 1, s} + A^{*} (s) ((θ_{i - 1} γ_{i - 1} + {\bar{θ}}_{i} γ_{i}))}{ϑ_{i} - s}, & if ϑ_{i} \neq s, \\ - (ϕ ζ_{i, s}^{(1)} + ϑ_{i + 1} χ_{i + 1, s}^{(1)} + A^{* (1)} (s) ((θ_{i - 1} γ_{i - 1} + {\bar{θ}}_{i} γ_{i}))), & if ϑ_{i} = s, \end{matrix} \end{matrix}$

where

$\begin{matrix} ρ_{i, s}^{(l)} & = & \{\begin{matrix} \frac{ϑ_{i + 1} ρ_{i + 1, s}^{(l)} + A^{* (l)} (s) (θ_{i - 1} η_{i - 1} + {\bar{θ}}_{i} η_{i}) + l ρ_{i, s}^{(l - 1)}}{(ϑ_{i} - s)}, & if ϑ_{i} \neq s, \\ - \frac{ϑ_{i + 1} ρ_{i + 1, s}^{(l + 1)} + A^{* (l + 1)} (s) (θ_{i - 1} η_{i - 1} + {\bar{θ}}_{i} η_{i})}{l + 1}, & if ϑ_{i} = s . \end{matrix} \end{matrix}$

$\begin{matrix} χ_{i, s}^{(l)} & = & \{\begin{matrix} \frac{ϕ ζ_{i, s}^{(l)} + ϑ_{i + 1} χ_{i + 1, s}^{(l)} + A^{* (l)} (s) (θ_{i - 1} γ_{i - 1} + {\bar{θ}}_{i} γ_{i}) + l χ_{i, s}^{(l - 1)}}{ϑ_{i} - s}, & if ϑ_{i} \neq s, \\ - \frac{ϕ ζ_{i, s}^{(l + 1)} + ϑ_{i + 1} χ_{i + 1, s}^{(l + 1)} + A^{* (l + 1)} (s) (θ_{i - 1} γ_{i - 1} + {\bar{θ}}_{i} γ_{i})}{l + 1}, & if ϑ_{i} = s . \end{matrix} \end{matrix}$

Taking

s = β μ_{1}

in Equation (11), we get

P_{0, 0} (0)

in terms of

P_{1, 0} (0)

P_{N, 1} (0)

, and

P_{N, 2} (0)

(39) $P_{0, 0} (0) = ▵_{0} P_{1, 0} (0) + ▵_{1} P_{N, 1} (0) + ▵_{2} P_{N, 2} (0),$

where

$▵_{0} = \frac{1}{A^{*} (β μ_{1})}, ▵_{1} = \frac{- β ν ζ_{2, β μ_{1}} - β μ_{2} χ_{2, β μ_{1}}}{A^{*} (β μ_{1})}, ▵_{2} = \frac{- β μ_{2} ρ_{2, β μ_{1}}}{A^{*} (β μ_{1})} .$

Taking

s = ϕ + β ν

in Equation (12), using the expressions of

P_{1, 1} (0)

and

P_{2, 1}^{*} (s),

we get

P_{1, 0} (0)

in terms of

P_{N, 1} (0)

as follows:

(40) $P_{1, 0} (0) = κ_{1} P_{N, 1} (0),$

where

$κ_{1} = ψ_{1} - β μ_{1} ζ_{2, ϕ + β ν} .$

Taking

β μ_{2} = s

into Equation (16) and using the expressions for

P_{1, 2} (0),

P_{1, 1}^{*} (s)

and

P_{2, 2}^{*} (s),

along with Equation (40), we obtain

P_{N, 2} (0)

in terms of

P_{N, 1} (0)

as follows:

(41) $P_{N, 2} (0) = κ_{2} P_{N, 1} (0),$

where

$κ_{2} = \frac{ϕ ζ_{1, β μ_{2}} + β μ_{1} χ_{2, β μ_{2}} + ϕ κ_{1} τ_{1, β μ_{2}} - γ_{1}}{η_{1} - β μ_{1} ρ_{2, β μ_{2}}} .$

From Equations (22), (28) and (30), and by applying Equation (40) into Equation (32), Equation (41) into Equations (35) and (37), and Equations (40) and (41) into Equation (38), we can express all

P_{i, j} (0)

in terms of

P_{N, 1} (0) .

The only unknown

P_{N, 1} (0)

can be determined from Equation (21).

$P_{N, 1} (0) = λ {(κ_{1} ▵_{0} + ▵_{1} + κ_{2} ▵_{2} + ψ_{1} + \sum_{i = 2}^{N} ψ_{i} + \sum_{i = 2}^{N} (γ_{i} + κ_{2} η_{i}))}^{- 1} .$

Now, the steady-state probabilities at a pre-arrival epoch

P_{i, j}^{-}

can be simply obtained from the rate probabilities

P_{i, j} (0)

through Equation (21). Finally, we formulate the steady-state probabilities at an arbitrary epoch in terms of the corresponding probabilities at a pre-arrival epoch. In particular, setting

s = 0

in Equations (10)–(18), and subsequently invoking relation (21) together with the requisite algebraic reductions, yields

$\begin{matrix} P_{N, 1} = \frac{λ}{δ_{N}} θ_{N - 1} P_{N - 1, 1}^{-}, \\ P_{i, 1} = (\frac{δ_{i + 1} - ϕ}{δ_{i}}) P_{i + 1, 1} + \frac{λ}{δ_{i}} (θ_{i - 1} P_{i - 1, 1}^{-} - θ_{i} P_{i, 1}^{-}), i = N - 1, \dots, 3, \\ P_{2, 1} = (\frac{δ_{3} - ϕ}{δ_{2}}) P_{3, 1} + \frac{λ}{δ_{2}} (P_{1, 0}^{-} + P_{1, 1}^{-} - θ_{2} P_{2, 1}^{-}), \end{matrix}$

$\begin{matrix} P_{1, 1} = (\frac{β μ_{1}}{ϕ + β ν}) P_{2, 1} - (\frac{λ}{ϕ + β ν}) P_{1, 1}^{-}, \\ P_{N, 2} = \frac{ϕ}{ϑ_{N}} P_{N, 1} + \frac{λ}{ϑ_{N}} θ_{N - 1} P_{N - 1, 2}^{-}, \\ P_{i, 2} = (\frac{ϑ_{i + 1}}{ϑ_{i}}) P_{i + 1, 2} + \frac{ϕ}{ϑ_{i}} P_{i, 1} + \frac{λ}{ϑ_{i}} (θ_{i - 1} P_{i - 1, 2}^{-} - θ_{i} P_{i, 2}^{-}), i = N - 1, \dots, 2, \\ P_{1, 2} = \frac{μ_{1}}{μ_{2}} P_{2, 2} + \frac{ϕ}{β μ_{2}} P_{1, 1} - \frac{λ}{β μ_{2}} P_{1, 2}^{-}, \\ P_{1, 0} = \frac{ν}{μ_{1}} P_{2, 1} + \frac{μ_{2}}{μ_{1}} P_{2, 2} + \frac{λ}{β μ_{1}} (P_{0, 0}^{-} - P_{1, 0}^{-}), \end{matrix}$

where

P_{0, 0}

can be computed by using the normalization condition, that is,

$P_{0, 0} = 1 - P_{1, 0} - \sum_{i = 1}^{N} (P_{i, 1} + P_{i, 2}) .$

4. Performance Measures

To guide the design and optimization of the queuing system, we evaluate standard performance indices computed from the steady-state probabilities at arbitrary epochs obtained in the previous section. Throughout, let $λ$ denote the external arrival rate, N the system capacity, $α$ the per-customer reneging rate when applicable, and $θ_{i} \in [0, 1]$ the probability that an arriving customer decides to join upon seeing i customers in the system (so that ${\bar{θ}}_{i} : = 1 - θ_{i}$ is the balking probability). We write $P_{i, 0}$ for the probability of being in a state with i customers while the server is idle on a working vacation, $P_{i, 1}$ for a state with i customers and the server in a (reduced-rate) working–vacation service mode, and $P_{i, 2}$ for a state with i customers and the server in the regular busy mode. The loss probability upon arrival is denoted by

$Π_{L} : = P_{N, 1}^{-} + P_{N, 2}^{-},$

i.e., the probability that an arriving customer finds the system full (pre-arrival probabilities).

(i). Mean system size and mean sojourn time.

By definition and Little’s law applied to the effective throughput $λ (1 - Π_{L})$ , we have

(42) $\begin{matrix} L_{s} & = P_{1, 0} + \sum_{i = 1}^{N} i (P_{i, 1} + P_{i, 2}), \end{matrix}$

(43) $\begin{matrix} W_{s} & = \frac{L_{s}}{λ (1 - Π_{L})} . \end{matrix}$

(ii). Server-state occupancy probabilities.

The long-run fractions of time that the server is (a) idle during a working vacation, (b) busy during a working vacation, and (c) busy during a normal busy period are

(44) $Π_{I} = \sum_{i = 0}^{1} P_{i, 0}, Π_{W} = \sum_{i = 1}^{N} P_{i, 1}, Π_{B} = \sum_{i = 1}^{N} P_{i, 2} .$

(iii). Flow rates of joining, balking, and reneging.

Let $J_{s}$ denote the average rate of customers that actually enter service (i.e., join the system), $B_{r}$ the average balking rate, and $R_{r}$ the average reneging rate. Then,

(45) $\begin{matrix} J_{s} & = λ (P_{0, 0} + P_{1, 0} + P_{1, 1} + P_{1, 2}) + \sum_{i = 2}^{N} λ θ_{i} (P_{i, 1} + P_{i, 2}), \end{matrix}$

(46) $\begin{matrix} B_{r} & = \sum_{i = 2}^{N} λ {\bar{θ}}_{i} (P_{i, 1} + P_{i, 2}), \end{matrix}$

(47) $\begin{matrix} R_{r} & = \sum_{i = 2}^{N} (i - 2) α (P_{i, 1} + P_{i, 2}) . \end{matrix}$

The indices (42)–(43) together with $(Π_{I}, Π_{W}, Π_{B})$ and $(J_{s}, B_{r}, R_{r})$ provide a coherent basis for performance evaluation and, in particular, for selecting service-rate parameters that balance congestion, abandonment, and utilization under capacity constraints.

Remark 1.

Section 4 reports deterministic evaluations of the stationary distribution and associated indices via the exact recursions and transform-based formulae developed earlier; hence, no stochastic simulation or hardware-dependent averaging was involved. For completeness, we will specify the computational environment (standard workstation $(O (N))$ complexity per evaluation) and provide a clear rationale for input choices: service and vacation rates were selected to satisfy the structural ordering $(ν < μ_{2} < μ_{1})$ ; balking and reneging rates were set within ranges consistent with classical impatience models; and buffer sizes were varied to capture light, moderate, and heavy traffic regimes while ensuring stability.

5. Cost Model and Optimization Study

Building on the performance indices derived in the previous section, we now formulate the steady-state expected total cost per unit time for the proposed queuing system. The decision vector comprises three continuous, nonnegative service/vacation rates,

$(μ_{1}, μ_{2}, ν) \in R_{+}^{3},$

where

μ_{1}

and

μ_{2}

are the normal service rates of the two servers (with Server 1 faster than Server 2), and

ν

is the working-vacation service rate of Server 2. In practice, queuing managers aim to reduce operating costs subject to stability and operational constraints. Our objective is to determine optimal values

$(μ_{1}^{*}, μ_{2}^{*}, ν^{*}) \in arg min_{(μ_{1}, μ_{2}, ν)} F (μ_{1}, μ_{2}, ν),$

where F denotes the steady-state cost function defined below.

Unit cost elements.

Let $C_{1}, \dots, C_{9} \geq 0$ denote the following unit costs: $C_{1}$ :

Cost per unit time when Server 2 is idle during a working-vacation period;

$C_{2}$ :

Cost per unit time when Server 2 is busy during a working-vacation period;

$C_{3}$ :

Cost per unit time when Server 2 is busy during a normal (non-vacation) busy period;

$C_{4}$ :

Holding cost per unit time per customer present in the system;

$C_{5}$ :

Penalty cost per unit time due to balking or reneging;

$C_{6}$ :

Penalty cost per unit time per lost customer when the system is blocked;

$C_{7}$ :

Cost per unit of normal service effort;

$C_{8}$ :

Cost per unit of feedback service effort;

$C_{9}$ :

Fixed purchase (or capacity) cost per server unit.

Performance indices (given).

Let $Π_{I}, Π_{W}, Π_{B}, Π_{L} \in [0, 1]$ denote, respectively, the steady-state probabilities that Server 2 is idle on working vacation, and busy on working vacation, busy in normal operation, and the system is blocked; let $L_{s} \geq 0$ be the steady-state expected system size; let $B_{r}$ and $R_{r}$ be the steady-state rates of balking and reneging, respectively; let $λ > 0$ be the external arrival rate; and let $\bar{β} \in [0, 1]$ denote the mean feedback probability. These indices are determined by the stochastic model specified in the previous section.

Total expected cost.

The steady-state expected total cost per unit time is

(48) $\begin{matrix} F (μ_{1}, μ_{2}, ν) & = C_{1} Π_{I} + C_{2} Π_{W} + C_{3} Π_{B} + C_{4} L_{s} + C_{5} (B_{r} + R_{r}) \\ + C_{6} λ Π_{L} + (μ_{1} + μ_{2} + ν) (C_{7} + \bar{β} C_{8}) + 2 C_{9} . \end{matrix}$

Optimization problem.

We seek the cost-minimizing rates under the natural ordering and feasibility (stability) constraints:

(49) $\begin{matrix} (μ_{1}^{*}, μ_{2}^{*}, ν^{*}) \in & arg min_{(μ_{1}, μ_{2}, ν)} F (μ_{1}, μ_{2}, ν) \\ subject to & 0 < ν < μ_{2} < μ_{1}, \\ and the steady - state / stability conditions of the model . \end{matrix}$

Solution approach.

Because F in (48) is a highly nonlinear function of $(μ_{1}, μ_{2}, ν)$ through the embedded performance indices, closed-form optimality conditions are generally intractable. We therefore adopt a numerical approach based on the Bat algorithm, a metaheuristic inspired by echolocation, to efficiently explore the feasible region in (49) and obtain high-quality approximations to $(μ_{1}^{*}, μ_{2}^{*}, ν^{*})$ .

Comments

The selection of the decision variables $(μ_{1}, μ_{2}, ν)$ and the cost components $(C_{1}, \dots, C_{9})$ is motivated by both modeling fidelity and practical relevance. We briefly elaborate on their roles.

Service and vacation rates.

The triplet $(μ_{1}, μ_{2}, ν)$ forms the natural decision vector because service rates constitute the most direct levers for system performance. In particular,

$μ_{1}$ and $μ_{2}$ represent the normal operating capacities of two heterogeneous servers. Empirically, such asymmetry is common in practice, where one server is technologically superior or operated by a more skilled resource. Distinguishing them allows us to capture realistic differences in throughput and workload distribution.

$ν$ models the effective rate during a working-vacation period, a concept increasingly relevant in systems where partial service continues during maintenance, energy-saving modes, or off-peak operation. Including $ν$ in the optimization enables one to quantify the trade-off between reduced productivity during vacations and the associated cost savings.

Cost components.

The decomposition of the total cost into nine elements $(C_{1}, \dots, C_{9})$ ensures that all critical operational aspects are explicitly represented:

$C_{1}$ – $C_{3}$ differentiate between idle and busy states of Server 2 across normal and vacation regimes, thereby capturing the utilization-dependent expenditure structure.

$C_{4}$ – $C_{6}$ encode congestion and dissatisfaction costs: holding cost for customers in queue, penalties for impatience (balking/reneging), and loss penalties under blocking. This reflects the service quality dimension and its direct impact on customer retention.

$C_{7}$ – $C_{8}$ measure the marginal effort of normal versus feedback service, thereby penalizing overuse of resources and accounting for the additional burden of reprocessing tasks.

$C_{9}$ accounts for fixed capacity investments, ensuring that expansion or contraction of server capability is consistently weighed against operational benefits.

Modeling rationale.

This parameterization balances analytical tractability with operational richness. By separating performance-determining rates from cost-incurring states, the model provides a transparent mapping from managerial decisions to economic outcomes. The chosen structure also aligns with the literature on cost optimization in heterogeneous and vacation-based queuing systems, ensuring comparability and theoretical robustness.

Remark 2.

The proposed framework can, in principle, be extended to accommodate multi-class customer populations or priority-based service policies. In the multi-class setting, the state descriptor must be augmented to reflect both the number and the class of customers in service or waiting, leading to vector-valued probability generating functions and Laplace–Stieltjes transforms, together with class-dependent recursive relations. Similarly, under priority-based disciplines (preemptive or non-preemptive), the balance equations and boundary conditions must be reformulated to incorporate class-specific service rules, which naturally give rise to block-structured transition dynamics. Analytically, such extensions often require the use of matrix-analytic or quasi-birth–death techniques to preserve tractability. While feasible, these generalizations substantially increase the model’s dimensional and computational complexity, and their rigorous treatment lies beyond the scope of the present work.

6. Numerical Results

In this section, we provide a detailed numerical investigation of the multiple working vacation $G I / M / 2 / N$ queuing model. The inter-arrival process is examined under a range of canonical distributions, namely deterministic, exponential, and Erlang-2, so as to capture a spectrum of arrival-time variability. The numerical experiments, presented through carefully constructed graphs and tables, are intended to elucidate the influence of key queuing parameters on both the principal performance measures and the associated cost-optimization criteria. To this end, we have developed an R program, (R version 4.5.2), implemented by the authors, which computationally validates and illustrates the analytical formulas derived in the preceding sections. For the entire numerical analysis, representative values of the system parameters and cost coefficients are chosen in an illustrative yet systematic manner, ensuring that the essential qualitative behaviors of the model are faithfully exhibited.

In Table 2, we report the steady-state joint distributions observed at both pre-arrival and arbitrary epochs for three representative inter-arrival time laws—deterministic, exponential, and Erlang–2—all calibrated with rate parameter $λ = 0.9$ . The remaining system parameters are specified as $μ_{1} = 2.7$ , $μ_{2} = 2.3$ , $ν = 1.9$ , $\bar{β} = 0.2$ , $α = 0.6$ , and $ϕ = 1.0$ . Customer balking is modeled through the state-dependent function ${\bar{θ}}_{i} = i / N$ for $2 \leq i \leq N - 1$ , where the total system capacity is fixed at $N = 9$ . A noteworthy feature emerges in the exponential case: because of the memoryless property of the exponential distribution, the pre-arrival probabilities $P_{i, j}^{-}$ coincide exactly with the arbitrary-epoch probabilities $P_{i, j}$ , reflecting the fundamental Markovian nature of this setting.

6.1. Numerical Results of System Performance Measures

Table 3 and Figure 2, Figure 3 and Figure 4 report a sensitivity analysis of the principal performance measures with respect to the service rates during normal and working-vacation regimes, namely $μ_{1}$ , $μ_{2}$ , and $ν$ . Throughout, inter-arrival times are taken to be deterministic with mean $1 / λ$ . Unless otherwise stated, the baseline parameter vector is

$(λ, μ_{1}, μ_{2}, ν, \bar{β}, α, ϕ) = (1.2, 3.3, 2.5, 1.6, 0.2, 1.1, 5.0),$

and the balking function is

${\bar{θ}}_{i} = 1 - \frac{1}{i + 1}, 2 \leq i \leq N - 1, N = 7 .$

Main findings. The numerical results show that increasing any of the service rates $μ_{1}$ , $μ_{2}$ , or $ν$ leads to a monotone decrease in the congestion and impatience indicators:

$L_{s} ↓, W_{s} ↓, B_{r} ↓, R_{r} ↓,$

where

L_{s}

is the mean system size,

W_{s}

the mean sojourn time,

B_{r}

the average balking rate, and

R_{r}

the average reneging rate. In parallel, the effective joining rate

J_{s}

increases with faster service, reflecting improved responsiveness and reduced anticipated delay.

Regarding server-state utilization, the probability that Server 2 is busy during normal (non-vacation) periods, denoted $Π_{B}$ , decreases as each of $μ_{1}$ , $μ_{2}$ , and $ν$ increases. During working vacations, the busy probability of Server 2, denoted $Π_{W}$ , decreases with $μ_{1}$ and with $ν$ , but increases with $μ_{2}$ . This asymmetric response is consistent with the heterogeneity of the two servers and the division of work between normal and vacation regimes.

Effect of balking. Comparisons between systems with and without balking reveal that

$L_{s}, W_{s}, J_{s}, R_{r}, Π_{B} are all larger in the no-balking system,$

while

$B_{r} (trivially) and Π_{W} are larger when balking is present .$

These patterns are qualitatively in line with queuing intuition: allowing balking filters out customers facing long anticipated waits, thereby lowering congestion measures and shifting occupancy from normal periods to working vacations.

Secondly, in Table 4 and Figure 5, Figure 6 and Figure 7, we analyze the impact of feedback probability $\bar{β},$ arrival rate $λ,$ working vacation rate $ϕ$ , and impatience rate $α$ on various performance measures. Here, we assume that the inter-arrival times follow an Erlang-2 distribution. The balking function is taken as ${\bar{θ}}_{i} = i / N^{2},$ for $2 \leq i \leq N - 1$ , where $N = 7 .$ The parameters are taken as $λ = 1.0$ , $μ_{1} = 2.8$ , $μ_{2} = 2.4$ , $ν = 1.8$ , $α = 0.9$ , and $ϕ = 2.0 .$ From Table 4 and Figure 5, Figure 6 and Figure 7, we observe the following:

With fixed values of $\bar{β}$ , $ϕ$ , and $α,$ an increase in $λ$ leads to increases in $L_{s}$ and $W_{s},$ as expected. This in turn increases $J_{s}$ , $B_{r}$ , and $R_{r} .$ Consequently, the probability of customer loss due to system size limitation $Π_{L}$ significantly rises. Notably, the probability that Server 2 is idle during the working vacation period $Π_{I}$ decreases with $λ .$

With fixed values of $\bar{β}$ , $λ,$ and $α,$ when the working vacation rate $ϕ$ increases, Server 2 rapidly switches to the normal busy period at which the customers are served at a higher rate. Therefore, the characteristics $L_{s}$ , $R_{r}$ , $B_{r}$ , $W_{s}$ , and $Π_{L}$ decrease. This implies an increase in $J_{s},$ and in $Π_{I}$ because the mean working vacation time $1 / ϕ$ decreases. This trend matches absolutely with the realistic situation.

With fixed values of $\bar{β}$ , $λ,$ and $ϕ,$ when the impatience rate $α$ increases, the system characteristics $J_{s}$ , $R_{r}$ , and $Π_{I}$ increase, while $Π_{L}$ , $B_{r}$ , $L_{s}$ , and $W_{s}$ decrease. This relationship highlights that higher impatience rates lead to smaller system sizes on average and higher average reneging rates.

For fixed values of $α$ , $ϕ$ , and $λ$ , an increase in the feedback probability $\bar{β}$ leads to higher values of $L_{s}$ and $W_{s},$ as intuitively expected. Consequently, this increase results in elevated values of $B_{r}$ , $R_{r}$ , and the probability of customer loss due to system size limitation $Π_{L}$ . Conversely, $J_{s}$ and the probability that Server 2 is idle during the working vacation period $Π_{I}$ monotonically decrease.

Table 4

Effects of $λ$ , $ϕ$ , $α$ , and $\bar{β}$ on system characteristics.

	$\bar{β} = 0.3$		$\bar{β} = 0.5$		$\bar{β} = 0.7$
	$λ = 0.7$	$λ = 1.0$	$λ = 0.7$	$λ = 1.0$	$λ = 0.7$	$λ = 1.0$
$L_{s}$	0.3649478	0.5285393	0.5156073	0.7485630	0.8706453	1.2585121
$W_{s}$	0.5213540	0.5285393	0.7365818	0.7485630	1.2437789	1.2585122
$J_{s}$	0.6995574	0.9984455	0.6989834	0.9965871	0.6967290	0.9899939
$B_{r}$	0.0004426	0.0015545	0.0010166	0.0034129	0.0032710	0.0100061
$R_{r}$	0.0003104	0.0013975	0.0011433	0.0047919	0.0068955	0.0251020
	$\bar{β} = 0.3$		$\bar{β} = 0.5$		$\bar{β} = 0.7$
	$ϕ = 2.0$	$ϕ = 4.0$	$ϕ = 2.0$	$ϕ = 4.0$	$ϕ = 2.0$	$ϕ = 4.0$
$L_{s}$	0.5285393	0.5253591	0.7485630	0.7430044	1.2585121	1.2486714
$W_{s}$	0.5285393	0.5253591	0.7485630	0.7430044	1.2585122	1.2486714
$J_{s}$	0.9984455	0.9984870	0.9965871	0.9966743	0.9899939	0.9901889
$B_{r}$	0.0015545	0.0015130	0.0034129	0.0033257	0.0100061	0.0098111
$R_{r}$	0.0013975	0.0013191	0.0047919	0.0045670	0.0251020	0.0243441
	$\bar{β} = 0.3$		$\bar{β} = 0.5$		$\bar{β} = 0.7$
	$α = 0.6$	$α = 0.9$	$α = 0.6$	$α = 0.9$	$α = 0.6$	$α = 0.9$
$L_{s}$	0.5288857	0.5285393	0.7501344	0.7485630	1.2703945	1.2585121
$W_{s}$	0.5288857	0.5285393	0.7501344	0.7485630	1.2703948	1.2585122
$J_{s}$	0.9984379	0.9984455	0.9965519	0.9965871	0.9897164	0.9899939
$B_{r}$	0.0015621	0.0015545	0.0034481	0.0034129	0.0102836	0.0100061
$R_{r}$	0.0010108	0.0013975	0.0035746	0.0047919	0.0199789	0.0251020

Figure 5

Effects of $λ,$ $\bar{β}$ , and $ϕ$ on $Π_{L}$ and $Π_{I}$ .

[Figure omitted. See PDF]

Figure 6

Effects of $ϕ,$ $\bar{β}$ , and $α$ on $Π_{L}$ and $Π_{I}$ .

[Figure omitted. See PDF]

Figure 7

Effects of $α,$ $\bar{β}$ , and $λ$ on $Π_{L}$ and $Π_{I}$ .

[Figure omitted. See PDF]

6.2. Numerical Results of Cost Optimization

Our aim in this subsection is to determine the optimal service rates $μ_{1}$ , $μ_{2}$ , and $ν$ that minimize the total expected cost function $F (μ_{1}, μ_{2}, ν)$ under various operational scenarios. The unit cost coefficients are prescribed as

$C_{1} = 30, C_{2} = 20, C_{3} = 10, C_{4} = 15, C_{5} = 12, C_{6} = 10, C_{7} = 10, C_{8} = 4, C_{9} = 5 .$

Customer balking is modeled by the quadratic function

${\bar{θ}}_{i} = \frac{i^{2}}{N^{2}}, 2 \leq i \leq N - 1,$

with the system capacity fixed at

N = 5

. Unless otherwise specified, the parameters

$\bar{β} = 0.3, ϕ = 0.2, λ = 1.8, μ_{2} = 0.9, ν = 0.6$

remain constant throughout the analysis. For the numerical optimization, we implement the Bat Algorithm, a population-based metaheuristic, with a population size of

P = 20

and a maximum number of iterations

T_{max} = 30

In Figure 8 and Figure 9, the convexity and optimality of the total expected cost function $F (μ_{1}, μ_{2}, ν)$ are clearly observed. These figures illustrate that there exist specific values of the service rates $μ_{1},$ $μ_{2}$ , and $ν$ that minimize the total expected cost function for the selected model parameters.

According to Table 5, as the arrival rate $λ$ increases, we observe that the optimal service rates at the primary and secondary stations, $μ_{1}^{*}$ and $μ_{2}^{*}$ , increase, whereas the auxiliary-rate $ν^{*}$ decreases. Consequently, the optimal expected cost

$F (μ_{1}^{*}, μ_{2}^{*}, ν^{*})$

is strictly increasing in

λ

. This behavior is natural: higher inflow intensifies congestion, enlarging the expected system size and—holding the optimality conditions fixed—raising the total expected cost. By contrast, when the impatience (reneging) rate

α

increases, all three optimal rates

μ_{1}^{*}

μ_{2}^{*}

, and

ν^{*}

rise, and so does the objective

F (μ_{1}^{*}, μ_{2}^{*}, ν^{*})

. This accords with intuition: greater impatience penalizes delay more heavily, incentivizing faster service to mitigate reneging, but the higher service intensities themselves entail larger operating costs, thereby elevating the overall optimum. Moreover, the presence of balking amplifies these effects. In models with balking, the optimal rates

μ_{1}^{*}

μ_{2}^{*}

, and

ν^{*}

, together with the resulting optimal cost

F (μ_{1}^{*}, μ_{2}^{*}, ν^{*})

, exceed their counterparts in otherwise identical systems without balking. Thus, classical impatience behaviors—balking and reneging—unambiguously worsen total expected cost by compelling more aggressive (and costlier) service policies. The comparative statics reported in Table 5 are consistent with operational practice and provide actionable guidance for real-time implementations: appropriate tuning of the service rates

μ_{1}

μ_{2}

, and

ν

can hedge against higher arrival intensities and customer impatience, while explicitly accounting for balking to avoid systematic under-provision of capacity.

6.3. Discussions

Remark 3.

On the compact feasible set ${ν_{min} \leq ν \leq μ_{2} \leq μ_{1} \leq μ_{1, max}}$ , the mapping $(μ_{1}, μ_{2}, ν) \mapsto F$ is continuous (indeed locally Lipschitz) since it is a smooth composition of linear recursions in the stationary probabilities. Hence, a global minimizer exists. The terms $C_{1} Π_{I} + C_{2} Π_{W} + C_{3} Π_{B}$ price server phases (idle/WV/regular), $C_{4} L_{s}$ penalizes congestion, $C_{5} (B_{r} + R_{r})$ penalizes behaviorally lost work, $C_{6} λ Π_{L}$ accounts for blocked demand, and $(C_{7} + \bar{β} C_{8}) (μ_{1} + μ_{2} + ν)$ proxies energy/labor with rework overhead.

We initialize bat uniformly over the feasible cone $ν \leq μ_{2} \leq μ_{1}$ , project after each move to enforce ordering, and use logarithmic step sizes for rates to improve exploration across scales. A small amount of Gaussian jitter on the best incumbent accelerates escape from flat basins. We report the best feasible incumbent after T iterations and confirm robustness through five independent runs.

Remark 4.

We adapted the Bat algorithm in three essential ways to reflect the mathematical structure of our queuing system. First, the decision variables are constrained by the intrinsic service-rate ordering

$0 < ν < μ_{2} < μ_{1},$

and rather than introducing penalty terms, every candidate solution is projected directly onto this feasible cone (after proposing moves in log-rate space to preserve positivity and scale-invariance). Second, fitness evaluation is based not on simulation but on the exact recursion formulas derived from the Laplace–Stieltjes transform of the system, with special analytic limits at the singular points $δ_{i}, ϑ_{i}$ to maintain numerical stability and continuity. Third, because the resulting objective is deterministic, piecewise–smooth, and nonconvex, we modified the Bat dynamics to add controlled jitter around the incumbent best solution in order to escape plateaus, while keeping acceptance and loudness schedules less aggressive than in the stochastic case. These adaptations ensure that the search respects the system’s structural constraints, that each evaluation is mathematically exact and stable, and that the algorithm remains effective on the nonconvex landscape induced by the working-vacation queue with balking, reneging, and feedback.

Remark 5.

We would like to stress that the numerical results in Section 6 are based on exact recursive formulations of the stationary distribution and therefore are not subject to sampling variability as in simulation-based studies. Consequently, conventional statistical tests are not directly applicable. To ensure robustness, however, we performed several independent optimization runs with different initializations and obtained highly consistent solutions, indicating that the reported outcomes are stable with respect to the metaheuristic search. The performance measures themselves are deterministic functionals of the stationary probabilities, so their accuracy hinges on numerical stability rather than statistical fluctuations, which has been carefully verified in the implementation.

Remark 6.

The convergence of the proposed optimization procedure must be understood in a numerical rather than purely analytical sense. Because the underlying cost functional is evaluated exactly (through the backward recursions and analytic limits at singular points), the sequence of best solutions generated by the Bat algorithm is deterministic and exhibits monotone improvement, eventually stabilizing within a narrow tolerance. However, as is well known for population-based metaheuristics, global optimality cannot be guaranteed theoretically; what can be claimed is convergence to a numerically stable and high-quality local optimum, with robustness reinforced by projection into the ordered feasible cone and by independent replications of the search. From a computational perspective, the recursion structure makes each fitness evaluation linear in the buffer size, so the optimization remains tractable: in practice, for realistic system dimensions, the entire optimization completes within seconds to minutes on standard hardware. This ensures that the solution time is fully acceptable for decision-support and design purposes, even though a theoretical guarantee of global optimality is absent.

Remark 7.

In practical environments such as call centers, healthcare units, or flexible manufacturing, the model would be implemented by first calibrating its structural parameters—arrival intensity, patience (reneging) distributions, balking probabilities, service-time laws under regular versus working-vacation modes, and feedback frequencies—directly from operational data. Concretely, this requires access to high-resolution transactional logs (arrival timestamps, abandonment times, service completions, switch-overs to reduced-capacity states, and post-service return behavior). Standard statistical procedures (maximum likelihood, nonparametric hazard estimation, survival analysis) would then be used to fit the inter-arrival and service-time distributions as well as the patience and balking functions. With these empirically derived inputs, the queuing model can be parameterized and the optimization framework applied to support staffing, scheduling of reduced-capacity periods, and service-level decisions tailored to the specific application domain.

Remark 8.

The adoption of the Bat Algorithm was motivated by the structure of our optimization problem, which involves a deterministic yet piecewise-smooth objective derived from LST-based recursions with removable singularities. Such a landscape is nonconvex and may exhibit flat regions, making derivative-free methods particularly suitable. The Bat Algorithm, with its built-in balance of global exploration and adaptive local exploitation, proves efficient when coupled with our log-space parameterization and projection scheme enforcing $(ν < μ_{2} < μ_{1})$ . Although alternative metaheuristics (e.g., GA, PSO, DE) could also be applied and are expected to produce comparable optima, the difference would primarily concern convergence speed and robustness. Our choice therefore reflects adequacy and computational efficiency rather than exclusivity, and the methodology remains transferable to other metaheuristic frameworks.

7. Model Application in Practice

–. We consider a flexible manufacturing facility that must concurrently support make-to-order (MTO) jobs—customer-specific production released upon demand—and opportunistic make-to-stock (MTS) runs that replenish standard inventory. The shop comprises two heterogeneous machines (servers) operating under a first-in, first-out (FIFO) discipline. Both machines primarily process MTO work. Incoming MTO orders that find both machines busy enter a finite buffer of size N; an arrival that encounters a full buffer is blocked and lost; otherwise, it joins the queue.

–. To sustain high utilization during light-load periods while preserving responsiveness to MTO demand, one machine (Server 2) may enter a working vacation (WV) whenever it becomes idle. During a WV, Server 2 continues to process jobs at a deliberately reduced service rate $ν$ (for example, under a low-power or maintenance mode). The other machine (Server 1) remains fully available for MTO work at rate $μ_{1}$ . This mechanism reflects the widely used modeling idea that a server need not switch completely off; instead, it can operate at a lower productive rate rather than halt altogether [30,78,79]. At the completion of a WV, Server 2 adapts to the system state: if the system is empty, it immediately initiates a new WV, capturing extended low-load spells; otherwise, it returns to normal operation and serves at its regular MTO rate $μ_{2} > ν$ . This interruption/continuation logic, standard in the WV literature, captures practical rules that preserve responsiveness when backlog exists [78,79].

–. Customer impatience is represented through two complementary behaviors. Balking allows a newly arriving order to decline entry when congestion is visible; in a finite-capacity setting, this is modeled by a state-dependent joining probability $θ_{i}$ that depends on the observed system size i, generalizing the classical economics-of-balking perspective [80]. Reneging captures abandonment by customers who have joined but lose patience while waiting; we posit an exponential patience (hazard) rate $α$ while in the queue, following standard models of queues with abandonment [81,82,83]. Finally, we incorporate Bernoulli feedback at service completion: with probability $1 - β$ , a job is routed back to the queue for rework (e.g., remedial processing or quality corrections), whereas with probability $β$ , it departs permanently. This abstraction is classical for rework loops and repeated service attempts [84,85].

–. In aggregate, the model integrates four operational realities—finite capacity, working vacations with interruption, impatience (balking and reneging), and rework via Bernoulli feedback—within a unified and tractable framework. The purpose is to support design decisions such as sizing the buffer N, tuning the aggressiveness of WV operation (choice of $ν$ and the vacation policy), and quantifying the trade-offs among loss, abandonment, throughput, and cycle time in a two-server, heterogeneous MTO–MTS environment [86].

Remark 9.

While the present work is developed in a purely theoretical framework, the incorporation of empirical operational logs would offer a natural and valuable extension. Such data could serve to calibrate the impatience and feedback parameters—typically idealized in analytical models—by confronting them with observed inter-arrival, patience, and retry patterns. This would not only validate the adequacy of the exponential-type assumptions but also reinforce the practical relevance of the theoretical findings without altering their mathematical consistency.

Remark 10.

In large-scale systems, the cumulative cost of repeated fitness evaluations may pose challenges; however, this can be overcome through closed-form recursions, analytic treatment of singularities, and projection-based feasibility enforcement, which together ensure numerical stability and preserve scalability without simulation overhead.

8. Conclusions and Perspectives

We developed a steady-state analysis for a finite-buffer $G I / M / 2 / N$ queue with two heterogeneous servers, impatience in the form of balking and reneging, Bernoulli feedback, and multiple working vacations. Leveraging the supplementary-variable technique and a tailored recursion, we derived balance equations and system-size distributions at both pre-arrival and arbitrary epochs. From these, we obtained standard performance measures, including blocking and reneging probabilities, mean queue length, mean sojourn time, throughput, and server utilizations. On the managerial side, we formulated an economic objective that trades off congestion, abandonment, and capacity costs, and we optimized the decision variables using the Bat Algorithm. A numerical study implemented in R demonstrated how service heterogeneity, vacation parameters, buffer size, feedback probability, and impatience intensities jointly shape delay, loss, and total cost, thereby providing actionable guidance for the design of call centers, flexible manufacturing facilities, and telecommunication access networks.

Our analysis assumes exponential service times and a renewal input process. While these choices encompass a broad range of applications, richer service-time models and strongly correlated arrivals would require additional structure such as phase-type or Markov-modulated input. A further limitation is the assumption of FCFS with no routing asymmetry; incorporating priority or skill-based dispatching under heterogeneity would improve fidelity in many service contexts.

Perspectives. Several extensions are both impactful and technically tractable. One avenue replaces exponential service with phase-type distributions, yielding a $G I / P H / 2 / N$ counterpart amenable to matrix-analytic and matrix-geometric methods while retaining the WV/feedback/impatience structure [66,87]. Another introduces nonstationarity and burstiness through Markov-modulated or batch Markovian arrivals, or via time-varying staffing, following classical MAP/MMPP/BMAP frameworks and time-dependent staffing models [88,89,90,91]. Reliability considerations can be incorporated by adding breakdown and repair states coupled with WV control and impatience/feedback, connecting to network performance bounds and product-form insights [92,93]. Endogenizing dispatching and staffing with heterogeneous speeds—for example, via skill-based or priority routing—opens the door to the joint design of routing and WV control, building on established service-system modeling [94,95]. Multi-class demand with class-dependent balking, reneging, and feedback would enable analysis of fairness and service-level guarantees using scheduling and prioritization tools [95]. On the data side, renewal laws and impatience hazards can be estimated from operational logs, and feedback/WV parameters can be calibrated by likelihood or Bayesian approaches that align with the supplementary-variable framework [96,97]. Learning-augmented control—combining analytical models with online policy learning to adapt service rates, admission, or vacation durations under stability constraints—is a promising direction, with index policies and reinforcement learning furnishing principled mechanisms [98,99]. Risk-sensitive and robust design objectives, such as CVaR or chance constraints, would address tail performance and model uncertainty [100,101], while rare-event techniques and large deviations can quantify extreme blocking and delay beyond average-case metrics [102,103,104]. Finally, computational acceleration exploiting sparsity and Toeplitz structure in the balance equations—together with state-of-the-art matrix-analytic libraries and fast matrix-exponential routines—would scale the numerics to large buffers and richer phase structures [87,105].

Palm–time viewpoint. Across these extensions, maintaining an explicit Palm–time conversion under non-Poisson input is essential, since performance targets often mix arrival-based and time-based metrics. This separation prevents inadvertent PASTA assumptions and leads to correct comparisons of customer-experienced versus system-averaged quantities.

On the case $N = \infty$ . When the buffer is unbounded, positive recurrence is no longer guaranteed. A Foster–Lyapunov function $V (i, j) = i$ yields

$L V (i, j) \leq λ θ_{i} - (μ_{1} + μ_{2} 1_{{j = 2}} + ν 1_{{j = 1}}) + {(i - 2)}^{+} α,$

so stability obtains whether the effective service capacity—accounting for the WV mix and reneging—dominates the effective admission rate. Our focus on finite N ensures positive recurrence without additional assumptions.

Closing remark. By unifying heterogeneity, impatience, feedback, and working vacations within a renewal-input framework, this study provides analytical foundations and deployable levers for modern service operations. The outlined directions connect the model to established methodological pillars and point toward robust, data-informed, and learning-enabled designs.

Author Contributions

Conceptualization, A.G.; methodology, A.G.; validation, A.G. and S.B.; formal analysis, A.G.; investigation, A.G.; writing—original draft preparation, A.G. and S.B.; writing—review and editing, A.G. and S.B. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The original contributions presented in this study are included in the article. Further inquiries can be directed to the corresponding author.

Acknowledgments

The authors gratefully acknowledge the Editor-in-Chief, the Associate Editor, and the three anonymous referees for their insightful comments and valuable suggestions, which greatly enhanced both the clarity and the overall quality of this work.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Notations

Notation for the queuing model with working vacations, balking, and reneging.

Primitives and Parameters
$λ$	External arrival rate; mean inter-arrival time $1 / λ$ .
$A (u), a (u), A^{*} (s)$	Inter-arrival cdf, density, and Laplace–Stieltjes transform.
N	Total capacity (in service + waiting).
$θ_{i}, {\bar{θ}}_{i} = 1 - θ_{i}$	Join/balking probability on seeing i customers.
$α$	Exponential impatience (reneging) rate per waiting customer.
$β, \bar{β} = 1 - β$	Departure vs. feedback probability at service completion.
$μ_{1}, μ_{2}$	Regular service rates of Server 1 and Server 2 $(μ_{2} \leq μ_{1})$ .
$ν$	Server 2 service rate during a working vacation.
$ϕ$	Working-vacation termination rate.
State processes and stationary objects
$L (t)$	Total customers in system at time t (queue + in service).
$J (t) \in {0, 1, 2}$	Server 2 state: 0 = idle on WV; 1 = busy on WV; 2 = busy (regular).
$W (t)$	Remaining inter-arrival time at t.
$P_{i, j} (u, t)$	$P (L (t) = i, J (t) = j, u \leq W (t) < u + d u)$ (density in u).
$P_{i, j} (u)$	Stationary version of $P_{i, j} (u, t)$ as $t \to \infty$ .
$P_{i, j}^{*} (s)$	LST of $P_{i, j} (u)$ in u.
$P_{i, j} (0)$	Boundary (rate) probability at $W = 0$ (arrival instant).
$P_{i, j}$	Arbitrary-epoch stationary probability of $(L, J) = (i, j)$ .
$P_{i, j}^{-}$	Pre-arrival (Palm) probability; $P_{i, j}^{-} = λ^{- 1} P_{i, j} (0)$ .
$δ_{i}$	$β (μ_{1} + ν) + ϕ + (i - 2) α$ (WV balance parameter).
$ϑ_{i}$	$β (μ_{1} + μ_{2}) + (i - 2) α$ (regular balance parameter).
Performance measures
$Π_{I}, Π_{W}, Π_{B}$	Fractions of time: Server 2 idle on WV/busy on WV/busy regular.
$Π_{L}$	Loss probability on arrival: $P_{N, 1}^{-} + P_{N, 2}^{-}$ .
$L_{s}$	Mean system size: $P_{1, 0} + \sum_{i = 1}^{N} i (P_{i, 1} + P_{i, 2})$ .
$W_{s}$	Mean sojourn time: $W_{s} = \frac{L_{s}}{λ (1 - Π_{L})}$ .
$J_{s}$	Effective joining rate.
$B_{r}, R_{r}$	Average rates of balking and reneging.
Cost model
$C_{1}, \dots, C_{9}$	Unit cost coefficients used in the economic objective.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Schematic representation of the $G I / M / 2 / N$ queuing model.

Figure 2 Effects of $μ_{1},$ $μ_{2}$ , and $ν$ on $Π_{B}$ and $Π_{W}$ .

Figure 3 Effects of $μ_{2},$ $μ_{1}$ , and $ν$ on $Π_{B}$ and $Π_{W}$ .

Figure 4 Effects of $ν,$ $μ_{2}$ , and $μ_{1}$ on $Π_{B}$ and $Π_{W}$ .

Figure 8 Effects of $μ_{1}$ and $μ_{2}$ on the total expected cost function.

Figure 9 Effects of $μ_{1}$ and $ν$ on the total expected cost function.

Table 1

State-of-the-art summary.

© 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

Cost Optimization in a GI/M/2/N Queue with Heterogeneous Servers, Working Vacations, and Impatient Customers via the Bat Algorithm

Content area

Abstract

Full text

1. Introduction

2. Model Mathematical Description

3. Steady-State Solution

4. Performance Measures

5. Cost Model and Optimization Study

6. Numerical Results

6.1. Numerical Results of System Performance Measures

6.2. Numerical Results of Cost Optimization

6.3. Discussions

7. Model Application in Practice

8. Conclusions and Perspectives