Queuing-Inventory System with Catastrophes in the Warehouse: Case Rare Catastrophes

1. Introduction

To increase the adequacy of the developed mathematical models of a real queuing-inventory systems (QISs), it is necessary to take into account the possibility of various risks. Among these risks, catastrophes (disasters) that lead to the destruction of all (or part) of inventories are more important. Generally speaking, in real systems, catastrophes are rare events. However, the rate of catastrophes is not so rare that it should not be taken into account when developing mathematical models of the real QISs. In other words, catastrophes are rare events, but they must be taken into account, since they significantly affect the performance of the system. Unfortunately, this real risk is not taken into account in the vast majority of works devoted to QISs models.

Note that in the available literature, the issues of the impact of catastrophes on the operation of classical queuing systems have been studied in detail, see, for instance, [1,2,3,4,5,6,7]. These issues have been studied in connection with assessing the reliability of servers in queuing systems, since in them catastrophes are interpreted as failure of the system servers, see [8,9] and their references lists.

The catastrophes in a QISs warehouse gives rise to new problems, and these problems have been poorly studied in the available literature. Models of QMS with randomly catastrophes that destroy only one item in a warehouse were developed in works [10,11], in which catastrophes are interpreted as destructive customers. Recently, QISs models with the possibility of a catastrophe destroying all items in the warehouse have been considered in [12]. The vast majority of works, including the mentioned papers devoted to models with catastrophes, study QISs models with an infinite queue. To analyze such models, the well-established in classical book [13] matrix-geometric method (MGM) is often used, for more detailed describing and applications of MGM see recent two volume books [14,15] as well as book [16]. However, models with an infinite queue are not adequate to the real situation, since in real QISs, as a rule, there are finite waiting room. However, models of QISs with finite waiting room for c-customers have been poorly studied, see works [17,18,19,20,21,22,23].

This work is a continuation of the research begun in the indicated above work [12]. Namely, here we study a Markov QIS model with catastrophes that destroy the entire warehouse and with a finite buffer for waiting consumer customers (c-customers). Note that imposing a restriction on the length of the queue of c-customers leads to a fundamental change (compared to the infinite model) in the methodology for studying the system. In addition, new performance measures of the system appear here. We understand that the Poisson/exponential assumption is a rough approximation to the real situation, but, firstly, this paper is one of the first works in this direction, and secondly, here a methodology for studying these systems is proposed.

The main contributions of this work can be summarized as follows:

• Exact and approximate methods to study QIS model with finite waiting room are developed.
• High-accurate closed-form approximate formulas for calculating the steady-state probabilities, as well as performance measures of the investigated QIS in the case of rare catastrophes are developed.
• The developed approximate formulas make it possible to calculate the performance measures of large-scale QISs using closed-form formulas and minimize the expected total cost (ETC) by choosing the optimal value of reorder point.

This paper is organized as follows. In Section 2, we describe the model of investigated finite QIS. Exact and approximate methods to calculate the steady-state probabilities as well as performance measures are developed in Section 3. The results of numerical examples are shown in Section 4. Concluding remarks are given in Section 5.

2. The Model

Consider single server finite QIS in which warehouse has maximum capacity $S$. Input flow of homogeneous c-customers forms Poisson process with rate ${\lambda }^{+}$. Customer homogeneity means that each customer requires the same amount of inventory. The service times of the c-customers are independent identically distributed (i.i.d.) random variables with a common average value ${\mu }^{-1}$ and at the end of service the inventory level decreases by one unit. The waiting room for queuing c-customers has a finite size $R,R<\infty$. This means that if, when a c-customer arrives, there are already $R$ c-customers in the queue, then the arriving new c-customer is lost with probability (w.p.) 1; otherwise, the arriving c-customer will join the queue if the server is busy. A combined sales scheme is used, i.e., if upon arrival of c-customer inventory level is zero, then, in accordance to the Bernoulli trials, it either, w.p. ${\phi }_1$ joins the queue (backorder sale scheme), or w.p. ${\phi }_2=1-{\phi }_1$ leaves the system without inventory (lost sale scheme).

