Analysis of semi-open queueing networks using lost customers approximation with an application to robotic mobile fulfilment systems

Abstract

We consider a semi-open queueing network (SOQN), where one resource from a resource pool is needed to serve a customer. If on arrival of a customer some resource is available, the resource is forwarded to an inner network to complete the customer’s order. If no resource is available, the new customer waits in an external queue until one becomes available (“backordering”). When a resource exits the inner network, it is returned to the resource pool. We develop a new solution approach. In a first step we modify the system such that new arrivals are lost if the resource pool is empty (“lost customers”). We adjust the arrival rate of the modified system such that the throughputs in all nodes of the inner network are pairwise identical to those in the original network. Using queueing theoretical methods, in a second step we reduce this inner network to a two-station system including the resource pool. For this two-station systems, we invert the first step and obtain a standard SOQN which can be solved analytically. We apply our results to storage and delivering systems with robotic mobile fulfilment systems (RMFSs). Instead of sending pickers to the storage area to search for the ordered items and pick them, robots carry shelves with ordered items from the storage area to picking stations. We model the RMFS as an SOQN to determine the minimal number of robots.

1 Introduction

Queueing networks can be classified as follows according to (Chen and Yao 2001, p. 21ff.), Roy (2016), Azadeh et al. (2017b): open queueing networks (OQN), where customers arrive from the exterior, request service at several nodes and then leave the system; closed queueing networks (CQN), without external arrivals and departures, and a fixed number of customers; semi-open queueing networks (SOQN), which have characteristics of both OQNs and CQNs (Jia and Heragu 2009). An SOQN resembles an OQN because customers arrive from the exterior and leave the system after service. It resembles a CQN because there is an overall capacity constraint for the network, which is realized as follows: Any customer needs a resource from an associated resource pool for service. If on arrival of a customer a resource is available, the resource enters the network (which for definiteness will be termed “inner network”) to complete the customer’s order. If there is no resource available, the new customer has to wait in an external queue until one becomes available. When a resource exits the inner network, it returns to the resource pool.

SOQNs are adopted for performance analysis of manufacturing systems and service systems, e.g. in logistics, communication, warehousing, health care (Roy 2016).

We focus on SOQNs where the inner network consists of exponential single server nodes with infinite waiting room under either FCFS regime or processor sharing regime. The literature on SOQNs is overwhelming, so we point only to the most relevant sources for our present investigation. An overview of SOQNs and solution methods is available in Jia and Heragu (2009), Ekren et al. (2014) and Roy (2016). Roy (2016) also compares numerical accuracy of several methods. A recent article, not included in these reviews, is Kim et al. (2018). The most common solution approaches are matrix-geometric method, aggregation method, network decomposition approach, parametric decomposition method and performance bounds. For SOQNs with an inner network consisting of more than one node, closed-form expressions for steady-state distributions are not available.

We develop a new solution approach, visualised in Fig. 1: (i) We modify the system such that new arrivals are lost if the resource pool is empty (“lost customers”). For this modification, closed-form expressions for the steady-state distribution in product form are available. We adjust the arrival rate at the modified system such that the throughputs in all nodes of the resource network ($=$ inner network $+$ resource pool) are pairwise identical to those in the original network. We also prove that idling probabilities for nodes with constant service rates are pairwise identical. Moreover, we provide closed-form expressions for throughputs and for these idling probabilities. (ii) Using queueing theoretical methods (Norton’s Theorem), we reduce this inner network to a two-station system including the resource pool. (iii) For this system, we invert step (i) and obtain a standard SOQN with a two-station resource network and an external queue which can be solved analytically.

Step (ii) resembles the standard application of matrix-geometrical methods (MGM) as surveyed in Roy (2016). The inner network is converted into a two-station network by aggregation. “The challenge...lies in decomposing the original subnetwork...into two subnetworks” (Roy 2016, p. 1742), to be aggregated. Several recipes are provided. While the aggregations use heuristics, our use of Norton’s Theorem is exact because of the product form distribution behind, and we escape from the decomposition problem because we end up with a single node. The MGM as developed and used (e.g.) in Jia and Heragu (2009), Ekren et al. (2014), Lamballais et al. (2017) and Buitenhek et al. (2000) requires as final step a numerical approach to solve for steady-state distributions. Our final step uses steady-state distributions obtained in closed form for computing performance indices.

Convention To distinguish the versions of SOQNs which we elaborate on, we use terms borrowed from inventory theory, where unsatisfied demand can be either “backordered” or is called “lost sales” if backordering is not possible. SOQN-BO denotes a standard SOQN with external queue for waiting customers when the resource pool is empty (BO $=$ backordering), and SOQN-LC an SOQN with lost sales (LC $=$ lost customers). The first version is often called SOQN with infinite buffer, the second one SOQN with finite buffer.

Remark The approximation of BO-systems by LC-systems follows (Krenzler 2016, Section 2.1.4, p. 34ff.). In the literature, results for queueing networks with lost customers are available, e.g. lost sales models in Otten (2018), Schwarz et al. (2006) and Krishnamoorthy et al. (2011). For a review of systems with lost sales, we refer to Bijvank and Vis (2011).

We apply our new solution approach to the investigation of robotic mobile fulfilment systems (RMFSs), such as the Kiva System (nowadays Amazon Robotics), the GreyOrange Butler or the Swisslog CarryPick. RMFSs are a new type of warehousing system which have received attention recently, due to increasing growth in the e-commerce sector. The principle of RMFS (see Fig. 7) is, instead of sending pickers to the storage area to search for ordered items and pick them, to carry shelves—called pods—with ordered items by robots from the storage area to picking stations. At every picking station resides a picker, a person who takes items from the pods and packs them into boxes according to customers’ orders. When the picker does not need the pod any longer, the robot either transports the pod directly back to the storage area or first makes a stopover at a replenishment station. RMFSs pose many hard decision problems at strategic, tactical and operational levels. An overview is included in (Merschformann et al. 2019, Section 4). The literature on RMFS is, similar to that on SOQNs, already overwhelming, so we only point to some references closely related to our investigations. RMFSs are modelled as SOQNs in several articles. For an overview, we refer to (Azadeh et al. 2017a, Section 7.2 and Table 4) and (Azadeh et al. 2017b, Section 6). There articles are classified according to the decision problem of interest and the relevant methodology. Recent overviews are presented in Lamballais et al. (2019) and (Boysen et al. 2019, Section 7). Most articles analyse decision problems different from what we look at in this paper: Decisions on optimal number of robots. A similar problem is tackled in Yuan and Gong (2017) where two protocols are compared, pooled and dedicated robots using OQNs instead of SOQNs, for a model without replenishment station. In Zou et al. (2018), the optimal number of robots is determined for different battery recovery strategies using a nested SOQN. Different types of warehousing systems use also semi-open or open queueing models, e.g. Ekren and Akpunar (2021) and Ekren et al. (2013).

Main contributions of the paper

• We determine closed-form expressions for stability, throughputs and some idling probabilities in an SOQN.

• We develop a new approximation method for an SOQN to determine general performance metrics.

• We model an RMFS as an SOQN and apply our method to calculate the optimal number of robots. We develop a simulation tool for a queueing model of the RMFS to analyse the quality of our approximation method.

Structure of the paper

In Sect. 2, we describe a general SOQN (=SOQN-BO) and analyse its stability. In Sect. 2.3, we determine throughputs and idle probabilities in steady state. In Sect. 3, we introduce our new approximation method for SOQNs, the approach is depicted in Fig. 1. In Sect. 3.1, we analyse a modification of the SOQN-BO, where newly arriving customers are lost, if the resource pool is empty (SOQN-LC). In Sect. 3.2, we adjust the input of the SOQN-LC to meet the throughput of the original SOQN-BO. In Sect. 3.3, we calculate node-throughputs and idle probabilities in steady state. We reduce the complexity of the modified SOQN-LC using Norton’s theorem in Sect. 4.1, and we reinvent the external queue in Sect. 4.2. In Sect. 5, we model an RMFS as an SOQN and apply our results to determine the optimal number of robots. We formulate an algorithm to calculate the minimal number of robots for stability and present numerical examples. Furthermore, we discuss the limitations of the approximation. Section 6 summarises our conclusions.

Fig. 1

Overview of the models. We use the same colour for parts, which we change in a single approximation step

Notation and preliminaries

• ${\mathbb {N}}:=\left\{ 1,2,3,\ldots \right\}$, ${\mathbb {N}}_{0}:=\{0\}\cup {\mathbb {N}}$, ${\mathbb {R}}_{0}^{+}:=[0,\infty )$ and ${\mathbb {R}}^{+}:=(0,\infty )$.

• The vector ${\mathbf {0}}$ is a row vector of appropriate size with all entries equal to 0. The vector ${\mathbf {e}}$ is a column vector of appropriate size with all entries equal to 1. The vector ${\mathbf {e}}_{i}=(0,\ldots ,0,\underset{{i\text {-th element}}}{\underbrace{1}},0,\ldots ,0)$ is a vector of appropriate size.

• The function $1_{\left\{ expression\right\} }$ is 1 if expression is true and 0 otherwise.

• Empty sums are 0, and empty products are 1.

• We call a “generator” a matrix $M\in {\mathbb {R}}^{K\times K}$ with countable index set K if all its off-diagonal elements are non-negative and all its row sums are zero.

• “Markov process” means a continuous-time homogeneous strong Markov process with discrete state space, which is regular with cadlag paths (right-continuous, left limits everywhere). A Markov process is regular if it is non-explosive (i.e. the sequence of jump times of the process diverges almost surely) and its transition intensity matrix is a generator.

2 SOQN with backordering

2.1 Description of the model

An SOQN (SOQN-BO) consists of a queueing network (“inner network”), a resource pool, and an external queue (Fig. 2). Customers arrive according to a Poisson process with rate $\lambda _{\text {BO}}>0$. Every customer requires exactly one resource from the resource pool for service. If on arrival of a customer a resource is available, the resource enters the inner network to complete the customer’s order. If on arrival no resource is available, the new customer waits in the external queue under the first-come, first-served (FCFS) regime until a resource becomes available (backordering).

Fig. 2

An SOQN with backordering and external queue

When the resource exits the inner network, it returns to the resource pool (referred to as node 0) and waits for the next customer. Whenever the external queue is not empty and a resource item is returned to the resource pool, this item is instantaneously synchronised with the customer at the head of the line. The resources therefore move in a closed network, called resource network. The maximal number of resources in the resource pool is $N$.

The inner network consists of $J\ge 1$ numbered service stations (nodes), denoted by $\overline{J}:=\left\{ 1,\ldots ,J\right\}$. Each station j consists of a single server with infinite waiting room under FCFS regime or processor sharing regime. Customers in the network are indistinguishable. The service times are exponentially distributed random variables with mean 1. If there are $n_{j}>0$ customers present at node j, service at node j is provided with intensity $\nu _{j}(n_{j})>0$. All service and inter-arrival times constitute an independent family of random variables.

Movements of resources in the inner network are governed by a Markovian routing mechanism: After synchronisation with a customer, a resource visits node j with probability $r(0,j)\ge 0$. A resource, when leaving node i, selects with probability $r(i,j)\ge 0$ to visit node j next, and then enters node j immediately. It starts service if it finds the server idle, otherwise it joins the tail of the queue at node j. This resource can also leave the inner network with probability $r(i,0)\ge 0$. It holds $\sum _{j=0}^{J}r(i,j)=1$ with $r(0,0):=0$ for all $i\in \overline{J}_{0}:=\left\{ 0,1,\ldots ,J\right\}$. Given the departure node i, the resource’s routing decision is made independently of the network’s history. We assume that the routing matrix $\mathcal {R}:=\left( r(i,j):i,j\in \overline{J}_{0}\right)$ is irreducible.

To obtain a Markovian process description, we denote by $X_{\text {ex}}(t)$ the number of customers in the external queue at time $t\ge 0$, by $Y_{0}(t)$ the number of resources in the resource pool at time $t\ge 0$ and by $Y_{j}(t)$, $j\in \overline{J}$, the number of resources present at node j in the inner network at time $t\ge 0$, either waiting or in service. We call this $Y_{j}(t)$ queue length at node $j\in \overline{J}_{0}$ at time $t\ge 0$. Then $\mathbf {Y}(t):=\left( Y_{j}(t):j\in \overline{J}_{0}\right)$ is the queue length vector of the resource network at time $t\ge 0$. We define the joint queue length process of the semi-open network with backordering by $Z_{\text {BO}}:=\left( \left( X_{\text {ex}}(t),\mathbf {Y}(t)\right) :t\ge 0\right) .$ Due to the independence and memorylessness assumptions, $Z_{\text {BO}}$ is an irreducible Markov process with state space

$$\begin{aligned} E&:=\left\{\left( 0,n_{0},n_{1},\ldots ,n_{J}\right) :n_{j}\in \left\{ 0,\ldots ,N\right\} \,\forall j\in \overline{J}_{0},\sum _{j\in \overline{J}_{0}}n_{j}=N\right \}\\&\quad \cup \left \{\left( n_{\text {ex}},0,n_{1},\ldots ,n_{J}\right) :n_{\text {ex}}\in {\mathbb {N}},\,n_{j}\in \left\{ 0,\ldots ,N\right\} \,\forall j\in \overline{J},\sum _{j\in \overline{J}}n_{j}=N\right \}. \end{aligned}$$

2.2 Stability

In this section, we analyse the stability of the system and $Z_{\text {BO}}$, which has infinitesimal generator ${\mathbf {Q}}:=\left( q(z;{\tilde{z}}):z,{\tilde{z}}\in E\right)$ with the following transition rates for $\left( n_{\text {ex}},\mathbf {n}\right) ,\ \left( 0,\mathbf {n}\right) \in E$, where $\mathbf {n}:=\left( n_{j}:j\in \overline{J}_{0}\right)$:

$$\begin{aligned}&q\left( \left( n_{\text {ex}},\mathbf {n}\right) ;\left( n_{\text {ex}}+1,\mathbf {n}\right) \right) =\lambda _{\text {BO}}\cdot 1_{\left\{ n_{0}=0\right\} },\ n_{\text {ex}}\ge 0,\\&q\left( \left( 0,\mathbf {n}\right) ;\left( 0,\mathbf {n}-\mathbf {e}_{0}+\mathbf {e}_{i}\right) \right) =\lambda _{\text {BO}}\cdot r(0,i)\cdot 1_{\left\{ n_{0}>0\right\} },\ i\in \overline{J},\\&q\left( \left( n_{\text {ex}},\mathbf {n}\right) ;\left( n_{\text {ex}},\mathbf {n}-\mathbf {e}_{i}+\mathbf {e}_{j}\right) \right) =\nu _{i}(n_{i})\cdot r(i,j)\cdot 1_{\left\{ n_{i}>0\right\} },\ n_{\text {ex}}\ge 0,\ i,j\in \overline{J},\\&q\left( \left( 0,\mathbf {n}\right) ;\left( 0,\mathbf {n}-\mathbf {e}_{i}+\mathbf {e}_{0}\right) \right) =\nu _{i}(n_{i})\cdot r(i,0)\cdot 1_{\left\{ n_{i}>0\right\} },\ i\in \overline{J},\\&q\left( \left( n_{\text {ex}},\mathbf {n}\right) ;\left( n_{\text {ex}}-1,\mathbf {n}-\mathbf {e}_{i}+\mathbf {e}_{j}\right) \right) =\nu _{i}(n_{i})\cdot r(i,0)\cdot r(0,j)\cdot 1_{\left\{ n_{\text {ex}}>0\right\} }\cdot 1_{\left\{ n_{i}>0\right\} },\ \\&\quad n_{\text {ex}}>0,\ i\in \overline{J}. \end{aligned}$$

Furthermore, $q(z;{\tilde{z}})=0$ for any other pair $z\ne {\tilde{z}}$, and

$$\begin{aligned} q\left( z;z\right) =-\sum _{\begin{array}{c} {\tilde{z}}\in E,\\ {\tilde{z}}\ne z \end{array} }q\left( z;{\tilde{z}}\right) \qquad \forall z\in E. \end{aligned}$$

$Z_{\text {BO}}$ is a level-independent quasi-birth-and-death process in the sense of (Latouche and Ramaswami 1999, Def. 1.3.1, p. 12). The level is the length $n_{\text {ex}}$ of the external queue. The phase is the state of the resource network. For level zero, the phase space is

