On queueing-inventory-location problems

A network of queueing-inventory systems is considered where the inventories are replenished by a common central server. Travel times for transport of the items send out from the center and an adaptive dispatching regime to direct transportation towards the locations with queueing-inventory systems are incorporated in the model. The stationary distribution of the system is of product form and is explicitly given. An optimization procedure is developed to find for given locations of the queueing-inventory systems an optimal location for the replenishment center. Optimization is with respect to overall utilization of the resources measured in total throughput of the queueing-inventory systems.


Introduction
In standard OR-literature queueing theory and inventory theory usually are dealt with as different fields of research, see e.g. the volume "Stochastic Models" in the series "Handbooks in Operations Research and Management Science" Heyman and Sobel (1990) with separated sections on "Queueing Theory" Cooper (1990) and "Stochastic Inventory Theory" Porteus (1990).As the editors put it "Queueing theory is more descriptive than normative, inventory is the reverse,…" Heyman and Sobel (1990)[p. ix].
On the other side, planning and analysing complex systems in connection with nowadays large supply chains calls for models which encompass integration of production processes and inventory holding.Increasing complexity of supply chains requires integration of other areas of OR, e.g.reliability, logistics (transportation), and location analysis, see the discussion in Heckmann and Nickel (2019).
Realizing the need for models which integrate congestion phenomena (queueing) and inventory control, over the last 30 years queueing-inventory models of various structures and different level of complexity have been developed, see the review of selected literature below.
But observing the complexity of supply chains or production plants of today, it is obvious that integration of queueing and inventory aspects only will not meet all the features inherent in these systems.
The aim of this article is to integrate into queueing-inventory systems spatial aspects of distributed queueing-inventory systems.Starting this project it is immediately visible that aspects of transport have to be integrated as well.Integration of location decisions, logistics activities, production, and inventory management generates problems with decision making of different levels: Strategic, tactical, and operational decisions are mixed up.Heckmann and Nickel emphasize that "…making location decisions ignoring primary logistics activities like production or distribution may result in excessive costs incurred throughout the lifetime of the facilities supporting the logistics system."Heckmann and Nickel (2019)[p. 456].
In modelling supply chains the areas of queues, inventories, location analysis constitute important sub-problems which from the very definition of a supply chain are strongly interconnected.Nevertheless, due to the system's overall complexity often the related parts of production, storage, transportation, and decision for location of facilities are separated to obtain feasible (optimization-) problems where classical solution procedures from the respective fields could be applied to solve the isolated sub-problems.This separation is often justified by considering the location, production, inventory, and transport control problems as being problems that occur on different levels of decision making: Strategic, tactical, and operational scales of management.
On the other side, for several of the mentioned composite problems it has been shown recently that separation of the partial problems leads to sub-optimal solutions of the relevant global optimization problems, see for separated location-transportation control problems Salhi and Rand (1989).
Our aim is to construct a network model which integrates aspects of production with associated inventories, transport, and location of a common replenishment system (for to satisfy reorders of inventories) under random influences.Incorporating randomness excludes to consider many details which are part of a deterministic modeling process as it is presented in Cordeau et al. (2006) and Song and Wu (2022).
Main contributions of the paper.For the sketched network model • we compute the steady state distribution in explicit product form, i.e. the state process is separable, • we find an optimal location for the central replenishment server (supporting strategic decisions), • we determine the needed inventory to satisfy on a make-to-order basis the arriving demand (supporting tactical decisions), • we construct an adaptive dispatching regime for to send out items from the replenishment station (supporting operational decisions).
Structure of the paper.In Sect. 2 the general problem and its components are described in more detail.In Sect. 3 a short review of related literature is given.The formal model with all details is presented in Sect. 4. Section 4.1 is devoted to the analysis of the system in steady state with short sections on cost analysis in Sect.4.2 and generalizations in Sect.4.3.The location problem for the central unit is described in Sect. 5. A reduced model is developed in Sect.5.1 which enables to determine the optimal location of the center and the needed inventory under optimal location in Sect.5.2.Conclusions are summarized in Sect.6.In Appendix 1 technical necessities are reviewed.Notations and conventions.(x) + := max(x, 0) for x ∈ R. Empty sums are 0, empty products are 1.N 0 = {0, 1, 2, . . .}, N + = {1, 2, . . .}. J = {1, . . ., J }, and J + = {1, . . ., J , J + 1} for J ∈ N + .
Increasing≡non-decreasing, decreasing≡non-increasing. All random variables and processes which occur are defined on a common probability space ( , F , P). Expectations under P are expressed as E(•).
For a probability measures π on a discrete space we write shortly E π (•) for expectations under π.

