Multi-Server Queuing Production Inventory System with Emergency Replenishment

: We consider a multi-server production inventory system with an unlimited waiting line. Arrivals occur according to a non-homogeneous Poisson process and exponentially distributed service time. At the service completion epoch, one unit of an item in the on-hand inventory decreases with probability δ , and the customer leaves the system without taking the item with probability ( 1 − δ ) . The production inventory system adopts an ( s , S ) policy where the processing of inventory requires a positive random amount of time. The production time for a unit item is phase-type distributed. Furthermore, assume that an emergency replenishment of one item with zero lead time takes place when the on-hand inventory level decreases to zero. The emergency replenishment is incorporated in the system to ensure customer satisfaction. We derive the stationary distribution of the system and some main performance measures, such as the distribution of the production on/off time in a cycle and the mean emergency replenishment cycle time. Numerical experiments are conducted to illustrate the system performance. A cost function is constructed, and we examine the optimal number of servers to be employed. Furthermore, we numerically calculate the optimal values of the production starting level and maximum inventory level.


Introduction
Consider a c-server (c ≥ 2) queuing inventory model with an inventory-level dependent Poisson arrival of customers. The service provided by each server has a duration that is exponentially distributed with parameter µ. On the completion of a service, the customer is served an item with probability δ. With probability 1 − δ =δ, the customer is not served the item. The replenishment of inventory is through a production process that is controlled by the (s, S) policy. To produce one unit of an item, the time required is phase-type distributed with representation (β, T) of order m. This production time can be interpreted as the time until the underlying Markov chain {ζ(t), t ≥ 0}, with a finite state space {1, 2, . . . , m, m + 1}, reaches the single absorbing state m + 1, conditioned on the fact that the initial state of this process is selected as one of the states {1, 2, . . . , m} according to the initial probability vector β = (β 1 , . . . , β m ). The transition rates within the set {1, 2, . . . , m} are defined by the generator T, and the transition rates into the absorbing state are given by T 0 = −Te.
The mean production of the item is calculated by η ′ = −βT −1 e (see Neuts [1]). All distributions involved are assumed to be independent of each other. When the onhand inventory drops to zero, one unit of the same is purchased locally (emergency Mathematics 2022, 10, 3839 2 of 25 replenishment) at a higher cost, to avoid customer loss. The time required for this is assumed to be negligible. This system is analyzed in the equilibrium state (a condition whose existence is established in the sequel). The objective of the proposed model is to minimize the total cost by determining the optimal number of servers and the optimal (s, S) policy.
As an example of the model: consider customers arriving at a car show room; a few of them decide to buy the car at the end of the advisors description of the vehicle. The remaining customers leave without buying the car. Further, as a justification for our assumption for an inventory-level dependent customer arrival rate, we notice the tendency among customers of products to queue up for an item when it is scarce with the hope that they will receive the item served within the next few replenishments. On the other end, we can cite the following situation that is commonly seen: when an item has a good brand name and is available without much of a delay, customers throng the shop where it is available, even when the waiting line is long.
Another example for the model described: If the item is abundantly available, customers are confident in receiving the item even when the queue size is large. On the other hand, if the item is scarce, naturally customers will be discouraged from joining the system. This is the reason for our assumption of an "inventory-level-dependent arrival rate" rather than "arrival rate" depending on the customers already present. In addition, it may be noted that if the arrival rate depends on the number of customers already present, then the present analysis needs drastic changes-either one has to resort to a truncation procedure to achieve a LIQBD process or go for analysis of the asymptotic quasi-Toeplitz Markov chain.
The rest of the article is organized as follows. Section 2 is the literature review. Section 3 provides our notations, assumptions and system analysis. Some important performance measures, including the distribution of the production on (or off) time and mean emergency replenishment time are also given in Section 3. Section 4 provides a numerical example, sensitivity analysis and some graphical illustrations. Finally, our conclusions are given in Section 5.