$$\begin{aligned} {\widetilde{E}}_{0}:=\left \{\left( n_{0},n_{1},\ldots ,n_{J}\right) :\ n_{j}\in \left\{ 0,\ldots ,N\right\} \,\forall j\in \overline{J}_{0},\sum _{j\in \overline{J}_{0}}n_{j}=N\right \}. \end{aligned}$$

When level $n_{\text {ex}}> 0$, the resource pool is empty. Hence, for positive levels the phase space is

$$\begin{aligned} \widetilde{E}_{\text {+}}:=\left \{\left( 0,n_{1},\ldots ,n_{J}\right) :\,n_{j}\in \left\{ 0,\ldots ,N\right\} \,\forall j\in \overline{J},\,\sum _{j\in \overline{J}}n_{j}=N\right\}. \end{aligned}$$

Arranging states by level, the infinitesimal generator ${\mathbf {Q}}$ of $Z_{\text {BO}}$ is

$$\begin{aligned} {\mathbf {Q}}=\left( \begin{array}{ccccc} {\mathbf {B}}_{0} &{} {\mathbf {B}}_{1}\\ {\mathbf {B}}_{2} &{} {\mathbf {A}}_{0} &{} {\mathbf {A}}_{1}\\ &{} {\mathbf {A}}_{-1} &{} {\mathbf {A}}_{0} &{} {\mathbf {A}}_{1}\\ &{} &{} &{} \ddots \ddots &{} \ddots \end{array}\right) , \end{aligned}$$

where ${\mathbf {B}}_{0}\in {\mathbb {R}}^{{\widetilde{E}}_{0}\times {\widetilde{E}}_{0}}$, ${\mathbf {B}}_{1}\in {\mathbb {R}}^{{\widetilde{E}}_{0}\times \widetilde{E}_{\text {+}}}$, ${\mathbf {B}}_{2}\in {\mathbb {R}}^{\widetilde{E}_{\text {+}}\times {\widetilde{E}}_{0}}$ and ${\mathbf {A}}_{-1}$, ${\mathbf {A}}_{0}$, ${\mathbf {A}}_{1}\in {\mathbb {R}}^{\widetilde{E}_{\text {+}}\times \widetilde{E}_{\text {+}}}$ are matrices.

${\mathbf {A}}_{1}$ is a non-negative matrix with the following positive elements

$$\begin{aligned} a_{1}\left( \left( 0,n_{1},\ldots ,n_{J}\right) ;\left( 0,n_{1},\ldots ,n_{J}\right) \right) =\lambda _{\text {BO}}. \end{aligned}$$

${\mathbf {A}}_{-1}$ is a non-negative matrix with at most the following non-negative elements

$$\begin{aligned}&a_{-1}\left( \left( 0,n_{1},\ldots ,n_{J}\right) ;\left( 0,n_{1},\ldots ,n_{i}-1,\ldots ,n_{j}+1,\ldots ,n_{J}\right) \right) \nonumber \\&\quad =\nu _{i}(n_{i})\cdot r(i,0)\cdot r(0,j)\cdot 1_{\left\{ n_{i}>0\right\} },\quad i,j\in \overline{J}. \end{aligned}$$

${\mathbf {A}}_{0}$ has non-negative off-diagonal elements and strictly negative diagonals. The off-diagonal elements are

$$\begin{aligned}&a_{0}\left( \left( 0,n_{1},\ldots ,n_{J}\right) ;\left( 0,n_{1},\ldots ,n_{i}-1,\ldots ,n_{j}+1,\ldots ,n_{J}\right) \right) \\&\quad =\nu _{i}(n_{i})\cdot r(i,j)\cdot 1_{\left\{ n_{i}>0\right\} },\quad i,j\in \overline{J}. \end{aligned}$$

Let $\varvec{\pi }_{\text {BO}}:=\left( \pi _{\text {BO}}\left( n_{\text {ex}},\mathbf {n}\right) :\left( n_{\text {ex}},\mathbf {n}\right) \in E\right)$ be the steady-state distribution of the Markov process $Z_{\text {BO}}$. The global balance equations $\varvec{\pi }_{\text {BO}}\cdot \mathbf {Q=0}$ are:

For $n_{\text {ex}}=0$

$$\begin{aligned}&\pi _{\text {BO}}\left( 0,\mathbf {n}\right) \\&\qquad \cdot \Big (\lambda _{\text {BO}}+\sum _{i\in \overline{J}}\sum _{j\in \overline{J}\setminus \{i\}}\nu _{i}(n_{i})\cdot r(i,j)\cdot 1_{\left\{ n_{i}>0\right\} }+\sum _{i\in \overline{J}}\nu _{i}(n_{i})\cdot r(i,0)\cdot 1_{\left\{ n_{i}>0\right\} }\Big )\\&\quad =\sum _{i\in \overline{J}}\pi _{\text {BO}}\left( 0,\mathbf {n}+\mathbf {e}_{0}-\mathbf {e}_{i}\right) \cdot \lambda _{\text {BO}}\cdot r(0,i)\cdot 1_{\left\{ n_{i}>0\right\} }\\&\qquad +\sum _{i\in \overline{J}}\sum _{j\in \overline{J}\setminus \{i\}}\pi _{\text {BO}}\left( 0,\mathbf {n}+\mathbf {e}_{i}-\mathbf {e}_{j}\right) \cdot \nu _{i}(n_{i}+1)\cdot r(i,j)\cdot 1_{\left\{ n_{j}>0\right\} }\\&\qquad +\sum _{i\in \overline{J}}\pi _{\text {BO}}\left( 0,\mathbf {n}+\mathbf {e}_{i}-\mathbf {e}_{0}\right) \cdot \nu _{i}(n_{i}+1)\cdot r(i,0)\cdot 1_{\left\{ n_{0}>0\right\} }\\&\qquad +\sum _{i\in \overline{J}}\sum _{j\in \overline{J}\setminus \{i\}}\ \ \pi _{\text {BO}}\left( 1,\mathbf {n}+\mathbf {e}_{i}-\mathbf {e}_{j}\right) \cdot \nu _{i}(n_{i}+1)\cdot r(i,0)\cdot r(0,j)\cdot 1_{\left\{ n_{j}>0\right\} }\cdot 1_{\left\{ n_{0}=0\right\} }\\&\qquad +\sum _{i\in \overline{J}}\pi _{\text {BO}}\left( 1,\mathbf {n}\right) \cdot \nu _{i}(n_{i})\cdot r(i,0)\cdot r(0,i)\cdot 1_{\left\{ n_{i}>0\right\} }\cdot 1_{\left\{ n_{0}=0\right\} }. \end{aligned}$$

For $n_{\text {ex}}>0$, which implies $n_{0}=0,$

$$\begin{aligned}&\pi _{\text {BO}}\left( n_{\text {ex}},\mathbf {n}\right) \\&\qquad \cdot \Big (\lambda _{\text {BO}}+\sum _{i\in \overline{J}}\sum _{j\in \overline{J}\setminus \{i\}}\nu _{i}(n_{i})\cdot r(i,j)\cdot 1_{\left\{ n_{i}>0\right\} }+\sum _{i\in \overline{J}}\nu _{i}(n_{i})\cdot r(i,0)\cdot 1_{\left\{ n_{i}>0\right\} }\Big )\\&\quad =\pi _{\text {BO}}\left( n_{\text {ex}}-1,\mathbf {n}\right) \cdot \lambda _{\text {BO}}\\&\qquad +\sum _{j\in \overline{J}}\sum _{i\in \overline{J}\setminus \{i\}}\pi _{\text {BO}}\left( n_{\text {ex}},\mathbf {n}+\mathbf {e}_{i}-\mathbf {e}_{j}\right) \cdot \nu _{i}(n_{i}+1)\cdot r(i,j)\cdot 1_{\left\{ n_{j}>0\right\} }\\&\qquad +\sum _{j\in \overline{J}}\sum _{i\in \overline{J}\setminus \{i\}}\pi _{\text {BO}}\left( n_{\text {ex}}+1,\mathbf {n}+\mathbf {e}_{i}-\mathbf {e}_{j}\right) \cdot \nu _{i}(n_{i}+1)\cdot r(i,0)\cdot r(0,j)\cdot 1_{\left\{ n_{j}>0\right\} }\\&\qquad +\sum _{i\in \overline{J}}\pi _{\text {BO}}\left( n_{\text {ex}}+1,\mathbf {n}\right) \cdot \nu _{i}(n_{i})\cdot r(i,0)\cdot r(0,i)\cdot 1_{\left\{ n_{i}>0\right\} }. \end{aligned}$$

No closed-form expression is known for $\varvec{\pi }_{\text {BO}}$ in case of $J>1$. Latouche and Ramaswami developed a logarithmic reduction algorithm for level-independent quasi-birth-and-death processes to compute the steady-state distribution (Latouche and Ramaswami 1993, Latouche and Ramaswami 1999, Theorem 6.4.1 and Lemma 6.4.3, p. 142ff.). For $J=1$, we calculate a closed-form expression for the steady-state distribution, see Sect. 2.4.

To determine the stability condition of the system, we define traffic equations:

$$\begin{aligned} \eta _{j}=\sum _{i\in \overline{J}_{0}}\eta _{i}\cdot r(i,j),\qquad j\in \overline{J}_{0}. \end{aligned}$$

(1)

Denote by $\lambda _{\text {BO,max}}$ the throughput of node 0 in the closed network depicted in Fig. 3, which is obtained from the original SOQN when infinitely many customers reside in the external queue. In Lavenberg (1978) this network is called saturated—we will call this network stability network (“stb”). In this network, resources which enter node 0 spend zero time there and jump to the next node according to the branching vector $(r(0,j):j\in \overline{J})$. In Lavenberg (1978), it is proved that an SOQN (with backordering) is stable if $\lambda _{\text {BO}}<\lambda _{\text {BO,max}}$, and that $\lambda _{\text {BO}}>$ $\lambda _{\text {BO,max}}$ implies instability. We use matrix geometrical methods to show that $\lambda _{\text {BO}}<\lambda _{\text {BO,max}}$ is sufficient and necessary for stability. To simplify notation, we define

$$\begin{aligned} C^{\text {stb}}(\overline{J},N):=\sum _{\sum _{j\in \overline{J}}n_{j}=N}\prod _{j=1}^{J}\left( \prod _{i=1}^{n_{j}}\frac{\eta _{j}}{\nu _{j}(i)}\right) . \end{aligned}$$

Proposition 1

The system is stable if and only if $\lambda _{\text {BO}}<\lambda _{\text {BO,max}}$ with

$$\begin{aligned} \lambda _{\text {BO,max}}=\eta _{0}\cdot \frac{C^{\text {stb}}(\overline{J},N-1)}{C^{\text {stb}}(\overline{J},N)}. \end{aligned}$$

(2)

Proof

We apply matrix-geometric methods using (Latouche 2011, Theorem 1): Given the irreducible inter-level generator matrix ${\mathbf {A}}:={\mathbf {A}}_{-1}+{\mathbf {A}}_{0}+{\mathbf {A}}_{1}$ of $Z_{\text {BO}}$ and the stochastic solution $\mathbf {\varvec{\alpha }}:=(\alpha ({\widetilde{\mathbf {n}}}):{\widetilde{\mathbf {n}}}\in \widetilde{E}_{\text {+}})$ of $\mathbf {\varvec{\alpha }}\cdot {\mathbf {A}}={\mathbf {0}}$, the process $Z_{\text {BO}}$ is stable if and only if

$$\begin{aligned} \mathbf {\varvec{\alpha }}\cdot {\mathbf {A}}_{1}\cdot {\mathbf {e}}<\mathbf {\varvec{\alpha }}\cdot {\mathbf {A}}_{-1}\cdot {\mathbf {e}}. \end{aligned}$$

(3)

Because ${\mathbf {A}}_{1}$ is a diagonal matrix with $\lambda _{\text {BO}}$ on its diagonal and $\mathbf {\varvec{\alpha }}$ is a stochastic vector, the left-hand side of (3) is $\lambda _{\text {BO}}$. We define

$$\begin{aligned} \lambda _{\text {BO,max}}:=\mathbf {\varvec{\alpha }}\cdot {\mathbf {A}}_{-1}\cdot {\mathbf {e}}. \end{aligned}$$

To determine $\mathbf {\varvec{\alpha }}$, we note that the non-negative non-diagonal elements of the generator ${\mathbf {A}}$ are for $i\ne j$ of the form

$$\begin{aligned}&a\left( \left( 0,n_{1},\ldots ,n_{J}\right) ;\left( 0,n_{1},\ldots ,n_{i}-1,\ldots ,n_{j}+1,\ldots ,n_{J}\right) \right) \\&\quad =\nu _{i}(n_{i})\cdot \Big (r(i,j)+r(i,0)\cdot r(0,j)\Big )\cdot 1_{\left\{ n_{i}>0\right\} }. \end{aligned}$$

The diagonal elements are chosen to set row sums to zero.

We now solve for all ${\widetilde{\mathbf {n}}}:=\left( 0,n_{1},\ldots ,n_{J}\right) \in \widetilde{E}_{\text {+}}$

$$\begin{aligned}&\alpha \left( {\widetilde{\mathbf {n}}}\right) \cdot \sum _{i\in \overline{J}}\nu _{i}(n_{i})\cdot \sum _{j\in \overline{J}}\Big (r(i,j)+r(i,0)\cdot r(0,j)\Big )\cdot 1_{\left\{ n_{i}>0\right\} } \\&\quad =\sum _{j\in \overline{J}}\sum _{i\in \overline{J}}\alpha \left( {\widetilde{\mathbf {n}}}+\mathbf {e}_{i}-\mathbf {e}_{j}\right) \cdot \nu _{i}(n_{i}+1)\cdot \Big (r(i,j)+r(i,0)\cdot r(0,j)\Big )\cdot 1_{\left\{ n_{j}>0\right\} }.\nonumber \end{aligned}$$

(4)

Fig. 3

Stability network – A closed network described by the inter-level generator matrix ${\mathbf {A}}$ from the proof of Proposition 1. Interpretation I: a generalised Gordon–Newell network with zero service time at node 0. Interpretation II: a classical Gordon–Newell network obtained after rerouting at node 0

Equation (4) is the global balance equation of a generalised Gordon–Newell network with node set $\overline{J}_{0}$, $N$ customers and zero service time at node 0.

Equation (4) has another interpretation, which allows us to use standard algorithms for Gordon–Newell networks for performance analysis. Note that the status of node 0 is invariant ($=0$), therefore we define $\alpha '(n_{1},\ldots ,n_{J}) :=\alpha \left( 0,n_{1},\ldots ,n_{J}\right)$ and the routing matrix $\mathcal {R}':=\left( r'(i,j):i,j\in \overline{J}\right)$ with

$$\begin{aligned} r'(i,j):=r(i,j)+r(i,0)\cdot r(0,j). \end{aligned}$$

(5)

$\mathcal {R}'$ is obtained from $\mathcal {R}$ when every resource directed to 0 skips this node and jumps to the next node according to $(r(0,j):j\in \overline{J})$. So Eq. (4) can be written as

$$\begin{aligned}&\alpha '\left( n_{1},\ldots ,n_{J}\right) \cdot \sum _{i\in \overline{J}}\nu _{i}(n_{i})\sum _{j\in \overline{J}}r'(i,j)\cdot 1_{\left\{ n_{i}>0\right\} }\nonumber \\&\quad =\sum _{j\in \overline{J}}\sum _{i\in \overline{J}}\alpha '\left( (n_{1},\ldots ,n_{J})+\mathbf {e}_{i}-\mathbf {e}_{j}\right) \cdot \nu _{i}(n_{i}+1)\cdot r'(i,j)\cdot 1_{\left\{ n_{j}>0\right\} }. \end{aligned}$$

(6)

Equation (6) is the global balance equation of a Gordon–Newell network with nodes $\overline{J}=\left\{ 1,2,\ldots ,J\right\}$, $N$ customers, routing matrix $\mathcal {R}'$ and steady-state distribution

$$\begin{aligned} \alpha '(n_{1},\ldots ,n_{J})=\left[ C_{\text {BO}}^{\text {stb}'}(\overline{J},N)\right] ^{-1}\prod _{j=1}^{J}\left( \prod _{i=1}^{n_{j}}\frac{\eta '_{j}}{\nu _{j}(i)}\right) \end{aligned}$$

(7)

where ${\varvec{\eta }}':=\left( \eta '_{j}:j\in \overline{J}\right)$ is a solution of the traffic equation ${\varvec{\eta }}'\cdot \mathcal {R}'={\varvec{\eta }}'$ and

$$\begin{aligned} C_{\text {BO}}^{\text {stb}'}(\overline{J},N):=\sum _{\sum _{j\in \overline{J}}n_{j}=N}\prod _{j=1}^{J}\left( \prod _{i=1}^{n_{j}}\frac{\eta '_{j}}{\nu _{j}(i)}\right) \end{aligned}$$