Problem statement
We investigate a system which consists of a set of production facilities, each with an associated inventory (≡ queueing-inventory system).The positions of the locations with productioninventory system are known and fixed in the plane.We seek a (nearly) optimal position for a common replenishment facility according to a suitable optimization criterion which incorporates the effects of transportation, varying demand, and holding inventories.
To be more specific: We are given J ≥ 1 queueing-inventory systems, numbered 1, . . ., J at positions a j = (a j1 , a j2 ) in the plane.At each location the production unit is represented by a single server with ample waiting room under First-Come-First-Serve (FCFS) regime.External demand (represented by customers) at location j follows a Poisson-λ j stream.Production is according to customers' local demand on a make-to-order basis and to satisfy a customer's demand the production server needs one (single) item of raw material from the associated inventory.For the inventory a base-stock level is prescribed which is location specific.
If a customer's service is completed the customer departs from the system and the consumed item from the inventory is removed formally from the inventory at this instant.Concurrently an order for a new item of raw material for this location is placed which arrives at the replenishment center without time delay.So replenishment orders from the inventories follow local continuous review base-stock policies.
The items to be stored in the inventories are manufactured by a central replenishment unit which is represented by a single server with ample waiting room under FCFS.Whenever an item of raw material is produced it is send out to a location according to an adaptive dispatching regime which will be specified below.We assume that items are exchangeable.This implies that we can decide about the target location when the item is produced.
Customers' behaviour on arrival at the locations can be described as following a combined backordering/lost-sales principle: Customers who find on arrival an empty queueing system and an item on stock in the local inventory enter service immediately.If on arrival there is already a customer in service and the inventory is not empty these new customers enter the system and queue up (backordering).If the inventory is depleted every new arrival is rejected by the system (lost sales).
The local production-inventory systems with central supplier constitute a standard (integrated) queueing-inventory problem as investigated e.g. in Otten et al. (2016).We add to this model a location-transportation problem by integrating the transport of the manufactured items from the replenishment center to the locations and a decision for locating the center.
Transportation is by trucks which transport single items, starting when an item is send out.For modelling this feature we use infinite server systems which are standard devices of queueing theory to model transportation systems.We assume that there is ample transportation capacity, i.e. whenever an item is manufactured and ready to be send out a truck is available.Such an integrated system is termed "logistics network" in Cordeau et al. (2006)[p. 60]: "A logistics network is a set of suppliers, manufacturing plants and warehouses organized to manage the procurement of raw materials, their transformation into finished products, and the distribution of finished products to customers."In Cordeau et al. (2006) a deterministic one-period optimization problem is formulated and two approaches for solving the problem are described.Our approach is completely different.The main goals of our investigation are • to find specifications of the system which allow to determine explicitly the asymptotic and stationary behaviour of the system which is driven by stochastic processes, especially random arrival processes of demand, stochastic service and production times, random travel times for delivering the replenishment items, and • to find an optimal location for the central replenishment center which supplies the predetermined locations with items to be hold on stock for the ongoing production processes, and • to include an adaptive routing regime that decides automatically which orders for replenishment by the locations should be satisfied next, and • to include (different from e.g.Cordeau et al. (2006)) the time dimension in our optimization procedure in the sense of Tapiero (1971).We evaluate the time average of local total production (≡ satisfied external demand) over a long time horizon which due to stabilization of the system is computed as steady state throughput.
Summarizing, our model integrates strategic (facility location) and tactical/operational (allocation, scheduling, service, stock holding) aspects of decision making.These integrated logistics-location-service models occur e.g. in supply chain planning and operation, see e.g.Melo et al. (2009) and Heckmann and Nickel (2019).Clearly, there arise many difficulties with such integration procedures in combined location-routing-service problems.A typical example is discussed in Min et al. (1998)[p.10].

Literature review
The interplay of location analysis and decisions with aspects and problems from other areas of OR has been investigated and described in many research articles.We restrict our description to integrated models under random influences, i.e. the abundant literature around (one-period) deterministic optimization models will not be considered here.
Queueing-location problems, i.e., investigation and optimization of combined queueing and location models is a well established part of Operations Research.Larson (1974) initiated research on location problems within the scope of queueing systems, follow up papers in the context of discrete location problems are e.g. by Larson, Berman andcoauthors (e.g. Berman et al. (1985, 1987)).Relevant survey chapters in collections are Drezner and Hamacher (2004)[Chapter 11], Mirchandani and Francis (1990)[Chapter 13], (vehicle routing problems under stochastic side constraints).A recent very detailed survey is Berman and Krass (2019).Two main research directions are: (1) Clients are fixed, servers move.In Drezner et al. (1990) and Scott et al. (1999) the authors consider how to locate mobile servers in a plane to serve demands which occur as Poisson processes at fixed locations.The mobile servers are modelled as standard queueing systems.Travel times of the mobile servers to and from the clients are incorporated in the service times, a survey is Berman and Krass (2004).An investigation of determining locations in case of emerging queueing phenomena is Dan and Marcotte (2019).
(2) Servers are fixed, clients move.See e.g Berman and Drezner (2007), Aboolian et al. (2008Aboolian et al. ( , 2009) ) with additional references.Demand is generated locally according to renewal processes, the service is provided by M/M/k/∞ systems.A detailed review is Berman and Krass (2019).
Location-inventory problems (LIPs) are often summarized under the heading "location analysis", which means to find optimal locations for warehouses or central suppliers.LIPs aim "to integrate strategic supply chain decisions with tactical and operational inventory management decisions", see Farahani et al. (2015) with a survey of research on basic LIPs and more evolved variants.A basic LIP is described in Farahani et al. (2015)[Section2].
Location routing problems (LRPs) combine location analysis and vehicle routing decision.First principles of vehicle routing are described in Laporte (1988).A recent survey mainly with locations on networks is Albareda-Sambola and Rodriguez-Pereira (2019).Closer to location in the plane are Salhi and Nagy (2009) and Manzour-al-Ajdad et al. (2012).
Location-inventory-routing problems (LIRPs) as defined in Song and Wu (2022) constitute an integrated approach to model parts of a supply chain, including location of distribution centers, inventory holding, allocation decisions, and routing schedules for transportation.Predecessors of that investigation are summarized in Song and Wu (2022) [p. 3].
Integration of strategic and tactical/operational aspects of planning is a common topic of almost all of the mentioned work.In Salhi and Rand (1989) for LRPs it is shown that separating decisions on location and routing can lead to sub-optimal decisions.Under the heading "Why logistics matters in location modelling" this problem is discussed indepth in Heckmann and Nickel (2019) [Section 6.2], stating as "main conclusion …that making location decisions ignoring primary logistics activities …may result in excessive costs …." Literature on queueing-inventory problems is overwhelming as can be seen from the recent review Krishnamoorthy et al. (2021).Earlier reviews with additional comments on specific aspects of queueing-inventory systems are Krishnamoorthy et al. (2011) and the "Short review of known results" in Melikov et al. (2016).This points out that the integrated model which was introduced to the research community only 30 years ago independently of another in Melikov and Molchanov (1992) and Sigman and Simchi-Levi (1992) is of interest to theoretical research and practical applications.Instead of repeating the review in Krishnamoorthy et al. (2021) I refer only to fundamental methodological approaches here.
A first branch exploits approximation procedures, Melikov and Molchanov (1992) and Sigman and Simchi-Levi (1992) are typical examples.Realizing that the global balance equations for the dedicated composite models do not have an easy to find explicit solution the authors turn to suitable procedures to obtain approximate solutions for the steady state distribution.This technique reoccurs in many other articles, e.g. by Melikov, Molchanov, Ponomarenko, Koroliuk, Bagirova, et al. as demonstrated in Krishnamoorthy et al. (2021).
A second branch of methods to attack queueing-inventory systems are Markov decision processes and stochastic dynamic optimization which are standard methods to obtain optimal policies in classical inventory theory, see Puterman (1990).Especially, Berman and coworkers applied these techniques to obtain structural information on optimal policies to control queueing-inventory systems, see e.g.Berman and Kim (1999), Berman and Sapna (2001).
A third branch utilizes the state space structure of single queues with inventory as being a product N 0 × F with queue length in N 0 and with F counting for the inventory size (possibly plus some external environment).Any Markov process which fits into this class can be described as Quasi-birth-death process (QBD) and the steady state distribution (if it exists) can be obtained numerically with matrix-geometrical methods developed and described for this class of processes in Neuts (1981), Latouche and Ramaswami (1999).This applies especially if the arrival process of demands is not Poisson.Seemingly, most of the results published on queueing-inventory systems are obtained via matrix-geometrical methods, see Krishnamoorthy et al. (2021), to name a few authors we refer to the work of the group around Krishnamoorthy, Chakravarthy, Shajin, Narayanan, Lakshmy, et al.-for details see Krishnamoorthy et al. (2021).
A fourth branch of articles exploits internal structures of the describing Markov processes for the queueing-inventory systems.These processes exhibit in the internal structures of the state processes local or partial balance properties leading to surprising simple stationary distributions which are called product form equilibrium, which means that asymptotically and in equilibrium the joint distribution of the queueing and the inventory component decouple (≡ being stochastically independent).According to Krishnamoorthy et al. (2021) the number of articles where this characterization applies is minor compared to branches one and three.Contributions to this field started with Schwarz et al. (2006Schwarz et al. ( , 2007)), and are continued among others with Krishnamoorthy and Narayanan (2013), Shajin et al. (2018) and other members of the group around Krishnamoorthy, and with Saffari et al. (2013) and some further articles by Saffari, Haji, and their coworkers.
Remarkably, almost all work is on single servers with an attached inventory.The investigation of networks of queues with attached inventory is just in the beginning.The research presented in this article is related to investigation of network models in Otten et al. (2016Otten et al. ( , 2020)), Otten (2017)[Section 3.4].

