<figure><img src="image_1.png" /></figure>

HAL
open science

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

Abdelhak Guendouzi, Salim Bouzebda

► To cite this version:

Abdelhak Guendouzi, Salim Bouzebda. Cost Optimization in a GI/M/2/N Queue with Heterogeneous Servers, Working Vacations, and Impatient Customers via the Bat Algorithm. Mathematics, 2025, 13 (21), pp.3559. 10.3390/math13213559. hal-05352333

HAL Id: hal-05352333

https://hal.science/hal-05352333v1

Submitted on 7 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

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

Abdelhak Guendouzi <fnref n="1" />,† and Salim Bouzebda <fnref n="2" />,*,† ID

## Abstract

This paper analyzes a finite-capacity GI/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 (GI), 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 $(\mu_1, \mu_2, \nu, \alpha, \phi, \beta)$.

**Keywords:** GI/M/2/N queue; heterogeneous servers; multiple working vacations; impatient customers; Bernoulli feedback; Bat algorithm

**MSC:** 60K25; 68M20; 90B22

## 1. Introduction

Modern service systems routinely face finite buffers, rework loops, energy-aware operation, and behavioral customer responses. Treating these features jointly and without

<figure><img src="image_3.png" /></figure>

### check for updates

Academic Editor: Andrea Scozzari

Received: 19 September 2025

Revised: 7 October 2025

Accepted: 4 November 2025

Published: 6 November 2025

**Citation:** Guendouzi, A.; Bouzebda, S. Cost Optimization in a GI/M/2/N Queue with Heterogeneous Servers, Working Vacations, and Impatient Customers via the Bat Algorithm. Mathematics 2025, 13, 3559. https://doi.org/10.3390/math13213559

**Copyright:** © 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/).

<footnote n="1">¹ Institute of Science, University Center of Aflou, B. P. 603, Aflou 03001, Algeria</footnote>

<footnote n="2">² Université de Technologie de Compiègne, LMAC (Laboratory of Applied Mathematics of Compiègne), CS 60 319, 60 203 Compiègne, France</footnote>

<footnote>* Correspondence: salim.bouzebda@utc.fr</footnote>