is the normalisation constant. Because of the special structure (5) of $\mathcal {R}'$, $\eta '_{j}:=\eta _{j}$ for all $j\in J$ is a solution of ${\varvec{\eta }}'\cdot \mathcal {R}'={\varvec{\eta }}'$, see (Krenzler et al. 2016, Proposition 2.1). Consequently, $C_{\text {BO}}^{\text {stb}'}(\overline{J},N)=C^{\text {stb}}(\overline{J},N)$, and we can switch between both interpretations without recalculating ${\varvec{\eta }}'$ and $C_{\text {BO}}^{\text {stb}'}(\overline{J},N)$, and obtain, e.g.

$$\begin{aligned} \alpha '(n_{1},\ldots ,n_{J})=\left[ C^{\text {stb}}(\overline{J},N)\right] ^{-1}\prod _{j=1}^{J}\left( \prod _{i=1}^{n_{j}}\frac{\eta {}_{j}}{\nu _{j}(i)}\right) . \end{aligned}$$

We now calculate $\lambda _{\text {BO,max}}$ explicitly.

$$\begin{aligned}&\lambda _{\text {BO,max}}=\alpha \cdot {\mathbf {A}}_{-1}\cdot {\mathbf {e}} =\sum _{(0,m_{1},\ldots ,m_{J})\in \widetilde{E}_{\text {+}}}\big (\alpha \cdot {\mathbf {A}}_{-1}\big )\left( 0,m_{1},\ldots ,m_{J}\right) \\&\quad =\sum _{(0,m_{1},\ldots ,m_{J})\in \widetilde{E}_{\text {+}}}\Bigg [\sum _{\left( 0,n_{1},\ldots ,n_{J}\right) \in E_{1}}\alpha \left( 0,n_{1},\ldots ,n_{J}\right) \\&\qquad \cdot a_{-1}\left( \left( 0,n_{1},\ldots ,n_{J}\right) ;\left( 0,m_{1},\ldots ,m_{J}\right) \right) \Bigg ]\\&\quad =\sum _{\left( 0,n_{1},\ldots ,n_{J}\right) \in \widetilde{E}_{\text {+}}}\Bigg [\alpha \left( 0,n_{1},\ldots ,n_{J}\right) \\&\qquad \cdot \sum _{(0,m_{1},\ldots ,m_{J})\in \widetilde{E}_{\text {+}}} \quad a_{-1}\left( \left( 0,n_{1},\ldots ,n_{J}\right) ;\left( 0,m_{1},\ldots ,m_{J}\right) \right) \Bigg ]\\&\quad =\sum _{\left( 0,n_{1},\ldots ,n_{J}\right) \in \widetilde{E}_{\text {+}}}\alpha \left( 0,n_{1},\ldots ,n_{J}\right) \cdot \Bigg [\sum _{i=1}^{J}\sum _{j=1}^{J}\nu _{i}(n_{i})\cdot r(i,0)\cdot r(0,j)\cdot 1_{\left\{ n_{i}>0\right\} }\Bigg ]\\&\quad =\sum _{\left( 0,n_{1},\ldots ,n_{J}\right) \in \widetilde{E}_{\text {+}}}\underbrace{\alpha \left( 0,n_{1},\ldots ,n_{J}\right) }_{=\alpha '(n_{1},\ldots ,n_{J})}\cdot \Bigg [\sum _{i=1}^{J}\nu _{i}(n_{i})\cdot r(i,0)\cdot 1_{\left\{ n_{i}>0\right\} }\cdot \underbrace{\sum _{j=1}^{J}r(0,j)}_{=1}\Bigg ]\\&\quad =\sum _{\left( n_{1},\ldots ,n_{J}\right) \in \widetilde{E}_{\text {+}}}\alpha '(n_{1},\ldots ,n_{J})\cdot \Bigg [\sum _{i=1}^{J}\nu _{i}(n_{i})\cdot 1_{\left\{ n_{i}>0\right\} }\cdot r(i,0)\Bigg ]\\&\quad =\sum _{i=1}^{J}\Bigg [\underbrace{\sum _{n_{i}=0}^{N}\sum _{{\mathop {\sum _{j\in \overline{J}\backslash \left\{ i\right\} }n_{j}=N-n_{i}}\limits ^{n_{j}\in \left\{ 0,\ldots ,N\right\} ,\ j\in \overline{J}\backslash \left\{ i\right\} }}}\alpha '(n_{1},\ldots ,n_{J})\cdot \nu _{i}(n_{i})\cdot 1_{\left\{ n_{i}>0\right\} }}_{(*)}\Bigg ]\cdot r(i,0). \end{aligned}$$

The expression ($*$) is the throughput through node i in the Gordon–Newell network with routing matrix $\mathcal {R}'$:

$$\begin{aligned}&TH_{i}^{\text {stb}}(N):=\sum _{n_{i}=0}^{N}\sum _{{\mathop {\sum _{j\in \overline{J}\backslash \left\{ i\right\} }n_{j}=N-n_{i}}\limits ^{n_{j}\in \left\{ 0,\ldots ,N\right\} ,\ j\in \overline{J}\backslash \left\{ i\right\} }}}\alpha '(n_{1},\ldots ,n_{J})\cdot \nu _{i}(n_{i})\cdot 1_{\left\{ n_{i}>0\right\} },\qquad i\in \overline{J}. \end{aligned}$$

This yields

$$\lambda _{\text {BO,max}}=\sum _{i=1}^{J}TH_{i}^{\text {stb}}(N)\cdot r(i,0).$$

(8)

According to (Bolch et al. 1998, p. 374, (8.14)) it holds $TH_{i}^{\text {stb}}(N)=\eta _{i}\cdot \frac{C^{\text {stb}}(\overline{J},N-1)}{C^{\text {stb}}(\overline{J},N)}.$ Therefore,

$$\begin{aligned} \lambda _{\text {BO,max}}&=\sum _{i=1}^{J}\eta _{i}\cdot \frac{C^{\text {stb}}(\overline{J},N-1)}{C^{\text {stb}}(\overline{J},N)}\cdot r(i,0)=\frac{C^{\text {stb}}(\overline{J},N-1)}{C^{\text {stb}}(\overline{J},N)}\underbrace{\sum _{i=1}^{J}\eta _{i}\cdot r(i,0)}_{=\eta _{0}}.~~~ \end{aligned}$$

$\square$

Remark 1

The right-hand side of Eq. (8) is the throughput of node 0 in the stability network in Fig. 3. Formally, this is

$$\begin{aligned} \lambda _{\text {BO,max}}=TH_{0}^{\text {stb}}(N)\text { with }TH_{0}^{\text {stb}}(N):=\sum _{i=1}^{J}TH_{i}^{\text {stb}}(N)\cdot r(i,0). \end{aligned}$$

(9)

The advantage of representation (9) is that it uses throughputs $TH_{i}^{\text {stb}}(N)$, $i\in \overline{J}$, of a classical Gordon–Newell network with routing matrix $\mathcal {R}'$. We can calculate these throughputs efficiently with standard methods, e.g. mean value analysis (MVA). Using Eq. (2), another representation of $TH_{0}^{\text {stb}}(N)$ is

$$\begin{aligned} TH_{0}^{\text {stb}}(N)=\eta _{0}\cdot \frac{C^{\text {stb}}(\overline{J},N-1)}{C^{\text {stb}}(\overline{J},N)}. \end{aligned}$$

(10)

To calculate efficiently the constants on the right-hand side of Eq. (10), we can use, for example, the convolution algorithm. Both algorithms are illustrated in (Bolch et al. 1998, p. 371ff., Section 8.1 and p. 384ff., Section 8.2).

2.3 Throughputs and idle times

We consider an SOQN-BO in steady state. Let ${\varvec{\eta }}:=(\eta _{j}:j\in \overline{J}_{0})$ be a solution of Eq. (1). ${\varvec{\eta }}$ is unique up to a constant, which implies that the following formula (11) does not depend on that constant. Equation (11) occurs as an approximation in a slightly different setting in Dallery (1990) as (26).

Proposition 2

The local throughput at the nodes $j\in \overline{J}_{0}$ is

$$\begin{aligned} TH_{\text {BO},j}=\lambda _{\text {BO}}\cdot \frac{\eta _{j}}{\eta _{0}}. \end{aligned}$$

(12)

Proof

We define for $j\in \overline{J}_{0}$ in steady state:

• the mean number of departures from j per time unit is $D_{j}$, and

• the mean number of arrivals at j per time unit is $V_{j}$.

From the steady state assumption follows $TH_{\text {BO},j}=V_{j}=D_{j}$. For any $j\in \overline{J}_{0}$ it holds $V_{j}=\sum _{i\in \overline{J}_{0}}D_{i}\cdot r(i,j)$. Therefore, the vector $V_{}:=(V_{j}:j\in \overline{J}_{0})$ fulfils the set of equations $V_{j}=\sum _{i\in \overline{J}_{0}}V_{i}\cdot r(i,j)$, $j\in \overline{J}_{0}$, which is $V_{}=V_{}\cdot \mathcal {R}$. This implies that $V_{j}=\eta _{j}\cdot K$ for some constant $K>0$. Because of $\lambda _{\text {BO}}=V_{0}=\eta _{0}\cdot K$ we have $K=\frac{\lambda _{\text {BO}}}{\eta _{0}}$, and therefore, $V_{j}=\lambda _{\text {BO}}\cdot \frac{\eta _{j}}{\eta _{0}},\ j\in \overline{J}_{0}$. $\square$

Corollary 1

Let $Y_{\text {BO}}:=(Y_{\text {BO},j}:j\in \overline{J})$ denote a random vector which is distributed according to the stationary queue length at the nodes in $\overline{J}$ of the SOQN-BO. If the service rate at node j is independent of the queue length, i.e. $\nu _{j}(\cdot )=\nu _{j},j\in \overline{J}$, then the probability that node j is idling is

$$\begin{aligned} P(Y_{\text {BO},j}=0)=1-\lambda _{\text {BO}}\cdot \frac{\eta _{j}}{\eta _{0}}\cdot \nu _{j}^{-1}. \end{aligned}$$

This is also the proportion of time that node j is idling.

Proof

We define for $j\in \overline{J}$ in steady state:

• the mean number of customers in service is $B_{j}$,

• the mean service time is $S_{j}$, and

• the arrival intensity is $\lambda _{j}$.

According to Little’s formula, $B_{j}=\lambda _{j}\cdot S_{j}$ for every node j. In steady state, the arrival rate $\lambda _{j}$ at node j equals its throughput $TH_{\text {BO},j}$. Hence, from Proposition 2, we have $\lambda _{j}=TH_{\text {BO},j}=\lambda _{\text {BO}}\cdot \left( {\eta _{j}}/{\eta _{0}}\right)$. For node j, with constant service rate $\nu _{j}$, the mean service time $S_{j}$ is $\nu _{j}^{-1}$. Inserting these $\lambda _{j}$ and $S_{j}$ into Little’s formula yields $B_{j}=\lambda _{\text {BO}}\cdot \left( {\eta _{j}}/{\eta _{0}}\right) \cdot \nu _{j}^{-1}$ for the mean number of customers in service. Consequently, the probability that node j is idling is

$$\begin{aligned} P(Y_{\text {BO},j}=0)&= 1-P(Y_{\text {BO},j}>0)=1-E\left[ 1_{\left\{ Y_{\text {BO},j}>0\right\} }\right] =1-B_{j}\\&=1-\lambda _{\text {BO}}\cdot \frac{\eta _{j}}{\eta _{0}}\cdot \nu _{j}^{-1}. \end{aligned}$$

$\square$

2.4 Special case: $J=1$

We consider an SOQN-BO with inner network consisting of one node only.

Theorem 1

For $J=1$, the joint queue length process $Z_{\text {BO}}$ is stable if and only if $\lambda _{\text {BO}}<\nu _{1}(N)$. For stable $Z_{\text {BO}}$ the stationary distribution is $\varvec{\pi }_{\text {BO}}:=\left( \pi _{\text {BO}}\left( n_{\text {ex}},n_{0},n_{1}\right) :(n_{\text {ex}},n_{0},n_{1})\in E\right)$ with

$$\begin{aligned} \pi _{\text {BO}}(n_{\text {ex}},n_{0},n_{1})&=\pi _{\text {BO}}(n_{\text {ex}},N-n_{1},n_{1})\nonumber \\&=\left[ C_{\text {BO}}(\left\{ 1\right\} ,N)\right] ^{-1}\cdot \left( \frac{\lambda _{\text {BO}}}{\nu _{1}(N)}\right) ^{n_{\text {ex}}} \cdot \prod _{m=1}^{n_{1}}\frac{\lambda _{\text {BO}}}{\nu _{1}(m)} \end{aligned}$$

(12)

and normalisation constant

$$\begin{aligned} C_{\text {BO}}(\left\{ 1\right\} ,N){:}{=}\sum _{m=0}^{N-1}\prod _{\ell =1}^{m}\frac{\lambda _{\text {BO}}}{\nu _{1}(\ell )}+ \left( \frac{1}{1-\frac{\lambda _{\text {BO}}}{\nu _{1}(N)}}\right) \cdot \prod _{\ell =1}^{N}\frac{\lambda _{\text {BO}}}{\nu _{1}(\ell )}. \end{aligned}$$

(13)

Proof

Let $(\widehat{X}_{\text {ex}},\widehat{Y}_{0},\widehat{Y}_{1})$ denote a random vector with distribution $\varvec{\pi }_{\text {BO}}$. In (Avi-Itzhak and Heyman 1973, equations (20) and (21)) the authors calculated $P(\widehat{X}_{\text {ex}}+ \widehat{Y}_{1}=m)$. From this we obtain the probabilities for all queues, because for $m\le N$ it holds $\widehat{X}_{\text {ex}}+ \widehat{Y}_{1}=m\Leftrightarrow [\widehat{X}_{\text {ex}}=0\wedge \widehat{Y}_{1}=m\wedge \widehat{Y}_{0}=N-m]$ and for $m> N$ it holds $\widehat{X}_{\text {ex}}+\widehat{Y}_{1}=m \Leftrightarrow [\widehat{X}_{\text {ex}}=m-N\wedge \widehat{Y}_{1}=N \wedge \widehat{Y}_{0}=0]$. $\square$

Proposition 3

Let $(\widehat{X}_{\text {ex}},\widehat{Y}_{0},\widehat{Y}_{1})$ denote a random vector which is distributed according to the steady-state distribution of $Z_{\text {BO}}$ in the SOQN-BO with $J=1$.

(i) For the marginal distributions it holds:

$$\begin{aligned}&P(\widehat{X}_{\text {ex}}=0)= \left[ C_{\text {BO}}(\left\{ 1\right\} ,N)\right] ^{-1}\cdot \sum _{n_{1}=0}^{N}\prod _{m=1}^{n_{1}}\frac{\lambda _{\text {BO}}}{\nu _{1}(m)}, \end{aligned}$$

(14)

$$\begin{aligned}&P(\widehat{X}_{\text {ex}}=n_{\text {ex}})=\left[ C_{\text {BO}}(\left\{ 1\right\} ,N)\right] ^{-1} \left( \frac{\lambda _{\text {BO}}}{\nu _{1}(N)}\right) ^{n_{\text {ex}}} \prod _{m=1}^{N}\frac{\lambda _{\text {BO}}}{\nu _{1}(m)},\ n_{\text {ex}}>0, \end{aligned}$$

(15)

$$\begin{aligned}&P(\widehat{Y}_{0}=N-n_{1},\widehat{Y}_{1}=n_{1})=\left[ C_{\text {BO}}(\left\{ 1\right\} ,N)\right] ^{-1}\prod _{m=1}^{n_{1}}\frac{\lambda _{\text {BO}}}{\nu _{1}(m)}, \ 0\le n_{1}<N, \end{aligned}$$

(16)

$$\begin{aligned}&P(\widehat{Y}_{0}=0,\widehat{Y}_{1}=N)=\left[ C_{\text {BO}}(\left\{ 1\right\} ,N)\right] ^{-1} \frac{1}{1-\frac{\lambda _{\text {BO}}}{\nu _{1}(N)}} \prod _{m=1}^{N}\frac{\lambda _{\text {BO}}}{\nu _{1}(m)}. \end{aligned}$$

(17)

(ii) The average external queue length is

$$\begin{aligned} {\widehat{L_{\text {ex}}}}&=P(\widehat{Y}_{0}=0,\widehat{Y}_{1}=N)\cdot \frac{\lambda _{\text {BO}}}{\nu _{1}(N)-\lambda _{\text {BO}}} \end{aligned}$$

(18)

and the average waiting time of customers in the external queue is

$$\begin{aligned} {\widehat{W_{\text {ex}}}}=\frac{{\widehat{L_{\text {ex}}}}}{\lambda _{\text {BO}}}=P(\widehat{Y}_{0}=0,\widehat{Y}_{1}=N)\cdot \frac{1}{\nu _{1}(N)-\lambda _{\text {BO}}}. \end{aligned}$$

(19)

Proof

(i) All probabilities can be expressed in terms of Eq. (12), resulting in routine calculations.

(ii) The stationary average external queue length is obtained directly from the marginal distribution (15) in (i). Equation (19) follows by Little’s law, see e.g. Little and Graves (2008). $\square$

An important question about SOQN is whether more resources yield more throughput. Proposition 4 shows that even for $J=1$ this is not the case.

Proposition 4

Let $J=1$. Then the following are equivalent:

• (i) $\lambda _{\text {BO,max}}=TH_{0}^{\text {stb}}(\cdot )$ is non-decreasing on ${\mathbb {N}}$.

• (ii) $\nu _{1}$ is non-decreasing on ${\mathbb {N}}$.

Proof

Because of Eq. (2) and (9), (i) is equivalent to

$$\begin{aligned} \forall \;N\in {\mathbb {N}}:\;C^{\text {stb}}(\left\{ 1\right\} ,N-1)\cdot C^{\text {stb}}(\left\{ 1\right\} ,N+1)\le C^{\text {stb}}(\left\{ 1\right\} ,N)^{2}. \end{aligned}$$

(20)

We note

$$\begin{aligned} ~~~C^{\text {stb}}(\left\{ 1\right\} ,N+1)&=C^{\text {stb}}(\left\{ 1\right\} ,N)\cdot \frac{\eta _{1}}{\nu _{1}(N+1)},\\ C^{\text {stb}}(\left\{ 1\right\} ,N-1)&=C^{\text {stb}}(\left\{ 1\right\} ,N)\cdot \frac{\nu _{1}(N)}{\eta _{1}},\end{aligned}$$

implying

$$~~~(20) \Leftrightarrow \forall \;N\in {\mathbb {N}}:\;\frac{\nu _{1}(N)}{\nu _{1}(N+1)}\cdot C^{\text {stb}}(\left\{ 1\right\} ,N)^{2}\le C^{\text {stb}}(\left\{ 1\right\} ,N)^{2}$$

and by $C^{\text {stb}}(\left\{ 1\right\} ,N)^{2}>0$ this is equivalent to $\forall \;N\in {\mathbb {N}}:\;\nu _{1}(N)\le \nu _{1}(N+1),$ which is (ii). $\square$

Van der Wal (1989) proved that in general it holds (ii) $\Rightarrow$ (i).

Proposition 5

Let $J\ge 1$ and all service rates be non-decreasing in the number of customers, i.e. $\nu _{j}(n+1)\ge \nu _{j}(n)$, $n\in \{0,\ldots ,N-1\}$, for all $j\in \overline{J}$. Then $\lambda _{\text {BO,max}}=TH_{0}^{\text {stb}}(\cdot )$ is non-decreasing on ${\mathbb {N}}$.

3 Lost customer approximation of SOQN-BO

The steady-state distribution of an SOQN-BO is unknown for $J>1$. For a modified system, we get closed product-form results for $J\ge 1$.

3.1 SOQN with lost customers

We consider a modification of the SOQN-BO model from Sect. 2.1: Newly arriving customers are rejected and lost for the system if the resource pool is empty (Fig. 4). This property is termed, e.g. “lost customers”, “lost sales”, “lost arrivals” or “loss systems”. We denote this system SOQN-LC. Customers arrive in a Poisson stream with rate $\lambda _{\text {LC}}>0$. Because of loss of customers, the effective arrival rate $\lambda _{\text {eff}}(\lambda _{\text {LC}})$ is smaller.

Such SOQN-LC can be investigated using Gordon–Newell network theory for the resource network (after suitable modification), see (Chen and Yao 2001, p. 21).

Fig. 4

Transition from SOQN-BO (left) to SOQN-LC (right)

To obtain a Markovian process description, we denote by $Y_{0}(t)$ the number of resources in the resource pool at time $t\ge 0$ and by $Y_{j}(t)$, $j\in \overline{J}$, the number of resources at node j in the inner network at time $t\ge 0$, either waiting or in service (queue length at node $j\in \overline{J}$). Then $\mathbf {Y}(t):=\left( Y_{j}(t):j\in \overline{J}_{0}\right)$ is the local queue length vector of the resource network at time $t\ge 0$. We define the joint queue length process of this semi-open network by

$$\begin{aligned} Z_{\text {LC}}:=\left( \mathbf {Y}(t):t\ge 0\right) . \end{aligned}$$

Due to the usual independence and memorylessness assumptions (see Sect. 2.1), $Z_{\text {LC}}$ is an irreducible Markov process on state space

$$\begin{aligned} E_{LC}&:=\left \{\left( n_{0},n_{1},\ldots ,n_{J}\right) :\,n_{j}\in \left\{ 0,\ldots ,N\right\} \,\forall j\in \overline{J}_{0},\,\sum _{j\in \overline{J}_{0}}n_{j}=N\right \}. \end{aligned}$$

The stationary distribution $\varvec{\pi }_{\text {LC}}:=\left( \pi _{\text {LC}}\left( \mathbf {n}\right) :\mathbf {n}\in E_{LC}\right)$ of $Z_{\text {LC}}$ in product form is available (Chen and Yao 2001, p. 22, Theorem 2.5): For $\mathbf {n}:=\left( n_{j}:j\in \overline{J}_{0}\right) \in E_{LC}$

$$\begin{aligned} \pi _{\text {LC}}\left( \mathbf {n}\right) =\left[ C_{\text {LC}}(\overline{J}_{0},N)\right] ^{-1}\cdot \left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}\cdot \prod _{j=1}^{J}\left( \prod _{i=1}^{n_{j}} \frac{\eta _{j}}{\nu _{j}(i)}\right) \end{aligned}$$

(21)

with normalisation constant

$$\begin{aligned} C_{\text {LC}}(\overline{J}_{0},N)&{:}{=}\sum _{\sum _{j\in \overline{J}_{0}}n_{j}=N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}\cdot \prod _{j=1}^{J}\left( \prod _{i=1}^{n_{j}}\frac{\eta _{j}}{\nu _{j}(i)}\right) \nonumber \\&=\sum _{n_{0}=0}^{N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}\sum _{\sum _{j\in \overline{J}}n_{j}=N-n_{0}}\prod _{j=1}^{J}\left( \prod _{\ell =1}^{n_{j}}\frac{\eta _{j}}{\nu _{j}(\ell )}\right) . \end{aligned}$$

(22)

3.2 Adjustment

We use the modified system (SOQN-LC) to approximate the SOQN-BO. First, we ensure that both systems process in the mean the same number of customers, i.e. they have the same throughput at (synchronisation node) 0. Our main idea is: To compensate customer losses, we increase the input rate $\lambda _{\text {LC}}$ of the modified system until the desired throughput is reached. We will prove this in Theorem 2. Before, we calculate the throughput of both systems.

Lemma 1

The throughput of the SOQN-BO in steady state is $\lambda _{\text {BO}}$. The throughput of the SOQN-LC in steady state is

$$\begin{aligned} \lambda _{\text {eff}}(\lambda _{\text {LC}})=\lambda _{\text {LC}}\cdot \Bigg (1-\underbrace{\frac{C^{\text {stb}}(\overline{J},N)}{C_{\text {LC}}(\overline{J}_{0},N)}}_{=\pi _{\text {LC,0}}(0)}\Bigg ) \end{aligned}$$

where $\pi _{\text {LC,0}}(0)$ is the probability of an empty resource pool in the SOQN-LC.

Proof

Because all customers pass the SOQN-BO, in steady state its throughput is $\lambda _{\text {BO}}$. In the SOQN-LC, a portion $\pi _{\text {LC,0}}(0)$ of the arrivals is lost. From the steady-state distribution (21), we find the idling probability $\pi _{\text {LC,0}}(0)$ of the resource pool in the SOQN-LC as

$$\begin{aligned} \pi _{\text {LC,0}}(0)&:=\sum _{\sum _{j\in \overline{J}}n_{j}=N}\pi _{\text {LC}}\left( 0,n_{1},\ldots ,n_{J}\right) \overset{(21)}{=}\frac{C^{\text {stb}}(\overline{J},N)}{C_{\text {LC}}(\overline{J}_{0},N)}. \end{aligned}$$

(23)

Then the effective arrival rate $\lambda _{\text {eff}}(\lambda _{\text {LC}})$, which coincides with the throughput of the system, is $\lambda _{\text {eff}}(\lambda _{\text {LC}})=\lambda _{\text {LC}}\cdot \left( 1-\pi _{\text {LC,0}}(0)\right) .$ $\square$

We adjust $\lambda _{\text {LC}}$ so that both systems have the same throughput. We assume that both systems are stable. For the SOQN-BO, by Proposition 1 stability is equivalent to $\lambda _{\text {BO}}\in \left( 0,\lambda _{\text {BO,max}}\right)$. For the SOQN-LC, stability is granted for any arrival rate $\lambda _{\text {LC}}\in \left( 0,\infty \right)$, because the state space $E_{LC}$ is finite.

Theorem 2

For every stable SOQN-BO, there exists an SOQN-LC with arrival rate $\lambda _{\text {LC}}$ such that both systems have the same throughput in steady state. Formally, with $\lambda _{\text {BO,max}}$ from (2) this means

$$\begin{aligned} \text {For all }\lambda _{\text {BO}}\in \left( 0,\lambda _{\text {BO,max}}\right) \ \text {exists }\lambda _{\text {LC}}\in \left( 0,\infty \right) \ \text {with}\ \lambda _{\text {eff}}(\lambda _{\text {LC}})=\lambda _{\text {BO}}. \end{aligned}$$

Proof

We show that for any $\lambda _{\text {BO}}\in \left( 0,\lambda _{\text {BO,max}}\right)$, the function $\lambda _{\text {eff}}(\cdot )$ from Lemma 1 takes values larger and smaller than a prescribed $\lambda _{\text {BO}}$, and is continuous. By the intermediate value theorem, there exists $\lambda _{\text {LC}}$ with $\lambda _{\text {eff}}(\lambda _{\text {LC}})=\lambda _{\text {BO}}$.

(i) $\lambda _{\text {eff}}(\lambda _{\text {LC}})=\lambda _{\text {LC}}\cdot \left( 1-\frac{C^{\text {stb}}(\overline{J},N)}{C_{\text {LC}}(\overline{J}_{0},N)}\right)$ can be larger than any $\lambda _{\text {BO}}$. We have

$$\begin{aligned} C_{\text {LC}}(\overline{J}_{0},N)&=\sum _{\sum _{j\in \overline{J}_{0}}n_{j}=N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}\cdot \prod _{j=1}^{J}\left( \prod _{i=1}^{n_{j}}\frac{\eta _{j}}{\nu _{j}(i)}\right) \nonumber \\&=\sum _{n_{0}=0}^{N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}\cdot \sum _{\sum _{j\in \overline{J}}n_{j}=N-n_{0}}\prod _{j=1}^{J}\left( \prod _{i=1}^{n_{j}}\frac{\eta _{j}}{\nu _{j}(i)}\right) \nonumber \\&=\sum _{n_{0}=0}^{N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}\cdot C^{\text {stb}}(\overline{J},N-n_{0}). \end{aligned}$$

(24)

To simplify notation, we define $b(n_{0}):=C^{\text {stb}}(\overline{J},N-n_{0})$. Then

$$\begin{aligned}&\lambda _{\text {eff}}(\lambda _{\text {LC}}) =\lambda _{\text {LC}}\cdot \left( 1-\frac{C^{\text {stb}}(\overline{J},N)}{C_{\text {LC}}(\overline{J}_{0},N)}\right) \\&\quad=\lambda _{\text {LC}} \cdot \left( 1-\frac{b(0)}{\sum _{n_{0}=0}^{N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}b(n_{0})}\right) \\&\quad =\lambda _{\text {LC}}\cdot \left( \frac{\sum _{n_{0}=1}^{N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}b(n_{0})}{\sum _{n_{0}=0}^{N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}b(n_{0})}\right) \\&\quad =\eta _{0}\cdot \left( \frac{\sum _{n_{0}=1}^{N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}-1}b(n_{0})}{b(0)+\sum _{n_{0}=1}^{N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}b(n_{0})}\right) \\&\quad =\eta _{0}\cdot \left( \frac{b(1)+\sum _{n_{0}=2}^{N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}-1}b(n_{0})}{b(0)+\sum _{n_{0}=1}^{N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}b(n_{0})}\right) . \end{aligned}$$