The queueing-inventory-location model
We have a set of locations J := {1, . . ., J } located in the plain at a j = (a j1 , a j2 ), each equipped with an exponential single server as production server with infinite waiting room under First-Come-First-Served (FCFS) and an attached inventory of maximal size b j ≥ 1, j = 1, . . ., J .The service rates μ j (•) of the production servers depend on the (local) queue length.If at location j are n j > 0 customers present (in service or waiting) then μ j (n j ) > 0. If necessary we set μ j (0) = 0. We assume throughout that the μ j (n j ) are non-decreasing in n j ∈ N 0 .Recall that service of a customer (production of a demanded unit) at j needs one item from inventory which is consumed formally at the instant when service of this customer is finished.This customer leaves the system immediately and at the same time instant an order for manufacturing a new item arrives at the replenishment station (base stock policy).
The replenishment station (identified by number J + 1) is located at x = (x 1 , x 2 ) and is an exponential single server with infinite waiting room under FCFS with service rate ν > 0.
An item produced at J + 1 is dispatched to one of the stations, say j, which are not satisfied with items in the local inventory (k j items on stock at j) or on the way to j (m j items on transport to j) with m j + k j < b j .This restriction is due to the maximal capacity b j (= base-stock level) at j.We incorporate in the dispatching rules the condition that the total stock (≡ inventory position or system stock in standard inventory systems with backordering Porteus (1990)[p. 605]) at location j includes the items already on transport to j. (We do not include items which are ordered but not produced and on transport to the respective location because this would in case of base-stock policy under lost sales always sum up to b j .)Given an item is on leave from station J +1, station j is selected with branching probability r J +1, j (z) (to be defined below) if the global state of the queueing-inventory-transportation system is z ∈ E (≡ state space).
Transportation of an item by a truck from center J + 1 at x to location j at a j is modelled by serving the truck (as a customer) at an exponential infinite server (ample service capacity).For simplicity of presentation we assume that trucks travel with unit speed, i.e. if the distance between J + 1 and j is d j (x) := d(x, a j ), the mean travel (service) time is d j (x).Here d(•) is any suitable distance function (metric).So, we assume that travel times are exponentially distributed with rate d j (x) −1 .Modeling traffic on a lane by infinite server queueing system is a standard device from queueing theory in transportation theory, see e.g.Newell (1982)[Section 6].(In the algorithm to determine necessary inventory capacity to fulfill predetermined demand realistic speeds will be incorporated, see Algorithm 5.13.) We construct a multivariate Markov process Z = (Z (t) : t ≥ 0) which describes the evolution over time of this system.The state space of Z is A typical local state of location j at time t ≥ 0 is The meaning of (4.2) is: V j (t) = m j items are on transport with destination j, the size of the onhand inventory at j is Y j (t) = k j , and the queue length is X j (t) = n j at the production server at j.The queue length at the central replenishment server J +1 is implicitly determined as J j=1 (b j − m j − k j ).Imposing the usual (conditional) independence assumptions on the arrival streams, service and transport times, and routing (dispatching) decisions, it follows by standard arguments that is a strong Markov process with countable state space E from (4.1).Summarizing: The model encompasses a set of J exponential single server queues representing production which are coupled with a closed network representing support of the production.There are b = b 1 + • • • + b J customers cycling in the closed network according to problem specific rules.These customers represent throughout one cycle different objects (have alternating identities): Orders from an inventory waiting at the replenishment station, trucks carrying new items, and items stored in an inventory.
State dependent routing with blocking (here due to V j (t) + Y j (t) ≤ b j ) poses specific difficulties when separability is an aim of the modelling process.The following definition follows ideas of Towsley (1980) who constructed queueing networks with state dependent routing protocols which exhibit detailed balance and product form steady state in the sense of Gordon-Newell and Jackson networks.

Remark 4.2
The notation h j (•) is used for easier reading.In fact these functions are functions of the destination branch j and the base stock level b j : The extended notion will be used if necessary.Towsley (1980)[Theorem 3] showed that the above functions can be taken as In Daduna (1985) it was shown that in a simple star-like network C ∈ R may be allowed, which introduced a blocking scheme as needed for the more general network in the present investigation.