Along with c-customers, the system also receives a Poisson flow of negative customers (n-customers) with rate ${\lambda }^{-}$. Negative customers require no service or inventory, but upon the arrival of n-customers, one c-customer is pushed out from the system, if any. Procedure of pushing out the c-customer is managed as follows: (1) if there is a queue of c-customers, then only the c-customer is pushed out of the queue; (2) if there is no queue of c-customers but only c-customer is on server, then a n-customer force out of the system a c-customer which is on server (in these cases the inventory level does not change, since it is assumed that stocks are released after the completion of servicing a c-customer); (3) if there are no c-customers in the system (in queue or on server), then arrived n-customer does not affect the operation of the system.

Catastrophes form a Poisson flow with parameter κ, and when a catastrophe occurs, all inventory items are instantly destroyed. A catastrophe destroys even items that are in the status of release to the c-customer. In this case, the interrupted c-customer returns to the queue, i.e., the catastrophe only destroys the inventory and does not force c-customers out of the system. Catastrophes does not affect the operation of the system warehouse, if the inventory level is zero.

In order to be specific, here $\left(s,S\right)$ inventory replenishment policy is considered (sometimes this policy is called “Up to $S$ ” as well). This means that when the inventory level drops to the order replenishment point $s$, $0\le s<S,$ a replenishment order is placed and upon replenishment, the inventory level is restored to level $S$, regardless of how many items were in inventory.

The lead times of the replenishments i.i.d. random variables with a common average value ${\nu }^{-1}$.

The problem is to find the joint distribution of the number of c-customers in the system and the inventory level of the system, as well as to calculate the key performance measures: the average number of items in warehouse, the average order size, the average reorder rate, includes the average length of the queue and loss rate of c-customers.

3. Steady-State Analysis

In this section, the Markov model is developed and both exact and approximate methods to calculate the steady-state probabilities are proposed.

3.1. An Exact Approach

In this subsection, an exact method to calculate the steady-state probabilities as well as performance measures is proposed. As in Melikov et al (2023) [12], let $X_t$ be the number of с-customers at time $t$ and $Y_t$ be the inventory level at time $t$. So, the process $Z_t=\left\{\left(X_t,Y_t\right),t\ge 0\right\}$ forms a two-dimensional continuous time Markov chain (2D CTMC) with state spac

$$E={\bigcup }_{m=0}^S{E}_m,$$

(1)

where $E_m=\left\{\left(0,m\right),\left(1,m\right),\cdots ,\left(R,m\right)\right\}$ is the subset of states in which inventory level is equal to $m,m=\mathrm{0,1},\cdots ,S.$

The transition rate from state $\left(n_1,m_1\right)$ to state $\left(n_2,m_2\right)$ is denoted by $q\left(\left(n_1,m_1\right),\left(n_2,m_2\right)\right)$ . By taking into account the assumptions related to operating of the investigated QIS, we conclude that these transitions rates are calculated as follows