Hence, it holds

$$\lim _{\lambda _{\text {LC}}\rightarrow \infty }\lambda _{\text {eff}}(\lambda _{\text {LC}})=\eta _{0}\cdot \frac{b(1)}{b(0)}=\lambda _{\text {BO,max}}. $$

Therefore, $\lambda _{\text {eff}}(\lambda _{\text {LC}})$ can be larger than any arrival rate $\lambda _{\text {BO}}\in \left( 0,\lambda _{\text {BO,max}}\right)$.

(ii) $\lambda _{\text {eff}}(\lambda _{\text {LC}})$ can be smaller than any $\lambda _{\text {BO}}\in \left( 0,\lambda _{\text {BO,max}}\right)$, because

$$\begin{aligned} \lim _{\lambda _{\text {LC}}\rightarrow 0}\lambda _{\text {eff}}(\lambda _{\text {LC}})=\lim _{\lambda _{\text {LC}}\rightarrow 0}\lambda _{\text {LC}}\cdot \underbrace{\left( 1-\pi _{\text {LC,0}}(0)\right) }_{>0\ \text {and }<1}=0. \end{aligned}$$

(iii) From $\lambda _{\text {eff}}(\lambda _{\text {LC}})=\lambda _{\text {LC}}\cdot \left( 1-\frac{b(0)}{\sum _{n_{0}=0}^{N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}\cdot b(n_{0})}\right)$ follows that $\lambda _{\text {eff}}$ is a continuous function of $\lambda _{\text {LC}}\in \left( 0,\infty \right)$, which proves our claim by the intermediate value theorem. $\square$

Henceforth, we call $\lambda _{\text {LC}}$ with $\lambda _{\text {eff}}(\lambda _{\text {LC}})=\lambda _{\text {BO}}$ adjusted arrival rate for $\lambda _{\text {BO}}$.

Explicit results for adjusted $\lambda _{\text {LC}}$ for $N=1$ and $N=2$ can be found in Remark 3 in Appendix A. The proof of Proposition 6 is given in Appendix A.

Proposition 6

If the service rates $\nu _{j}(\cdot )$, $j\in \overline{J}$, are non-decreasing, $\lambda _{\text {LC}}$ in Theorem 2is unique.

3.3 Throughputs and idle times