Literature Review
Queuing inventory has its origin in Melikov et al. [2] and Sigman and Levy [3]. These were followed by a spate of articles by Berman et al. [4][5][6]. In none of these did the authors attempt to derive the product-form solution. The first attempt to obtain a productform solution was by Schwarz et al. [7] for the (s, S), (r, Q) and random reorder point control policies.
The crucial assumption leading to the product-form solution is that, on inventory level dropping to zero during the replenishment lead time, no new customer is admitted to the system until the next replenishment materializes. Safari et al. [8] extends Schwarz et al. [7] to arbitrarily distributed lead time. Baek and Moon [9] is another paper providing a productform solution. In their research paper, Krishnamoorthy and Raju (see [10]) introduced the concept of a local purchase in a (s, S) inventory system.
Yue et al. [11] investigated the results on a multi-server queuing inventory system with non-homogeneous Poisson arrivals, including the inventory-level-dependent arrival of customers. When the inventory level depletes to zero, one item is replenished immediately with zero lead time. They obtained a stability condition and computed the stationary distribution of the system. Furthermore, they derived the total cost function and illustrated numerical examples of an optimal (s, S) policy.
Haughton and Isotupa [12] analyzed an inventory system with lost sales with regular and emergency orders. They assumed that the lead time of regular orders and emergency orders were, respectively, stochastic and deterministic. There were higher costs per item for the emergency order. They compared the total costs for this system to a system without emergency order. Recently Yue and Qin [13] analyzed a production system with service time and a production vacation of random duration.
Refer the recent survey paper by Krishnamoorthy et al. [14] for more detailed reports on inventory with a positive service time and production-inventory system.

Notations
In the sequel, we need the following notations and abbreviations:

Assumptions
Consider an (s, S) production inventory system with c-servers. Demands occur according to a non-homogeneous Poisson process with rate λ j , which is dependent on the on-hand inventory level j, 1 ≤ j ≤ S and λ 1 ≤ λ 2 ≤ · · · ≤ λ S . It takes an exponentially distributed time with parameter µ, to serve each customer. The on-hand inventory decreases by one unit with probability δ at the moment of a service completion and with probabilityδ = (1 − δ), the customer leaves the system without taking the item at the service completion epoch.
When the on-hand inventory depletes to s, the production process is switched on. Each production is of one unit and production process is kept in the on mode until the inventory level becomes S. It takes phase-type distributed with representation (β, T) of order m to produce one item. The mean production of the item is calculated by Further, each server serves independently of the others. The arrival, service and production processes are independent of each other.
In order to prevent the loss of demands when the inventory level equal to zero, an emergency replenishment policy is adopted in the system. When the on-hand inventory depletes to zero, one item is replenished immediately with zero lead time. Thus, there is no loss of demand in this model with the emergency replenishment policy in place. The production system is such that it takes negligible time for the item produced to reach the retail shop. Furthermore, we assume that c < s. The warehouse has limited capacity as we assume that S has an upper bound-say, M. Hence, c < s < S ≤ M.
The structure of the system under study is given in Figure 1.

Analysis
Let N(t), I(t) and J(t) denote, respectively, the number of customers in the system, the number of items in the inventory, phase of production at time t. For any time t, define C(t) = 0, if production is off 1 production is on. The production process is in on mode where 1 ≤ I(t) ≤ s and is in off mode if I(t) = S. However, when the inventory level lies between s + 1 and S − 1, the production can be in on or off mode.
The behavior of the system under study can be described by the irreducible continuous time Markov chain (CTMC) Ω = {(N(t), I(t), C(t), J(t)), t ≥ 0}, which is a levelindependent quasi-birth and death process (LIQBD) with state space To define CTMC Ω, it is necessary to write down, for any pair of above states, the intensity of the transitions between those states. Thus, the infinitesimal generator of the process Ω is given by Now, consider for 1 ≤ n ≤ c, and for n > c, we have

Stability Condition
Let ϕ = (ϕ 1 , ϕ 2 , . . . , ϕ S ) be the steady state probability vector of The Markov chain is stable if and only if (see Neuts [1]) the left drift rate exceeds the right drift rate. Thus, From the relation ϕA = 0, we have From the relation (3), we have the stability condition as with

Steady State Probability Vector
Assuming that (3) is satisfied, we briefly outline the computation of the steady state system state probability. Let y = (y 0 , y 1 , y 2 , . . . ) be the steady-state probability vector of Q. Then, yQ = 0, ye = 1.
We see that y, under the assumption that the stability condition (3) holds, is obtained as (see Theorem 3.1.1 of Neuts [1]) where R is the minimal non-negative solution to the matrix quadratic equation and the boundary equations are given by The normalizing condition ye = 1 gives where