Definition 4.3
The infinitesimal generator (intensity matrix) of Z is denoted by Q = (q(z, z ) : z, z ∈ E) with the following positive rates.(For readability we only depict the relevant local states for the respective transitions, e.g.we write (. . .m j , k j , n j ; . . . ) instead of (m 1 , k 1 , n 1 ; . . .m j , k j , n j ; . . .m J , k J , n J ).)For j = 1, . . ., J we have q((. . .m j , k j , n j ; . . .), (. . . The diagonal elements q(z, z) of Q are selected in a way that row sums are zero.All other rates q(z, z ), z = z , are zero.

Steady-state analysis
Performance evaluation for the queueing-inventory-location model is either based on long time evaluation with averaging over time or based on observing the system in steady state (if the system is stable).For a system described by a Markov process these methods yield the same performance metrics if the process is ergodic.

Theorem 4.4 The Markov process Z is ergodic if and only if for all j
In a first step we show that a non-normalized solution of this system of global equilibrium equations is For j ∈ {1, . . ., J } we obtain in (4.11) (inserting the definition of the routing probabilities) which after simple manipulation is (b).
which coincides with (a) in this case.This finishes the first part of the proof.Because irreducibility of Z follows directly from the definition of E and routing R in (4.4) summability of (x(z) : z ∈ E) guarantees ergodicity of Z .Direct computation separates terms associated to queue lengths, i.e.J j=1 n j i=1 λ j μ j (i) in (4.9).Because H (b 1 , . . .b J ; J ) is a finite sum this finishes the second part of the proof.
Consider the system from Theorem 4.4 with transportation times = 0.This is a classical starlike system of queues with additional attached inventories at the branches.For this system Otten (2017) [Theorem 3.4.4]has found the stationary distribution in product form under the branching regime (termed "weak priorities") from Daduna (1985).
the stationary distribution of an exponential single server queue with Poisson-λ arrival stream and state dependent service rates μ (•).Then the result of Theorem 4.4 for where is a probability measure on state space Corollary 4.5 says that the stationary distribution of Z is separable with factors π , = 1, . . ., J , and θ .Separability opens the path for successful investigation of a system because several performance metrics are then directly accessible.We shall henceforth denote by a random vector which is distributed according to the stationary distribution (4.8), respectively (4.17).
Possibly, the most important performance metric for production systems is the throughput which quantifies the total output of the system.It is defined either as a time average over a long time horizon or as a steady state measure per time unit.Roughly, the throughput quantifies the utilization of the system's resources.More formally, we define as the total amount of served demand at location j within time horizon [0, T ] with base stock levels b 1 , . . ., b J and the center being located at x.This definition incorporates (contrasting to the standard one-period optimization procedures in deterministic location-inventory problems) a time dimension in the sense of Tapiero (1971) which is adequate for strategic decisions which are concerned e.g. with decision making for the replenishment center.
The relevant performance criterion for the composite system is then the asymptotic time average of total satisfied demand Because we consider a stable (ergodic) system described by a Markov process the asymptotic time average is the stationary amount of served demand over all locations per time unit which is For readability we shall often abbreviate T H := T H(b 1 , . . ., b J ; x) and similarly other terms.
Theorem 4.6 The steady state mean number of served customers (demand) per time unit at location j ∈ {1, . . ., J } is the local throughput at j: Proof To simplify notation we consider the case j = J and use the abbreviation d (x) =: η −1 again.The product form of the steady state distribution applies in ( ) below to cancel the queue lengths related terms together with 1 (n J >0) • μ J (n J ) and yields Using the extended notion for the routing probabilities from Definition 4.1 we have We insert this into the sum (and common factor) which is relevant for location J and obtain where we reintroduced the shorthand notation h This finishes the proof.
T H is not the total throughput (= overall expected number of departures in the network per time unit) of the closed system consisting of the replenishment server, the transportation and inventory nodes with stationary distribution θ .Nevertheless, the T H j resemble the expressions of local throughputs at specified nodes in a Gordon-Newell network with these exponential nodes.
The state dependent routing probabilities from Definition 4.1 determine an adaptive routing control for sending out produced items from the replenishment server.r J +1,(•) (m 1 , k 1 ; . . .m J , k J ) selects the destination for a finished item according to the free stock capacity, including items on transport.Conditioned on state (m 1 , k 1 , n 1 ; . . .m J , k J , n J ) ∈ E the probability for sending an item to location j is proportional to (b j − (m j + k j )) + .Being interested in global performance of the queueing-inventory-location system it is natural to ask whether it is possible to compute the steady state expected routing probabilities.By ergodicity these expectations would be the long time averages of the routing process as well.We determine the "expected virtual routing probabilities" in the sense of distinguishing virtual waiting times from actual waiting times in a queueing system: We compute in steady state the probability that at time t ≥ 0 an item would be send out to location j if a departure from node J + 1 would occur at time t, j = 1, . . ., J , see e.g.Wolff (1989) [Section 5.13] .
Proposition 4.7 The expected routing probability to location to j, j = 1, . . ., J , is Proof We demonstrate the case j → J and abbreviate again d (x) =: η −1 .Using separability at ( ) below we obtain For the terms referring to J we obtain for fixed (m 1 , k 1 ; . . .m J −1 , k J −1 ) (using the extended notation for the routing probabilities and noticing h A bit surprising with the result of Proposition 4.7 is that r (J + 1, j) is proportional to T H j .From the hindsight, a little reflection indicates that this should be the case.

Remark 4.8
The state-dependent routing protocol of Definition 4.1 due to Towsley (1980) has found successors in the literature which often requires additional structural properties of the networks, e.g.local balance properties or quasi-reversibility of the nodes considered in isolation, for reviews see e.g.Huisman and Boucherie (2011) or Balsamo (2000).
Although the stationary distribution (4.8) has an appealing explicit form it usually poses difficulties for computational evaluation of, e.g., throughputs.Similarly to standard closed queueing networks computational algorithms can be developed, see Sauer (1983) for a Convolution Algorithm to compute normalization constants and Mean Value Analysis for determining more detailed information on expected queue lengths, expected passage times, etc. Explicit formulas are provided there for branches of length 1, i.e. consisting of a single queue.

Cost analysis, revenue
The explicit form of the stationary distribution (4.8) allows to compute cost for running the system over a long time horizon, respectively in stationary state per time unit, under the prescribed system layout.Because revenue is obtained from successfully served demand we have immediately from Theorem 4.6 the Proposition 4.9 Let e j > 0 be the profit obtained from serving one unit of demand at location j.Then the total expected revenue R per time unit is The costs for running the system include (assumed all to be > 0) • w( j): waiting cost for a customer (demand) at the server of location j (waiting or in service), • i( j): holding cost for an item on stock in the inventory at location j, • t( j): transport cost for an item (on a truck) on the way from replenishment server J + 1 to location j, • s( j): shortage cost for each demand rejected and lost at location j, • c( j): capacity cost per time unit for providing the inventory at location j (e.g.rent, insurance), • w(J + 1): waiting cost for an order at the replenishment server at location J + 1 (waiting or in service).
The cost function is The time average costs over long time horizon are (asymptotically) determined by ergodicity of Z .
Proposition 4.10 Denote by (V 1 , Y 1 , X 1 ; . . .V J , Y J , X J ) a vector which is distributed according to the stationary distribution π of the queueing-inventory-location system.The stationary costs per time unit of the system are The expectations E π j (n j ) can be directly evaluated, while the E θ (m j ) and E θ (k j ) can be computed using mean value analysis, and λ j − T H j can be computed using the convolution algorithm following Sauer (1983)[Section 3].

Proof
The expectation (4.23) is where ( ) follows from separability of the system.To transform the last expression ) ) we note that (4.21) in the proof of Theorem 4.6 can for any j be expressed as

Generalizations and modifications
We have presented our composite model in a simple form which to a certain extend included assumptions which oversimplify the system with the aim to provide easier access to our results.In this section we first sketch some easy generalizations.Thereafter we discuss a class of related problems which seemingly need much more effort to approach explicit solutions.
Replenishment server The single exponential server with state-independent service rate ν can be substituted • by a node with general service discipline Daduna (2001)[Definition 9.1] with exponential service time request or a node with symmetric service discipline Daduna (2001)[Definition 9.5] according to Kelly (1979)[Section 3.3], or • by a replenishment network as in Otten et al. (2020).A necessary restriction is that departure from the replenishment network into the transportation subsystem is from a unique dedicated "departure node".
The proofs of both of these generalizations follow from Towsley (1980).
Modeling transport times For simplicity of presentation we assumed that for given distance d j (x) = d(x, a j ) the mean travel time (with travelling speed 1) is d j (x), and we applied an exponential-d −1 j (x) distribution for the service time at the infinite server which models travelling the lane.As a result of the so-called "insensitivity property of infinite servers" we can apply other travel time distributions with the same mean and obtain an explicit stationary distribution for a "supplemented Markov process" with a state space where the infinite server components (here the Y j (t)) carry as additional information the residual service time (travel time) or the obtained service time as supplementary variable, for more details see Towsley (1980) and Chandy et al. (1977).The most relevant consequence for our project is obvious: We are in a position to apply our findings to more realistic situations, especially to non-deterministic travel times with means d j (x) and a small variance which allows some overtaking.
Control policies for inventories The inventory control policies for the inventories in the queueing-inventory systems in this article are base-stock policies.These policies have been investigateded in many previous research papers related to the subject of our investigation.Other inventory control policies of interest are of type (r , Q) and (r , S).We sketch here only research on queueing-inventory systems which provided as an outcome the stationary distribution in explicit expressions.Research which is relevant for networks of queueing-inventory systems with a common replenishment system and base-stock contol of the inventories is communicated in Otten (2017), Otten et al. (2016Otten et al. ( , 2020)).
Employing the popular (r , Q) and (r , S) policies seemingly poses severe difficulties in case of networks of queueing-inventory systems with a common replenishment system.In Otten (2017)[Section 7] for different classes of networks of queueing-inventory systems with stochastic non-zero lead times the local inventory control is by location specific (r j , S j ) policies which means that whenever at location j the stock size decreases to the reorder point r j a replenishment order is sent out.When the replenishment arrives at location j the inventory is filled up to maximal stock size S j .For special network structures (with transportation times = 0) for r j ∈ {0, 1} an explicit product form of the stationary distribution is obtained.Seemingly all other research articles which provided explicit expressions of the stationary distributions under (r , Q) and (r , S) policies are on single queueing-inventory systems.Examples are Krishnamoorthy and Narayanan (2013), Baek and Moon (2014), Saffari et al. (2013), Schwarz et al. (2006).More details are provided in the survey Krishnamoorthy et al. (2021) [Section 5.2.1.1].
Results aside from these streams of research are provided in Schwarz et al. (2007).In a classical Jackson network, considered as a Flexible Manufacturing System (FMS), to some servers an inventory is attached and manufacturing at this server needs exactly one item from the associated inventory.The inventory control is by either (r j , Q j ) or (r j , S j ) policy with node-specific control limits.

Location analysis
The overall productivity of a system is usually measured by its throughput, i.e. the steady state mean number of produced items per time unit.Similarly, for the queueing-inventorylocation system the most important single performance measure is the satisfied overall long time external demand (≡ total production over all locations).This is for systems running under stability conditions the expected number of produced units (demanded from the exterior), i.e. the overall throughput as defined in (4.19) in Sect.4.1.Theorem 4.6 enables us to compute this directly for a given system layout as dealt with in Theorem 4.4, i.e. with fixed locations, base stock levels, and local production functions (the μ j (•)) and an adaptive routing function r J +1, j (•).
In this section we discuss the possibility to find for given locations a j = (a j1 , a j2 ), j = 1, . . ., J , a location x = (x 1 , x 2 ) for the replenishment center such that throughput is maximized.We therefore indicate in the overall throughput T H(x) the location x = (x 1 , x 2 ) as decision variable.
where E (x) (•) refers to the expectation under process distribution determined by the central location placed at x.We write T H j (x) for local throughput of the server at location j with center located at x. Typically this optimization problem is considered as a strategic decision problem.The classical problem is to determine for given a j and some distance measure d(•, •) a location x such that for certain problem specific weights v j ≥ 0, j = 1, . . ., J , the total weighted distance is minimized.This is the well known "Generalized Weber Problem".To find suitable weights v j usually is a non trivial problem.Information on the Weber problem can be found in Wesolowsky (1993) or Drezner et al. (2004).
The problem we are faced with is that we are not only interested in minimizing (weighted total) distances, but in maximizing the overall satisfied demand which is a global revenue.Additionally, the (strategic) decision for an optimal location has to pay attention for longtime development of the system.Especially one needs a prediction of the future demand.Consequently, we start our investigation with setting in force the Assumption 5.1 The long time expected demand λ j > 0 per time unit at the locations j = 1, . . ., J , is (approximately) known.
This assumption is reasonable if the amount of satisfied demand is a decision criterion.This assumption imposes a necessary condition on the system's specification if we are interested in running a system which stabilizes over time.The proof is a direct consequence of Theorem 4.4.
Corollary 5.2 For stabilizing the queueing-inventory system the service rates μ j (•) (≡ capacity functions) of the production servers at the locations must be (asymptotically) larger than the arrival rates λ j , j = 1, . . ., J , in the sense that holds for all j = 1, . . ., J .For state independent service rates μ j this means λ j < μ j .
Consequently, we set in force Assumption 5.3 The service rates can (and will) be selected in a way that (5.3) is fulfilled for all j = 1, . . ., J .(But the μ j (•) need not be fixed in advance according to the next corollary.) A somewhat surprising but important consequence for the global optimization problem to find the location for the center is the following corollary, the proof of which is direct from Theorem 4.6 under Assumption 5.3.

Corollary 5.4
The production rates μ j (•), j = 1, . . ., J , at the locations are not relevant for placing the central location such that the overall throughput T H (x) is maximized.
So, we have already proved, although in a stylized (but rather detailed) model, that the strategic decision for location of the center and the tactical decision for installing service capacities at the locations can indeed be separated.This is different from other integrated problem settings, e.g. with locating depots and transportation decisions in Salhi and Rand (1989).Nevertheless, the rates μ j (•), j = 1, . . ., J , will reoccur below.
The main problem is now to maximize with respect to location x in the plane with H (b 1 , . . ., b j , . . ., b J , {1, . . ., J }) given in Theorem 4.4.

The reduced model
We did not succeed in solving this maximization problem for the detailed queueing-inventorylocation model investigated so far.The main difficulties are to incorporate the specified base stock levels b j and the adaptive routing for departures from the replenishment server J + 1 which evaluates online the local total stock sizes (onhand inventory + items on transport).We therefore reduce the complexity of the model to obtain tractable processes.The main components of the reduced model are • the central replenishment server where items are produced in a make-to-order regime, • transportation of items from the center to the locations, • the locations with inventories where items are delivered with a rate that enables the server to meet with the location's throughput the external demand of the original model (in the reduced model), and • reorders which are send without time delay from the inventories to the replenishment servers when an item is delivered, and • an overall limit on the total stock capacity.
The restrictions imposed by the local base-stock levels are weakened.Instead of prescribing local restrictions b j , we fix the total number of "orders send out + items on transport + items as onhand-inventory" = b 1 + • • • + b J =: b.Because we assumed for the local base-stock levels that b j ≥ 1, j = 1, . . ., J , holds, b ≥ J holds for the original system.For easier reading we allow in this section b ≥ 1.In the final Algorithm 5.16 for determining the base-stock levels the restrictions b j ≥ 1 will be re-invented.
This reduced model is a closed starlike network of Gordon-Newell type with 2 • J + 1 nodes and b customers cycling (having different interpretation at the different nodes of the network).The central node is the exponential-ν replenishment single server and there are J branches numbered 1, . . ., J .Each branch consists of a tandem with (i) an exponentiald j (x) −1 infinite server (which represents transport), and (ii) a location with state dependent exponential-μ j (•) single server, j = 1, . . ., J .An item which is served (delivered at a location) departs immediately from that location and is sent to the center, occurring there as an order.
The most difficult problem is to select suitably the routing probabilities for items on leave from the central server.Recall that the predicted loads for the locations are λ j , j = 1, . . ., J (Assumption 5.1).
A simple policy is to distribute the total amount of items produced by the replenishment server (≡ the network's total throughput!) proportional to these predicted loads.Therefore the local throughputs should be proportional to the local demands, which implies fair sharing of the available capacity according to the locations' demand.Selection of this dispatching policy is supported by the following observation: From Proposition 4.7 we know that in the detailed model from Sect. 4 the expected routing probabilities are proportional to the local throughputs.
This suggests to apply in the reduced model a routing scheme R = ( r (J + 1, j), j = 1, . . ., J ) according to (5.4) Summarizing this construction, we describe the evolution over time of this reduced model by Definition 5.5 Let for j = 1, . . ., J , denote V j (t) the number of items on the way to location j at time t ≥ 0, Y j (t) the size of the inventory at j at time t ≥ 0, and Y J +1 (t) the number of orders waiting or in service at the central replenishment server J + 1 at time t ≥ 0. The network process Z = ( Z (t) : t ≥ 0) with ) is an ergodic Markov process with state space The intensity matrix (infinitesimal generator) Q = ( q(z, z ) : z, z ∈ E) of Z has the following positive rates.For readability we only depict the relevant local states for the respective transitions, e.g.we write (k J +1 ; . . .m j , k j ; . . . ) instead of (k J +1 ; m 1 , k 1 ; . . .m j , k j ; . . .m J , k J ).For j = 1, . . ., J we have q((k J +1 ; . . .m j , k j ; . . .), (k J +1 − 1; . . . 123 The diagonal elements q(z, z) of Q are selected in a way that row sums are zero.All other rates q(z, z ), z = z , are zero.
The intensity matrix Q determines the steady state distribution of the Gordon-Newell network process Z and the asymptotic throughput (respectively the steady state throughput) of the network.We summarize without proofs the facts necessary to evaluate the performance of the reduced model.These are standard results in Gordon-Newell network theory, see e.g.Daduna (2001).
Proposition 5.6 The Markov process Z is ergodic.If the center J + 1 is located at x = (x 1 , x 2 ), the stationary distribution ψ = ψ(b; x) is with G(b, {1, . . ., J , J + 1}; x) as normalization constant (5.5) . ., J , are obtained as the probability solution of the standard traffic equation.Note that α j occurs twice in (5.5) for each j = 1, . . ., J .The network's total throughput, i.e. the overall departure rate from all nodes is (5.6) The throughput from location j, i.e. the departure rate from the single server node at location j is (5.7) The throughput from all locations, i.e. the sum of the departure rates from all single server nodes at the locations 1, . . ., J is Before proceeding, a short discussion is necessary.A first important observation about the predicted demand rates λ j occurring implicitly in (5.5) and (5.6) is that within the factors (α j • d j (x)) and (α j /μ j (•)) of the stationary distribution only the relative demands The second observation is that in taking instead of (b 1 , . . ., b J ) the sum b := b 1 +• • •+b J , we have "averaged" the base stock levels.
This leads to the third observation that the stationary distribution (5.5) and the throughput (5.8) depend on the decision variable x ∈ R 2 for locating the center and on the potential decision variable for providing the overall inventory capacity b.Both variables will turn out to be of fundamental importance for optimal responding to the total predicted demands because • x will be the result of a strategic decision for locating the replenishment center, • b will be a result of tactical decisions if the sum of the demands λ 1 + • • • + λ J are given.
The main problem is that these decision variables are intertwined via (5.6) and optimization of T H(b; x) has at a first glance to be concurrently in b and x according to T H loc (b; x) subject to T H loc, j (b; x) ≥ λ j , j = 1, . . ., J .
Note that the side constraints guarantee that the total produced and sent-out items are delivered in such a way that the local demands are satisfied.It will turn out that similar to the solution procedure in Lange and Daduna (2022)  T H loc (b; x) .
If the manufacturing capacity ν of the replenisment server J + 1 is too small compared with the demand that has to be delivered there might be no solution of the Optimization Problem 5.7 because there emerge bottlenecks in the network.Nevertheless, all optimization sub-problems in Optimization Problem 5.8 have solutions.We assume (if necessary) that the capacities guarantee that a feasible solution of the Optimization Problem 5.7 exists, i.e. the side constraints in the Optimization Problem 5.7 can be satisfied.The astonishing fact is that the solutions (in fact there will be only one common solution for all b) to all sub-problems of Optimization Problem 5.8 satisfy automatically the side constraints in Optimization Problem 5.7, whenever the ν and the μ j (•) guarantee the existence of a feasible solution.
The Optimization Problems 5.7 and 5.8 will be solved in a way similar to the location problem in "logistics and services networks under congestion" in Lange and Daduna (2022).It turns out that the present problem can be attacked almost in the same way.We therefore will only sketch proofs and refer for details to Lange and Daduna (2022).In our solution procedure we will exploit the standard "Weber Problem".For an introduction with general remarks on extensions and on the history of the problem see Wesolowsky (1993) and Drezner et al. (2004).The next theorem is a modification of Theorem 4.4 in Lange and Daduna (2022) [Section 4.2].
Theorem 5.9 Consider locations a j = (a j1 , a j2 ) ∈ R 2 , j = 1, . . ., J , in the plane and associated weights v j := λ j /(λ 1 + • • • + λ J ).Let x * ∈ R 2 be a solution of the standard Weber problem with weighted distances: (5.9) Then x * is a solution of the Optimization Problem 5.8 for any b ≥ 1 as well.
The proof relies on the observation that the reduced model for the queueing-inventory-location network can be described in terms of a closed queueing network of Gordon-Newell type as sketched above.The proof of the next theorem will be given implicitly by proving correctness of Algorithm 5.13 below.
Theorem 5.10 If a solution of the Optimization Problem 5.7 exists for the capacity ν, the service intensities μ j (•), and the demands λ j , j = 1, . . ., J , the solution for the optimal value of b ≥ 1 is uniquely determined if x * is given.
The results of Theorems 5.9 and 5.10 are striking and comments are necessary.(i) If the Optimization Problem 5.7 has a solution, i.e., the side constraints are satified with capacities ν and μ j (•), j = 1, . . ., J , these capacities do not matter for optimizing the overall throughput with respect to the location of the center.Similarly, the number b of total inventory is not relevant for finding the optimal location x * .The relevant information for the location decision only comprises for j = 1, . . ., J , • distances d j (x) := d(a j , x), which determine travel times, and of items to be dispatched to location j.
(ii) It is intuitive that increasing the loading capacity at the center increases the throughput at any location.However, the following is less intuitive: Fix the capacity at the replenishment center and the service rates at all but one dedicated location and increase the service rate at the dedicated location, then the throughput at all locations increases.Both facts are consequences of Theorem 14.B.13 of Shaked and Shanthikumar (1994).Theorem 5.9 guarantees that under changes of capacities at the locations the decision for the optimal location has not to be renewed, as long as the solution of the Optimization Problem 5.7 exists.
(iii) Theorem 5.9 does not propose that throughput at the locations is independent of local properties.Details about functional dependencies will be given below.Moreover, it is not clear in advance whether a prescribed overall throughput can be met with sufficiently high total inventory capacity b by a given set of parameter values.If this is the case, Theorem 5.10 (implicitly in connection with the main result of van der Wal (1989)) guarantees that by successively adding inventory we can increase the throughput until the total requirements from the locations can be satisfied.Otherwise, if capacities do not suffice, bottlenecks occur.The relevant proofs show, that we can increase the throughput by increasing the speed of production ν at the replenishment center (capacity).Theorem 5.9 states that in any case the selected position for the center remains optimal.
The next lemma indicates how the "Weber Problem" enters our solutions.

Proof
The proof is along the lines of the proof of Theorem 4.7 in Lange and Daduna (2022).We sketch some details to fix ideas.
We finish the proof by observing that for g ∈ {0, . . ., b} holds The representation (5.11) of G(b, {1, . . ., J , J + 1}; x) shows that it is the same normalization constant as that of a Gordon-Newell network with b customers, a single server with rate ν, J stations with the service rates μ j (•), and a single infinite server with visit ratio 1/3 and exponentially distributed service time with mean ).The consequences of Lemma 5.11 are manifold.The main benefit is helping to prove Theorem 5.9.It separates the sum of the weighted travel times from the other ingredients of the network characteristics.This weighted sum is exactly what the Weber problem is about: For any solution x * of the Weber problem (5.2) with v j = λ j /(λ 1 + minimizes this weighted sum over x ∈ R 2 .The proof of Theorem 5.9 follows the lines of the proof of Theorem 4.4 in Lange and Daduna (2022), which is rather technical and tedious.A short path to the main idea is as follows.
We have to investigate throughputs i.e.quotients of normalization constants.According to Lemma 5.11 these can be split into a term which can be interpreted as coming from a single infinite server node, and a second term which reminds normalization constants of a single server Gordon-Newell network.Realizing that throughput of a closed network with only one infinite server node increases with the speed of service (≡ speed of transport which is maximized by the solution of the Weber problem) we have to show that the attached Gordon-Newell parts, represented by the C g (b, J + ), do not disturb this isotonicity.A main technical ingredient is the monotonicity property of throughputs in Gordon-Newell networks with non decreasing service rates of van der Wal (1989) plus some combinatorial identities.
A second benefit of Lemma 5.11 is that it enables us to construct a rather simple algorithm to evaluate the optimized throughputs successively in the number b ≥ 1 of customers (items) in the system.This enables us to find the minimum total size of needed inventory to realize the required throughputs.
Remark 5.12 It is possible that the location of the center coincides with one of the locations.If the center's location is, say x = a j , this results in distance to this node d j (x) = 0 and the service times at the infinite server (travel times) in branch j are zero.In this case we always have m j = 0.

Determining needed inventory sizes
We demonstrate the power of Lemma 5.11 by showing how to determine efficiently (i) the total inventory needed to fulfill the overall demand, and (ii) the size of the local inventories.We assume that the center's location x and capacities μ j (•) and ν are fixed and sufficiently high to satisfy demands λ j eventually.We assumed up to now that trucks travel from replenishment station to the locations with unit speed (1 km/hour), which implies that d j (x) is exactly the mean time to travel distance d j (x).For the present demonstration we allow general speed S > 0 for trucks.The mean time for travelling distance d j (x) is then the mean service time d j (x)/S at the infinite server j, j = 1, . . ., J .
Compute with μ j (k) from Assumption 5.3, j = 1, . . ., J , and replenishment rate ν Output is the minimal total inventor y needed to guarantee the required total stock size for satisfying demand λ Proof Following the discussion before Definition 5.5, G(b, J + ; x) from Proposition 5.6 can be interpreted as normalization constant in a Gordon-Newell network with nodes 1, . . ., J +1 having visit ratios α j and service rates μ j (n j ) at stations j = 1, . . ., J , ν at station J +1, and an additional infinite server node with visiting ratio 1/3 and exponential-( are nondecreasing in n j and n, from van der Wal's theorem in van der Wal (1989) it follows that the throughput of this artificial Gordon-Newell network is nondecreasing in b.As can be seen from the proof in van der Wal (1989) the throughput is strictly increasing in b.This guarantees that the algorithm stops after a finite number of iterations, because we assumed that the capacities are high enough to satisfy all the demands eventually.
Taking b ← b + J guarantees that minimal base stock level 1 can be guaranteed for all locations.
A direct consequenc is that with Algorithm 5.13 we obtain automatically that the total throughput, which is ≥ λ, is distributed to the locations as desired.
Corollary 5.14 Assume that for the reduced model (see Proposition 5.6) the optimal location for the center at x * is found and the necessary total inventory to satisfy the overall demand It remains to determine the distribution of the total inventory over the locations which is needed to re-transform the results from the reduced model to the model from Sect. 4. Because we did not restrict the number of items present at, say location j, in the reduced model all b items can "queue up" at j.The problem is to restrict the inventory places positioned at j in a way that in the long run (≡ in steady state) the mean overflow at j is minimized for all j ∈ {1, . . ., J }.We need the following conditional probability of the vector of joint inventory positions.The proof is immediate from Proposition 5.6.Lemma 5.15 Let Z = ( Y J +1 , V 1 , Y 1 , . . ., V J , Y J ) denote a vector which is distributed according to the stationary distribution ψ = ψ(b; x * ) from Proposition 5.6 (see (5.5)) with b from Algorithm 5.13.The conditional distribution of ( Y 1 , . . ., Y J ) given ( Y J +1 = 0, V 1 = 0, . . ., V J = 0) is with (5.12) The result of the Lemma 5.15 is surprising because P( Y 1 = k 1 , . . ., Y J = k J | Y J +1 = 0, V 1 = 0, . . ., V J = 0) is independent of x, respectively x * .This is a consequence of the product form stationary distribution.Therefore a tempting heuristic to distribute the total inventory is as follows.Remark.The generalzations described in Sect.4.3 apply to the reduced model directly as a consequence of insensitivity theory for queueing networks, see Daduna (2001)[Section 9].

Conclusion
We have developed a model for to integrate in a single Markovian system the features of manufacturing at various locations, inventory holding and control, replenishment production, transportation with dispatching according to system state dependent adaptive regimes.This enabled an optimization procedure with respect to an overall performance measure for the system.We developed two models, which together enabled us to explicitly compute the stationary distribution of the respective system in product form and to solve the integrated location problem.Possibly the most important benefit is to single out the decision for the location of the central manufacturing unit in terms of a generalized Weber problem.
According to the classification in Min et al. (1998) or Nagy and Salhi (2007) the problem tackled in this article is related to location-allocation problem because when the location of the replenishment center is fixed it is assumed that there are only radial trips from the center to the locations.The problem is related to location-routing problems because it deals with "location planning with tour planning aspects taken into account, see Nagy and Salhi (2007)[p.1]",and because we have to decide about subsequent radial tours of the trucks to visit different locations by applying the congestion dependent scheduling regimes.The relation to transportation-location problems in the sense of Cooper (1972Cooper ( , 1976) is obvious because we consider both of these problems, especially that the center has limitations on its capacity to manufacture and ship the items Cooper (1972)[p. 94].Addressing all these sub-problems, the investigation in this article is related to location-inventory-routing problems as considered in Song and Wu (2022).Our optimization criterion encompasses furthermore a time dimension as it is considered in Tapiero (1971) for the classical setting of Cooper.Aspects considered in the present article, which are not included in the above standard problem classes, are: • The effect of limited manufacturing capacities at the locations and the central production facility.These limitations result in congestion, queueing problems, and blocking, respectively rejection of demand.• The effect of including "time dimension of transportation-location-allocation problems" in the optimization criterion as discussed in Tapiero (1971)[p. 383].In Tapiero (1971)[Section 5] questions concerning delivery time lags are sketched.
and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

Appendix: Buzen's Algorithm
We present the version of Buzen's algorithm (which is often termed "convolution algorithm") using notation which is directly applicable to prove the Algorithm for determining the minimal number of overall inventory.We have introduced the stationary distribution of the reduced model as the stationary distribution of a Gordon-Newell network with probability solution of the traffic equation.Buzen's algorithm, see Bruell and Balbo (1980)[Section 2.2.1], applies without assuming this restriction.Because our application is to compute the C g (b, J + ) from Lemma 5.11 we have J +1 j=1 α j = 2/3.