Theorem 2 only guarantees that for a stable SOQN-BO with arrival rate $\lambda _{\text {BO}}$ an SOQN-LC with the same resource network and adjusted arrival rate exists such that the resource pools have the same throughput. Because an SOQN-LC can be investigated with standard Gordon–Newell network techniques, all local throughputs can be computed. We shall prove, that these local throughputs are the same as those in the SOQN-BO. This suggests to use local performance characteristics of the queues in the resource network of the SOQN-LC as approximation for performance measures of the SOQN-BO.

Proposition 7

The local throughput $TH_{\text {LC},j}$ at nodes $j\in \overline{J}_{0}$ in the SOQN-LC with adjusted arrival rate is pairwise the same as that of the respective nodes in the SOQN-BO given in Proposition 2. With

$$\begin{aligned}&C_{\text {LC}}(\overline{J}_{0},N)=\sum _{n_{0}=0}^{N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}\sum _{\sum _{j\in \overline{J}}n_{j}=N-n_{0}}\prod _{j=1}^{J}\left( \prod _{\ell =1}^{n_{j}}\frac{\eta _{j}}{\nu _{j}(\ell )}\right) \end{aligned}$$

it holds:

$$TH_{\text {LC},j}=\eta _{j}\cdot \frac{C_{\text {LC}}(\overline{J}_{0},N-1)}{C_{\text {LC}}(\overline{J}_{0},N)}=\lambda _{\text {BO}}\cdot \frac{\eta _{j}}{\eta _{0}},\quad j\in \overline{J}_{0}.$$

Proof

It was shown in the proof of Theorem 2 that with $\lambda _{\text {LC}}$ as adjusted arrival rate for $\lambda _{\text {BO}}$ it holds

$$\begin{aligned} TH_{\text {LC},0}=\lambda _{\text {eff}}(\lambda _{\text {LC}})=\eta _{0}\cdot \frac{C_{\text {LC}}(\overline{J}_{0},N-1)}{C_{\text {LC}}(\overline{J}_{0},N)}. \end{aligned}$$

It is well known (Chen and Yao 2001, Section 2.3), that the joint queue length vector $Z_{\text {LC}}:= \left( \left( Y_{j}(t):j\in \overline{J}_{0}\right) : t\ge 0\right)$ behaves stochastically as the joint queue length vector of a Gordon–Newell network consisting of nodes in $\overline{J}_{0}$, where node 0 has service rate $\lambda _{\text {LC}}$ and nodes in $\overline{J}$ have the same characteristics as in the SOQN-LC. Because $\lambda _{\text {eff}}(\lambda _{\text {LC}})= \lambda _{\text {BO}}$ we have $\lambda _{\text {BO}}=\eta _{0}\cdot \frac{C_{\text {LC}}(\overline{J}_{0},N-1)}{C_{\text {LC}}(\overline{J}_{0},N)}$ and from the formula for throughputs in Gordon–Newell networks it follows

$$\begin{aligned} TH_{\text {BO},j}=\lambda _{\text {BO}}\cdot \frac{\eta _{j}}{\eta _{0}}=\eta _{j}\cdot \frac{C_{\text {LC}}(\overline{J}_{0},N-1)}{C_{\text {LC}}(\overline{J}_{0},N)}=TH_{\text {LC},j},\quad j\in \overline{J}_{0}. \end{aligned}$$

$\square$

The explicit formulas for the throughputs in Proposition 7 allow to determine efficiently the steady-state marginal distribution of the queue length at every node $j\in \overline{J}$ without knowing the adjusted value $\lambda _{\text {eff}}(\lambda _{\text {LC}})$ of $\lambda _{\text {LC}}$. This leads especially to

Proposition 8

Let $Y_{\text {LC}}:=(Y_{\text {LC},j}:j\in \overline{J})$ denote a random vector which is distributed according to the stationary queue length at the nodes in $\overline{J}$ of the SOQN-LC with adjusted arrival rate.

If the service rate at node j does not depend on the queue length, i.e. $\nu _{j}(\cdot )=\nu _{j},j\in \overline{J}$, then the probabilities that the nodes $j\in \overline{J}_{0}$ in the SOQN-LC with adjusted arrival rate are idling are pairwise the same as those of the respective nodes in the SOQN-BO given in Corollary 1:

$$\begin{aligned} P(Y_{\text {LC},j}=0)=1-\lambda _{\text {BO}}\cdot \frac{\eta _{j}}{\eta _{0}}\cdot \nu _{j}^{-1}. \end{aligned}$$

Proof

According to Proposition 7, the throughputs are equal. The rest of the proof is the same as the proof of Corollary 1. $\square$

4 Approximation of the external queue

We have shown that, after adjusting $\lambda _{\text {LC}}$, the behaviour of the resource network of the original SOQN-BO can be approximated well by the behaviour of the resource network of the modified SOQN-LC. The behaviour of the external queue of the SOQN-BO is represented only by the modified arrival intensity of the associated SOQN-LC. On the other side characteristics of the external queue are important performance measures of the original system. We solve this problem in a two-step approach to approximate the external queue:

Step 1:

In Sect. 4.1, we reduce the modified system to a simpler SOQN-LC.

Step 2:

In Sect. 4.2, we combine the results from Sects. 4.1 and  2.4, obtaining a simple SOQN-BO to approximate the external queue.

4.1 Reduced SOQN with lost customers

Because the joint queue length vector $Z_{\text {LC}}:= \left( \left( Y_{j}(t):j\in \overline{J}_{0}\right) : t\ge 0\right)$ of the SOQN-LC from Sect. 3.1 can be studied via a Gordon–Newell network, we can reduce complexity further by applying Norton’s theorem (Chandy et al. 1975) to construct a two-node Gordon–Newell network as shown in Fig. 5, with the same throughput.

Fig. 5

Step 1: Reduction of complexity

The inner network is replaced by a composite node ($\overline{J}:=\left\{ 1\right\}$): A single exponential-1-server with infinite waiting room under FCFS regime with a queue-length-dependent service intensity $\varphi (\cdot )$ (Chandy et al. 1975, p. 39, eq. (20)), with

$$\begin{aligned} \varphi (0) =0,\qquad \varphi (m) =\eta _{0}\cdot \frac{C^{\text {stb}}(\overline{J},m-1)}{C^{\text {stb}}(\overline{J},m)},\quad m\in \left\{ 1,\ldots ,N\right\} . \end{aligned}$$

(25)

Remarkably, $\varphi (N)$ equals $\lambda _{\text {BO,max}}$ in (2). We deduce from (10) that

$$\begin{aligned} \varphi (m)=TH_{0}^{\text {stb}}(m),\quad m\in \left\{ 1,\ldots ,N\right\} , \end{aligned}$$

(26)

which is independent of $\lambda _{\text {LC}}$. The normalisation constants $C^{\text {stb}}(\overline{J},m)$, $m\in \left\{ 0,\ldots ,N\right\}$, can be calculated by the convolution algorithm or mean value analysis (MVA), see (Bolch et al. 1998, p. 371ff., Section 8.1 and p. 384ff., Section 8.2).

4.2 Back to backordering

Recall that Theorem 2 guarantees, that for every stable SOQN-BO there exists an SOQN-LC with adjusted arrival rate $\lambda _{\text {LC}}$ such that both systems have the same throughput in steady state. Our next step is to apply an “inversion” of that construction to the reduced SOQN-LC on the right side of Fig. 5. This results in removing the lost-customer property to get again the backordering property as shown in Fig. 6. The recipe is simple: We reopen the external queue and reduce the arrival intensity $\lambda _{\text {LC}}$ to $\lambda _{\text {BO}}$ while the service rate at the single queue of the inner network is $\nu _{1}(m):=\varphi (m)=TH_{0}^{\text {stb}}(m)$ from Eq. (26). The final result is: The external queue of the large SOQN-BO with $J>1$ approximated by a reduced SOQN-BO with $J=1$. This yields the following approximating performance characteristics for the external queue of the SOQN-BO of Sect. 2.

Fig. 6

Step 2: Transition from reduced SOQN-LC to reduced SOQN-BO

Denote by $\left( X_{\text {ex}},\mathbf {Y}\right)$ a vector distributed according to the stationary distribution of $Z_{\text {BO}}$, $L_{\text {ex}}$ the mean external queue length, and $W_{\text {ex}}$ the mean waiting time at the external queue length in steady state. These characteristics are approximated by the respective quantities from the system at the right side of Fig. 6, given in Proposition 3.

$$\begin{aligned} P(X_{\text {ex}}=0)\approx P(\widehat{X}_{\text {ex}}=0)\overset{(14)}{=} \left[ C_{\text {BO}}(\left\{ 1\right\} ,N)\right] ^{-1}\cdot \sum _{n_{1}=0}^{N}\prod _{m=1}^{n_{1}}\frac{\lambda _{\text {BO}}}{TH_{0}^{\text {stb}}(m)} \end{aligned}$$

and for $n_{\text {ex}}>0 :$

$$\begin{aligned} P(X_{\text {ex}}=0)&\approx P(\widehat{X}_{\text {ex}}=n_{\text {ex}})\\&\overset{(15)}{=}\left[ C_{\text {BO}}(\left\{ 1\right\} ,N)\right] ^{-1} \cdot \left( \frac{\lambda _{\text {BO}}}{TH_{0}^{\text {stb}}(N)}\right) ^{n_{\text {ex}}} \cdot \prod _{m=1}^{N}\frac{\lambda _{\text {BO}}}{TH_{0}^{\text {stb}}(m)} \end{aligned}$$

with

$$C_{\text {BO}}(\left\{ 1\right\} ,N)\overset{(13)}{=}\sum _{n_{1}=0}^{N-1}\prod _{m=1}^{n_{1}}\frac{\lambda _{\text {BO}}}{TH_{0}^{\text {stb}}(m)}+ \frac{1}{1-\frac{\lambda _{\text {BO}}}{TH_{0}^{\text {stb}}(N)}}\cdot \prod _{m=1}^{N}\frac{\lambda _{\text {BO}}}{TH_{0}^{\text {stb}}(m)}.$$

The mean external queue length is approximated with Eq. (17), (18) by

$$\begin{aligned}&L_{\text {ex}}\approx {\widehat{L_{\text {ex}}}}\nonumber \\&\quad =\frac{1}{C_{\text {BO}}(\left\{ 1\right\} ,N)} \cdot \frac{1}{1-\frac{\lambda _{\text {BO}}}{TH_{0}^{\text {stb}}(N)}}\cdot \prod _{m=1}^{N}\frac{\lambda _{\text {BO}}}{TH_{0}^{\text {stb}}(m)} \cdot \frac{\lambda _{\text {BO}}}{TH_{0}^{\text {stb}}(N)-\lambda _{\text {BO}}}. \end{aligned}$$

(27)

With our approximation we arrive at the same formula for $L_{\text {ex}}$ as (Dallery 1990, eq. (22)) obtained by aggregation technique in (Dallery 1990, Section 6). Using (19) the average waiting time of customers in the external queue is approximated as

$$\begin{aligned} W_{\text {ex}}\approx \widehat{W_{\text {ex}}} = \frac{\widehat{L_{\text {ex}}}}{\lambda _{\text {BO}}}. \end{aligned}$$

(28)

Remark: We expect that the results are close to the true values, but at present we do not have strict error bounds.

5 Application to RMFS

We evaluate analytically the performance of a robotic mobile fulfilment system (RMFS) modelled as an SOQN-BO. In an RMFS, robots are expensive resources. Therefore, we determine the minimal number of robots needed to stabilize the system or to maintain a required quality of service. From Sect. 2, an SOQN-BO can be described by a level-independent quasi-birth-death-process, so numerical schemes are at hand. Due to the large state space direct application of these matrix-geometric methods is not practical. For example, for 10 robots, we need to calculate ca. $\left( 9\cdot 10^{4}\right) ^{2}$ entries of a special matrix. Therefore, approximative matrix-geometric methods are developed.

5.1 Description of RMFS

The components and the order fulfilment processes in an RMFS with an illustrated example are depicted in Fig. 7. Central components are:

• movable shelves, called pods, on which items are stored,

• storage area—the area where the pods are stored,

• workstations, where

  • the items are picked from pods by pickers (picking stations) or

  • the items are stored to pods (replenishment stations),

• mobile robots, which can move underneath pods and carry them to workstations.

Figure 7 illustrates order fulfilment processes in an RMFS. On the upper left hand we have three customers’ orders which contain different items, distinguished by colours. To fulfil orders, we send them or parts of them to picking stations. To these stations we send pods with the necessary items. Each pod is carried by a robot. In this way, customers’ orders generate tasks for robots. The robots, with their pods, queue up in front of the picking stations. A picker takes all the necessary items from the pod at the head of the queue. Then he sends the pod with its robot back to the storage area. As soon as the customer’s order or part of it is fulfilled, we remove it from the picking station. The order in which we send the customers’ orders, how we split them apart, and which pod we send, is a complex topic, see Xie et al. (2021). In the present paper, we focus on the generated robots’ tasks, which we will call just tasks.

In the example, each customer’s order is split into two parts. Three parts are sent to picking station 1, and three to picking station 2. To fulfil these partial orders, a robot transports one pod to picking station 1, and another robot transports one pod to picking station 2, from the storage area.

From time to time, we need to refill pods. To do this, we send these pods to the replenishment station. There, employees refill the pods and send them back to the storage area. In the example, after picking, pod 2 is sent to the replenishment station to refill it with blue items.

Fig. 7

Order fulfilment processes in RMFS. The circled numbers refer to the processes in Fig. 8

5.2 Modelling as SOQN

An SOQN-BO as model for the RMFS from Fig. 7 is depicted in Fig. 8. The RMFS is open with respect to tasks and closed with respect to robots, which are the resources. It has two picking and one replenishment station.

Fig. 8

RMFS modelled as an SOQN-BO. The circled numbers refer to the processes in Fig. 7

Customers’ orders arrive at the RMFS with rate $\lambda _{\text {CO}}$ and generate tasks. The number of tasks, which a single customer's orders can generate, depends on many parameters. In particular, it depends on the efficiency of the algorithm which tries to find an optimal match between customers’ orders and pods. The matching problem is NP-hard, and, to the best of our knowledge, there are no formulas known to determine how many pods an order will need. Therefore, we assume that there exists an average pod/order ratio $\sigma _{\text {pod/order}}$ which we find empirically for an RMFS. We assume this ratio depends only on the pods’ contents and customers’ order contents and is independent of the number of robots.

The matching algorithm adds some delay in order to assign pods to orders. We assume that this delay depends only on the pods’ contents, customers’ order contents and order input rate and is independent of the number of robots. We assume that we can find the average delay $W_{\text {alg}}$ empirically for our RMFS.

Thus, the customers’ orders generate a stream of “bring a pod to a picking station” tasks with rate $\lambda _{\text {BO}}=\lambda _{\text {CO}}\cdot \sigma _{\text {pod/order}}>0$. The delay, introduced by the matching algorithm, does not change this rate.

We reduce the complexity of creating a task and model the task stream as a Poisson stream with rate $\lambda _{\text {BO}}=\lambda _{\text {CO}}\cdot \sigma _{\text {pod/order}}$. To be processed ($=$ to enter the inner network), each such task requires exactly one idle robot from the robot pool (resource pool), which is henceforth referred to as node 0. If there is no idle robot available, the new task has to wait in an external queue until a robot becomes available (“backordering”). The maximal number of robots in the resource pool is $N$. The inner network in the example in Fig. 8 consists of 11 nodes, denoted by

$$\begin{aligned} \overline{J}:=\left\{ \text {sp},\text {pp}_{1},\text {pp}_{2},\text {p}_{1},\text {p}_{2},\text {p}_{1}\text {s},\text {p}_{2}\text {s},\text {p}_{1}\text {r},\text {p}_{2}\text {r},\text {r},\text {rs}\right\} . \end{aligned}$$

The meaning and notations of nodes are given in Table 1. The robot with assigned task moves through the network. The following processes occur from the perspective of a robot:

• The idle robot waits to be assigned to a task (bring a particular pod).

• The robot moves with the assigned task to a pod.

• With this pod the robot moves with probability $q_{\text {pp}_{1}}\in \left( 0,1\right)$ to picking station 1 and with probability $q_{\text {pp}_{2}}\in \left( 0,1\right)$ to station 2, $q_{\text {pp}_{1}}+q_{\text {pp}_{2}}=1$.