<footnote>† These authors contributed equally to this work.</footnote>
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–10] and Little’s law [11]). For recent references, see [12–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–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–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–35], multi-server configurations M/M/c [36–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 GI/·/· capture such features parsimoniously. For $GI/M/1$ and finite-capacity variants with WV, notable contributions include explicit sojourn-time characterizations, interruption mechanisms, and phase-type vacation structures [31,42–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 sup-plementary 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, ex-
plicit probability-generating functions across vacation and busy periods have been ob-
tained in [47], while finite buffers, policy comparisons, and cost optimization under impa-
tience/WV interactions appear in [48–55]. A related and practically pervasive mechanism
is feedback: after service completion, a customer returns to the queue with fixed proba-
bility, modeling rework loops or packet retransmissions. Originating with Takács [56],
feedback has been analyzed in Markovian and renewal settings [57–64]. Analytically, feed-
back 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 au-
tomated systems but frequently violated in human-service or mixed-technology envi-
ronments. Heterogeneity ($\mu_1 \neq \mu_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–69], with focused developments for
$M/M/2$ (including WV) and transient analyses under impatience [70–72]. In contrast,
renewal-input $GI/M/2/N$ queues with heterogeneous servers remain comparatively un-
derexplored; 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 $GI/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 \ge 1}$ be i.i.d. inter-arrival times with cdf $A(u)$, density $a(u)$, and mean $1/\lambda \in (0, \infty)$; service times are exponential with heterogeneous rates $\mu_1, \mu_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{\theta}(i) \in [0, 1]$ when $i$ customers are present upon arrival (with $\bar{\theta}(N) = 1$ under finite capacity), while reneging is modeled via exponential patience with hazard $\alpha(i)$ that may depend on the queue length. Feedback is Bernoulli with probability $\bar{\beta} \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.
<table><caption>**Table 1.** State-of-the-art summary.</caption><thead><tr><th>Reference Block</th><th>Arr.</th><th>Fin N</th><th>2H</th><th>WV</th><th>Balk.</th><th>Ren.</th><th>FB</th><th>Palm/Time</th><th>Scope/Method</th></tr></thead><tbody><tr><td>M/M/1 [30]</td><td>P</td><td>•</td><td>○</td><td>•</td><td>○</td><td>○</td><td>○</td><td>○</td><td>Canonical working-vacation idea; single server baseline.</td></tr><tr><td>WV, multi-server (homog.) [35–38,40,41]</td><td>P/GI</td><td>•</td><td>•</td><td>•</td><td>○</td><td>○</td><td>○</td><td>○</td><td>Homogeneous c; PH/Erlang/discrete-time variants; some finite buffers.</td></tr><tr><td>GI input with WV (mostly GI/M/1) [42,45,46]</td><td>GI</td><td>•</td><td>○</td><td>•</td><td>○</td><td>○</td><td>○</td><td>•</td><td>Renewal via supplementary variable/PH; mainly single-server.</td></tr><tr><td>WV + impatience (Markovian) [48,50,52–54]</td><td>P</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>○</td><td>○</td><td>Balking &amp; reneging with WV; chiefly M/M/1 or homogeneous servers.</td></tr><tr><td>Feedback queues (no WV) [56–59,61,63]</td><td>P/GI</td><td>•</td><td>•</td><td>○</td><td>○</td><td>○</td><td>•</td><td>○</td><td>Classical rework/feedback; little coupling with WV.</td></tr><tr><td>Heterogeneous M/M/2 [65–67,70–72]</td><td>P</td><td>•</td><td>•</td><td>○</td><td>○</td><td>•</td><td>○</td><td>○</td><td>Two unequal servers; limited impatience/WV interaction.</td></tr><tr><td>GI/M/2/N (heterogeneous) [73]</td><td>GI</td><td>•</td><td>•</td><td>○</td><td>•</td><td>•</td><td>○</td><td>○</td><td>Finite capacity under renewal input; no WV/feedback/cost pipeline.</td></tr><tr><td>Present work</td><td>P/GI</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>Unified GI/M/2/N with heterogeneity, WV, state-dep. balking, exp. reneging, Bernoulli feedback; explicit Palm vs. time; supplementary variable + LST recursion (O(N)); full performance; economic cost + Bat optimization [74–76] Sections 3–5.</td></tr></tbody></table>

**Legend:** Arr. (P = Poisson; GI = renewal). Fin N = finite buffer. 2H = two heterogeneous servers. WV = working vacations. Balk./Ren. = balking/reneging. FB = feedback/rework. “Palm/Time” indicates an explicit separation of pre-arrival vs. time-stationary observables. Marks: ● addressed; ● partially addressed; ○ not addressed.

**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 $\mu_1$ and $\mu_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 $(\mu_1, \mu_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 GI 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 $\mu_1 \neq \mu_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 GI 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 $\pi^{(A)}$ denote the pre-arrival stationary law on the augmented state space and $\pi$ 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 $\pi^{(A)}$ and $\pi$ 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 $\pi^{(A)}$ and $\pi$ 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 $\mu_1, \mu_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–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 GI/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 \ge 0$), Laplace-Stieltjes transform

$$A^*(s) = \int_{[0, \infty)} e^{-su} dA(u), \quad s \ge 0,$$

and mean inter-arrival time

$$\mathbb{E}[U] = -\frac{d}{ds}A^*(s)\Big|_{s=0} = \frac{1}{\lambda} \in (0, \infty).$$
*   **Finite capacity.** The total system capacity is $N \in \mathbb{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 $\theta_i$ and balks with probability $\bar{\theta}_i := 1 - \theta_i$. We assume

$$ \theta_0 = \theta_1 = 1, \quad 0 \le \theta_{i+1} \le \theta_i \le 1 (2 \le i \le N-1), \quad \theta_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 $(\theta_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 $\mu_1$ and $\mu_2$, respectively, with $\mu_2 \le \mu_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 $\phi$. During a working vacation, Server 2 remains active but at a reduced exponential service rate $\nu \le \mu_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 \ge 2$ customers, the queue length equals $i - 2$. Each waiting customer runs an independent exponential impatience timer $T \sim \text{Exp}(\alpha)$. 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 $\beta$ and, with complementary probability $\bar{\beta} := 1 - \beta$, 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.

<figure><img src="image_4.png" /><figcaption>Figure 1. Schematic representation of the GI/M/2/N queuing model.</figcaption></figure>
### 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 \ge 0\}$, with the state space defined as

$$ \Omega = \{(0,0,u) \cup (1,0,u) \cup (i,1,u) \cup (i,2,u) : 1 \le i \le N, u \ge 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{cases} 0, & \text{if Server 2 is idle during the working vacation period,} \\ 1, & \text{if Server 2 is busy during the working vacation period,} \\ 2, & \text{if Server 2 is busy during the regular service period.} \end{cases} $$

* $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:

$$ P_{i,0}(u,t) du = \mathbb{P}(L(t) = i, u \le W(t) < u + du, J(t) = 0), \quad u \ge 0, i = 0,1, $$

$$ P_{i,j}(u,t) du = \mathbb{P}(L(t) = i, u \le W(t) < u + du, J(t) = j), \quad u \ge 0, j = 1,2, 1 \le i \le N. $$

In steady-state, as $t \to \infty$, we have

$$ P_{i,0}(u) = \lim_{t \to \infty} P_{i,0}(u,t), \quad i = 0,1, $$

$$ P_{i,j}(u) = \lim_{t \to \infty} P_{i,j}(u,t), \quad j = 1,2, 1 \le i \le N. $$

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 + dt$. Taking $\lim_{t \to \infty}$ and after some simplification, we have the following set of differential difference equations.

$$ -P_{0,0}^{(1)}(u) = \beta[\mu_1 P_{1,0}(u) + \nu P_{1,1}(u) + \mu_2 P_{2,1}(u)], \quad (1) $$

$$ -P_{1,0}^{(1)}(u) = -\beta\mu_1 P_{1,0}(u) + \beta[\nu P_{2,1}(u) + \mu_2 P_{2,2}(u)] + a(u)P_{0,0}(0), \quad (2) $$

$$ -P_{1,1}^{(1)}(u) = -(\phi + \beta\nu)P_{1,1}(u) + \beta\mu_1 P_{2,1}(u), \quad (3) $$

$$ -P_{2,1}^{(1)}(u) = -[\beta(\mu_1 + \nu) + \phi]P_{2,1}(u) + [\beta(\mu_1 + \nu) + \alpha]P_{3,1}(u) \\ +a(u)[P_{1,0}(0) + P_{1,1}(0)], \quad (4) $$

$$ -P_{i,1}^{(1)}(u) = -[\beta(\mu_1 + \nu) + \phi + (i-2)\alpha]P_{i,1}(u) + [\beta(\mu_1 + \nu) + (i-1)\alpha]P_{i+1,1}(u) \\ +a(u)[\theta_{i-1}P_{i-1,1}(0) + \bar{\theta}_i P_{i,1}(0)], \quad 3 \le i \le N-1, \quad (5) $$

$$ -P_{N,1}^{(1)}(u) = -[\beta(\mu_1 + \nu) + \phi + (N-2)\alpha]P_{N,1}(u) + a(u)[\theta_{N-1}P_{N-1,1}(0) + P_{N,1}(0)], \quad (6) $$

$$ -P_{1,2}^{(1)}(u) = -\beta\mu_2 P_{1,2}(u) + \phi P_{1,1}(u) + \beta\mu_1 P_{2,2}(u), \quad (7) $$

$$ -P_{i,2}^{(1)}(u) = -[\beta(\mu_1 + \mu_2) + (i-2)\alpha]P_{i,2}(u) + [\beta(\mu_1 + \mu_2) + (i-1)\alpha]P_{i+1,2}(u) \\ +\phi P_{i,1}(u) + a(u)[\theta_{i-1}P_{i-1,2}(0) + \bar{\theta}_i P_{i,2}(0)], \quad 2 \le i \le N-1, \quad (8) $$
$$ -P_{N,2}^{(1)}(u) = -[\beta(\mu_1 + \mu_2) + (N-2)\alpha)P_{N,2}(u) + \phi P_{N,1}(u) + a(u)[\theta_{N-1}P_{N-1,2}(0) + P_{N,2}(0)], \quad (9) $$

where $P_{i,0}(0), i = 0, 1$ and $P_{i,j}(0), j = 1, 2, 1 \le i \le 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{aligned} P_{i,0}^*(s) &= \int_0^\infty e^{-su} P_{i,0}(u) du, \quad i = 0, 1, \\ P_{i,j}^*(s) &= \int_0^\infty e^{-su} P_{i,j}(u) du, \quad j = 1, 2, \quad 1 \le i \le N. \end{aligned} $$

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 \le i \le 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 \le i \le N$, define

$$ \delta_i := \beta(\mu_1 + \nu) + \phi + (i - 2)\alpha, \quad \vartheta_i := \beta(\mu_1 + \mu_2) + (i - 2)\alpha. $$

Multiplying Equations (1)–(9) by $e^{-su}$ and integrating with respect to $u$ over $[0, \infty)$, we obtain the transformed relations

$$ -s P_{0,0}^*(s) = -P_{0,0}(0) + \beta[\mu_1 P_{1,0}^*(s) + \nu P_{1,1}^*(s) + \mu_2 P_{1,2}^*(s)]. \quad (10) $$

$$ (\beta\mu_1 - s)P_{1,0}^*(s) = -P_{1,0}(0) + \beta[\nu P_{2,1}^*(s) + \mu_2 P_{2,2}^*(s)] + A^*(s)P_{0,0}(0). \quad (11) $$

$$ (\phi + \beta\nu - s)P_{1,1}^*(s) = P_{1,1}(0) + \beta\mu_1 P_{2,1}^*(s). \quad (12) $$

$$ (\delta_2 - s)P_{2,1}^*(s) = -P_{2,1}(0) + (\delta_3 - \phi)P_{3,1}^*(s) + A^*(s)[P_{1,0}(0) + P_{1,1}(0) + \bar{\theta}_2 P_{2,1}(0)]. \quad (13) $$

$$ (\delta_i - s)P_{i,1}^*(s) = -P_{i,1}(0) + (\delta_{i+1} - \phi)P_{i+1,1}^*(s) + A^*(s)[\theta_{i-1}P_{i-1,1}(0) + \bar{\theta}_i P_{i,1}(0)], \quad 3 \le i \le N-1. \quad (14) $$

$$ (\delta_N - s)P_{N,1}^*(s) = -P_{N,1}(0) + A^*(s)[\theta_{N-1}P_{N-1,1}(0) + P_{N,1}(0)]. \quad (15) $$

$$ (\beta\mu_2 - s)P_{1,2}^*(s) = -P_{1,2}(0) + \phi P_{1,1}^*(s) + \beta\mu_1 P_{2,2}^*(s). \quad (16) $$

$$ (\vartheta_i - s)P_{i,2}^*(s) = -P_{i,2}(0) + \phi P_{i,1}^*(s) + \vartheta_{i+1} P_{i+1,2}^*(s) + A^*(s)[\theta_{i-1}P_{i-1,2}(0) + \bar{\theta}_i P_{i,2}(0)], \quad 2 \le i \le N-1. \quad (17) $$

$$ (\vartheta_N - s)P_{N,2}^*(s) = -P_{N,2}(0) + \phi P_{N,1}^*(s) + A^*(s)[\theta_{N-1}P_{N-1,2}(0) + P_{N,2}(0)]. \quad (18) $$

Summing (10)–(18) and simplifying yields

$$ -A^*(s)\left(\sum_{i=0}^{1} P_{i,0}(0) + \sum_{i=1}^{N} (P_{i,1}(0) + P_{i,2}(0))\right) = s\left(\sum_{i=0}^{1} P_{i,0}^*(s) + \sum_{i=1}^{N} (P_{i,1}^*(s) + P_{i,2}^*(s))\right). \quad (19) $$

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

$$ \sum_{i=0}^{1} P_{i,0}(0) + \sum_{i=1}^{N} (P_{i,1}(0) + P_{i,2}(0)) = \lambda. \quad (20) $$
Equation (20) expresses conservation of flow: the long-run average inflow to the system per unit time equals the mean arrival rate $\lambda$. 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 \le i \le 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{aligned}
P_{i,j}^{-} &= \lim_{t \to \infty} \mathbb{P}(L(t) = i, S(t) = j | U(t) = 0) \\
&= \lim_{t \to \infty} \frac{\mathbb{P}(L(t) = i, S(t) = j, U(t) = 0)}{\mathbb{P}(U(t) = 0)}.
\end{aligned}
$$

Consequently,

$$ P_{i,j}^{-} = \frac{1}{\lambda} P_{i,j}(0), \quad \begin{array}{l} j = 0, i \in \{0, 1\}, \\ j \in \{1, 2\}, 1 \le i \le N, \end{array} \qquad (21) $$

where $\lambda$ 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{aligned}
P_{i,0}(0) (i = 0, 1), & \quad P_{i,j}(0) (j = 1, 2, 1 \le i \le N), \\
P_{i,0}^{*}(s) (i = 0, 1), & \quad P_{i,j}^{*}(s) (j = 1, 2, 1 \le i \le N),
\end{aligned}
$$

recursively in terms of $P_{N,0}(0)$ (details below). Throughout, $A^{*}(\cdot)$ denotes the relevant transform and $A^{*(k)}$ its $k$th derivative.

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

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

$$ P_{N-1,1}(0) = \psi_{N-1} P_{N,1}(0), \quad \psi_{N-1} := \frac{1 - A^*(\delta_N)}{A^*(\delta_N) \theta_{N-1}}. \qquad (22) $$

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

Substituting (22) back into (15) gives

$$ (\delta_N - s) P_{N,1}^{*}(s) = \left( A^{*}(s)(\theta_{N-1}\psi_{N-1} + \psi_N) - \psi_N \right) P_{N,1}(0), \qquad (23) $$

where we set $\psi_N \equiv 1$. Hence, for $s \ne \delta_N$,

$$ P_{N,1}^{*}(s) = \frac{A^{*}(s)(\theta_{N-1}\psi_{N-1} + \psi_N) - \psi_N}{\delta_N - s} P_{N,1}(0). \qquad (24) $$

Step 3: The value at $s = \delta_N$ and higher derivatives.

Differentiating (23) with respect to $s$ and then setting $s = \delta_N$ yields

$$ P_{N,1}^{*}(\delta_N) = -A^{*(1)}(\delta_N)(\theta_{N-1}\psi_{N-1} + \psi_N)P_{N,1}(0). \qquad (25) $$

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

$$ (\delta_N - s) P_{N,1}^{*(l)}(s) - l P_{N,1}^{*(l-1)}(s) = A^{*(l)}(s)(\theta_{N-1}\psi_{N-1} + \psi_N) P_{N,1}(0). \qquad (26) $$
Step 4: Compact representation.

Combining (24)–(26), we obtain

$$P_{N,1}^*(s) = \zeta_{N,s} P_{N,1}(0), \qquad (27)$$

where

$$\zeta_{N,s} = \begin{cases} \frac{A^*(s)(\theta_{N-1}\psi_{N-1} + \psi_N) - \psi_N}{\delta_N - s}, & s \neq \delta_N, \\ -A^{*(1)}(s)(\theta_{N-1}\psi_{N-1} + \psi_N), & s = \delta_N, \end{cases}$$

and

$$\zeta_{N,s}^{(l)} = \begin{cases} \frac{A^{*(l)}(s)(\theta_{N-1}\psi_{N-1} + \psi_N) + l\zeta_{N,s}^{(l-1)}}{\delta_N - s}, & s \neq \delta_N, \\ -\frac{A^{*(l+1)}(s)(\theta_{N-1}\psi_{N-1} + \psi_N)}{l+1}, & s = \delta_N, \end{cases}$$

with $\zeta_{N,s}^{(l)}$ denoting the $l$th derivative of $\zeta_{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 = \delta_{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:

$$P_{N-2,1}(0) = \psi_{N-2}P_{N,1}(0), \qquad (28)$$

with

$$\psi_{N-2} = \frac{\psi_{N-1} - (\delta_N - \phi)\zeta_{N,\delta_{N-1}} - A^*(\delta_{N-1})\bar{\theta}_{N-1}\psi_{N-1}}{A^*(\delta_{N-1})\theta_{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:

$$P_{N-1,1}^*(s) = \zeta_{N-1,s} P_{N,1}(0), \qquad (29)$$

where

$$\zeta_{N-1,s} = \begin{cases} \frac{A^*(s)(\theta_{N-2}\psi_{N-2} + \bar{\theta}_{N-1}\psi_{N-1}) + (\delta_N - \phi)\zeta_{N,s} - \psi_{N-1}}{(\delta_{N-1} - s)}, & \text{if } s \neq \delta_{N-1}, \\ -(A^{*(1)}(s)(\theta_{N-2}\psi_{N-2} + \bar{\theta}_{N-1}\psi_{N-1}) + (\delta_N - \phi)\zeta_{N,s}^{(1)}), & \text{if } s = \delta_{N-1}, \end{cases}$$

where

$$\zeta_{N-1,s}^{(l)} = \begin{cases} \frac{A^{*(l)}(s)(\theta_{N-2}\psi_{N-2} + \bar{\theta}_{N-1}\psi_{N-1}) + (\delta_N - \phi)\zeta_{N,s}^{(l)} + l\zeta_{N-1,s}^{(l-1)}}{(\delta_{N-1} - s)}, & \text{if } s \neq \delta_{N-1}, \\ -\frac{A^{*(l+1)}(s)(\theta_{N-2}\psi_{N-2} + \bar{\theta}_{N-1}\psi_{N-1}) + (\delta_N - \phi)\zeta_{N,s}^{(l+1)}}{l+1}, & \text{if } s = \delta_{N-1}. \end{cases}$$

For $i = N - 2, N - 3, \dots, 3$ in Equation (14), $P_{i,1}(0)$, for $2 \le i \le N - 3$, and $P_{i,1}^*(s)$, for $3 \le i \le N - 2$, are expressed in terms of $P_{N,1}(0)$ as follows:

$$P_{i-1,1}(0) = \psi_{i-1}P_{N,1}(0), \quad i = N - 2, N - 3, \dots, 3, \qquad (30)$$
where

$$ \psi_{i-1} = \frac{\psi_i - (\delta_{i+1} - \phi)\zeta_{i+1,\delta_i} - A^*(\delta_i)\bar{\theta}_i\psi_i}{A^*(\delta_i)\theta_{i-1}} \quad i = N-2, N-3, \dots, 3, $$

and

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

where

$$ \zeta_{i,s} = \begin{cases} \frac{A^*(s)(\theta_{i-1}\psi_{i-1} + \bar{\theta}_i\psi_i) + (\delta_{i+1} - \phi)\zeta_{i+1,s} - \psi_i}{(\delta_i - s)}, & \text{if } s \neq \delta_i, \\ -(A^{*(1)}(s)(\theta_{i-1}\psi_{i-1} + \bar{\theta}_i\psi_i) + (\delta_{i+1} - \phi)\zeta_{i+1,s}^{(1)}), & \text{if } s = \delta_i, \end{cases} $$

where

$$ \zeta_{i,s}^{(l)} = \begin{cases} \frac{A^{*(l)}(s)(\theta_{i-1}\psi_{i-1} + \bar{\theta}_i\psi_i) + (\delta_{i+1} - \phi)\zeta_{i+1,s}^{(l)} - l\zeta_{i,s}^{(l-1)}}{(\delta_i - s)}, & \text{if } s \neq \delta_i, \\ -\frac{(\delta_{i+1} - \phi)\zeta_{i+1,\delta_i}^{(l+1)} + A^{*(l+1)}(s)(\theta_{i-1}\psi_{i-1} + \bar{\theta}_i\psi_i)}{l+1}, & \text{if } s = \delta_i. \end{cases} $$

Substituting $s = \delta_2$ into Equation (13) yields $P_{1,1}(0)$ in terms of $P_{N,1}(0)$ and $P_{1,0}(0)$ as follows:

$$ P_{1,1}(0) = \psi_1 P_{N,1}(0) + \omega P_{1,0}(0), \quad (32) $$

where

$$ \psi_1 = \frac{\psi_2 - (\delta_3 - \phi)\zeta_{3,\delta_2} - A^*(\delta_2)\bar{\theta}_2\psi_2}{A^*(\delta_2)}, \quad \text{and} \quad \omega = -\frac{A^*(\delta_2)}{A^*(\delta_2)} = -1. $$

Using Equation (32) in Equation (13), we obtain $P_{2,1}^*(s)$ in terms of $P_{N,1}(0)$ as follows:

$$ P_{2,1}^*(s) = \zeta_{2,s} P_{N,1}(0), \quad (33) $$

where

$$ \zeta_{2,s} = \begin{cases} \frac{-\psi_2 + (\delta_3 - \phi)\zeta_{3,s} + A^*(s)(\psi_1 + \bar{\theta}_2\psi_2)}{(\delta_2 - s)}, & \text{if } s \neq \delta_2, \\ -((\delta_3 - \phi)\zeta_{3,s}^{(1)} + A^{*(1)}(s)(\psi_1 + \bar{\theta}_2\psi_2)), & \text{if } s = \delta_2, \end{cases} $$

and

$$ \zeta_{2,s}^{(l)} = \begin{cases} \frac{(\delta_3 - \phi)\zeta_{3,s}^{(l)} + A^{*(l)}(s)(\psi_1 + \bar{\theta}_2\psi_2) - l\zeta_{2,s}^{(l-1)}}{(\delta_2 - s)}, & \text{if } s \neq \delta_2, \\ -\frac{(\delta_3 - \phi)\zeta_{3,s}^{(l+1)} + A^{*(l+1)}(s)(\psi_1 + \bar{\theta}_2\psi_2)}{l+1}, & \text{if } s = \delta_2. \end{cases} $$

From Equation (12), we have

$$ P_{1,1}^*(s) = \zeta_{1,s} P_{N,1}(0) + \tau_{1,s} P_{1,0}(0), \quad (34) $$
where

$$
\begin{align*}
\zeta_{1,s} &= \begin{cases}
\frac{\beta\mu_1\zeta_{2,s} - \psi_1}{(\phi + \beta\nu - s)}, & \text{if } s \neq \phi + \beta\nu, \\
-\beta\mu_1\zeta_{2,s}^{(1)}, & \text{if } s = \phi + \beta\nu.
\end{cases} \\
\tau_{1,s} &= \begin{cases}
\frac{\omega}{(\phi + \beta\nu - s)}, & \text{if } s \neq \phi + \beta\nu, \\
0, & \text{if } s = \phi + \beta\nu,
\end{cases}
\end{align*}
$$

where

$$
\begin{align*}
\zeta_{1,s}^{(l)} &= \begin{cases}
\frac{\beta\mu_1\zeta_{2,s}^{(l)} - l\zeta_{1,s}^{(l-1)}}{(\phi + \beta\nu - s)}, & \text{if } s \neq \phi + \beta\nu, \\
-\frac{\beta\mu_1\zeta_{2,s}^{(l+1)}}{l+1}, & \text{if } s = \phi + \beta\nu.
\end{cases} \\
\tau_{1,s}^{(l)} &= \begin{cases}
\frac{l\tau_{1,s}^{(l-1)}}{(\phi + \beta\nu - s)}, & \text{if } s \neq \phi + \beta\nu, \\
0, & \text{if } s = \phi + \beta\nu.
\end{cases}
\end{align*}
$$

Putting $s = \vartheta_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:

$$
P_{N-1,2}(0) = \eta_{N-1}P_{N,2}(0) + \gamma_{N-1}P_{N,1}(0), \quad (35)
$$

where

$$
\eta_{N-1} = \frac{1 - A^*(\vartheta_N)}{A^*(\vartheta_N)\vartheta_{N-1}} \quad \text{and} \quad \gamma_{N-1} = -\frac{\phi \zeta_{N,\vartheta_N}}{A^*(\vartheta_N)\vartheta_{N-1}}.
$$

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

$$
P_{N,2}^*(s) = \rho_{N,s} P_{N,2}(0) + \chi_{N,s} P_{N,1}(0), \quad (36)
$$

where

$$
\begin{align*}
\rho_{N,s} &= \begin{cases} \frac{A^*(s)(\theta_{N-1}\eta_{N-1} + \eta_N) - \eta_N}{\vartheta_N - s}, & \text{if } \vartheta_N \neq s, \\ -A^{*(1)}(s)(\theta_{N-1}\eta_{N-1} + \eta_N), & \text{if } \vartheta_N = s. \end{cases} \\
\chi_{N,s} &= \begin{cases} \frac{\phi\zeta_{N,s} + A^*(s)\theta_{N-1}\gamma_{N-1}}{\vartheta_N - s}, & \text{if } \vartheta_N \neq s, \\ -(\phi\zeta_{N,s}^{(1)} + A^{*(1)}(s)\theta_{N-1}\gamma_{N-1}), & \text{if } \vartheta_N = s, \end{cases}
\end{align*}
$$

where
$$ \rho_{N,s}^{(l)} = \begin{cases} \frac{A^{*(l)}(s)(\theta_{N-1}\eta_{N-1} + \eta_N) + l\rho_{N,s}^{(l-1)}}{(\vartheta_N - s)}, & \text{if } \vartheta_N \neq s, \\ -\frac{A^{*(l+1)}(\vartheta_N)(\theta_{N-1}\eta_{N-1} + \eta_N)}{l+1}, & \text{if } \vartheta_N = s. \end{cases} $$

$$ \chi_{N,s}^{(l)} = \begin{cases} \frac{\phi\zeta_{N,s}^{(l)} + A^{*(l)}(s)\theta_{N-1}\gamma_{N-1} + l\chi_{N,s}^{(l-1)}}{\vartheta_N - s}, & \text{if } \vartheta_N \neq s, \\ -\frac{\phi\zeta_{N,\vartheta_N}^{(l+1)} + A^{*(l+1)}(\vartheta_N)\theta_{N-1}\gamma_{N-1}}{l+1}, & \text{if } \vartheta_N = s, \end{cases} $$

where $\eta_N = 1$ and $\gamma_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 \le i \le N-2$, and $P_{i,2}^*(s)$, $2 \le i \le N-1$, can be obtained from Equation (17) in terms of $P_{N,1}(0)$ and $P_{N,2}(0)$ as

$$ P_{i-1,2}(0) = \eta_{i-1}P_{N,2}(0) + \gamma_{i-1}P_{N,1}(0), \quad 2 \le i \le N-1, \qquad (37) $$

where

$$ \begin{aligned} \eta_{i-1} &= \frac{\eta_i - \vartheta_{i+1}\rho_{i+1,\vartheta_i} - A^*(\vartheta_i)\bar{\theta}_i\eta_i}{A^*(\vartheta_i)\theta_{i-1}}, \\ \gamma_{i-1} &= \frac{\gamma_i - \vartheta_{i+1}\chi_{i+1,\vartheta_i} - A^*(\vartheta_i)\bar{\theta}_i\gamma_i - \phi\zeta_{i,\vartheta_i}}{A^*(\vartheta_i)\theta_{i-1}}. \end{aligned} $$

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

$$ P_{i,2}^*(s) = \rho_{i,s} P_{N,2}(0) + \chi_{i,s} P_{N,1}(0), \qquad (38) $$

where

$$ \begin{aligned} \rho_{i,s} &= \begin{cases} \frac{-\eta_i + \vartheta_{i+1}\rho_{i+1,s} + A^*(s)(\theta_{i-1}\eta_{i-1} + \bar{\theta}_i\eta_i)}{(\vartheta_i - s)}, & \text{if } \vartheta_i \neq s, \\ -(\vartheta_{i+1}\rho_{i+1,s}^{(1)} + A^{*(1)}(s)(\theta_{i-1}\eta_{i-1} + \bar{\theta}_i\eta_i)), & \text{if } \vartheta_i = s. \end{cases} \\ \chi_{i,s} &= \begin{cases} \frac{-\gamma_i + \phi\zeta_{i,s} + \vartheta_{i+1}\chi_{i+1,s} + A^*(s)((\theta_{i-1}\gamma_{i-1} + \bar{\theta}_i\gamma_i))}{\vartheta_i - s}, & \text{if } \vartheta_i \neq s, \\ -(\phi\zeta_{i,s}^{(1)} + \vartheta_{i+1}\chi_{i+1,s}^{(1)} + A^{*(1)}(s)((\theta_{i-1}\gamma_{i-1} + \bar{\theta}_i\gamma_i))), & \text{if } \vartheta_i = s, \end{cases} \end{aligned} $$

where

$$ \rho_{i,s}^{(l)} = \begin{cases} \frac{\vartheta_{i+1}\rho_{i+1,s}^{(l)} + A^{*(l)}(s)(\theta_{i-1}\eta_{i-1} + \bar{\theta}_i\eta_i) + l\rho_{i,s}^{(l-1)}}{(\vartheta_i - s)}, & \text{if } \vartheta_i \neq s, \\ -\frac{\vartheta_{i+1}\rho_{i+1,s}^{(l+1)} + A^{*(l+1)}(s)(\theta_{i-1}\eta_{i-1} + \bar{\theta}_i\eta_i)}{l+1}, & \text{if } \vartheta_i = s. \end{cases} $$
$$ \chi_{i,s}^{(l)} = \begin{cases} \frac{\phi\zeta_{i,s}^{(l)} + \vartheta_{i+1}\chi_{i+1,s}^{(l)} + A^{*(l)}(s)(\theta_{i-1}\gamma_{i-1} + \bar{\theta}_i\gamma_i) + l\chi_{i,s}^{(l-1)}}{\vartheta_i - s}, & \text{if } \vartheta_i \neq s, \\ -\frac{\phi\zeta_{i,s}^{(l+1)} + \vartheta_{i+1}\chi_{i+1,s}^{(l+1)} + A^{*(l+1)}(s)(\theta_{i-1}\gamma_{i-1} + \bar{\theta}_i\gamma_i)}{l+1}, & \text{if } \vartheta_i = s. \end{cases} $$

Taking $s = \beta\mu_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)$ as

$$ P_{0,0}(0) = \Delta_0 P_{1,0}(0) + \Delta_1 P_{N,1}(0) + \Delta_2 P_{N,2}(0), \quad (39) $$

where

$$ \Delta_0 = \frac{1}{A^*(\beta\mu_1)}, \quad \Delta_1 = \frac{-\beta\nu\zeta_{2,\beta\mu_1} - \beta\mu_2\chi_{2,\beta\mu_1}}{A^*(\beta\mu_1)}, \quad \Delta_2 = \frac{-\beta\mu_2\rho_{2,\beta\mu_1}}{A^*(\beta\mu_1)}. $$

Taking $s = \phi + \beta\nu$ 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:

$$ P_{1,0}(0) = \kappa_1 P_{N,1}(0), \quad (40) $$

where

$$ \kappa_1 = \psi_1 - \beta\mu_1\zeta_{2,\phi+\beta\nu}. $$

Taking $\beta\mu_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:

$$ P_{N,2}(0) = \kappa_2 P_{N,1}(0), \quad (41) $$

where

$$ \kappa_2 = \frac{\phi\zeta_{1,\beta\mu_2} + \beta\mu_1\chi_{2,\beta\mu_2} + \phi\kappa_1\tau_{1,\beta\mu_2} - \gamma_1}{\eta_1 - \beta\mu_1\rho_{2,\beta\mu_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) = \lambda \left( \kappa_1 \Delta_0 + \Delta_1 + \kappa_2 \Delta_2 + \psi_1 + \sum_{i=2}^{N} \psi_i + \sum_{i=2}^{N} (\gamma_i + \kappa_2 \eta_i) \right)^{-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

$$ P_{N,1} = \frac{\lambda}{\delta_N} \theta_{N-1} P_{N-1,1} $$

$$ P_{i,1} = \left( \frac{\delta_{i+1} - \phi}{\delta_i} \right) P_{i+1,1} + \frac{\lambda}{\delta_i} \left( \theta_{i-1} P_{i-1,1}^{-} - \theta_i P_{i,1}^{-} \right), \quad i = N-1, \dots, 3, $$

$$ P_{2,1} = \left( \frac{\delta_3 - \phi}{\delta_2} \right) P_{3,1} + \frac{\lambda}{\delta_2} \left( P_{1,0}^{-} + P_{1,1}^{-} - \theta_2 P_{2,1}^{-} \right), $$
$$P_{1,1} = \left( \frac{\beta\mu_1}{\phi + \beta\nu} \right) P_{2,1} - \left( \frac{\lambda}{\phi + \beta\nu} \right) P_{1,1}^{-},$$

$$P_{N,2} = \frac{\phi}{\theta_N} P_{N,1} + \frac{\lambda}{\theta_N} \theta_{N-1} P_{N-1,2}^{-},$$

$$P_{i,2} = \left( \frac{\theta_{i+1}}{\theta_i} \right) P_{i+1,2} + \frac{\phi}{\theta_i} P_{i,1} + \frac{\lambda}{\theta_i} \left( \theta_{i-1} P_{i-1,2}^{-} - \theta_i P_{i,2}^{-} \right), \quad i = N-1, \dots, 2,$$

$$P_{1,2} = \frac{\mu_1}{\mu_2} P_{2,2} + \frac{\phi}{\beta\mu_2} P_{1,1} - \frac{\lambda}{\beta\mu_2} P_{1,2}^{-},$$

$$P_{1,0} = \frac{\nu}{\mu_1} P_{2,1} + \frac{\mu_2}{\mu_1} P_{2,2} + \frac{\lambda}{\beta\mu_1} \left( P_{0,0}^{-} - P_{1,0}^{-} \right),$$

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 $\lambda$ denote the external arrival rate, $N$ the system capacity, $\alpha$ the per-customer reneging rate when applicable, and $\theta_i \in [0, 1]$ the probability that an arriving customer decides to join upon seeing $i$ customers in the system (so that $\bar{\theta}_i := 1 - \theta_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

$$\Pi_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 $\lambda(1 - \Pi_L)$, we have

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

$$W_s = \frac{L_s}{\lambda (1 - \Pi_L)}. \qquad (43)$$

(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

$$\Pi_I = \sum_{i=0}^{1} P_{i,0}, \quad \Pi_W = \sum_{i=1}^{N} P_{i,1}, \quad \Pi_B = \sum_{i=1}^{N} P_{i,2}. \qquad (44)$$
(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,

$$J_s = \lambda (P_{0,0} + P_{1,0} + P_{1,1} + P_{1,2}) + \sum_{i=2}^{N} \lambda \theta_i (P_{i,1} + P_{i,2}), \quad (45)$$

$$B_r = \sum_{i=2}^{N} \lambda \bar{\theta}_i (P_{i,1} + P_{i,2}), \quad (46)$$

$$R_r = \sum_{i=2}^{N} (i-2) \alpha (P_{i,1} + P_{i,2}). \quad (47)$$

The indices (42)–(43) together with $(\Pi_I, \Pi_W, \Pi_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 ($v < \mu_2 < \mu_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,

$$(\mu_1, \mu_2, v) \in \mathbb{R}_+^3,$$

where $\mu_1$ and $\mu_2$ are the normal service rates of the two servers (with Server 1 faster than Server 2), and $v$ 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

$$(\mu_1^*, \mu_2^*, v^*) \in \arg \min_{(\mu_1, \mu_2, v)} F(\mu_1, \mu_2, v),$$

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

Unit cost elements.

Let $C_1, \dots, C_9 \ge 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;
C9: Fixed purchase (or capacity) cost per server unit.

Performance indices (given).

Let $\Pi_I, \Pi_W, \Pi_B, \Pi_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 \ge 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 $\lambda > 0$ be the external arrival rate; and let $\bar{\beta} \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

$$F(\mu_1, \mu_2, \nu) = C_1 \Pi_I + C_2 \Pi_W + C_3 \Pi_B + C_4 L_s + C_5 (B_r + R_r) + C_6 \lambda \Pi_L + (\mu_1 + \mu_2 + \nu) (C_7 + \bar{\beta} C_8) + 2 C_9. \quad (48)$$

Optimization problem.

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

$$\begin{array}{l} (\mu_1^*, \mu_2^*, \nu^*) \in \arg \min_{(\mu_1, \mu_2, \nu)} F(\mu_1, \mu_2, \nu) \\ \text{subject to} \quad 0 < \nu < \mu_2 < \mu_1, \end{array} \quad (49)$$

and the steady-state/stability conditions of the model.

Solution approach.

Because $F$ in (48) is a highly nonlinear function of $(\mu_1, \mu_2, \nu)$ 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 $(\mu_1^*, \mu_2^*, \nu^*)$.

Comments

The selection of the decision variables $(\mu_1, \mu_2, \nu)$ 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 $(\mu_1, \mu_2, \nu)$ forms the natural decision vector because service rates constitute the most direct levers for system performance. In particular,

* $\mu_1$ and $\mu_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.

* $\nu$ 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 $\nu$ 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 GI/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, expo-
nential, and Erlang-2—all calibrated with rate parameter $\lambda = 0.9$. The remaining system
parameters are specified as $\mu_1 = 2.7$, $\mu_2 = 2.3$, $\nu = 1.9$, $\bar{\beta} = 0.2$, $\alpha = 0.6$, and $\phi = 1.0$. Cus-
tomer balking is modeled through the state-dependent function $\bar{\theta}_i = i/N$ for $2 \le i \le 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.

<table><caption>**Table 2.** Steady-state distributions at pre-arrival epoch and arbitrary epoch.</caption><thead><tr><th>i</th><th>P<sub>i,0</sub><sup>-</sup></th><th>P<sub>i,1</sub><sup>-</sup></th><th>P<sub>i,2</sub><sup>-</sup></th><th>P<sub>i,0</sub></th><th>P<sub>i,1</sub></th><th>P<sub>i,2</sub></th></tr></thead><tbody><tr><td colspan="7" style="text-align: center; font-weight: bold;">Deterministic inter-arrival time distribution</td></tr><tr><td>0</td><td>0.3589552</td><td></td><td></td><td>0.7055414</td><td></td><td></td></tr><tr><td>1</td><td>0.2661013</td><td>0.06782457</td><td>0.12072168</td><td>0.1208911</td><td>0.02869207</td><td>0.00995806</td></tr><tr><td>2</td><td></td><td>0.07912866</td><td>0.07143672</td><td></td><td>0.06173432</td><td>0.04550012</td></tr><tr><td>3</td><td></td><td>0.01260821</td><td>0.01819160</td><td></td><td>0.01022743</td><td>0.01339617</td></tr><tr><td>4</td><td></td><td>0.00151615</td><td>0.00300192</td><td></td><td>0.00126551</td><td>0.00236619</td></tr><tr><td>5</td><td></td><td>0.00013569</td><td>0.00034105</td><td></td><td>0.00011576</td><td>0.00028012</td></tr><tr><td>6</td><td></td><td>0.00000878</td><td>0.00002671</td><td></td><td>0.00000762</td><td>0.00002256</td></tr><tr><td>7</td><td></td><td>0.00000038</td><td>0.00000139</td><td></td><td>0.00000034</td><td>0.00000119</td></tr><tr><td>8</td><td></td><td>0.00000001</td><td>0.00000003</td><td></td><td>0.00000002</td><td>0.00000004</td></tr><tr><td>9</td><td></td><td>0.00000000</td><td>0.00000000</td><td></td><td>0.00000000</td><td>0.00000000</td></tr><tr><td>Sum</td><td>0.6250565</td><td>0.16122245</td><td>0.21372110</td><td>0.8264325</td><td>0.10204305</td><td>0.07152445</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">Exponential inter-arrival time distribution</td></tr><tr><td>0</td><td>0.6413016</td><td></td><td></td><td>0.6413016</td><td></td><td></td></tr><tr><td>1</td><td>0.2238212</td><td>0.02989575</td><td>0.02623701</td><td>0.2238212</td><td>0.02989575</td><td>0.02623701</td></tr><tr><td>2</td><td></td><td>0.04733494</td><td>0.01944150</td><td></td><td>0.04733494</td><td>0.01944150</td></tr><tr><td>3</td><td></td><td>0.00614876</td><td>0.00444061</td><td></td><td>0.00614876</td><td>0.00444061</td></tr><tr><td>4</td><td></td><td>0.00061891</td><td>0.00064103</td><td></td><td>0.00061891</td><td>0.00064103</td></tr><tr><td>5</td><td></td><td>0.00004733</td><td>0.00006390</td><td></td><td>0.00004733</td><td>0.00006390</td></tr><tr><td>6</td><td></td><td>0.00000266</td><td>0.00000443</td><td></td><td>0.00000266</td><td>0.00000443</td></tr><tr><td>7</td><td></td><td>0.00000010</td><td>0.00000020</td><td></td><td>0.00000010</td><td>0.00000020</td></tr><tr><td>8</td><td></td><td>0.00000000</td><td>0.00000001</td><td></td><td>0.00000000</td><td>0.00000001</td></tr><tr><td>9</td><td></td><td>0.00000000</td><td>0.00000000</td><td></td><td>0.00000000</td><td>0.00000000</td></tr><tr><td>Sum</td><td>0.8651229</td><td>0.08404845</td><td>0.05082868</td><td>0.8651229</td><td>0.08404845</td><td>0.05082868</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">Erlang-2 inter-arrival time distribution</td></tr><tr><td>0</td><td>0.90496984</td><td></td><td></td><td>0.5950249</td><td></td><td></td></tr><tr><td>1</td><td>0.08299433</td><td>0.00385205</td><td>0.00451111</td><td>0.3586692</td><td>0.01288227</td><td>0.01095993</td></tr><tr><td>2</td><td></td><td>0.00233648</td><td>0.00125206</td><td></td><td>0.01663433</td><td>0.00525185</td></tr><tr><td>3</td><td></td><td>0.00003961</td><td>0.00004322</td><td></td><td>0.00030899</td><td>0.00025859</td></tr><tr><td>4</td><td></td><td>0.00000048</td><td>0.00000082</td><td></td><td>0.00000404</td><td>0.00000577</td></tr><tr><td>5</td><td></td><td>0.00000000</td><td>0.00000001</td><td></td><td>0.00000004</td><td>0.00000008</td></tr><tr><td>6</td><td></td><td>0.00000000</td><td>0.00000000</td><td></td><td>0.00000000</td><td>0.00000000</td></tr><tr><td>7</td><td></td><td>0.00000000</td><td>0.00000000</td><td></td><td>0.00000000</td><td>0.00000000</td></tr><tr><td>8</td><td></td><td>0.00000000</td><td>0.00000000</td><td></td><td>0.00000000</td><td>0.00000000</td></tr><tr><td>9</td><td></td><td>0.00000000</td><td>0.00000000</td><td></td><td>0.00000000</td><td>0.00000000</td></tr><tr><td>Sum</td><td>0.98796417</td><td>0.00622862</td><td>0.00580721</td><td>0.9536941</td><td>0.02982966</td><td>0.01647621</td></tr></tbody></table>

## 6.1. Numerical Results of System Performance Measures

Table 3 and Figures 2–4 report a sensitivity analysis of the principal performance measures with respect to the service rates during normal and working-vacation regimes,
namely $\mu_1, \mu_2$, and $v$. Throughout, inter-arrival times are taken to be deterministic with mean $1/\lambda$. Unless otherwise stated, the baseline parameter vector is

$$ (\lambda, \mu_1, \mu_2, v, \bar{\beta}, \alpha, \phi) = (1.2, 3.3, 2.5, 1.6, 0.2, 1.1, 5.0), $$

and the balking function is

$$ \bar{\theta}_i = 1 - \frac{1}{i+1}, \quad 2 \le i \le N-1, \quad N = 7. $$

<table><caption>**Table 3.** Effect of $\mu_1, \mu_2$, and $v$ on system characteristics.</caption><thead><tr><th></th><th colspan="3">μ<sub>1</sub> variation (with balking)</th><th colspan="3">μ<sub>1</sub> variation (with balking)</th></tr><tr><th></th><th>μ<sub>1</sub> = 3.0</th><th>μ<sub>1</sub> = 3.3</th><th>μ<sub>1</sub> = 3.6</th><th>μ<sub>1</sub> = 3.0</th><th>μ<sub>1</sub> = 3.3</th><th>μ<sub>1</sub> = 3.6</th></tr></thead><tbody><tr><th>L<sub>s</sub></th><td>0.4755609</td><td>0.4556073</td><td>0.4355588</td><td>0.5991304</td><td>0.5583294</td><td>0.5217978</td></tr><tr><th>W<sub>s</sub></th><td>0.3963007</td><td>0.3796727</td><td>0.3629657</td><td>0.4993088</td><td>0.4652965</td><td>0.4348463</td></tr><tr><th>J<sub>s</sub></th><td>1.0585704</td><td>1.0736609</td><td>1.0862906</td><td>1.1999195</td><td>1.1999433</td><td>1.1999592</td></tr><tr><th>B<sub>r</sub></th><td>0.1414296</td><td>0.1263391</td><td>0.1137094</td><td>0.0000805</td><td>0.0000567</td><td>0.0000408</td></tr><tr><th>R<sub>r</sub></th><td>0.0206344</td><td>0.0172738</td><td>0.0146332</td><td>0.1027059</td><td>0.0828194</td><td>0.0679238</td></tr></tbody><thead><tr><th></th><th colspan="3">μ<sub>2</sub> variation (with balking)</th><th colspan="3">μ<sub>2</sub> variation (without balking)</th></tr><tr><th></th><th>μ<sub>2</sub> = 2.5</th><th>μ<sub>2</sub> = 2.8</th><th>μ<sub>2</sub> = 3.1</th><th>μ<sub>2</sub> = 2.5</th><th>μ<sub>2</sub> = 2.8</th><th>μ<sub>2</sub> = 3.1</th></tr></thead><tbody><tr><th>L<sub>s</sub></th><td>0.4556073</td><td>0.4499070</td><td>0.4431046</td><td>0.5583294</td><td>0.5342463</td><td>0.5140029</td></tr><tr><th>W<sub>s</sub></th><td>0.3796727</td><td>0.3749225</td><td>0.3692539</td><td>0.4652965</td><td>0.4452205</td><td>0.4283468</td></tr><tr><th>J<sub>s</sub></th><td>1.0736609</td><td>1.0849904</td><td>1.0929703</td><td>1.1999433</td><td>1.1999590</td><td>1.1999692</td></tr><tr><th>B<sub>r</sub></th><td>0.1263391</td><td>0.1150096</td><td>0.1070297</td><td>0.0000567</td><td>0.0000410</td><td>0.0000308</td></tr><tr><th>R<sub>r</sub></th><td>0.0172738</td><td>0.0149741</td><td>0.0133365</td><td>0.0828194</td><td>0.0679570</td><td>0.0578190</td></tr></tbody><thead><tr><th></th><th colspan="3">v variation (with balking)</th><th colspan="3">v variation (without balking)</th></tr><tr><th></th><th>v = 1.6</th><th>v = 1.9</th><th>v = 2.2</th><th>v = 1.6</th><th>v = 1.9</th><th>v = 2.2</th></tr></thead><tbody><tr><th>L<sub>s</sub></th><td>0.4556073</td><td>0.4513631</td><td>0.4472588</td><td>0.5583294</td><td>0.5540824</td><td>0.5499365</td></tr><tr><th>W<sub>s</sub></th><td>0.3796727</td><td>0.3761359</td><td>0.3727157</td><td>0.4652965</td><td>0.4617568</td><td>0.4583014</td></tr><tr><th>J<sub>s</sub></th><td>1.0736609</td><td>1.0763325</td><td>1.0788873</td><td>1.1999433</td><td>1.1999443</td><td>1.1999452</td></tr><tr><th>B<sub>r</sub></th><td>0.1263391</td><td>0.1236675</td><td>0.1211127</td><td>0.0000567</td><td>0.0000557</td><td>0.0000548</td></tr><tr><th>R<sub>r</sub></th><td>0.0172738</td><td>0.0167911</td><td>0.0163336</td><td>0.0828194</td><td>0.0814966</td><td>0.0802134</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Effects of $\mu_1, \mu_2$, and $v$ on $\Pi_B$ and $\Pi_W$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 3. Effects of $\mu_2$, $\mu_1$, and $v$ on $\Pi_B$ and $\Pi_W$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4. Effects of $v$, $\mu_2$, and $\mu_1$ on $\Pi_B$ and $\Pi_W$.</figcaption></figure>

Main findings. The numerical results show that increasing any of the service rates $\mu_1$, $\mu_2$, or $v$ leads to a monotone decrease in the congestion and impatience indicators:

$$L_s \downarrow, W_s \downarrow, B_r \downarrow, R_r \downarrow,$$

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 $\Pi_B$, decreases as each of $\mu_1$, $\mu_2$, and $v$ increases. During working vacations, the busy probability of Server 2, denoted $\Pi_W$, decreases with $\mu_1$ and with $v$, but increases with $\mu_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, \Pi_B$ are all larger in the no-balking system,
while

$B_r$ (trivially) and $\Pi_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 Figures 5–7, we analyze the impact of feedback probability $\bar{\beta}$, arrival rate $\lambda$, working vacation rate $\phi$, and impatience rate $\alpha$ on various performance measures. Here, we assume that the inter-arrival times follow an Erlang-2 distribution. The balking function is taken as $\bar{\theta}_i = i/N^2$, for $2 \le i \le N-1$, where $N = 7$. The parameters are taken as $\lambda = 1.0$, $\mu_1 = 2.8$, $\mu_2 = 2.4$, $\nu = 1.8$, $\alpha = 0.9$, and $\phi = 2.0$. From Table 4 and Figures 5–7, we observe the following:

* With fixed values of $\bar{\beta}$, $\phi$, and $\alpha$, an increase in $\lambda$ 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 $\Pi_L$ significantly rises. Notably, the probability that Server 2 is idle during the working vacation period $\Pi_I$ decreases with $\lambda$.

* With fixed values of $\bar{\beta}$, $\lambda$, and $\alpha$, when the working vacation rate $\phi$ 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 $\Pi_L$ decrease. This implies an increase in $J_s$, and in $\Pi_I$ because the mean working vacation time $1/\phi$ decreases. This trend matches absolutely with the realistic situation.

* With fixed values of $\bar{\beta}$, $\lambda$, and $\phi$, when the impatience rate $\alpha$ increases, the system characteristics $J_s$, $R_r$, and $\Pi_I$ increase, while $\Pi_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 $\alpha$, $\phi$, and $\lambda$, an increase in the feedback probability $\bar{\beta}$ 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 $\Pi_L$. Conversely, $J_s$ and the probability that Server 2 is idle during the working vacation period $\Pi_I$ monotonically decrease.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Effects of $\lambda$, $\bar{\beta}$, and $\phi$ on $\Pi_L$ and $\Pi_I$.</figcaption></figure>
<table><caption>**Table 4.** Effects of $\lambda$, $\phi$, $\alpha$, and $\bar{\beta}$ on system characteristics.</caption><thead><tr><th></th><th colspan="2">β̄ = 0.3</th><th colspan="2">β̄ = 0.5</th><th colspan="2">β̄ = 0.7</th></tr><tr><th></th><th>λ = 0.7</th><th>λ = 1.0</th><th>λ = 0.7</th><th>λ = 1.0</th><th>λ = 0.7</th><th>λ = 1.0</th></tr></thead><tbody><tr><th>L<sub>s</sub></th><td>0.3649478</td><td>0.5285393</td><td>0.5156073</td><td>0.7485630</td><td>0.8706453</td><td>1.2585121</td></tr><tr><th>W<sub>s</sub></th><td>0.5213540</td><td>0.5285393</td><td>0.7365818</td><td>0.7485630</td><td>1.2437789</td><td>1.2585122</td></tr><tr><th>J<sub>s</sub></th><td>0.6995574</td><td>0.9984455</td><td>0.6989834</td><td>0.9965871</td><td>0.6967290</td><td>0.9899939</td></tr><tr><th>B<sub>r</sub></th><td>0.0004426</td><td>0.0015545</td><td>0.0010166</td><td>0.0034129</td><td>0.0032710</td><td>0.0100061</td></tr><tr><th>R<sub>r</sub></th><td>0.0003104</td><td>0.0013975</td><td>0.0011433</td><td>0.0047919</td><td>0.0068955</td><td>0.0251020</td></tr></tbody><thead><tr><th></th><th colspan="2">β̄ = 0.3</th><th colspan="2">β̄ = 0.5</th><th colspan="2">β̄ = 0.7</th></tr><tr><th></th><th>φ = 2.0</th><th>φ = 4.0</th><th>φ = 2.0</th><th>φ = 4.0</th><th>φ = 2.0</th><th>φ = 4.0</th></tr></thead><tbody><tr><th>L<sub>s</sub></th><td>0.5285393</td><td>0.5253591</td><td>0.7485630</td><td>0.7430044</td><td>1.2585121</td><td>1.2486714</td></tr><tr><th>W<sub>s</sub></th><td>0.5285393</td><td>0.5253591</td><td>0.7485630</td><td>0.7430044</td><td>1.2585122</td><td>1.2486714</td></tr><tr><th>J<sub>s</sub></th><td>0.9984455</td><td>0.9984870</td><td>0.9965871</td><td>0.9966743</td><td>0.9899939</td><td>0.9901889</td></tr><tr><th>B<sub>r</sub></th><td>0.0015545</td><td>0.0015130</td><td>0.0034129</td><td>0.0033257</td><td>0.0100061</td><td>0.0098111</td></tr><tr><th>R<sub>r</sub></th><td>0.0013975</td><td>0.0013191</td><td>0.0047919</td><td>0.0045670</td><td>0.0251020</td><td>0.0243441</td></tr></tbody><thead><tr><th></th><th colspan="2">β̄ = 0.3</th><th colspan="2">β̄ = 0.5</th><th colspan="2">β̄ = 0.7</th></tr><tr><th></th><th>α = 0.6</th><th>α = 0.9</th><th>α = 0.6</th><th>α = 0.9</th><th>α = 0.6</th><th>α = 0.9</th></tr></thead><tbody><tr><th>L<sub>s</sub></th><td>0.5288857</td><td>0.5285393</td><td>0.7501344</td><td>0.7485630</td><td>1.2703945</td><td>1.2585121</td></tr><tr><th>W<sub>s</sub></th><td>0.5288857</td><td>0.5285393</td><td>0.7501344</td><td>0.7485630</td><td>1.2703948</td><td>1.2585122</td></tr><tr><th>J<sub>s</sub></th><td>0.9984379</td><td>0.9984455</td><td>0.9965519</td><td>0.9965871</td><td>0.9897164</td><td>0.9899939</td></tr><tr><th>B<sub>r</sub></th><td>0.0015621</td><td>0.0015545</td><td>0.0034481</td><td>0.0034129</td><td>0.0102836</td><td>0.0100061</td></tr><tr><th>R<sub>r</sub></th><td>0.0010108</td><td>0.0013975</td><td>0.0035746</td><td>0.0047919</td><td>0.0199789</td><td>0.0251020</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Effects of $\phi$, $\bar{\beta}$, and $\alpha$ on $\Pi_L$ and $\Pi_I$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Effects of $\alpha$, $\bar{\beta}$, and $\lambda$ on $\Pi_L$ and $\Pi_I$.</figcaption></figure>
## 6.2. Numerical Results of Cost Optimization

Our aim in this subsection is to determine the optimal service rates $\mu_1, \mu_2$, and $\nu$ that minimize the total expected cost function $F(\mu_1, \mu_2, \nu)$ 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{\theta}_i = \frac{i^2}{N^2}, \quad 2 \le i \le N-1,$$

with the system capacity fixed at $N = 5$. Unless otherwise specified, the parameters

$$\bar{\beta} = 0.3, \quad \phi = 0.2, \quad \lambda = 1.8, \quad \mu_2 = 0.9, \quad \nu = 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 Figures 8 and 9, the convexity and optimality of the total expected cost function $F(\mu_1, \mu_2, \nu)$ are clearly observed. These figures illustrate that there exist specific values of the service rates $\mu_1, \mu_2$, and $\nu$ that minimize the total expected cost function for the selected model parameters.

According to Table 5, as the arrival rate $\lambda$ increases, we observe that the optimal service rates at the primary and secondary stations, $\mu_1^*$ and $\mu_2^*$, increase, whereas the auxiliary-rate $\nu^*$ decreases. Consequently, the optimal expected cost

$$F(\mu_1^*, \mu_2^*, \nu^*)$$

is strictly increasing in $\lambda$. 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 $\alpha$ increases, all three optimal rates $\mu_1^*, \mu_2^*$, and $\nu^*$ rise, and so does the objective $F(\mu_1^*, \mu_2^*, \nu^*)$. 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 $\mu_1^*, \mu_2^*$, and $\nu^*$, together with the resulting optimal cost $F(\mu_1^*, \mu_2^*, \nu^*)$, 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 $\mu_1, \mu_2$, and $\nu$ can hedge against higher arrival intensities and customer impatience, while explicitly accounting for balking to avoid systematic under-provision of capacity.

<table><caption>**Table 5.** Optimal values $\mu_1^*, \mu_2^*, \nu^*$, and $F(\mu_1^*, \mu_2^*, \nu^*)$ for different $\lambda$ and $\alpha$.</caption><thead><tr><th rowspan="2">λ</th><th rowspan="2"></th><th colspan="2">With Balking</th><th colspan="2">Without Balking</th></tr><tr><th>α = 0.5</th><th>α = 1.5</th><th>α = 0.5</th><th>α = 1.5</th></tr></thead><tbody><tr><td rowspan="4">1.2</td><td>μ<sub>1</sub><sup>*</sup></td><td>1.6220068</td><td>1.9243682</td><td>1.2092605</td><td>1.3692950</td></tr><tr><td>μ<sub>2</sub><sup>*</sup></td><td>0.5429882</td><td>0.5702495</td><td>0.5018050</td><td>0.5213838</td></tr><tr><td>ν<sup>*</sup></td><td>0.5206338</td><td>0.5651892</td><td>0.4405885</td><td>0.4677076</td></tr><tr><td>F(μ<sub>1</sub><sup>*</sup>, μ<sub>2</sub><sup>*</sup>, ν<sup>*</sup>)</td><td>85.2775457</td><td>88.1958150</td><td>82.1100803</td><td>83.8265865</td></tr></tbody></table>
<table><caption>Table 5. *Cont.*</caption><thead><tr><th rowspan="2">λ</th><th rowspan="2"></th><th colspan="2">With Balking</th><th colspan="2">Without Balking</th></tr><tr><th>α = 0.5</th><th>α = 1.5</th><th>α = 0.5</th><th>α = 1.5</th></tr></thead><tbody><tr><td rowspan="4">1.5</td><td>μ<sub>1</sub><sup>*</sup></td><td>1.8395036</td><td>2.2576245</td><td>1.2879081</td><td>1.5292813</td></tr><tr><td>μ<sub>2</sub><sup>*</sup></td><td>0.5598252</td><td>0.5966068</td><td>0.5057864</td><td>0.5313989</td></tr><tr><td>ν<sup>*</sup></td><td>0.4590795</td><td>0.5475410</td><td>0.3285328</td><td>0.3848557</td></tr><tr><td>F(μ<sub>1</sub><sup>*</sup>, μ<sub>2</sub><sup>*</sup>, ν<sup>*</sup>)</td><td>93.3784508</td><td>97.3268385</td><td>89.4158168</td><td>91.8879399</td></tr><tr><td rowspan="4">1.8</td><td>μ<sub>1</sub><sup>*</sup></td><td>2.0449603</td><td>2.5894097</td><td>1.3117024</td><td>1.6659500</td></tr><tr><td>μ<sub>2</sub><sup>*</sup></td><td>0.5619796</td><td>0.6001498</td><td>0.5079521</td><td>0.5444188</td></tr><tr><td>ν<sup>*</sup></td><td>0.3384183</td><td>0.4609499</td><td>0.1290864</td><td>0.1587587</td></tr><tr><td>F(μ<sub>1</sub><sup>*</sup>, μ<sub>2</sub><sup>*</sup>, ν<sup>*</sup>)</td><td>100.7312278</td><td>105.7525288</td><td>95.9237311</td><td>99.2935054</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>Figure 8. Effects of $μ_1$ and $μ_2$ on the total expected cost function.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 9. Effects of $μ_1$ and $ν$ on the total expected cost function.</figcaption></figure>
## 6.3. Discussions

**Remark 3.** On the compact feasible set $\{v_{\min} \le v \le \mu_2 \le \mu_1 \le \mu_{1,\max}\}$, the mapping $(\mu_1, \mu_2, v) \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\Pi_I + C_2\Pi_W + C_3\Pi_B$ price server phases (idle/WV/regular), $C_4L_s$ penalizes congestion, $C_5(B_r + R_r)$ penalizes behaviorally lost work, $C_6\Lambda\Pi_L$ accounts for blocked demand, and $(C_7 + \beta C_8)(\mu_1 + \mu_2 + v)$ proxies energy/labor with rework overhead.

We initialize bat uniformly over the feasible cone $v \le \mu_2 \le \mu_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 < v < \mu_2 < \mu_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 $\delta_i, \theta_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 ($v < \mu_2 < \mu_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 $v$ (for example, under a low-power or maintenance mode). The other machine (Server 1) remains fully available for MTO work at rate $\mu_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 $\mu_2 > v$. 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 $\theta_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 $\alpha$ while in the queue, following standard models of queues with abandonment [81–83]. Finally, we incorporate *Bernoulli feedback* at service completion: with probability $1 - \beta$, a job is routed back to the queue for rework (e.g., remedial processing or quality corrections), whereas with probability $\beta$, 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 $v$ 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 $GI/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 $GI/PH/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–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–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

$$ \mathcal{L}V(i, j) \leq \lambda \theta_i - (\mu_1 + \mu_2 \mathbb{1}_{\{j=2\}} + \nu \mathbb{1}_{\{j=1\}}) + (i - 2)^+ \alpha, $$

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.

**Funding:** This work was supported by ongoing institutional funding. No additional grants were obtained to carry out or direct this particular research.

**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

<table><tr><td>λ</td><td>External arrival rate; mean inter-arrival time 1/λ.</td></tr><tr><td>A(u), a(u), A<sup>*</sup>(s)</td><td>Inter-arrival cdf, density, and Laplace–Stieltjes transform.</td></tr><tr><td>N</td><td>Total capacity (in service + waiting).</td></tr><tr><td>θ<sub>i</sub>, θ&#773;<sub>i</sub> = 1 - θ<sub>i</sub></td><td>Join/balking probability on seeing i customers.</td></tr><tr><td>α</td><td>Exponential impatience (reneging) rate per waiting customer.</td></tr><tr><td>β, β&#773; = 1 - β</td><td>Departure vs. feedback probability at service completion.</td></tr><tr><td>μ<sub>1</sub>, μ<sub>2</sub></td><td>Regular service rates of Server 1 and Server 2 (μ<sub>2</sub> ≤ μ<sub>1</sub>).</td></tr><tr><td>ν</td><td>Server 2 service rate during a working vacation.</td></tr><tr><td>φ</td><td>Working-vacation termination rate.</td></tr></table>

### State processes and stationary objects

<table><tr><td>L(t)</td><td>Total customers in system at time t (queue + in service).</td></tr><tr><td>J(t) ∈ {0, 1, 2}</td><td>Server 2 state: 0 = idle on WV; 1 = busy on WV; 2 = busy (regular).</td></tr><tr><td>W(t)</td><td>Remaining inter-arrival time at t.</td></tr><tr><td>P<sub>i,j</sub>(u, t)</td><td>P(L(t) = i, J(t) = j, u ≤ W(t) < u + du) (density in u).</td></tr><tr><td>P<sub>i,j</sub>(u)</td><td>Stationary version of P<sub>i,j</sub>(u, t) as t → ∞.</td></tr><tr><td>P<sup>*</sup><sub>i,j</sub>(s)</td><td>LST of P<sub>i,j</sub>(u) in u.</td></tr><tr><td>P<sub>i,j</sub>(0)</td><td>Boundary (rate) probability at W = 0 (arrival instant).</td></tr><tr><td>P<sub>i,j</sub></td><td>Arbitrary-epoch stationary probability of (L, J) = (i, j).</td></tr><tr><td>P&#773;<sub>i,j</sub></td><td>Pre-arrival (Palm) probability; P&#773;<sub>i,j</sub> = λ<sup>-1</sup>P<sub>i,j</sub>(0).</td></tr><tr><td>δ<sub>i</sub></td><td>β(μ<sub>1</sub> + ν) + φ + (i - 2)α (WV balance parameter).</td></tr><tr><td>θ<sub>i</sub></td><td>β(μ<sub>1</sub> + μ<sub>2</sub>) + (i - 2)α (regular balance parameter).</td></tr></table>

### Performance measures

<table><tr><td>Π<sub>I</sub>, Π<sub>W</sub>, Π<sub>B</sub></td><td>Fractions of time: Server 2 idle on WV / busy on WV / busy regular.</td></tr><tr><td>Π<sub>L</sub></td><td>Loss probability on arrival: P&#773;<sub>N,1</sub> + P&#773;<sub>N,2</sub>.</td></tr><tr><td>L<sub>s</sub></td><td>Mean system size: P<sub>1,0</sub> + Σ<sup>N</sup><sub>i=1</sub> i (P<sub>i,1</sub> + P<sub>i,2</sub>).</td></tr><tr><td>W<sub>s</sub></td><td>Mean sojourn time: W<sub>s</sub> = λ(1 - Π<sub>L</sub>).</td></tr><tr><td>J<sub>s</sub></td><td>Effective joining rate.</td></tr><tr><td>B<sub>r</sub>, R<sub>r</sub></td><td>Average rates of balking and reneging.</td></tr></table>

### Cost model

<table><tr><td>C<sub>1</sub>, ..., C<sub>9</sub></td><td>Unit cost coefficients used in the economic objective.</td></tr></table>

## References

1. Erlang, A.K. Sandsynlighedsregning og Telefonsamtaler. *Nyt Tidss. Math.* **1909**, 20, 33–39.

2. Kingman, J.F.C. The first Erlang century—and the next. *Queueing Syst.* **2009**, 63, 3–12. [CrossRef]

3. Kleinrock, L. *Queueing Systems, Volume II: Computer Applications*; Wiley: New York, NY, USA, 1976.

4. Wolff, R.W. *Stochastic Modeling and the Theory of Queues*; Prentice-Hall International, Inc.: Englewood Cliffs, NJ, USA, 1989.

5. Asmussen, S. *Applied Probability and Queues*, 2nd ed.; Springer: New York, NY, USA, 2003; Volume 51. [CrossRef]

6. Baccelli, F.; Brémaud, P. *Elements of Queueing Theory: Palm Martingale Calculus and Stochastic Recurrences*, 2nd ed.; Springer: New York, NY, USA, 2003; Volume 26. [CrossRef]

7. Takagi, H. *Queueing Analysis, Volume 1: Vacation and Priority Systems*; North-Holland: Amsterdam, The Netherlands, 1991.
8. Shortle, J.F.; Thompson, J.M.; Gross, D.; Harris, C.M. *Fundamentals of Queueing Theory*, 5th ed.; John Wiley & Sons: Hoboken, NJ, USA, 2018. [CrossRef]

9. Borovkov, A.A. *Stochastic Processes in Queueing Theory*; Applications of Mathematics, No. 4; Springer: New York, NY, USA, 1976; p. xi+280; Translated from the Russian by Kenneth Wickwire. [CrossRef]

10. Khintchine, A.Y. *Mathematical Methods in the Theory of Queueing*; Griffin's Statistical Monographs & Courses, No. 7; Andrews, D.M., Quenouille, M.H., Translators; Hafner Publishing Co.: New York, NY, USA, 1960; p. 120.

11. Little, J.D.C. A Proof for the Queueing Formula $L = \lambda W$. *Oper. Res.* **1961**, *9*, 383–387. [CrossRef]

12. Bouchentouf, A.A.; Laxmi, P.V.; Ziad, I. M/M/c/DWV Queueing Model with Bernoulli Schedule Working Vacation Interruption, Waiting Servers, and Impatience as a Model of a Call Center. *Am. J. Math. Manag. Sci.* **2025**, *43*, 173–192. [CrossRef]

13. Yang, X.; Zhang, Y.; Wang, B.; Li, X.J. Performance Analysis and Cost Optimization of the M/M/1/N Queueing System with Working Vacation and Working Breakdown. *Mathematics* **2025**, *13*, 2980. [CrossRef]

14. Singh, M.; Jain, M.; Singh, V.P. Analytical and Computational Study of Bi-level Vacation Policy for Transient Markovian Queue with State-Dependent Rates, Waiting Server, and Discouragement. *Commun. Appl. Math. Comput.* **2025**, *1–26*. [CrossRef]

15. Dhibar, S.; Jain, M. Metaheuristics and strategic behavior of markovian retrial queue under breakdown, vacation and bernoulli feedback. *Appl. Intell.* **2025**, *55*, 273. [CrossRef]

16. Tamuli, A.; Das, D.; Choudhury, A.; Kushvaha, B. Optimal service design for heterogeneous queueing system with reverse balking and reneging. *Oper. Res.* **2025**, *25*, 33. [CrossRef]

17. Wu, C.H.; Yang, D.Y.; Lin, W.L. Modeling and cost optimization of controlling customer arrivals in queueing systems with two waiting areas. *J. Ind. Manag. Optim.* **2025**, *21*, 1251–1268. [CrossRef]

18. Mathavavisakan, N.M.; Indhira, K.; Ammar, S.I. Cost optimization and analysis of retrial queues with hybrid vacation policies and fault correction: A metaheuristic approach. *RAIRO-Oper. Res.* **2025**, *59*, 2683–2720. [CrossRef]

19. Ke, J.C.; Chang, F.M.; Liu, T.H. Bi-objective optimization of a retrial model with synchronous working vacation interruption. *Qual. Technol. Quant. Manag.* **2025**, *22*, 659–682. [CrossRef]

20. Bouchentouf, A.A.; Kumar, K.; Chahal, P.K. Particle Swarm Optimization for a redundant repairable machining system with working vacations and impatience in a multi-phase random environment. *Swarm Evol. Comput.* **2024**, *90*, 101688. [CrossRef]

21. Chahal, P.K.; Kumar, K.; Soodan, B.S. Grey wolf algorithm for cost optimization of cloud computing repairable system with N-policy, discouragement and two-level Bernoulli feedback. *Math. Comput. Simul.* **2024**, *225*, 545–569. [CrossRef]

22. Haghighi, A.M.; Mishev, D.P. *Queueing Models in Industry and Business*; Nova Publishers: New York, NY, USA, 2008; LCCN 2007043239.

23. Newell, G.F. *Applications of Queueing Theory*; Monographs on Applied Probability and Statistics; Chapman and Hall, Ltd.: London, UK, 1971; p. x+148.

24. Haghighi, A.M.; Mishev, D.P. *Delayed and Network Queues*; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2016; p. xiii+400.

25. Ancker, C.J., Jr.; Gafarian, A.V. Some queuing problems with balking and reneging. II. *Oper. Res.* **1963**, *11*, 928–937. [CrossRef]

26. Ancker, C.J., Jr.; Gafarian, A.V. Some queuing problems with balking and reneging. I. *Oper. Res.* **1963**, *11*, 88–100. [CrossRef]

27. Abou-El-Ata, M. The state-dependent queue: M/M/1/N with reneging and general balk functions. *Microelectron. Reliab.* **1991**, *31*, 1001–1007. [CrossRef]

28. Wang, K.H.; Chang, Y.C. Cost analysis of a finite M/M/R queueing system with balking, reneging, and server breakdowns. *Math. Methods Oper. Res.* **2002**, *56*, 169–180. [CrossRef]

29. Ammar, S.I.; Helan, M.M.; Al Amri, F.T. The busy period of an M/M/1 queue with balking and reneging. *Appl. Math. Model.* **2013**, *37*, 9223–9229. [CrossRef]

30. Servi, L.; Finn, S. M/M/1 queues with working vacations (M/M/1/WV). *Perform. Eval.* **2002**, *50*, 41–52. [CrossRef]

31. Xu, X.; Liu, M.; Zhao, X. The bulk input $M^X/M/1$ queue with working vacations. *J. Syst. Sci. Syst. Eng.* **2009**, *18*, 358–368. [CrossRef]

32. Baba, Y. The $M^X/M/1$ Queue with Multiple Working Vacation. *Am. J. Oper. Res.* **2012**, *2*, 217–224. [CrossRef]

33. Jain, M.; Agrawal, P.K. $M/E_k/1$ queueing system with working vacation. *Qual. Technol. Quant. Manag.* **2007**, *4*, 455–470. [CrossRef]

34. Li, J.h.; Tian, N.s.; Zhang, Z.G.; Luh, H.P. Analysis of the $M/G/1$ queue with exponentially working vacations—A matrix analytic approach. *Queueing Syst.* **2009**, *61*, 139–166. [CrossRef]

35. Zhang, M.; Hou, Z. $M/G/1$ queue with single working vacation. *J. Appl. Math. Comput.* **2012**, *39*, 221–234. [CrossRef]

36. Lin, C.H.; Ke, J.C. Multi-server system with single working vacation. *Appl. Math. Model.* **2009**, *33*, 2967–2977. [CrossRef]

37. Xu, X.; Geng, J.; Liu, M.; Guo, H. Stationary analysis for the fluid model driven by the $M/M/c$ working vacation queue. *J. Math. Anal. Appl.* **2013**, *403*, 423–433. [CrossRef]

38. Wu, C.H.; Wang, K.H.; Ke, J.C.; Ke, J.B. A heuristic algorithm for the optimization of M/M/s queue with multiple working vacations. *J. Ind. Manag. Optim.* **2012**, *8*, 1–17. [CrossRef]
39. Majid, S.; Manoharan, P. Analysis of a $M/M/c$ queue with single and multiple synchronous working vacations. *Appl. Appl. Math.* **2017**, *12*, 671–694.

40. Li, J.H.; Tian, N.S.; Liu, W.Y. Discrete-time GI/Geo/1 queue with multiple working vacations. *Queueing Syst.* **2007**, *56*, 53–63.
[CrossRef]

41. Li, J.h.; Tian, N.s. Analysis of the discrete time Geo/Geo/1 queue with single working vacation. *Qual. Technol. Quant. Manag.* **2008**, *5*, 77–89. [CrossRef]

42. Baba, Y. Analysis of a GI/M/1 queue with multiple working vacations. *Oper. Res. Lett.* **2005**, *33*, 201–209. [CrossRef]

43. Banik, A.; Gupta, U.; Pathak, S. On the GI/M/1/N queue with multiple working vacations—Analytic analysis and computation. *Appl. Math. Model.* **2007**, *31*, 1701–1710. [CrossRef]

44. Li, J.H.; Tian, N.S.; Ma, Z.Y. Performance analysis of GI/M/1 queue with working vacations and vacation interruption. *Appl. Math. Model.* **2008**, *32*, 2715–2730. [CrossRef]

45. Banik, A.D. Analysis of single working vacation in $GI/M/1/N$ and $GI/M/1/∞$ queueing systems. *Int. J. Oper. Res.* **2010**, *7*, 314–333. [CrossRef]

46. Zhang, M.; Hou, Z. Steady state analysis of the GI/M/1/N queue with a variant of multiple working vacations. *Comput. Ind. Eng.* **2011**, *61*, 1296–1301. [CrossRef]

47. Yue, D.; Yue, W.; Xu, G. Analysis of customers' impatience in an $M/M/1$ queue with working vacations. *J. Ind. Manag. Optim.* **2012**, *8*, 895–908. [CrossRef]

48. Selvaraju, N.; Goswami, C. Impatient customers in an M/M/1 queue with single and multiple working vacations. *Comput. Ind. Eng.* **2013**, *65*, 207–215. [CrossRef]

49. Laxmi, P.V.; Goswami, V.; Jyothsna, K. Analysis of finite buffer Markovian queue with balking, reneging and working vacations. *Int. J. Strateg. Decis. Sci.* **2013**, *4*, 1–24. [CrossRef]

50. Bouchentouf, A.A.; Guendouzi, A.; Majid, S. On impatience in Markovian $M/M/1/N/DWV$ queue with vacation interruption. *Croat. Oper. Res. Rev. CRORR* **2020**, *11*, 21–37. [CrossRef]

51. Goswami, V. A discrete-time queue with balking, reneging, and working vacations. *Int. J. Stoch. Anal.* **2014**, *8*, 358529. [CrossRef]

52. Laxmi, P.V.; Rajesh, P. Analysis of variant working vacations queue with customer impatience. *Int. J. Manag. Sci. Eng. Manag.* **2017**, *12*, 186–195. [CrossRef]

53. Bouchentouf, A.A.; Guendouzi, A.; Meriem, H.; Shakir, M. Analysis of a single server queue in a multi-phase random environment with working vacations and customers' impatience. *Oper. Res. Decis.* **2022**, *32*, 16–33. [CrossRef]

54. Bouchentouf, A.A.; Houalef, M.; Guendouzi, A. Analysis of a batch arrival multi-server queueing system with waiting servers, synchronous working vacations and impatient customers. *Acta Univ. Sapientiae Math.* **2022**, *14*, 37–60. [CrossRef]

55. Majid, S.; Bouchentouf, A.A.; Guendouzi, A. Analysis and optimisation of a $M/M/1/WV$ queue with Bernoulli schedule vacation interruption and customer's impatience. *Acta Univ. Sapientiae Math.* **2021**, *13*, 367–395. [CrossRef]

56. Takacs, L. A single-server queue with feedback. *Bell Syst. Tech. J.* **1963**, *42*, 505–519. [CrossRef]

57. Thangaraj, V.; Vanitha, S. On the analysis of $M/M/1$ feedback queue with catastrophes using continued fractions. *Int. J. Pure Appl. Math.* **2009**, *53*, 133–151.

58. Kumar, R.; Jain, N.K.; Som, B.K. Optimization of an $M/M/1/N$ feedback queue with retention of reneged customers. *Oper. Res. Decis.* **2014**, *24*, 45–58. [CrossRef]

59. Ayyappan, G.; Shyamala, S. Transient solution of an $M^{[X]}/G/1$ queueing model with feedback, random breakdowns, Bernoulli schedule server vacation and random setup time. *Int. J. Oper. Res.* **2016**, *25*, 196–211. [CrossRef]

60. Bouchentouf, A.A.; Guendouzi, A. Sensitivity analysis of feedback multiple vacation queueing system with differentiated vacations, vacation interruptions and impatient customers. *Int. J. Appl. Math. Stat.* **2018**, *57*, 104–121.

61. Bouchentouf, A.A.; Guendouzi, A. Cost optimization analysis for an $M^X/M/c$ vacation queueing system with waiting servers and impatient customers. *SeMA J.* **2019**, *76*, 309–341. [CrossRef]

62. Bouchentouf, A.A.; Guendouzi, A. The $M^X/M/c$ Bernoulli feedback queue with variant multiple working vacations and impatient customers: Performance and economic analysis. *Arab. J. Math.* **2020**, *9*, 309–327. [CrossRef]

63. Bouchentouf, A.A.; Guendouzi, A. Single server batch arrival Bernoulli feedback queueing system with waiting server, K-variant vacations and impatient customers. *Oper. Res. Forum* **2021**, *2*, 23. [CrossRef]

64. Guendouzi, A.; Saadaoui, A. Customer impatience in the $GI/M(k)/1/N$ queue with working vacation: Economic and performance analysis. *Croat. Oper. Res. Rev.* **2025**, *17*, 2026.0002. [CrossRef]

65. Singh, V.P. Two-server Markovian queues with balking: Heterogeneous vs. homogeneous servers. *Oper. Res.* **1970**, *18*, 145–159. [CrossRef]

66. Neuts, M.F.; Takahashi, Y. Asymptotic behavior of the stationary distributions in the $GI/PH/c$ queue with heterogeneous servers. *Z. Wahrsch. Verw. Gebiete* **1981**, *57*, 441–452. [CrossRef]

67. Shawky, A.I. On Poisson machine interference model: $M/M/2/k/N$ with balking, reneging and heterogeneous repairmen. *Commun. Fac. Sci. Univ. Ank. Sér. A1 Math. Stat.* **2004**, *53*, 7–12. [CrossRef]
68. Krishna Kumar, B.; Pavai Madheswari, S. An M/M/2 queueing system with heterogeneous servers and multiple vacations. *Math. Comput. Model.* **2005**, *41*, 1415–1429. [CrossRef]

69. Bouchentouf, A.A.; Messabihi, A. Heterogeneous two-server queueing system with reverse balking and reneging. *Opsearch* **2018**, *55*, 251–267. [CrossRef]

70. Krishnamoorthy, A.; Sreenivasan, C. An M/M/2 queueing system with heterogeneous servers including one with working vacation. *Int. J. Stoch. Anal.* **2012**, *16*, 145867. [CrossRef]

71. Ammar, S.I. Transient analysis of a two-heterogeneous servers queue with impatient behavior. *J. Egypt. Math. Soc.* **2014**, *22*, 90–95. [CrossRef]

72. Kumar, R.; Sharma, S. Transient analysis of a Markovian queuing model with multiple-heterogeneous servers, and customers' impatience. *Opsearch* **2021**, *58*, 540–556. [CrossRef]

73. Vijaya Laxmi, P.; Jyothsna, K. Balking and reneging multiple working vacations queue with heterogeneous servers. *J. Math. Model. Algorithms Oper. Res.* **2015**, *14*, 267–285. [CrossRef]

74. Yang, X. A New Metaheuristic Bat-Inspired Algorithm. In *Nature Inspired Cooperative Strategies for Optimization (NICSO 2010)*; González, J.R., Pelta, D.A., Cruz, C., Terrazas, G., Krasnogor, N., Eds.; Studies in Computational Intelligence; Springer: Berlin/Heidelberg, Germany, 2010; Volume 284, pp. 65–74. [CrossRef]

75. Shekhar, C.; Varshney, S.; Kumar, A. Optimal control of a service system with emergency vacation using bat algorithm. *J. Comput. Appl. Math.* **2020**, *364*, 112332. [CrossRef]

76. Jyothsna, K.; Vijaya Kumar, P.; Vijaya Laxmi, P. Bat algorithm for optimizing a working vacation queue with impatient clients and secondary service. *Opsearch* **2025**, *62*, 777–796. [CrossRef]

77. Cox, D.R. The analysis of non-Markovian stochastic processes by the inclusion of supplementary variables. *Proc. Camb. Philos. Soc.* **1955**, *51*, 433–441. [CrossRef]

78. Li, J.; Tian, N. The M/M/1 queue with working vacations and vacation interruptions. *J. Syst. Sci. Syst. Eng.* **2007**, *16*, 121–127. [CrossRef]

79. Zhang, M.; Hou, Z. Performance analysis of M/G/1 queue with working vacations and vacation interruption. *J. Comput. Appl. Math.* **2010**, *234*, 2977–2985. [CrossRef]

80. Naor, P. The regulation of queue size by levying tolls. *Econometrica* **1969**, *37*, 15–24. [CrossRef]

81. Stanford, R.E. Reneging phenomena in single channel queues. *Math. Oper. Res.* **1979**, *4*, 162–178. [CrossRef]

82. Ward, A.R.; Glynn, P.W. A diffusion approximation for a Markovian queue with reneging. *Queueing Syst.* **2003**, *43*, 103–128. [CrossRef]

83. Ward, A.R.; Glynn, P.W. A diffusion approximation for a GI/GI/1 queue with balking or reneging. *Queueing Syst.* **2005**, *50*, 371–400. [CrossRef]

84. Choi, B.D.; Kim, B.; Choi, S.H. On the M/G/1 Bernoulli feedback queue with multi-class customers. *Comput. Oper. Res.* **2000**, *27*, 269–286. [CrossRef]

85. Choi, B.D.; Kim, B.; Choi, S.H. An M/G/1 queue with multiple types of feedback, gated vacations and FCFS policy. *Comput. Oper. Res.* **2003**, *30*, 1289–1309. [CrossRef]

86. Gross, D.; Shortle, J.F.; Thompson, J.M.; Harris, C.M. *Fundamentals of Queueing Theory*, 4th ed.; Wiley Series in Probability and Statistics; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2008; p. xiv+500. [CrossRef]

87. Latouche, G.; Ramaswami, V. *Introduction to Matrix Analytic Methods in Stochastic Modeling*; ASA-SIAM Series on Statistics and Applied Probability, Society for Industrial and Applied Mathematics (SIAM); American Statistical Association: Alexandria, VA, USA, 1999; p. xiv+334. [CrossRef]

88. Neuts, M.F. A versatile Markovian point process. *J. Appl. Probab.* **1979**, *16*, 764–779. [CrossRef]

89. Heffes, H.; Lucantoni, D.M. A Markov Modulated Characterization of Packetized Voice and Data Traffic. *IEEE J. Sel. Areas Commun.* **1986**, *4*, 856–868. [CrossRef]

90. Lucantoni, D.M. New Results on the Single Server Queue with a Batch Markovian Arrival Process. *Commun. Stat. Stoch. Model.* **1991**, *7*, 1–46. [CrossRef]

91. Green, L.V.; Kolesar, P.J.; Whitt, W. Coping with Time-Varying Demand when Setting Staffing Requirements for a Service System. *Prod. Oper. Manag.* **2007**, *16*, 13–39. [CrossRef]

92. Kelly, F.P. *Reversibility and Stochastic Networks. with a New Preface*, 1979th ed.; Cambridge University Press: Cambridge, UK, 2011.

93. Boucherie, R.J.; van Dijk, N.M. (Eds.) *Queueing Networks. A Fundamental Approach*; Springer: New York, NY, USA, 2011; Volume 154. [CrossRef]

94. Gans, N.; Koole, G.; Mandelbaum, A. Telephone Call Centers: Tutorial, Review, and Research Prospects. *Manuf. Serv. Oper. Manag.* **2003**, *5*, 79–141. [CrossRef]

95. Harchol-Balter, M. *Performance Modeling and Design of Computer Systems. Queueing Theory in Action*; Cambridge University Press: Cambridge, UK, 2013. [CrossRef]
96. Brown, L.; Gans, N.; Mandelbaum, A.; Sakov, A.; Shen, H.; Zeltyn, S.; Zhao, L. Statistical analysis of a telephone call center: A queueing-science perspective. *J. Am. Stat. Assoc.* **2005**, *100*, 36–50. [CrossRef]

97. Garnett, O.; Mandelbaum, A.; Reiman, M.I. Designing a Call Center with Impatient Customers. *Manuf. Serv. Oper. Manag.* **2002**, *4*, 208–227. [CrossRef]

98. Whittle, P. *Restless Bandits: Activity Allocation in a Changing World*; Cambridge University Press: Cambridge, UK, 1988. [CrossRef]

99. Sutton, R.S.; Barto, A.G. *Reinforcement Learning—An Introduction*, 2nd ed.; MIT Press: Cambridge, MA, USA, 2018.

100. Rockafellar, R.T.; Uryasev, S. Optimization of Conditional Value-at-Risk. *J. Risk* **2000**, *2*, 21–41. [CrossRef]

101. Ben-Tal, A.; El Ghaoui, L.; Nemirovski, A. *Robust Optimization*; Princeton Series in Applied Mathematics; Princeton University Press: Princeton, NJ, USA, 2009; p. xxii+542. [CrossRef]

102. Shwartz, A.; Weiss, A. *Large Deviations for Performance Analysis. Queues, Communications, and Computing*, 1995th ed.; Routledge Revivals; Routledge: Milton, FL, USA, 2018. [CrossRef]

103.asmussen, S.r.; Glynn, P.W. *Stochastic Simulation: Algorithms and Analysis*; Stochastic Modelling and Applied Probability; Springer: New York, NY, USA, 2007; Volume 57, p. xiv+476. [CrossRef]

104. Juneja, S.; Shahabuddin, P. Rare-Event Simulation Techniques: An Introduction and Recent Advances. In *Handbooks in Operations Research and Management Science*; Henderson, S.G., Nelson, B.L., Eds.; Elsevier: Amsterdam, The Netherlands, 2006; Volume 13, pp. 291–350. [CrossRef]

105. Sidje, R.B. EXPOKIT: A Software Package for Computing Matrix Exponentials. *ACM Trans. Math. Softw.* **1998**, *24*, 130–156. [CrossRef]

**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.