$$q\left(\left(n_1,m_1\right),\left(n_2,m_2\right)\right)=\left\{\begin{array}{cc}{\lambda }^{+}{\phi }_1,& m_2=m_1=0,n_2=n_1+1,\\ {\lambda }^{+},& m_2=m_1>0,n_2=n_1+1,\\ {\lambda }^{-},& m_2=m_1,n_2=n_1-1,\\ \mu ,& m_2=m_1-1,n_2=n_1-1,\\ \mathrm{\kappa },& m_1>0,m_2=0,n_2=n_1,\\ \mathrm{\nu },& m_1\le s,m_2=S,n_2=n_1.\end{array}\right.$$

(2)

From relations (2) we conclude that each state of the constructed 2D CTMC can be reached from any other state through a finite number of transitions, i.e. considered chain is an irreducible one. In other words, for any positive values of the initial parameters, a stationary mode is existing. Let us denote by 𝑝(𝑛, 𝑚) the probability of the state (𝑛, 𝑚) ∈ 𝐸. These probabilities can be obtained by solving the system of balance equations (SBE) which is constructed by using relations (2).

$$\begin{array}{c}\mathrm{Case}1:\mathrm{When}\left(n,0\right)\in E_0:\\ \left({\lambda }^{+}{\phi }_1\mathrm{\chi }\left(n<R\right)+{\lambda }^{-}\mathrm{\chi }\left(n>0\right)+\mathrm{\nu }\right)p\left(n,0\right)={\lambda }^{+}{\phi }_{1}p\left(n-\mathrm{1,0}\right)\left(n>0\right)\\ +{\lambda }^{-}p\left(n+\mathrm{1,0}\right)\mathrm{\chi }\left(n<R\right)+\mu p\left(n+\mathrm{1,1}\right)\mathrm{\chi }\left(n<R\right)+\mathrm{\kappa }{\sum }_{m=1}^{S}p\left(n,m\right).\end{array}$$

(3)

$$\begin{array}{c}\mathrm{Case}2:\mathrm{When}\left(n,m\right)\in E_m,0<m\le s:\\ \left({\lambda }^{+}\mathrm{\chi }\left(n<R\right)+{\lambda }^{-}\mathrm{\chi }\left(n>0\right)+\mathrm{\nu }+\mu +\mathrm{\kappa }\right)p\left(n,m\right)={\lambda }^{+}p\left(n-1,m\right)\mathrm{\chi }\left(n>0\right)\\ +{\lambda }^{-}p\left(n+1,m\right)\mathrm{\chi }\left(n<R\right)+\mu p\left(n+1,m+1\right)\mathrm{\chi }\left(n<R\right).\end{array}$$

(4)

$$\begin{array}{c}\mathrm{Case}3:\mathrm{When}\left(n,m\right)\in E_m,s<m<S:\\ \left({\lambda }^{+}\mathrm{\chi }\left(n<R\right)+{\lambda }^{-}\mathrm{\chi }\left(n>0\right)+\mu +\mathrm{\kappa }\right)p\left(n,m\right)={\lambda }^{+}p\left(n-1,m\right)\mathrm{\chi }\left(n>0\right)\\ +{\lambda }^{-}p\left(n+1,m\right)\mathrm{\chi }\left(n<R\right)+\mu p\left(n+1,m+1\right)\mathrm{\chi }\left(n<R\right).\end{array}$$

(5)

$$\begin{array}{c}\mathrm{Case}4:\mathrm{When}\left(n,S\right)\in E_S:\\ \left({\lambda }^{+}\mathrm{\chi }\left(n<R\right)+{\lambda }^{-}\mathrm{\chi }\left(n>0\right)+\mu +\mathrm{\kappa }\right)p\left(n,S\right)={\lambda }^{+}p\left(n-1,S\right)\mathrm{\chi }\left(n>0\right)\\ +{\lambda }^{-}p\left(n+1,S\right)\mathrm{\chi }\left(n<R\right)+\mu p\left(n+1,m+1\right)\mathrm{\chi }\left(n<R\right)+\mathrm{\nu }{\sum }_{m=0}^{s}p\left(n,m\right).\end{array}$$

(6)

Here and below, χ$\left(A\right)$ denotes the indicator function of the event $A$. To SBE (3)-(6), a normalization condition should be added, i.e.

$${\sum }_{\left(n,m\right)\in E}p\left(n,m\right)=1.$$

(7)

The constructed SBE (3)-(7) is a system of linear algebraic equations of dimension $\left(R+1\right)·\left(S+1\right)$ and it can be solved numerically by using well-known software in case the moderate values of the parameters $R$ and $S$.

After calculating the steady-state probabilities, key performance measures of the QIS under study can be determined using a standard technique. The performance measures are divided into two groups: (1) inventory-related performance measures and (2) queuing-related performance measures. The first group of performance measures includes the average number of items in warehouse $\left(S_{av}\right)$, the average order size $\left(V_{av}\right)$ and the average reorder rate $\left(RR\right)$.

• The average number of items in warehouse (i.e. the average inventory level) is calculated as mathematical expectation of the appropriate random variable and is given by

  $$S_{av}={\sum }_{m=1}^{S}m{\sum }_{n=0}^{R}p\left(n,m\right).$$

  (8)
• Similar to (8), the average order size (i.e. the average size of replenished items from external source) is calculated as mathematical expectation of the appropriate random variable and is calculates as follows

  $$V_{av}={\sum }_{m=S-s}^{S}m{\sum }_{n=0}^{R}p\left(n,S-m\right).$$

  (9)
• An inventory order is placed in two cases: (1) when the inventory level drops to the reorder point $s$ after completing customer service in states $\left(n,s+1\right)\in E_{s+1}$, and (2) when catastrophes occur the in states $\left(n,m\right)\in E_m,m>0.$ Therefore, the average reorder rate is calculated as follows

  $$RR=\mu {\sum }_{n=1}^{R}p\left(n,s+1\right)+\kappa \left(1-{\sum }_{n=0}^{R}p\left(n,0\right)\right).$$

  (10)

The second group of performance measures includes the average length of the queue $\left(L_{av}\right)$ and loss rate of c-customers $\left(LR\right)$.

• The average length of the queue is calculated as mathematical expectation of the appropriate random variable and is given by

  $$L_{av}={\sum }_{n=1}^{R}n{\sum }_{m=0}^{S}p\left(n,m\right).$$

  (11)
• Losing c-customers occurs in three cases: (1) if at the time the c-customer arrives the waiting room is full (with probability 1), i.e. when the system are in one of the states $\left(R,m\right)\in E_m$,$m=\mathrm{0,1},\cdots ,S,$ (2) if at the time the c-customer arrives, the inventory level is zero and waiting room is not full (with probability ${\phi }_2$), i.e, when system are in one of the states $\left(n,0\right)\in E_0$,$n<R,$(3) when n-customer arrived, it displaced one c-customer. Therefore, the loss rate of c-customers is calculated as follows

  $$LR={\lambda }^{+}{\sum }_{m=0}^{S}p\left(R,m\right)+{\lambda }^{+}{\phi }_2{\sum }_{n=0}^{R-1}p\left(n,0\right)+{\lambda }^{-}\left(1-{\sum }_{m=0}^{S}p\left(0,m\right)\right).$$

  (12)

As it was mentioned above, the proposed exact approach can be used for investigation the QIS model with moderate state space. Computational difficulties arise for large-scale models, and an approximate method for calculating steady-state probabilities is then required. Below we develop an approximate method for solving this problem, which can be used in cases of rare catastrophes.

3.2. An Approximate Approach

In this subsection, we derive the closed-form approximate solution for steady-state probabilities of the investigated 2D CTMC, by using space merging approach, see [24]. This approach is highly accurate for systems with rare catastrophes, i.e. we will assume that $\kappa \ll min\left({\lambda }^{+},{\lambda }^{-},\mu \right).$ Note that the last assumption is not extraordinary, since in the opposite case (i.e. when the level of catastrophes is close to the rate of c-customers, the speed of their service, the rate of n-customers) the QIS under consideration is generally not effective.

In the case where the above assumption is fulfilled, the basic requirement for an adequate application of space merging method is satisfied. In this case, transition rates between states in each subset $E_m$ (see (1)) are much great than the transition rates between states from different subsets. So, in accordance to the space merging algorithm, each subset of states $E_m$ in (1) is combined into one merged state$<m>$, and the following merging function is determined in the initial state space (1)$:U\left(n,m\right)=<m>,\left(n,m\right)\in E$. The merged states form the set $\widehat{E}=\left\{<m>:m=\mathrm{0,1},\dots ,S\right\}.$ The approximate values of steady-state probabilities, $\widehat{p}\left(m,n\right),$ are calculate as follows

$$\widehat{p}\left(n,m\right)\approx {\rho }_m\left(n\right)\pi \left(<m>\right)$$

(13)

where ${\rho }_m\left(n\right)$ denotes probability of state $\left(n,m\right)$ within subset $E_m$ and $\pi \left(<m>\right)$ denotes probability of merged state $<m>\in \widehat{E}$.

From relations (2) we conclude that the steady-state probabilities ${\rho }_0\left(n\right),n=\mathrm{0,1},\cdots ,R,$within a split model with the state space $E_0$ coincide with the stationary distribution of a finite one-dimensional birth-death process in which birth rate is ${\lambda }^{+}{\phi }_1$ while death rate is ${\lambda }^{-}$. By the same way, from relations (2) we conclude that the steady-state probabilities ${\rho }_m\left(n\right),m>0,n=\mathrm{0,1},\cdots ,R,$within a split model with the state space $E_m$ are independent on $m,$and coincide with the stationary distribution of a finite one-dimensional birth-death process in which birth rate is ${\lambda }^{+}$while death rate is ${\lambda }^{-}$. In other words, steady-state probabilities within split models are calculated as

$${\rho }_m\left(n\right)=\left\{\begin{array}{cc}{\theta }^n\frac{1-\theta }{1-{\theta }^{R+1}},& m>0,n=\mathrm{0,1},\dots ,N,\\ {\theta }_0^n\frac{1-{\theta }_0}{1-{\theta }_0^{R+1}},& m=0,n=\mathrm{0,1},\dots ,N,\end{array}\right.$$

(14)

where ${\theta }_0={\lambda }^{+}{\phi }_1/{\lambda }^{-},\theta ={\theta }_0/{\phi }_1.$

Note 1. To simplify the notation, below the subscript $m$ for cases $m>0$ is omitted in state probabilities ${\rho }_m\left(n\right).$ In cases $\theta =1$ and/or${\theta }_0=1$, all state probabilities ${\rho }_m\left(n\right)=1/\left(R+1\right)$ for any $n,n=\mathrm{0,1},\dots ,R,$ and $m,m=\mathrm{0,1},\dots ,S$.

Let us denote the rate of transition from the merged state $<m_1>$to the merged state $<m_2>$by $q\left(<m_1>,<m_2>\right)$. Then, taking into account relations (2) and (14), we obtain the following formulas for calculating the indicated rates (all other transition rates are zero):

$$\begin{array}{c}\mathrm{Case}0\le m\le s:\\ q\left(<m>,<S>\right)=\mathrm{\nu }{\sum }_{n=0}^R{\rho }_m\left(n\right)=\mathrm{\nu }\end{array}$$

(15)

$$\begin{array}{c}\mathrm{Case}m>0:\\ q\left(<m>,<0>\right)=\mathrm{\kappa }{\sum }_{n=0}^R{\rho }_m\left(n\right)=\mathrm{\kappa };\end{array}$$

(16)

$$q\left(<m>,<m-1>\right)=\mu {\sum }_{n=1}^R{\rho }_m\left(n\right)=\mu \left(1-\rho \left(0\right)\right).$$

(17)

In other words, merged model represents one-dimensional Markov chain with state space $\widehat{E}$ where transition rates between states are calculated via formulas (15)-(17). Using the approach proposed in [12], we obtain the following closed-form formulas for calculating the probabilities of merged states:

$$\pi \left(0\right)=\frac{1+bc}{1+dc},$$

(18)

$$\pi \left(1\right)=d\pi \left(0\right)-b,$$

(19)

$$\pi \left(m\right)=a_m\pi \left(1\right),2\le m\le S,$$

(20)

where

$$d=\frac{\nu +\kappa }{\mu \left(1-\rho \left(0\right)\right)},b=\frac{\kappa }{\mu \left(1-\rho \left(0\right)\right)},$$

$$c={\sum }_{m=1}^S{a}_m,a_m=\left\{\begin{array}{cc}{\left(1+d\right)}^{m-1},& \mathrm{if}1\le m\le s+1,\\ {\left(1+d\right)}^s{\left(1+b\right)}^{m-s-1},& \mathrm{if}s+1<m\le S.\end{array}\right.$$

Eventually, taking into account formulas (13), (14), (18)-(20), we conclude that an approximate values of performance measures (8)-(12) can be calculated using the following explicit formulas

$$S_{av}={\sum }_{m=1}^{S}m\pi \left(m\right);$$

(21)

$$V_{av}={\sum }_{m=S-s}^{S}m\pi \left(S-m\right);$$

(22)

$$RR=\mu \left(1-\rho \left(0\right)\right)\pi \left(s+1\right)+\kappa \left(1-\pi \left(0\right)\right);$$

(23)

$$L_{av}={\sum }_{n=1}^{R}n\left({\rho }_0\left(n\right)\pi \left(0\right)+\rho \left(n\right)\left(1-\pi \left(0\right)\right)\right);$$

(24)

$$\begin{array}{c}LR={\lambda }^{+}{\phi }_2\pi \left(0\right)\left(1-{\rho }_0\left(0\right)\right)+{\lambda }^{+}{\rho }_0\left(R\right)\pi \left(0\right)+{\lambda }^{+}\left(1-\pi \left(0\right)\right)\rho \left(R\right)+\\ +{\lambda }^{-}\left(\pi \left(0\right)\left(1-{\rho }_0\left(0\right)\right)+\left(1-\pi \left(0\right)\right)\left(1-\rho \left(0\right)\right)\right)\end{array}$$

(25)

4. Numerical Experiments

Based on the purpose of the paper, we conducted numerical experiments that have three goals: (1) estimation of the accuracy of the proposed approximate formulas; (2) study the behavior of the performance measures versus initial parameters and (3) solve the optimization problem.

4.1. Accuracy of the Developed Approximate Formulas

The accuracy of the developed approximate formulas is investigated via numerical experiments. For this purpose, exact values of the steady-state probabilities (SSP) are determined from SBE (3)-(7) for the QIS with maximum capacity of warehouse $S=50$ and buffer size $R=30,$dimension of SBE is equal to 1581. The accuracy of the developed approximate formulas can be estimated using several norms, e.g. cosine similarity, Euclidean distance, Jaccard norm etc. To be specific, here we use a simple norm, that is, the maximum errors when calculating SSPs. Results of numerical experiments are shown in Table 1. In this table, along with an indication of the accuracy of calculating the SSPs, results are given that indicate the accuracy of calculating the performance measures (8)-(12). From this table, we conclude that the accuracy of the developed formulas to calculation of the SSPs as well as performance measures are sufficiently. From this table it is also clear that the accuracy of calculating the SSPs is greater than the accuracy of calculating performance indicators. This was to be expected, since the performance indicators are calculated through SSPs using operations of multiplication by large numbers, see formulas (8)-(12) and (21)-(25). We conducted a large number of experiments and demonstrated only a small part of them here. An interesting result of these experiments is that increasing the size of the system (i.e. increasing $S$ and $R$) leads to an increase in the accuracy of the approximate results.

4.2. Behavior of Performance Measures Versus Reorder Point

Performance measures are computed and numerically illustrated for different sets of values for the input parameters. In all experiments the values of following parameters are fixed: $S=50,N=30,{\phi }_1=0.4,\kappa =0.1.$ So, Figure 1 shows the behavior of performance measures as a functions of reorder point $s$ for the three different values of ${\lambda }^{+}$. From these plots, we conclude that inventory-related performance measures are almost independent of the rate of c-customers, but significantly dependent on the reorder point, see Figure 1, (a), (b) and (c). Note that all inventory-related performance measures are increasing functions versus reorder point. Measure$S_{av}$ is an increasing function versus $s$, and this fact was expected, see Figure 1 (a). At first glance, the increasing of $V_{av}$ versus $s$ seems unexpected, see Figure 1 (b). However, this fact has the following explanation: as $s$ increases, the reorder rate increases (see Figure 1 (c)), and as a result, the average order size increases. The $RR$ is an increasing function with respect to $s$ because as $s$ increases, the probability that the inventory level is positive, also increases and hence $RR$ becomes an increasing function, see formula (10) as well. For the selected data the rate of its increase becomes very high at large (possible) values of the reorder point, see Figure 1 (c). Opposite, the queuing-related performance measures are almost independent of the reorder point, see Figure 1, (d) and (e). For selected initial data the average number of customers in queue $\left(L_{av}\right)$for all values of ${\lambda }^{+}$ are very close to buffer size ($R=30$), see Figure 1 (d). Therefore, loss rate $\left(LR\right)$ is very close to ${\lambda }^{+}$ , see Figure 1 (e).

The dependence of performance measures on $s$ and ${\lambda }^{-}$ is shown in Figure 2. It is interesting to note that here behavior, even absolute values of the inventory-related performance measures (see Figure 2, (a)-(c)), are same as in Figure 1 (a)-(c). From Figure 2 (d) we conclude that $L_{av}$ is almost independent of $s$ and the rate of its decreasing versus ${\lambda }^{-}$ is very small, i.e., increasing ${\lambda }^{-}$ even by ten times leads to a change in the value of $L_{av}$ in the second digit after the decimal point. Similarly, from Figure 2 (e) we conclude that $LR$ is also almost independent of $s$, but here the rate of its decrease compared to ${\lambda }^{-}$ is noticeable.

The dependence of performance measures on $s$ and $\mu$ is shown in Figure 3. An increase in $S_{av}$ with respect to $s$ is obvious, since an increase in $s$ leads to an increase in the rate of filling the warehouse to its maximum size, see Figure 3 (a); on the other hand, $S_{av}$ is a decreasing function versus $\mu$ because an increase in $\mu$ results in an increase in the rate at which inventory is sold. Since $S_{av}$is a decreasing function with respect to $\mu$, therefore $V_{av}$ is increasing with respect to $\mu$, see Figure 3 (b); an increase in $s$ leads to increasing the probability that the inventory level drops to zero due to catastrophes, i.e., average order size is also increasing function versus $s$. Measure $RR$ is an increasing function with respect to $s$, since increasing of $s$ leads to increasing the probability that the stock level is greater than zero, i.e., the rate of replenishment of stocks due to catastrophes also increases, see Figure 3(c); this measure is also an increasing function with respect to $\mu$, since an increase in $\mu$ leads to increasing the probability that the inventory level drops to the re-order point $s$, i.e., the rate of replenishment increases (see formula (9) also). From Figure 3(d) we conclude that $L_{av}$ is an increasing function versus $s$ and a decreasing function versus $\mu$. And for large values of $s$, i.e., $s>25,$ the value of $L_{av}$ is practically independent of $\mu$. These facts are expected. Measure $LR$ is decreasing function versus on both $s$ and $\mu$, see Figure 3 (e). The decrease of this function with respect to $s$ is obvious, but its decrease relative to $\mu$ is not evident at first glance. The last fact has the following explanation: an increase in $\mu$ leads to increasing the probability that the inventory level is zero, i.e., the loss probability of arriving customers also increases and hence the measure $LR$ becomes increasing. However, the rate of increasing of $LR$ versus $\mu$ is very small, i.e., increasing $\mu$ even by three times leads to a change in the value of $LR$ in the second digit after the decimal point.

The dependence of performance measures on $s$ and ν is shown in Figure 4. Measure $S_{av}$ is an increasing one with respect to both $s$ and ν, see Figure 4 (a). Indeed, an increase in $s$ and ν leads to an increase in the rate of filling the warehouse to its maximum size. Measure $V_{av}$ is increasing versus $s$, while it is decreasing function versus ν. Here $V_{av}$ is decreasing versus ν, since $S_{av}$ is a increasing function with respect to ν, see Figure 4 (b); the increase in $V_{av}$ with respect to $s$ is explained as above, i.e., as s increases, the probability of inventory levels falling to zero due to catastrophes also increases, so the average order size is also an increasing function with respect to $s$. Measure $RR$ is an increasing one with respect to both $s$ and ν, since an increase in both $s$ and ν leads to increasing the probability that the stock level is greater than zero, hence the rate of replenishment of inventory due to catastrophes also increases, see Figure 4(c). From Figure 4 (d) we conclude that the $L_{av}$ is almost constant (does not decrease) depending on $s$ and ν; only for small values of $s$, i.e. s<10, we observe insignificant differences between the values of $L_{av}$ for different values of ν. Measure $LR$ is a decreasing function with respect to both $s$ and ν, see Figure 4 (e). Such behavior of this measure is expected, since an increase in both $s$ and ν leads to an increase in the average inventory level and, as a consequence, to a decrease in $LR$.

4.5. Optimization Problem

Third goal of performing numerical experiments is solving the optimization problem. To be specific, here minimization of Expected Total Cost (ETC) is considered. In this problem it is assumed that all load parameters as well as structural parameters of the QIS are fixed and only controllable parameter is reorder point. Similar to Melikov et al (2023) [12], ETC is defined as follows:

$$ETC\left(s\right)=\left(K+c_r·V_{av}\right)·RR+c_h·S_{av}+c_{ps}·\kappa ·S_{av}+c_l·LR+c_w{·L}_{av},$$

(26)

where $K$ is the fixed price of one order, $c_r$ is the unit price of the order size, $c_h$ is the unit inventory storage price per unit of time, $c_{ps}$ is the price of unit inventory destruction, $c_l$ is the cost for a single c-customer loss, $c_w$ is the price per unit time of queuing delay for a single c-customer.

The problem is to find a value (optimal) of $s$ that minimizes (26). For any values of initial parameters this problem has solution since admissible set for values of $s$ is finite and discrete, i.e., $0\le s\le S-1.$ Coefficients in (26) for hypothetical model are selected as $K=10,c_r=15,c_h=10,c_l=450,c_w=400,c_{ps}=15$. Some results of the minimization of (26) are demonstrated in Table 2. Here we assume that $N=30,{\phi }_1=0.4,{\lambda }^{+}=15,{\lambda }^{-}=1,\kappa =0.1,\mu =2,\nu =1.$ Optimal solution for indicated values of $S$ is $s*=0$. For completeness, Table 2 shows the values of the performance measures in the optimal solution, as well as the minimum value of the expected total cost that denoted by $ETC*.$

5. Conclusions

The model of single-server QIS with catastrophes in the warehouse and finite waiting room for consumer customers was proposed. When a catastrophe occurs, the entire inventory, as well as items that were in the status of release to the consumer, are instantly destroyed. However, catastrophes do not force consumers out of the system. It is assumed that if the inventory level upon arrival of consumer customer is zero, then, in accordance to the Bernoulli trials, it either joined the queue of infinite length, or left the system unserved. In the system “Up to S” replenishment policy is applied. In addition to c-customers, negative customers also enter the system. When a negative customer arrives, one of the consumers, if any, is pushed out. The mathematical model of the investigated QIS is constructed as two-dimensional Markov chain. Both exact and approximate methods for calculating steady-state probabilities and performance measures of the QIS under study are proposed. The exact method is based on balance equations while the approximate method is based on the space merging approach. It is noted that the proposed approximate method has high accuracy in the case of rare catastrophes. Closed-form formulas for calculating the performance measures were proposed. The results of numerical experiments are demonstrated.

A direction for further research should be the study of similar models with MAP flows of c- customers and/or n-customers, as well as with the PH distribution of service times for c-customers and/or lead times.