• The robot queues with the pod at the picking stations.

• After picking at picking station 1 resp. picking station 2, the robot:

  • either

    • carries the pod directly back to the storage area with probability $q_{\text {p}_{1}\text {s}}\in \left( 0,1\right)$ resp. $q_{\text {p}_{2}\text {s}}\in \left( 0,1\right)$, or

    • moves to the replenishment station with probability $q_{\text {p}_{1}\text {r}}\in \left( 0,1\right)$ resp. $q_{\text {p}_{2}\text {r}}\in \left( 0,1\right)$, whereby $q_{\text {p}_{1}\text {s}}+q_{\text {p}_{1}\text {r}}=1$ resp. $q_{\text {p}_{2}\text {s}}+q_{\text {p}_{2}\text {r}}=1$,

  • queues at the replenishment station, and

  • carries the pod back to the storage area and waits for the next task.

Each of these processes is modelled as a queue. Movements of the robots are modelled by processor-sharing nodes with exponential service times with intensities $\nu _{j}(n_{j}):=\mu _{j}\cdot \phi _{j}(n_{j})$, $j\in \overline{J}\setminus \left\{ \text {p}_{1},\text {p}_{2},\text {r}\right\}$, presented in Table 1.

The two picking stations and the replenishment station, which are referred to as node $\text {p}_{1}$, node $\text {p}_{2}$ resp. node $\text {r}$, consist of a single server with waiting room under the FCFS regime. The picking times and the replenishment times are exponentially distributed with rates $\nu _{\text {p}_{1}}$, $\nu _{\text {p}_{2}}$ resp. $\nu _{r}$.

Table 1 Overview of the nodes in the network

The robots travel among the nodes following an irreducible routing matrix $\mathcal {R}:=\left( r(i,j):i,j\in \overline{J}_{0}\right)$, whereby $\overline{J}_{0}:=\left\{ 0\right\} \cup \overline{J}$, which is given by

$$\begin{aligned} \mathcal {R}=\left( \begin{array}{c|cccccccccccc} &{} 0&{} \text {sp}&{} \text {pp}_{1}&{} \text {pp}_{2}&{} \text {p}_{1}&{} \text {p}_{2}&{} \text {p}_{1}\text {s}&{} \text {p}_{2}\text {s}&{} \text {p}_{1}\text {r}&{} \text {p}_{2}\text {r}&{} \text {r}&{} \text {rs}\\ \hline 0&{} &{} 1\\ \text {sp}&{} &{} &{} q_{\text {pp}_{1}}&{} q_{\text {pp}_{2}}\\ \text {pp}_{1}&{} &{} &{} &{} &{} 1\\ \text {pp}_{2}&{} &{} &{} &{} &{} &{} 1\\ \text {p}_{1}&{} &{} &{} &{} &{} &{} &{} q_{\text {p}_{1}\text {s}}&{} &{} q_{\text {p}_{1}\text {r}}\\ \text {p}_{2}&{} &{} &{} &{} &{} &{} &{} &{} q_{\text {p}_{2}\text {s}}&{} &{} q_{\text {p}_{2}\text {r}}\\ \text {p}_{1}\text {s}&{} 1\\ \text {p}_{2}\text {s}&{} 1\\ \text {p}_{1}\text {r}&{} &{} &{} &{} &{} &{} &{} &{} &{} &{} &{} 1\\ \text {p}_{2}\text {r}&{} &{} &{} &{} &{} &{} &{} &{} &{} &{} &{} 1\\ \text {r}&{} &{} &{} &{} &{} &{} &{} &{} &{} &{} &{} &{} 1\\ \text {rs}&{} 1 \end{array}\right) . \end{aligned}$$

We define the joint stochastic process Z of this system by

$$\begin{aligned} Z&:=\Big (\Big (X_{\text {ex}}(t),Y_{0}(t),Y_{\text {sp}}(t),Y_{\text {pp}_{1}}(t),Y_{\text {pp}_{2}}(t),Y_{\text {p}_{1}}(t),Y_{\text {p}_{2}}(t),Y_{\text {p}_{1}\text {s}}(t),Y_{\text {p}_{2}\text {s}}(t),\\&\quad \ \ Y_{\text {p}_{1}\text {r}}(t),Y_{\text {p}_{2}\text {r}}(t),Y_{\text {r}}(t),Y_{\text {rs}}(t)\Big ):t\ge 0\Big ). \end{aligned}$$

Due to the usual independence and memoryless assumptions, Z is an irreducible Markov process with state space

$$\begin{aligned} E&:=\Bigg \{\left( 0,k_{\text {idle robots}},n_{\text {sp}},n_{\text {pp}_{1}},n_{\text {pp}_{2}},n_{\text {p}_{1}},n_{\text {p}_{2}},n_{\text {p}_{1}\text {s}},n_{\text {p}_{2}\text {s}},n_{\text {p}_{1}\text {r}},n_{\text {p}_{2}\text {r}},n_{\text {r}},n_{\text {rs}}\right) :\\&\qquad n_{j}\in \left\{ 0,\ldots ,N\right\} \ \forall j\in \overline{J}_{0},\ \sum _{j\in \overline{J}_{0}}n_{j}=N\Bigg \}\\&\quad \cup \Bigg \{\left( n_{\text {ex}},0,n_{\text {sp}},n_{\text {pp}_{1}},n_{\text {pp}_{2}},n_{\text {p}_{1}},n_{\text {p}_{2}},n_{\text {p}_{1}\text {s}},n_{\text {p}_{2}\text {s}},n_{\text {p}_{1}\text {r}},n_{\text {p}_{2}\text {r}},n_{\text {r}},n_{\text {rs}}\right) :\\&\qquad n_{\text {ex}}\in \mathbb {N},\,n_{j}\in \left\{ 0,\ldots ,N\right\} \ \forall j\in \overline{J},\ \sum _{j\in \overline{J}}n_{j}=N\Bigg \}. \end{aligned}$$

5.3 Determine the minimal number of robots

The throughput $TH_{0}^{\text {stb}}(N)$ in Eq. (9) depends on the number of robots $N$. To find the minimal number of robots that stabilises the system, we check the stability criterion from Proposition 1. The maximal number of robots $N^{\max }$ is the number of pods or is determined by financial restrictions. Algorithm 1 determines the set ${\overline{N}}^{*}$of feasible numbers of robots for a stable system.

Remark 2

If $TH_{0}^{\text {stb}}(\cdot )$ is non-decreasing in ${\mathbb {N}}$, the algorithm can be shortened: Starting with one robot and adding a new robot in each step until the stability criterion is satisfied for the first time. Sufficient conditions for non-decreasing $TH_{0}^{\text {stb}}(\cdot )=\lambda _{\text {BO,max}}$ are given in Propositions 5 and 4 .

Stability does not guarantee acceptable turnover times of orders. Therefore, we consider additionally the turnover time of customers’ order. The turnover time of a customer’s order can be split into three main parts:

• 1. Waiting time until the matching algorithm has assigned all required pods to that order. By assumption, this time does not depend on the number of robots and is on average $W_{\text {alg}}>0$.

• 2. Waiting time of the first matched pod for an idle robot, time for transport to the picking station, waiting time for the picker at the picking station. During all these times, the order is coupled with at least one task. We call this turnover time for the task $TO_{\text {task}}(\lambda _{\text {LC}},N)$.

• 3. Time of an order between start of picking and its completion. This time depends on many factors, for example: How many orders can a picker complete with the same pod? Will all completed orders wait until a pod leaves? Is the order’s content in multiple pods? Will these pods arrive right after each other, or will there be many pods for other orders in between? Is there any complex merging procedure outside of the picking station? In our model, we use a simplifying assumption that the order needs on average $W_{\text {assembled}}>0$ from the time its first pod arrives at the picking station until the time picking for this order is completed.

With these assumptions, we can assume that the turnover time of an order is

$$\begin{aligned} TO_{\text {order}}(\lambda _{\text {LC}},N):=W_{\text {alg}}+TO_{\text {task}}(\lambda _{\text {LC}},N)+W_{\text {assembled}}. \end{aligned}$$

Even if $W_{\text {alg}}$ and $W_{\text {assembled}}$ are unknown, we can still use $TO_{\text {task}}(\lambda _{\text {LC}},N)$ as a lower bound for $TO_{\text {order}}(\lambda _{\text {LC}},N)$.

Because of the simplifying assumption about $W_{\text {alg}}$ and $W_{\text {assembled}}$, only the turnover time $TO_{\text {task}}(\lambda _{\text {LC}},N)$ of a task depends on $N$, and for the minimal number of robots we can focus on this.

The turnover time $TO_{\text {task}}(\lambda _{\text {LC}},N)$ of a task is measured from the time the task is received to the time the picker starts to process it:

$$\begin{aligned} TO_{\text {task}}(\lambda _{\text {LC}},N):=W_{\text {ex}}(N)+W_{\text {in}}(\lambda _{\text {LC}},N). \end{aligned}$$

$W_{\text {ex}}(N)$ is the average time a task spends waiting in the external queue until it enters the inner network. We can calculate it with Eq. (28). $W_{\text {in}}(\lambda _{\text {LC}},N)$ is the average time a task spends in the inner network until a picker starts to process it at one of the picking stations. Given the average waiting times $W_{j}(\lambda _{\text {LC}},N)$ at nodes $j\in \overline{J}=\left\{ \text {sp},\text {pp}_{1},\text {pp}_{2},\text {p}_{1},\text {p}_{2},\text {p}_{1}\text {s},\text {p}_{2}\text {s},\text {p}_{1}\text {r},\text {p}_{2}\text {r},\text {r},\text {rs}\right\}$ from arrival until service completion, and constant service rates $\nu _{j}$ at nodes $j\in \{\text {p}_{1},\text {p}_{2}$}, then

$$\begin{aligned} W_{\text {in}}(\lambda _{\text {LC}},N)&:=W_{\text {sp}}(\lambda _{\text {LC}},N)\\&\quad +r(\text {sp},\text {pp}_{1})\cdot \left( W_{\text {pp}_{1}}(\lambda _{\text {LC}},N)+W_{\text {p}_{1}}(\lambda _{\text {LC}},N)-1/\nu _{\text {p}_{1}}\right) \\&\quad +r(\text {sp},\text {pp}_{2})\cdot \left( W_{\text {pp}_{2}}(\lambda _{\text {LC}},N)+W_{\text {p}_{2}}(\lambda _{\text {LC}},N)-1/\nu _{\text {p}_{2}}\right) . \end{aligned}$$

We calculate $W_{j}(\lambda _{\text {LC}},N)$, $j\in \overline{J}$, with MVA.

In Algorithm 2, we determine the minimal number of robots for an acceptable turnover time of a task. We will call this time $TO_{task}^{max}$.

5.4 Numerical experiments

In this section, we show by example (with data taken from the literature) how to apply our approximation for SOQNs. We compare the results with simulations of the original SOQN, considering more details of the network. We investigate the performance of the system under prescribed utilisation levels for robots and discuss the quality of the approximation. Eventually, we investigate the effect of interarrival time variability on waiting times at some stations of the SOQN.

In our experiments, we take parameters from (Lamballais et al. 2019, Table 5.3 and Table 5.4). The maximal number of pods is $N^{\max }=550$, arrival rate of tasks is $468\,\frac{{\text {tasks}}}{{\text {h}}}=0.13\,\frac{{\text {tasks}}}{{\text {s}}}$ (in (Lamballais et al. 2019, Table 5.3 and Table 5.4), arrival rates are given in [order/hour]. Because every order generates one task, we use [task/hour] directly). Average travel time at node $\text {sp}$: $\mu _{\text {sp}}^{-1}={18.4}\text { s}$, at node $\text {pp}_{1}$: $\mu _{\text {pp}_{1}}^{-1}={34.5}\text { s}$, at node $\text {pp}_{2}$: $\mu _{\text {pp}_{2}}^{-1}={34.5}\text { s}$, at node $\text {p}_{1}\text {s}$: $\mu _{\text {p}_{1}\text {s}}^{-1}={34.5}\text { s}$, at node $\text {p}_{2}\text {s}$: $\mu _{\text {p}_{2}\text {s}}^{-1}={34.5}\text { s}$, at node $\text {p}_{2}\text {r}$: $\mu _{\text {p}_{2}\text {r}}^{-1}={34.5}\text { s}$, at node $\text {p}_{1}\text {r}$: $\mu _{\text {p}_{1}\text {r}}^{-1}={34.5}\text { s}$, at node $\text {rs}$: $\mu _{\text {rs}}^{-1}={34.5}\text { s}$. Average picking time of picking station 1: $\nu _{\text {p}_{1}}^{-1}={10}\text { s}$, of picking station 2: $\nu _{\text {p}_{2}}^{-1}={10}\text { s}$, average replenishment time (node $\text {r}$): $\nu _{r}^{-1}={30}\text { s}$.

We assume that moving robots do not interfere. Hence, our processor-sharing queues are infinite server queues, i.e. $\phi _{j}(n_{j})=n_{j}$ for all $j\in \overline{J}\setminus \left\{ \text {p}_{1},\text {p}_{2},\text {r}\right\}$.

We implemented our algorithm in R and used the package queueing, see (Canadilla 2017). In the worst case scenario—when we need to try all ($550-18+1$) of the robots—our implementation needs on average 83 seconds on a notebook with an i7-7600U CPU processor, 2.80GHz and 16GB RAM. We plotted important system parameters in Figs. 9, 10, 11 and 12. For better readability, we plotted data for a limited number of robots (until stabilising of curves).

Figure 9 shows maximal arrival rates $\lambda _{\text {BO}}$ for given numbers $N$ of robots to keep the system stable. In particular, the minimal number of robots for the system to be stable is 18 and for more than 40 robots, additional robots do not allow significantly higher arrival rates.

Figure 10 shows the throughputs of the nodes. By Proposition 2, these throughputs do not depend on $N$, and are pairwise the same for the original SOQN-BO and the adjusted approximation SOQN-LC. Idling probabilities of nodes $\text {p}_{1}$, $\text {p}_{2}$ and $\text {r}$ are 0.35, 0.35, and 0.22, obtained with Corollary 1.

Figure 11 shows adjusted arrival rates $\lambda _{\text {LC}}$ for an SOQN-LC to obtain an effective arrival rate $\lambda _{\text {eff}}(\lambda _{\text {LC}})=\lambda _{\text {BO}}=468\, \text {h}^{-1}$. Because in an SOQN-LC with many robots the probability of an empty resource pool is $\approx 0$, only few customers are lost: $\lambda _{\text {LC}}$ has to be adjusted only slightly, so $\lim _{N\rightarrow \infty }\lambda _{\text {LC}}\approx \lambda _{\text {BO}}$.

Figure 12 shows average turnover times. From arrival of an order (= one task) at the system until it is completed at a picking station. With only 18 robots, turnover times are extremely large, although the system is stable. Remarkably, turnover times decrease dramatically with only one additional robot.

Fig. 9

Maximal arrival rates $\lambda _{\text {BO}}$ for given numbers of robots to keep the system stable

Fig. 10

Throughputs for each node of the RMFS example

Fig. 11

Adjusted arrival rate $\lambda _{\text {LC}}$ for a system with lost customers such that the effective arrival rate is $\lambda _{\text {BO}}= 468\, \text {h}^{-1}$

Fig. 12

Average turnover time $TO_{\text {task}}(\lambda _{\text {LC}},N)$ of a task, which is the average delay of a task until it is completed

Comparison with simulation. We simulated the SOQN-BO model of the RMFS for 365 days 20 times for each number of robots using SimPy 3.0. Figure 13 shows mean waiting times in the external queue starting with 19 robots. The simulated SOQN-BO and the SOQN-LC approximation show the same qualitative behaviour. Although the approximation under-estimates the “true values” (obtained by simulation) in the region of 19-25 robots, the approximation answers the question “how many robots do we need to obtain a target turnover time” quite well. From 26 robots on, the approximation reflects the behaviour of the original system well. We omitted results for 18 robots because then the system operates on the edge of instability (utilisation $\approx 1$) with extremely large mean value and standard deviation. We ran 200 additional simulations for this setting. Figure 14 shows the large variability of the average waiting time. So, operating a system under such conditions cannot be recommended. More details on utilisation for the system are shown in Fig. 15.