Performance Measures
In this section, we obtain expressions for a few additional system performance measures.
• The expected number of customers in the system: • The expected number of items in the inventory: • The expected production rate: where η ′ is given in (1). • The expected rate at which the production process is switched on: • The expected rate at which the production process is switched off: Mathematics 2022, 10, 3839 8 of 25 • The expected rate at which the emergency replenishment occurs: • The mean number of busy servers:

Special Case
We investigate the problem discussed in the previous section for the case of exponentially distributed production time. The purpose of simplifying the assumptions is to investigate the system.
When the on-hand inventory depletes to s, the production process is switched on. Each production is of one unit, and the production process is kept in the on mode until the inventory level becomes S. It takes an exponentially distributed with parameter η to produce one item. The remainder of the assumptions are the same as in Section 3.2.
Thus, the CTMC Ω = {(N(t), I(t), C(t)), t ≥ 0} is a level independent quasi-birth and death process (LIQBD) with state space The infinitesimal generator of the process Ω is given in (2) A 0 are square matrices of order 2S − (s + 1) and Now, consider for 1 ≤ n ≤ c, and for n > c, we have

From the relation π
Solving the above set of equations, we obtain π(i, 1) for 1 ≤ i ≤ S − 1 and π(i, 0) for The unknown probability π(S, 0) can be found from the normalizing condition The Markov chain is stable if and only if (see Theorem 3.1.1 of Neuts [1]) the left drift rate exceeds the right drift rate. Thus, which implies the inequality (17).
Under the assumption that the stability condition (17) holds, where R is the minimal non-negative solution to the matrix quadratic equation and the boundary equations are given by The normalizing condition xe = 1 gives where Next, we proceed to derive certain system characteristics that shed light on its performance. Given ϵ > 0, there exists an N sufficiently large such that Except for heavy traffic, the above approximation is near to the exact value.