Figure 16 shows that our approximation works well for turnover times. To assess the quality of approximation better, recall that turnover times consist of transportation and waiting times for a picker. Average transportation times are equal in the original system and in the approximation. We calculate them from service times at appropriate nodes: $\mu _{\text {sp}}^{-1}+r(\text {sp},\text {pp}_{1})\cdot \mu _{\text {pp}_{1}}^{-1}+r(\text {sp},\text {pp}_{2})\cdot \mu _{\text {pp}_{1}}^{-1}$ $=0.0147\text {h}$. The hard part is to estimate average waiting times for pickers. Figure 17 shows that the approximation is good.

Although mean waiting times for replenishments are not needed for our optimisation, Fig. 18 demonstrates that our algorithm estimates further network characteristics well. (We admit that less impressive results may occur.)

Fig. 13

Average waiting times in the external queue, simulation vs. approximation

Fig. 14

Distribution of waiting times in the external queue for a system with 18 robots, obtained with 200 simulations

Fig. 15

Robots’ utilisation: Proportion of time robots do not wait for assigned tasks, simulation vs. approximation

Fig. 16

Average turnover times, simulation vs. approximation. The graphs for approximation and simulation coincide. The large fraction of turnover times for transportation is the same in both systems by construction

Fig. 17

Average waiting times for pickers, simulation vs. approximation

Fig. 18

Average waiting times for replenishment, simulation vs. approximation

Utilisation-dependent comparisons. Figures 13 and 15 emphasize the importance of systems’ utilisation for a good fit of our lost-customers approximation of SOQNs. In this section, we study the goodness-of-fit for the approximation under prescribed utilisations of the robots.

The network’s structure (inner network and resource pool) is the same as before. We consider robot utilisations 0.95, 0.90, 0.80, 0.70 and 0.60 and determine for given number of robots $N=1,2,\dots ,50$ the external arrival rate to the SOQN-BO which yields the required utilisation. Because simulation of the SOQN-BO is (especially for high utilisations) extremely time consuming we used an iteration based on the approximation procedure described in Sects. 2 and 3. For a guess of $\lambda _{BO}$ we determine in the associated lost customer approximation with adjusted arrival rate the robots’ utilisation. This is iterated until we meet the prescribed utilisation. So, by construction the utilisation of the approximation is in any case precisely of the required size. We checked whether the utilisation of the SOQN-BO is the same with that $\lambda _{BO}$ (depending on N) for $N=1,2,\dots ,50$. The results are presented in Fig. 19 and confirm that deviations are moderate: For $N\le 20$ the results fit well, and even in the worst cases (for $N=50$) the relative deviation is below $10\%$.

The deviation of the approximation of the average turnover times from the respective simulations are shown in Fig. 20 and are almost neglectable.

Fig. 19

Robots’ utilisation, simulation vs. approximation

Fig. 20

Average turnover times, simulation vs. approximation

Fig. 21

Average waiting times in the external queue, simulation vs. approximation

For approximations of the external queue the situation is different. In Fig. 21 we see that our method reproduces the behaviour of the system qualitatively pretty well: The shapes of the average waiting time curves for fixed utilisation are nearly the same. But the differences between the respective values of the simulated and approximated average waiting times for fixed N are in parts large, especially for high utilisations and high numbers of robots. (We presented only results for average waiting times up to $0.07\text {h}^{-1}$. Behaviour of curves for larger delay is similar.)

An explanation for the deviations of the approximated external mean waiting times from the simulated values may be as follows. Consider the resource network (inner network + resource pool) as a black box server with complicated service time structure for an approximately Poissonian arrival stream (Fig. 2). From the Pollaćek-Khintchin formula for mean waiting times we know that even in a simple $M/G/1/\infty$ queue (with state-independent service rate) the variability of the service time is crucial for the size of the mean waiting time at the system. We conjecture that our last approximation step which reduces the complex inner network to a single server queue with state dependent service rates does not encapsulate the variability of the passage time through the inner network sufficiently precise. We studied this problem theoretically in a simple model which fits into the class of systems described in Sect. 2. Surprisingly, we found all cases: Increasing, decreasing, and invariant variability of passage time variance in the reduction step, see the comment at the end of Appendix A.

Arrival streams We assumed so far that the arrival streams to the original SOQN and to the approximating system are Poissonian. If only mean interarrival time $\lambda ^{-1}>0$ is available, this is the classical assumption based on entropy arguments: Exponential distribution has maximal entropy among continuous distributions on $[0,\infty )$ with this mean. On the other side the deterministic interarrival variable with mean $\lambda ^{-1}>0$ has entropy 0. We performed experiments comparing SOQNs with internal structure as in the introductory example in the beginning of Sect. 5.4, now with deterministic and exponential interarrival times, comparing this with our approximation using Poissonian arrivals.

As worst case we observed results for average waiting times at the picker, see Fig. 22, with largest difference between simulation of the SOQN-BO with deterministic interarrival times and approximation using Poisson arrivals. While up to 25 robots the deviation of the approximation is below $\approx 10\%$, for 50 robots the deviation of the approximation from the simulation result is $\approx 20\%$. We included in Fig. 22 simulation results for the SOQN-BO with Poisson arrivals which fit very well with the approximation.

Figure 23 shows average waiting times at the external queue similar to Fig. 13. The approximation provides results in the midst between the simulation results with Poisson and deterministic arrival streams. We conclude that approximation of the external queue works well in both cases.

This in-between-property holds for the average waiting times at the picker as well, see Fig. 22, but it is not a universal property as can be seen from Fig. 24. We investigated waiting times at the replenisher. In this case for less than 32 robots the approximation under-estimates the average waiting time at the replenisher under deterministic arrivals, while for larger number of robots this average waiting time is over-estimated, although in the worst case for 50 robots for less than $\approx 10\,\%$, which seems to be satisfying.

Fig. 22

Average waiting time for the picker. Simulation with Poisson and deterministic arrival times and approximation

Fig. 23

Average waiting times in the external queue. Simulation with Poisson and deterministic arrivals vs. approximation

Fig. 24

Average waiting times for the replenisher. Simulation with Poisson and deterministic arrivals vs. approximation

6 Conclusion

Our contribution to modelling and performance evaluation of SOQNs is based on an interplay of exact procedures, partly originating from queueing network theory, and heuristic transformation of SOQNs with infinite external queue into SOQNs with finite external queue and vice-versa. Before realising this transformation we have investigated performance indices of the original SOQN to obtain exact stability conditions and some directly accessible mean values (throughputs) and steady-state probabilities. The advantage of the procedure is to obtain a closed form steady-state distribution for the external queue and the total population size in the inner network.

The subsequent application of our results to performance analysis of an RMFS validates the applicability of the obtained performance evaluation methods. An advantage of our procedure is the possibility to use in applications well-established performance algorithms from queueing network theory.

Proof of Proposition 6 and further details

Proof of Proposition 6

We show strict isotonicity of $\lambda _{\text {eff}}(\lambda _{\text {LC}})$ by induction in $N$. We write $\lambda _{\text {eff}}^{(N)}(\lambda _{\text {LC}})$ in a system with $N$ resources. Recall from Eq. (22) for $L=1,2,\ldots ,N$

$$\begin{aligned} C_{\text {LC}}(\overline{J}_{0},L)=\sum _{n_{0}=0}^{L}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}\sum _{\sum _{j\in \overline{J}}n_{j}=L-n_{0}}\prod _{j=1}^{J}\left( \prod _{i=1}^{n_{j}}\frac{\eta _{j}}{\nu _{j}(i)}\right) . \end{aligned}$$

We set $C_{\lambda _{\text {LC}}}(\overline{J}_{0},L):=C_{\text {LC}}(\overline{J}_{0},L)$ to indicate functional dependence on $\lambda _{\text {LC}}$. It holds

$$\begin{aligned} \lambda _{\text {eff}}^{(N)}(\lambda _{\text {LC}})&=\lambda _{\text {LC}}\cdot \left( 1-\frac{C^{\text {stb}}(\overline{J},N)}{C_{\text {LC}}(\overline{J}_{0},N)}\right) \\&=\lambda _{\text {LC}}\cdot \left( 1-\frac{\sum _{\sum _{j\in \overline{J}}n_{j}=N}\prod _{j=1}^{J}\left( \prod _{i=1}^{n_{j}}\frac{\eta _{j}}{\nu _{j}(i)}\right) }{\sum _{\sum _{j\in \overline{J}_{0}}n_{j}=N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}\cdot \prod _{j=1}^{J}\left( \prod _{i=1}^{n_{j}}\frac{\eta _{j}}{\nu _{j}(i)}\right) }\right) \\&=\eta _{0}\cdot \frac{\sum _{n_{0}=1}^{N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}-1}\sum _{\sum _{j\in \overline{J}}n_{j}=N-n_{0}}\prod _{j=1}^{J}\left( \prod _{i=1}^{n_{j}}\frac{\eta _{j}}{\nu _{j}(i)}\right) }{\sum _{n_{0}=0}^{N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}\sum _{\sum _{j\in \overline{J}}n_{j}=N-n_{0}}\prod _{j=1}^{J}\left( \prod _{i=1}^{n_{j}}\frac{\eta _{j}}{\nu _{j}(i)}\right) }\\&=\eta _{0}\cdot \frac{\sum _{n_{0}=0}^{N-1}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}\sum _{\sum _{j\in \overline{J}}n_{j}=N-1-n_{0}}\prod _{j=1}^{J}\left( \prod _{i=1}^{n_{j}}\frac{\eta _{j}}{\nu _{j}(i)}\right) }{\sum _{n_{0}=0}^{N}\left( \frac{\eta _{0}}{\lambda _{\text {LC}}}\right) ^{n_{0}}\sum _{\sum _{j\in \overline{J}}n_{j}=N-n_{0}}\prod _{j=1}^{J}\left( \prod _{i=1}^{n_{j}}\frac{\eta _{j}}{\nu _{j}(i)}\right) }\\&=\eta _{0}\cdot \frac{C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)}{C_{\lambda _{\text {LC}}}(\overline{J}_{0},N)}. \end{aligned}$$

We further define for

$$\begin{aligned}\text {}~L=1,2,\ldots ,N:~~~ C(\overline{J},L):=\sum _{\sum _{j\in \overline{J}}n_{j}=L}\prod _{j=1}^{J}\left( \prod _{i=1}^{n_{j}}\frac{\eta _{j}}{\nu _{j}(i)}\right) \qquad \qquad \end{aligned}$$

and obtain

$$\begin{aligned} ~~~ C_{\lambda _{\text {LC}}}(\overline{J}_{0},L)=C(\overline{J},L)+\frac{\eta _{0}}{\lambda _{\text {LC}}}\cdot C_{\lambda _{\text {LC}}}(\overline{J}_{0},L-1),\quad L=1,2,\ldots ,N.\quad \end{aligned}$$

(29)

We prove by induction

$$\begin{aligned} \lambda _{\text {eff}}^{(N)}(\lambda _{\text {LC}}+\varepsilon )>\lambda _{\text {eff}}^{(N)}(\lambda _{\text {LC}})\qquad \forall \lambda _{\text {LC}}>0,\ \varepsilon >0,\ \forall N=1,2,\ldots \end{aligned}$$

Base step: For $N=1$:

$$\begin{aligned} \lambda _{\text {eff}}^{(1)}(\lambda _{\text {LC}}+\varepsilon )-\lambda _{\text {eff}}^{(1)}(\lambda _{\text {LC}})=\frac{\eta _{0}}{\frac{\eta _{0}}{\lambda _{\text {LC}}+\varepsilon }+\sum _{j=1}^{J}\frac{\eta _{j}}{\nu _{j}(1)}}-\frac{\eta _{0}}{\frac{\eta _{0}}{\lambda _{\text {LC}}}+\sum _{j=1}^{J}\frac{\eta _{j}}{\nu _{j}(1)}}>0. \end{aligned}$$

Induction step: Assume the inequality holds for $1\le L\le N-1$. Then for $L=N$ holds

$$\begin{aligned}&\frac{1}{\eta _{0}}\cdot \left( \lambda _{\text {eff}}^{(N)}(\lambda _{\text {LC}}+\varepsilon )-\lambda _{\text {eff}}^{(N)}(\lambda _{\text {LC}})\right) =\frac{C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-1)}{C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N)}-\frac{C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)}{C_{\lambda _{\text {LC}}}(\overline{J}_{0},N)}\\&\quad =\frac{C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-1)\cdot C_{\lambda _{\text {LC}}}(\overline{J}_{0},N)-C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)\cdot C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N)}{C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N)\cdot C_{\lambda _{\text {LC}}}(\overline{J}_{0},N)}. \end{aligned}$$

Because the denominator is strictly positive, it suffices to show that the numerator is strictly positive. Using the induction assumption, we obtain from Eq. (29):

$$\begin{aligned}&C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-1)\cdot C_{\lambda _{\text {LC}}}(\overline{J}_{0},N)-C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)\cdot C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N)\\&\quad {\mathop {=}\limits ^{}}\left( C(\overline{J},N-1)+\frac{\eta _{0}}{\lambda _{\text {LC}}+\varepsilon }\cdot C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-2)\right) \\&\qquad \cdot \left( C(\overline{J},N)+\frac{\eta _{0}}{\lambda _{\text {LC}}}\cdot C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)\right) \\&\qquad -\left( C(\overline{J},N-1)+\frac{\eta _{0}}{\lambda _{\text {LC}}}\cdot C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-2)\right) \\&\qquad \cdot \left( C(\overline{J},N)+\frac{\eta _{0}}{\lambda _{\text {LC}}+\varepsilon }\cdot C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-1)\right) \\&\quad =C(\overline{J},N-1)\cdot C(\overline{J},N)+C(\overline{J},N-1)\cdot \frac{\eta _{0}}{\lambda _{\text {LC}}}\cdot C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)\\&\qquad +\frac{\eta _{0}}{\lambda _{\text {LC}}+\varepsilon }\cdot C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-2)\cdot C(\overline{J},N)\\&\qquad +\frac{\eta _{0}}{\lambda _{\text {LC}}+\varepsilon }\cdot C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-2)\cdot \frac{\eta _{0}}{\lambda _{\text {LC}}}\cdot C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)\\&\qquad -C(\overline{J},N-1)\cdot C(\overline{J},N)-C(\overline{J},N-1)\cdot \frac{\eta _{0}}{\lambda _{\text {LC}}+\varepsilon }\cdot C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-1)\\&\qquad -\frac{\eta _{0}}{\lambda _{\text {LC}}}\cdot C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-2)\cdot C(\overline{J},N)\\&\qquad -\frac{\eta _{0}}{\lambda _{\text {LC}}}\cdot C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-2)\cdot \frac{\eta _{0}}{\lambda _{\text {LC}}+\varepsilon }\cdot C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-1)\\&\quad =\frac{\eta _{0}}{\lambda _{\text {LC}}+\varepsilon }\cdot \frac{\eta _{0}}{\lambda _{\text {LC}}}\cdot \Big [C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-2)\cdot C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)\\&\qquad \underbrace{\qquad \qquad \qquad \qquad -C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-2)\cdot C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-1)\Big ]}_{=:A_{N}}\\&\qquad +\frac{\eta _{0}}{\lambda _{\text {LC}}}\cdot C(\overline{J},N-1)\cdot C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)+\frac{\eta _{0}}{\lambda _{\text {LC}}+\varepsilon }\cdot C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-2)\\&\qquad \cdot C(\overline{J},N)\\&\qquad -\frac{\eta _{0}}{\lambda _{\text {LC}}+\varepsilon }\cdot C(\overline{J},N-1)\cdot C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-1)-\frac{\eta _{0}}{\lambda _{\text {LC}}}\cdot C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-2)\\&\qquad \cdot C(\overline{J},N)\\&\quad =A_{N}+\frac{\eta _{0}}{\lambda _{\text {LC}}}\cdot \left[ C(\overline{J},N-1)\cdot C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)-C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-2)\cdot C(\overline{J},N)\right] \\&\qquad +\frac{\eta _{0}}{\lambda _{\text {LC}}+\varepsilon }\cdot \left[ C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-2)\cdot C(\overline{J},N)-C(\overline{J},N-1)\right. \\&\qquad \left. \cdot C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-1)\right] . \end{aligned}$$