Analysis of the Emergency Replenishment Cycle
The emergency cycle time T is the time length between two adjacent emergency replenishment epochs. Once the last remaining item in the inventory is taken away, an emergency replenishment will occur. Thus, a new item will be instantly purchased and brought to the inventory because it is assumed that there is no lead time for an emergency replenishment. Therefore, the state I(t) = 0 can be seen as a virtual (instantaneous) state. Now, consider the Markov chain {(N(t), I(t), C(t)), t ≥ 0} with state space {(n, i, 1), denotes the absorbing state. Thus, an emergency replenishment cycle time T is the first passage time from the state (n, 1, 1), 0 ≤ n ≤ N to the absorbing state {0}. Let This can be taken as the initial probability vector of an emergency replenishment cycle time T. Thus, there is only one item at the beginning of an emergency replenishment cycle time T. The first passage time from the state (n, 1, 1) to {0} follows a phase-type distribution with representation (α N , W N ). The corresponding infinitesimal generator is of the form W N W 0 N 0 0 where W N is a matrix of order (N + 1)(2S − (s + 1)) and is given by with Remaining matrices are given in Section 3.2. From page 46 of Neuts [1], the expectation of the emergency replenishment cycle time is given by The production process is switched on at a service completion epoch T 0 , which started with s + 1 items in the inventory and the production process in off mode. The production process, once turned on, is turned off only at an epoch T 1 at which the number of items in the inventory reaches the maximum level S. A production on cycle starts with the switching on of the production process as the inventory level reaches s from above and terminates with the inventory reaching S.
We analyze the length T(= T 1 − T 0 ) of the production on cycle as the time until absorption in a Markov chain {(N(t), I(t)), t ≥ 0} with state space {(n, i), 0 ≤ n ≤ N, 1 ≤ i ≤ S − 1} {S} where {S} denotes the absorbing state, which represents switching the production process off. Thus, the production on cycle time from the state (n, s) to the absorbing state {S} follows a phase-type distribution with representation (β N , x n (s + 1, 0) . Thus, the infinitesimal generator where Y N is a matrix of order (N + 1)(S − 1) and is given by where Thus, the expected length of production on time in a cycle is given by

Distribution of the Production off Time in a Cycle
The production process is switched off at an epoch T 0 when the inventory level touches S. The production process, once turned off, is turned on only at future epoch T 1 at which the number of items in the inventory decreases to s. We analyze the length T(= T 1 − T 0 ) of the production off cycle as the time until absorption in a Markov chain {(N(t), I(t)), t ≥ 0} with state space {(n, i), 0 ≤ n ≤ N, s + 1 ≤ i ≤ S} {s} where {s} denotes the absorbing state, which represents switching the production process on. The infinitesimal generator of this Markov chain is of the form and is given by where The production off cycle time T follows a phase-type distribution with representation Thus, the expected length of production off cycle T is given by

Additional Performance Measures
In this section, we obtain expressions for a few additional system performance measures.
• The expected number of customers in the system: • The expected number of items in the inventory: • The expected production rate: • The expected rate at which the production process is switched on: • The expected rate at which the production process is switched off: • The average emergency replenishment rate: where E T is defined in Section 3.7.3. • The mean number of busy servers:

An Optimization Problem
Based on the above performance measures, we constructed a cost function for determining the optimality of s, S and c.
Consider the following cost parameters: where E N , E I , E PR , E on , E ER and E BS are given in Sections 3.6 and 3.7.6.

Numerical Illustrations
Next, we proceed to a few numerical examples in order to bring out the system behavior with respect to certain parameters. Here, demands occur according to a nonhomogeneous Poisson process with rate λ j , which depends on the on-hand inventory level j; 1 ≤ j ≤ S and λ 1 ≤ λ 2 ≤ · · · ≤ λ S . We use where 0 ≤ γ ≤ 1 (see Alfares [15]). This is an increasing function of the on-hand inventory level. If γ = 0, then the demand rate is homogeneous Poisson process with rate λ. Choose γ = 0.1 in the following examples.

Example: Effect of λ
From Table 1, for the indicated input parameters, we note that the expected number of customers in the system increases with the increase in the arrival rate. This results in the fast depletion of inventory since the service rate is high. In turn, this results in an increased production rate since the production unit has to be on for a longer time duration. As a consequence, the emergency replenishment rate E ER increases, and the mean number of busy servers also increases.

Example: Effect of µ
A faster service rate ensures a reduced queue size and also a reduced inventory level (see Table 2). This, in turn, makes the production unit work over a longer duration in a cycle resulting in an increased production rate. Thus, emergency ordering rate drastically decreases as well as the mean number of busy servers.  Table 3: An increased production rate ensures an increase in the availability of inventory and hence a faster disposal of demands. It also results in an increased rate of production and reduced 'switching on' of the production mechanism. However, reduced 'switching on' results in an increase in emergency replenishment. This, in turn, ensures a reduction in the number of busy servers.  Table 4: δ is the probability of a customer being served one unit of inventoried item at the end of their service. Interestingly, an increase in the value of δ initially results in a decrease in the expected number of customers and then starts increasing with an increase in the value of δ. Nevertheless, the expected inventory level continues decreasing drastically with the increase in the value of δ. This is no surprise since the probability of a customer served the inventoried item results in a decrease in the expected items held. The expected production rate (E PR ) has to increase. However, E on decreases first and then increases because the production process is over a longer duration. A reasonably high production rate ensures reduced emergency purchases. The mean number of busy servers continues increasing.  Table 5 indicates that, by increasing the maximum level of inventory, the expected number of customers can be brought down. This is because an increase in the inventory level results in a larger number of services. On the other hand, an increase in the value of S results in an increase in the number of items held in the inventory (see Table 5, Figures 2 and 3).    Table 6 are indications of the fact that the production rate increases with an increase in the value of s. This is the consequences of the early start of production. However, the same thing brings down the duration of the 'production on' status (see Table 6). For a fixed value of S and varying s, E PR (see Figure 4) increases and E on (see Figure 5) decreases. But for a fixed value of s and varying S, both E PR (see Figure 4) and E on (see Figure 5) decrease. These are consequences of the delayed 'switching on' of production.    Table 7 provides the expected number of busy servers E BS (see Figure 6) with respect to s and S and the expected emergency replenishment values E ER (see Figure 7). For fixed s and varying (increasing) S, E BS and E ER , respectively, increase/decrease, which is a consequence of the greater availability of inventoried items.

Example: Optimal Cost
In order to study the variation in different parameters on the cost function described in Section 3.8, we first take the values (K, h cust , h inv , C PR , C ER , C idle , C busy ) = ($2000, $0.5, $5, $100, $250, $40, $5). Now, we numerically compute the optimal (s, S) pair (denoted by (s * , S * )) such that 4.6.1. Example (a) From Table 8, the optimal pair (s * , S * ) = (6, 16) and optimum cost F (c, s * , S * ) = $606.0135 for the input parameters γ = 0.1, c = 3, λ = 4, µ = 7, δ = 0.8 and T =  Table 8 provides the optimal running cost of the system for the indicated input parameters. The reason for the increase in cost for pairs (s, S) for fixed S while s is varied is evident: the increased holding cost. However, with S increasing, we notice that the system cost decreases drastically. This is attributed to reduced emergency purchases, reduced production cost per unit time and, to some extent, reduced holding costs. From Table 9, the optimal pair (s * , S * ) = (10, 16) and optimum cost F (c, s * , S * ) = $419.8059 for the indicated input parameters.  Table 9 provides the optimal running cost of the system for the indicated input parameters. The reason for the increase in cost for pairs (s, S) for fixed S while s is varied is evident: the increased holding cost. However, with S increasing, we notice that the system cost decreases drastically. This is attributed to reduced emergency purchases, reduced production cost per unit time and, to some extent, reduced holding costs. Now, we compute the optimal (s * , S * ) and optimal server c * such that At the optimal pair (s * , S * ) = (10, 16), the minimum cost F (c * , s * , S * ) = $390.7507 is obtained numerically when c * = 2.

Example: The relation between δ and Cost
In order to study the variation in different parameters on cost function, we first take the values (K, h cust , h inv , C PR , C ER , C idle , C busy ) = ($5000, $0.05, $5, $100, $200, $30, $5). Table 10 provides the effect of the probability of an item being served to a customer on the cost function for the input parameters S = 13, s = 6, γ = 0.1, λ = 5, µ = 6, T =  Figure 8, we found that the increase in the value of δ results initially decreases F (c, s, S) and then starts increasing with the increase in the value of δ. We obtained the minimum cost when δ = 0.2. on the cost function for the input parameters S = 13, s = 6, γ = 0.1, λ = 5, µ = 6, T = 263 Table 11 provides the optimal value of servers to be employed for the given input 264 parameters. We look at the effect of the probability of an item being served to a customer 265 on the number of servers to be employed. It is seen that with probability of item 266 being served to the customer at least equal to half the optimal number of servers to be 267 employed gets stabilized and the system cost also gets drastically cut.

Example: The Optimal Number of Servers and Minimum Cost
In order to study the variation in different parameters on the cost function, we first take the values (K, h cust , h inv , C PR , C ER , C idle , C busy ) = ($2000, $0.5, $5, $100, $125, $3, $2). Table 11 provides the optimal value of the servers to be employed for the given input parameters. We look at the effects of the probability of an item being served to a customer on the number of servers to be employed. When the probability of an item being served to the customer is equal to at least half, then the optimal number of servers to be employed is stabilized, and the system cost is also drastically cut. Table 11. Optimal server c * and corresponding minimum cost for (S, s, η, λ, µ) = (30, 12, 2.6, 5, 7). In this example, we give a comparison between the constant arrival case (λ) and inventory-level-dependent arrival case (λ i , 1 ≤ i ≤ S) in Figure 9 and found the two cases to be showing entirely different behaviors. Interestingly, an increase in the value of δ initially results in a decrease in the expected number of customers and then starts increasing with an increase in the value of δ in the non-homogeneous arrival case. However, in the homogeneous arrival case, as δ increases initially, the expected number of customers remains constant and then starts increasing.
These observations are dependent on the input values of the parameters.

Conclusions
We studied a multi-server production inventory model that uses emergency replenishment when a delay in the production process leads to an imminent stock out. Customers arrive according to a non-homogeneous Poisson process. The service time duration at all servers is independent with identically distributed exponential random variables. A customer receives one item from the inventory at the end of their service with the probability δ. When the on-hand inventory reaches zero, a local purchase (emergency replenishment) to bring the inventory level to 1 is made. Due to this, the inventory level never reaches zero. Under the stability condition, we established a long-run system state distribution. A cost function involving these control variables was constructed, and we numerically investigated the optimal values. The effects of various parameters on the system performance measures were also evaluated.
In a follow-up paper, we will extend the present work to the case of Markovian arrival process MAP and phase-type PH service time.