$A_{N}$ is positive: $\frac{\eta _{0}}{\lambda _{\text {LC}}+\varepsilon }\cdot \frac{\eta _{0}}{\lambda _{\text {LC}}}$ multiplied by the numerator of $\lambda _{\text {eff}}^{(N-1)}(\lambda _{\text {LC}}+\varepsilon )-\lambda _{\text {eff}}^{(N-1)}(\lambda _{\text {LC}})$. Therefore, $A_{N}>0$ by induction assumption. Thus, it suffices to prove

$$\begin{aligned}&\frac{\eta _{0}}{\lambda _{\text {LC}}}\cdot \left[ C(\overline{J},N-1)\cdot C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)-C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-2)\cdot C(\overline{J},N)\right] \\&\quad +\frac{\eta _{0}}{\lambda _{\text {LC}}+\varepsilon }\cdot \left[ C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-2)\cdot C(\overline{J},N)-C(\overline{J},N-1)\right. \\&\quad \left. \cdot C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-1)\right] \ge 0, \end{aligned}$$

which is equivalent to

$$\begin{aligned}&\frac{\eta _{0}}{\lambda _{\text {LC}}}\cdot \left[ \frac{C(\overline{J},N-1)}{C(\overline{J},N)}-\frac{C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-2)}{C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)}\right] \cdot C(\overline{J},N)\cdot C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)\\&\qquad -\frac{\eta _{0}}{\lambda _{\text {LC}}+\varepsilon }\cdot \left[ \frac{C(\overline{J},N-1)}{C(\overline{J},N)}-\frac{C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-2)}{C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-1)}\right] \cdot C(\overline{J},N)\\&\qquad \cdot C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-1)\\&\quad =\Bigg \{\underbrace{\left[ \frac{C(\overline{J},N-1)}{C(\overline{J},N)}-\frac{C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-2)}{C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)}\right] }_{=:B}\cdot \underbrace{\frac{\eta _{0}}{\lambda _{\text {LC}}}\cdot C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)}_{=:D} \\&\qquad -\underbrace{\left[ \frac{C(\overline{J},N-1)}{C(\overline{J},N)}-\frac{C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-2)}{C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-1)}\right] }_{=:C}\cdot \underbrace{\frac{\eta _{0}}{\lambda _{\text {LC}}+\varepsilon }\cdot C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-1)}_{=:E} \Bigg \}\\&\qquad \cdot C(\overline{J},N) \ge 0. \end{aligned}$$

(i) We show that $B\ge 0$, $C\ge 0$. Because increasing population increases throughput (Van der Wal 1989)

$$\begin{aligned}&\frac{C(\overline{J},N-1)}{C(\overline{J},N)}\ge \frac{C(\overline{J},N-2)}{C(\overline{J},N-1)}.\end{aligned}$$

Next, we show:

$$\frac{C(\overline{J},N-2)}{C(\overline{J},N-1)}\ge \frac{C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-2)}{C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)}. $$

We consider the right-hand side as throughput of a cyclic Gordon–Newell network with node set $\overline{J}_{0}:=\left\{ 0,1,\ldots ,J\right\}$, service rates $\mu _{j}(n):=\frac{\nu _{j}(i)}{\eta _{j}}$, $j=1,\ldots ,J$, $n=0,1,\ldots ,N$ and $\mu _{0}(n):=\frac{\lambda _{\text {LC}}}{\eta _{0}}$ and solution $\underbrace{(1,1,\ldots ,1)}_{J+1\text {-times}}$ of the routing matrix for the cycle. $\frac{C(\overline{J},N-2)}{C(\overline{J},N-1)}$ and $\frac{C_{\lambda _{\text {LC}}}(\overline{J},N-2)}{C_{\lambda _{\text {LC}}}(\overline{J},N-1)}$ as given in our derivations are the (average) throughput of a cycle following (Daduna et al. 2008, Definition 2.6).

The left-hand side is the throughput of a cyclic Gordon–Newell network which is obtained from the first cycle by deleting node 0 and skipping the gap by cycling customers.

Lemma 2.8 in Daduna et al. (2008) states that deleting a node of a cycle with skipping the gap increases the throughput. Consequently, $B\ge 0$ in a two-step conclusion, and similarly $C\ge 0$.

(ii) From the definitions of D and E it follows by direct comparison $D>E$.

(iii) The proof will be finished if we can show $B\ge C$. It is sufficient to prove with Shanthikumar and Yao (1986)

$$\begin{aligned} \frac{C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-2)}{C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)}\le \frac{C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-2)}{C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-1)}. \end{aligned}$$

(30)

The left-hand side is the throughput of the Gordon–Newell network according to Definition 2.2 in Shanthikumar and Yao (1986). The right-hand side is the throughput of the Gordon–Newell network with service rate at node 0 increased to $\lambda _{\text {LC}}+\varepsilon$. Corollary 3.1(i) in Shanthikumar and Yao (1986) states

$$\begin{aligned}\qquad \qquad \eta _{0}\cdot \frac{C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-2)}{C_{\lambda _{\text {LC}}}(\overline{J}_{0},N-1)}\le \eta _{0}\cdot \frac{C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-2)}{C_{\lambda _{\text {LC}}+\varepsilon }(\overline{J}_{0},N-1)},\end{aligned}$$

which verifies (30). $\square$

Remark 3

For $\lambda _{\text {LC}}$ with $\lambda _{\text {eff}}(\lambda _{\text {LC}})=\lambda _{\text {BO}}\in \left( 0,\lambda _{\text {BO,max}}\right) =\left( 0,\eta _{0}\cdot \tfrac{b(1)}{b(0)}\right)$ it holds

$$\begin{aligned} \lambda _{\text {LC}}={\left\{ \begin{array}{ll} \frac{\eta _{0}\cdot \lambda _{\text {BO}}}{\eta _{0}-\lambda _{\text {BO}}\cdot b(0)}, &{} N=1,\\ -\frac{\eta _{0}}{2\cdot \big (\eta _{0}\cdot b(1)-\lambda _{\text {BO}}\cdot b(0)\big )}\\ \cdot \left( \eta _{0}-\lambda _{\text {BO}}\cdot b(1)-\sqrt{\big (\eta _{0}+\lambda _{\text {BO}}\cdot b(1)\big ){}^{2}-4\cdot \lambda _{\text {BO}}^{2}\cdot b(0)}\right) , &{} N=2. \end{array}\right. } \end{aligned}$$

Proof

Due to Theorem 2 we have for $\lambda _{\text {LC}}\in (0,\infty )$

$$\begin{aligned}&\lambda _{\text {eff}}(\lambda _{\text {LC}})=\lambda _{\text {LC}}\cdot \left( 1-\frac{b(0)}{\sum _{n=0}^{N}\left( \frac{\lambda _{\text {LC}}}{\eta _{0}}\right) ^{n}\cdot b(n)}\right) \\&\quad =\lambda _{\text {LC}}\cdot \left( \frac{\sum _{n=1}^{N}b(n)\cdot \lambda _{\text {LC}}\cdot \left( \frac{\lambda _{\text {LC}}}{\eta _{0}}\right) ^{N-n}}{\sum _{n=0}^{N}b(n)\cdot \left( \frac{\lambda _{\text {LC}}}{\eta _{0}}\right) ^{N-n}}\right) \end{aligned}$$

and

$$\begin{aligned}& b(N)=1,\quad b(N-1)=\sum _{j=1}^{J}\frac{\eta _{j}}{\nu _{j}(1)}, \\& b(N-2)=\sum _{j=1}^{J}\frac{\eta _{j}}{\nu _{j}(1)\cdot \nu _{j}(2)} \quad +\sum _{j=1}^{J-1}\sum _{k=j+1}^{J}\frac{\eta _{j}\cdot \eta _{k}}{\nu _{j}(1)\cdot \nu _{k}(1)}.\end{aligned}$$

Let $\lambda _{\text {BO}}\in \left( 0,\eta _{0}\cdot \tfrac{b(1)}{b(0)}\right)$. First, we note that

$$\begin{aligned}&\eta _{0}\cdot b(1)-\lambda _{\text {BO}}\cdot b(0)>\eta _{0}\cdot b(1)-\eta _{0}\cdot \frac{b(1)}{b(0)}\cdot b(0)=0. \end{aligned}$$

(31)

$\lambda _{\text {eff}}(\lambda _{\text {LC}})=\lambda _{\text {BO}}$ is equivalent to

$$\begin{aligned}\frac{\sum _{n=1}^{N}\left( \frac{\lambda _{\text {LC}}}{\eta _{0}}\right) ^{N+1-n}\cdot \eta _{0}\cdot b(n)}{\sum _{n=0}^{N}\left( \frac{\lambda _{\text {LC}}}{\eta _{0}}\right) ^{N-n}\cdot b(n)}=\lambda _{\text {BO}},\end{aligned}$$

which is equivalent to:

$$\sum _{n=1}^{N}\big (\eta _{0}\cdot b(n)-\lambda _{\text {BO}}\cdot b(n-1)\big )\cdot \left( \frac{\lambda _{\text {LC}}}{\eta _{0}}\right) ^{N+1-n}-\lambda _{\text {BO}}\cdot b(N)=0.$$

(32)

If $N=1$, then Eq. (32) is equivalent to

$$\begin{aligned} \lambda _{\text {LC}}=\frac{\eta _{0}\cdot \lambda _{\text {BO}}\cdot b(1)}{\eta _{0}\cdot b(1)-\lambda _{\text {BO}}\cdot b(0)}=\frac{\eta _{0}\cdot \lambda _{\text {BO}}}{\eta _{0}-\lambda _{\text {BO}}\cdot b(0)}. \end{aligned}$$

If $N=2$, then Eq. (32) is equivalent to

$$\begin{aligned}&\big (\eta _{0}\cdot b(1)-\lambda _{\text {BO}}\cdot b(0)\big )\cdot \left( \frac{\lambda _{\text {LC}}}{\eta _{0}}\right) ^{2}+\big (\eta _{0}\cdot b(2)-\lambda _{\text {BO}}\cdot b(1)\big )\cdot \left( \frac{\lambda _{\text {LC}}}{\eta _{0}}\right) -b(2)\cdot \lambda _{\text {BO}}\\&\quad =0. \end{aligned}$$

Since the discriminant fulfils

$$\begin{aligned} \big (\eta _{0}\cdot b(2)-\lambda _{\text {BO}}\cdot b(1)\big )^{2}+4\cdot \big (\eta _{0}\cdot b(1)-\lambda _{\text {BO}}\cdot b(0)\big )\cdot b(2)\cdot \lambda _{\text {BO}}>0 \end{aligned}$$

by Eq. (31), it follows that

$$\begin{aligned} \frac{\lambda _{\text {LC}}}{\eta _{0}}&=-\frac{1}{2\cdot \big (\eta _{0}\cdot b(1)-\lambda _{\text {BO}}\cdot b(0)\big )} \cdot \Bigl (\eta _{0}\cdot b(2)-\lambda _{\text {BO}}\cdot b(1)\\&\quad \pm \sqrt{\big (\eta _{0}\cdot b(2)-\lambda _{\text {BO}}\cdot b(1)\big )^{2}+4\cdot \big (\eta _{0}\cdot b(1)-\lambda _{\text {BO}}\cdot b(0)\big )\cdot b(2)\cdot \lambda _{\text {BO}}}\Bigr ). \end{aligned}$$

Due to Eq. (31), the solution $\lambda _{\text {LC}}$ is positive only if

$$\begin{aligned}&\eta _{0}\cdot b(2)-\lambda _{\text {BO}}\cdot b(1)\\&\quad \pm \sqrt{\big (\eta _{0}\cdot b(2)-\lambda _{\text {BO}}\cdot b(1)\big )^{2}+4\cdot \big (\eta _{0}\cdot b(1)-\lambda _{\text {BO}}\cdot b(0)\big )\cdot b(2)\cdot \lambda _{\text {BO}}}<0. \end{aligned}$$

Let $c:=\eta _{0}\cdot b(2)-\lambda _{\text {BO}}\cdot b(1)$ and $d:=4\cdot \big (\eta _{0}\cdot b(1)-\lambda _{\text {BO}}\cdot b(0)\big )\cdot b(2)\cdot \lambda _{\text {BO}}$, which implies $d>0$. If $c\ge 0$, then $c+\sqrt{c^{2}+d}>0$. If $c<0$, then $c+\sqrt{c^{2}+d}>c+\sqrt{c^{2}}=c+|c|=c-c=0.$ Hence, the only positive solution (existence guaranteed by Theorem 2) is

$$\begin{aligned}&\lambda _{\text {LC}}=-\frac{\eta _{0}}{2\cdot \big (\eta _{0}\cdot b(1)-\lambda _{\text {BO}}\cdot b(0)\big )} \cdot \Bigl (\eta _{0}\cdot b(2)-\lambda _{\text {BO}}\cdot b(1)\\&\qquad -\sqrt{\big (\eta _{0}\cdot b(2)-\lambda _{\text {BO}}\cdot b(1)\big )^{2}+4\cdot \big (\eta _{0}\cdot b(1)-\lambda _{\text {BO}}\cdot b(0)\big )\cdot b(2)\cdot \lambda _{\text {BO}}}\Bigr )\\&\quad =\frac{-\eta _{0}}{2\cdot \big (\eta _{0}\cdot b(1)-\lambda _{\text {BO}}\cdot b(0)\big )}\\&\qquad \cdot \left( \eta _{0}-\lambda _{\text {BO}}\cdot b(1)-\sqrt{\big (\eta _{0}+\lambda _{\text {BO}}\cdot b(1)\big )^{2}-4\cdot \lambda _{\text {BO}}^{2}\cdot b(0)}\right) . \end{aligned}$$

$\square$

Transformation of passage time through the inner network. We consider the SOQN-BO from Sect. 2 and its approximation by an SOQN-LC in Sect. 3. Within the SOQN-LC the resource network behaves stochastically as a Gordon–Newell network with $N\ge 1$ customers. We take the inner network as a complex black-box server, which is able to serve several customers in parallel. It is well known that the distribution of customers’ passage time through that black box is in general not known. Only mean values are accessible, although complicated as well. At a first glance one expects to have passage times, say Y, with high variability, measured e.g. by the squared coefficient of variation $C^2(Y)=Var(Y)/E^2(Y)$. (Note, even Var(Y) is in general not known.) The problem is to evaluate the effect of the transformation which reduces the inner network to a single substitute station with state dependent service rates exploiting Norton’s Theorem. We have not been able to resolve the problem of transforming the variability in general, but we can show by examples why there will be no general answer to that question.

We consider three cases of SOQN-BO with Poisson-$\lambda$ arrival stream and only one resource. Then at most one customer can be served by the inner network, and the passage time through the inner network has a Phase-type distribution (PH). So the system behaves like a single server $M/PH/1/\infty$ queue and for simple PH-distributions we can compute e.g. the mean waiting time at the external queue explicitly. We fix in advance a parameter $1/\mu >0$.

(i) The inner system consists of only one single server station with service rate $\mu$. Then the reduced system is trivially the same and the coefficient of variation of the passage time through the inner network and the substitute is 1.

(ii) The inner system consists of $k>1$ single server stations in line (tandem queue) each with service rate $k\cdot \mu$. Then the passage time through the inner network is a k-stage Erlang distribution with mean $1/\mu$ and coefficient of variation $1/k <1$. The reduced system according to Norton’s Theorem is a single server station with service rate $\mu$ and coefficient of variation 1.

(iii) Let $C^2>1$. The inner system consists of two parallel single server stations with service rates $\mu _{1}$ and $\mu _{2}$. A customer traversing the inner network chooses station 1 with probability $\alpha \in (0,1)$, station 2 with probability $1-\alpha$. We take $\mu _1:=2\alpha \mu$ and $\mu _2:=2(1-\alpha ) \mu$, and $\alpha := (1/2)\cdot [1-\sqrt{(C^2-1)/(C^2+1)}]$. This yields a hyperexponential distribution with mean $1/\mu$ and coefficient of variation $C^2 >1$. The reduced system according to Norton’s Theorem is a single server station with service rate $\mu$ and coefficient of variation 1.

The conclusion is that the variability of passage time through the inner network under the transformation can be (i) maintained, (ii) decreased, (iii) increased.