*Article* **Markovian Demands on Two Commodity Inventory System with Queue-Dependent Services and an Optional Retrial Facility**

**K. Jeganathan 1, M. Abdul Reiyas 2, S. Selvakumar 1, N. Anbazhagan <sup>3</sup> , S. Amutha 4, Gyanendra Prasad Joshi 5,\*, Duckjoong Jeon <sup>6</sup> and Changho Seo 6,\***


**Abstract:** The use of a Markovian inventory system is a critical part of inventory management. The purpose of this study is to examine the demand for two commodities in a Markovian inventory system, one of which is designated as a major item (Commodity-I) and the other as a complimentary item (Commodity-II). Demand arrives according to a Poisson process, and service time is exponential at a queue-dependent rate. We investigate a strategy of (*s*, *Q*) type control for commodity-I with a random lead time but instantaneous replenishment for commodity-II. If the waiting hall reaches its maximum capacity of *N*, any arriving primary client may enter an infinite capacity orbit with a specified ratio. For orbiting consumers, the classical retrial policy is used. In a steady-state setting, the joint probability distributions for commodities and the number of demands in the queue and the orbit, are derived. From this, we derive a waiting time analysis and a variety of system performance metrics in the steady-state. Additionally, the physical properties of various performance measures are evaluated using various numerical assumptions associated with diverse stochastic behaviours.

**Keywords:** classical retrial policy; queue dependent service rate; waiting time analysis; infinite orbit

**MSC:** 60K20; 60K25

#### **1. Introduction**

The function or usage of one product may be dependent on another product in general. The first product is the major commodity and the second one is the complimentary of the first product, such as mobile phone with memory card, bike with helmet, printer with ink cartridge, computer with software, torch with battery, etc. From the production point of view, both commodities are correlated with each other. According to the demand of both commodities, a firm will sell them abundantly to the targeted population who are economically benefited with the purchase.

Elaborately if the cost of one product increases, the customer demand for the corresponding complimentary product decreases. So the customers' interest towards the product will be changed. Furthermore, it will spoil the existing quantity of the product and hence, the company may encounter lose in their sales. In order to maintain the goodwill of business, the company should sell their product along with its complimentary product. Furthermore any company must plan to introduce a new product as a compliment of

**Citation:** Jeganathan, K.; Reiyas, M.A.; Selvakumar, S.; Anbazhagan, N.; Amutha, S.; Joshi, G.P.; Jeon, D.; Seo, C. Markovian Demands on Two Commodity Inventory System with Queue-Dependent Services and an Optional Retrial Facility. *Mathematics* **2022**, *10*, 2046. https://doi.org/ 10.3390/math10122046

Academic Editor: Andreas C. Georgiou

Received: 24 April 2022 Accepted: 8 June 2022 Published: 13 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

another product that will gain customers interests towards the new product with positive feedback. It can be applied in any industry right from software to dairy products. These impacts make us analyze the economic strategy of multi inventory system which sells both products with an affordable cost.

In reality, many servers can adjust the speed of service according to customer satisfaction which impacts brand reputation. In a queuing system, some authors considered queue length-dependent service times (see Abolnikov [1], Dshalalow [2], Fakinos [3], Harris [4] and Ivnitskiy et al. [5]). Furthermore, if an arriving customer finds the waiting room full, the customer decides to reattempt to get the entry with certain proportion. Under such real conditions, the proposed model deals with two commodity perishable inventory system that is the novelty of this paper. The next section elaborates the Review Literature of this model. System description is presented in Section 3. Under stability conditions, the mathematical and waiting time of the model is studied in Sections 4 and 5, respectively. System characteristics and numerical analysis are presented in Sections 6 and 7, respectively. Conclusion of the deal is presented in the final section.

#### *Related Works*

The relevant features of the present study are discussed in various modeling of queuing inventory system. Before making procurement of some products in an inventory system, a customer should need some demonstration about the product. It will take some positive time for a server showing its demonstration. In that duration, any arriving customer may wait for his/her turn. An inventory system with arbitrary service time was first studied by Sigman and Simchi-levi [6]. Schwarz et al. [7] first developed an inventory queuing system for which lead time, service time and the time between each arrival into the system are all considered to be independent exponential distributions. This was discussed by Berman et al. [8]. Recently, regarding service facilities one can refer [9].

Gebhard [10] considered a M/M/1 queuing system with two rate of service policy. Sangeetha and Sivakumar [11] studied a perishable inventory system for MAP arrival and Phase type service distribution environment and found an optimal policies of service rates in order to minimize the expected total cost. Under the steady state approximation, the probability of empty state and mean queue length are also obtained. Furthermore Doo Il Choi et al. [12] individually studied the queue dependent service time with finite capacity and infinite capacity queue. In these models, one rate of service is fixed up to a certain level in the queue and another rate of service is provided after finding the queue size beyond that level. Jeganathan et al. [13] considered an inventory queuing system which provides k independent phases of services with each service rate depending upon two discrete states of queue length.

Keerthana et al. [14] considered the postponed and renewal demands on their inventory system in which they follow arbitrary demand distribution for the inter-arrival times. The economical advantages of an inventory system are analyzed with three kinds of retrial mode by Krishnamoorthy and Jose [15]. Amirthakodi and Sivakumar [16] investigated the feedback mechanism on their inventory system with an orbital search policy. Dhanya Shajin and Krishnamoorthy [17] explored the advanced reservation, cancellation, overbooking with an impatient customers. Furthermore, more queries regarding inventory and production inventory discussion, can be referred through [18–21].

Paul Manuel et al. [22] studied a perishable finite capacity retrial inventory system with service facilities. The inter-retrial time and life time of a stored item for each model discussed here are exponentially distributed. Furthermore, Paul Manuel et al. [23] counted a constant rate of retrial from infinite space orbit is dealt with a finite queuing perishable inventory. Kathiresan et al. [24] considered an inventory system with finite capacity waiting hall. In this model,any arriving customer can also use a finite capacity orbit whenever there is no vacant in waiting hall and there is a vacant in the orbit with constant rate of retrial. Most of the times, the nature of the retrial policy is dependent on the size of the customers in the orbit which is discussed as a classical retrial policy (CRP). Artaljeo et al. [25] and Ushakumari [26] both elaborately discussed CRP in an inventory system.

The merits of multi item service facilities are pointed out in some publications. The optimum ordering policy of multi-item inventory system under some different constraints of total cost function was studied by Veinott and Wagner [27] and Wagner et al. [28]. Alscher and Schneider [29] determined a multi-item inventory control model undertaking with varying costs of multi items. Kalpakam and Arivarignan [30] considered a joint (s,S) reordering policy of a multi-item inventory. In this model, S denotes the aggregate of maximum stock level of each commodity and whenever the aggregate stock level reaches to s, the quantity of each commodity are ordered up to its maximum stock level and the replenishment time of any joint new order is zero. Anbazhagan and Arivarignan [31–33] analyzed the different ordering policies of two commodities inventory system under the assumption of various random conditions. Sivakumar et al. [34] studied two commodity inventory queuing system where any one of two item is replaceable when a demanded item is not available. Instead of a finite waiting hall, Sivakumar [35] studied the same with retrial demand.

Yadavalli et al. [36] studied a two commodity coordinated inventory system in which demand of each commodity is sold with another commodity under a distinct Bernoulli schedule and arrival pattern of a customer is poisson. Furthermore, Yadavalli et al. [36] thoroughly investigated two perishable commodity inventory system with three types of customers where the demand of one commodity is always bulk. Anbazhagan et al. [37] introduced the gift item for a customer whenever the quantity of the demand of the main product is beyond certain level. Under a base stock ordering policy, a two commodity inventory system with a finite capacity waiting hall was analyzed by Gomathi et al. [38]. Anbazhagan and Jeganathan [39] independently studied the ordering policies of primary product and gift item, but compliment item may also be sold when a customer do not make any demand of primary product.

#### **2. System Description**

A single server two commodity inventory system with a limited queue of size *N* and an optional retrial facility of indefinite size is considered in this model, in which one commodity is designated as a significant item and the other as a complimentary item is taken into consideration. When a customer first enters the system, they purchase a large item and depart with a complementary item after the service is completed. The interval between the arrival times of any two customers is exponentially distributed with rate *λ* between them. The single server delivers queue-dependent service at a rate of *xμ*, where *x* denotes the number of consumers currently waiting in line at the time of the request. If there is no available space in the queue, any new customers who join in an endless retrial orbit with a rate of *qλ*. After some exponential time has passed, the customer can be reintroduced into the system to satisfy his or her demand, with the rate of reintroduction being *uλr*, where *u* represents the number of orbital clients present at the time.

For each commodity, a continuous review ordering policy is considered: (*s*, *Q*) and (0, *S*2)(instantaneous) are the ordering policies for the major item and complimentary item, respectively; further, anytime the level of stock for a major item falls below *s*, the system immediately places an order for quantity *Q* of the major item, and the time it takes for the ordered quantity of major item to arrive is exponential with rate *β*, as shown in the following diagram. However, whenever the stock of the complimentary item reaches zero, the order for quantity *S*<sup>2</sup> is placed immediately. The deteriorating times of both items are random and exponential in nature, with intensities of *γ*<sup>1</sup> and *γ*<sup>2</sup> for the major item and *γ*<sup>2</sup> for the complimentary item, respectively.

#### **3. Analysis of the Model**

Let *X*1(*t*), *X*2(*t*), *X*3(*t*), *W*(*t*) described number of demands in the orbit, current inventory level of first commodity (CIL1), current inventory level of second commodity (CIL2), number of demands in the waiting hall (queue), respectively. Assumption developed on the birth death process make a stochastic process *Y*(*t*) = {(*X*1(*t*), *X*2(*t*), *X*3, *W*(*t*)), *t* ≥ 0} and it is also said to be a continuous-time stochastic process (CTSP) having the state

space *E* such that *E* = {(℘1, ℘2, ℘3, ℘4) : ℘<sup>1</sup> = 0, 1, ··· ; ℘<sup>2</sup> = 0, 1, ··· , *S*1; ℘<sup>3</sup> = 1, ··· , *S*2; ℘<sup>4</sup> = 0, 1, ··· , *N*}.

#### *3.1. Construction of Infinitesimal Generator Matrix*

Indicating to the discrete state space and continuous time Markov chain, the transition matrix of *Y*(*t*) having the structure as follows:

$$B\_{\mathbf{a}} = \begin{pmatrix} \mathbb{B}\_{00} & \mathbb{B}\_{01} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots \\ \mathbb{B}\_{10} & \mathbb{B}\_{11} & \mathbb{B}\_{01} & \mathbf{0} & \mathbf{0} & \cdots \\ \mathbf{0} & \mathbb{B}\_{20} & \mathbb{B}\_{21} & \mathbb{B}\_{01} & \mathbf{0} & \cdots \\ \vdots & \vdots & \ddots & \ddots & \ddots & \ddots \end{pmatrix} / \mathbf{0}$$

where

$$\mathbb{B}\_{01} = \left\{ \begin{array}{ll} q\lambda & \wp\_1' = \wp\_{1\prime} & \wp\_1 \in \{0, 1, 2, \cdot \cdot \cdot \} \\ & \wp\_2' = \wp\_{2\prime} & \wp\_2 \in \{0, 1, 2, \cdot \cdot \cdot \cdot \} \\ & \wp\_3' = \wp\_{3\prime} & \wp\_3 \in \{1, 2, \cdot \cdot \cdot \cdot \cdot \} \\ & \wp\_4' = \wp\_{4\prime} & \wp\_4 = N \\ 0, & \text{otherwise}. \end{array} \right\}$$

$$
\mathbb{B}\_{\wp\_1 0} = \begin{cases}
\wp\_1 \lambda\_{\nu \prime} & \wp\_1^{\prime} = \wp\_1 - 1, \quad \wp\_1 \in \{1, 2, \cdots, \cdot\} \\
& \wp\_2^{\prime} = \wp\_2, \quad \wp\_2 \in \{0, 1, 2, \cdots, \cdot\} \mathcal{S}\_1 \\
& \wp\_3^{\prime} = \wp\_3, \quad \wp\_3 \in \{1, 2, \cdots, \cdot\} \mathcal{S}\_2 \\
& \wp\_4^{\prime} = \wp\_4, \quad \wp\_4 \in \{0, 1, 2, \cdots, \cdot\} \mathcal{N} - 1 \\
0, & \text{otherwise}.
\end{cases}
$$

and

<sup>B</sup>℘1,℘ <sup>1</sup> <sup>=</sup> ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ *β* ℘ <sup>1</sup> = ℘1, ℘<sup>1</sup> ∈ {0, 1, 2, . . .} ℘ <sup>2</sup> = ℘<sup>2</sup> + *Q*, ℘<sup>2</sup> ∈ {0, 1, 2, . . . ,*s*} ℘ <sup>3</sup> = ℘3, ℘<sup>3</sup> ∈ {1, 2, . . . , *S*2} ℘ <sup>4</sup> = ℘4, ℘<sup>4</sup> ∈ {0, 1, . . . , *N*} *λ* ℘ <sup>1</sup> = ℘1, ℘<sup>1</sup> ∈ {0, 1, 2, . . .} ℘ <sup>2</sup> = ℘2, ℘<sup>2</sup> ∈ {0, 1, 2, . . . , *S*1} ℘ <sup>3</sup> = ℘3, ℘<sup>3</sup> ∈ {1, 2, . . . , *S*2} ℘ <sup>4</sup> = ℘4, ℘<sup>4</sup> ∈ {0, 1, . . . , *N* − 1} ℘2*γ*<sup>1</sup> ℘ <sup>1</sup> = ℘1, ℘∈{0, 1, 2, . . .} ℘ <sup>2</sup> = ℘<sup>2</sup> − 1, ℘<sup>2</sup> ∈ {1, 2, . . . , *S*1} ℘ <sup>3</sup> = ℘3, ℘<sup>3</sup> ∈ {1, 2, . . . , *S*2} ℘ <sup>4</sup> = ℘4, ℘<sup>4</sup> ∈ {0, 1, . . . , *N*} ℘3*γ*<sup>2</sup> ℘ <sup>1</sup> = ℘1, ℘<sup>1</sup> ∈ {0, 1, 2, . . .} ℘ <sup>2</sup> = ℘2, ℘<sup>2</sup> ∈ {0, 1, 2, . . . , *S*<sup>1</sup> − 1} ℘ <sup>3</sup> = *S*2, ℘<sup>3</sup> = 1 ℘ <sup>3</sup> = ℘<sup>3</sup> − 1, ℘<sup>3</sup> ∈ {2, 3, . . . , *S*2} ℘ <sup>4</sup> = ℘4, ℘<sup>4</sup> ∈ {0, 1, . . . , *N*} ℘4*μ* ℘ <sup>1</sup> = ℘1, ℘<sup>1</sup> ∈ {0, 1, 2, . . .} ℘ <sup>2</sup> = ℘<sup>2</sup> − 1, ℘<sup>2</sup> ∈ {1, 2, . . . , *S*1} ℘ <sup>3</sup> = *S*2, ℘<sup>3</sup> = 1 ℘ <sup>3</sup> = ℘<sup>3</sup> − 1, ℘<sup>3</sup> ∈ {2, 3, . . . , *S*2} ℘ <sup>4</sup> = ℘<sup>4</sup> − 1, ℘<sup>4</sup> ∈ {1, 2, . . . , *N*} −[*H*(*s* − ℘2)*β* + ℘2*γ*<sup>1</sup> ℘ <sup>1</sup> = ℘1, ℘<sup>1</sup> ∈ {0, 1, 2, . . .} +℘3*γ*<sup>2</sup> + *v*<sup>2</sup> ¯ *<sup>δ</sup>*0℘2℘4*<sup>μ</sup>* ℘ <sup>2</sup> = ℘2, ℘<sup>2</sup> ∈ {0, 1, 2, . . . , *S*1} *<sup>δ</sup>N*℘<sup>4</sup> *<sup>q</sup><sup>λ</sup>* + ¯ *<sup>δ</sup>N*℘4℘1*λr*] ℘ <sup>3</sup> = ℘3, ℘<sup>3</sup> ∈ {1, 2, . . . , *S*2} ℘ <sup>4</sup> = ℘4, ℘<sup>4</sup> ∈ {0, 1, . . . , *N*} 0, otherwise.

The infinitesimal generator matrix, *B*, is to be obtained by the transitions as follows:


### *3.2. Matrix Geometric Approximation*

#### Steady State Analysis

Consider the point at which this truncation procedure stops for the matrix-geometric approximation to be *K*. In order to identify the steady state of the considered system using Neut's Rao truncation approach, we make the assumptions that *Bi*<sup>0</sup> = *BK*<sup>0</sup> and *Bi*<sup>1</sup> = *BK*<sup>1</sup> for all *i* ≥ *K*. In addition, the updated generator matrix for the truncated system with the following structure is created.

$$
\hat{B}\_{\;i} = \begin{pmatrix}
\mathbb{B}\_{00} & \mathbb{B}\_{01} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots \\
\mathbb{B}\_{10} & \mathbb{B}\_{11} & \mathbb{B}\_{01} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots \\
\mathbf{0} & \mathbb{B}\_{20} & \mathbb{B}\_{21} & \mathbb{B}\_{01} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots \\
\vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \\
\mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbb{B}\_{K0} & \mathbb{B}\_{K1} & \mathbb{B}\_{01} & \mathbf{0} & \mathbf{0} & \cdots \\
\mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbb{B}\_{K0} & \mathbb{B}\_{K1} & \mathbb{B}\_{01} & \mathbf{0} & \cdots \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \ddots \\
\end{pmatrix},
$$

**Theorem 1.** *The steady-state probability vector, χ, where*

$$\begin{array}{c} \chi = (\chi^{(0)}, \chi^{(1)}, \dots, \chi^{(\mathfrak{S}\_1)}),\\ \chi^{(\wp\_2)} = \chi^{(\wp\_2,0)}, \chi^{(\wp\_2,1)}, \dots, \chi^{(\wp\_2,\mathfrak{S}\_2)} \wp\_2 = 0, 1, \dots, S\_1 \end{array}$$

$$\chi^{(\wp\_2,\wp\_3)} = \chi^{(\wp\_2,\wp\_3,0)}, \chi^{(\wp\_2,\wp\_3,1)}, \dots, \chi^{(\wp\_2,\wp\_3,N)} \wp\_2 = 0, 1, \dots, S\_1 \\ \text{and} \\ \wp\_3 = 0, 1, \dots, S\_2.$$

*that corresponds to the generator matrix is denoted by* B*<sup>K</sup> where* B*<sup>K</sup>* = B*K*<sup>0</sup> + B*K*<sup>1</sup> + B<sup>01</sup> *is given by*

$$\chi^{(i)}\;=\;\chi^{(Q)}\;e\_{i\prime}\quad i=0,1,\ldots,\mathrm{S}\_1\dots$$

*where*

$$e\_i = \begin{cases} (-1)^{Q-i} F\_Q E\_{Q-1}^{-1} F\_{Q-1} \cdots F\_{i+1} E\_i^{-1}, & i = 0, 1, \cdots, Q-1 \\ I\_\prime & i = Q \\ (-1)^{2Q-i+1} \sum\_{j=0}^{S\_1 - i} [(F\_Q E\_{Q-1}^{-1} F\_{Q-1} \cdots F\_{s+1-j} E\_{s-j}^{-1}) \times \\ \text{ } G E\_{S\_1 - j}^{-1} (F\_{S\_1 - j} E\_{S\_1 - j-1}^{-1} F\_{S\_1 - j-1} \cdots F\_{i+1} E\_i^{-1})], & i = Q+1, Q+2, \cdots, S\_1 \end{cases}$$

*and χ*(*Q*) *is obtained by solving*

$$\begin{aligned} &\chi^{(Q)}[(-1)^{Q}\sum\_{j=0}^{s-1}[(F\_{Q}E\_{Q-1}^{-1}F\_{Q-1}\cdots F\_{s+1-j}E\_{s-j}^{-1})GE\_{S\_1-j}^{-1}(F\_{S\_1-j}E\_{S\_1-j-1}^{-1}F\_{S\_1-j-1}\cdots F\_{i+1}E\_i^{-1})] \\ &F\_{Q+1}+E\_Q+(-1)^{Q}F\_QE\_{Q-1}^{-1}F\_{Q-1}\cdots F\_1E\_0^{-1}G) = \mathbf{0} \\ &\text{and }\sum\_{j=1}^{S\_1}\chi^{(i)}\mathbf{e}=1. \end{aligned}$$

**Proof.** We have

$$
\chi \mathbb{B}\_K = \mathbf{0} \text{ and } \chi \mathbf{e} = 1 \tag{1}
$$

where

$$\begin{array}{rcl} [\mathbb{B}\_K]\_{ij} &=& \begin{cases} \begin{array}{rcl} E\_i & j = i, & \quad i = 0, 1, 2, \cdots, \text{'}, S\_1; \\ F\_i & j = i - 1, & \quad i = 1, 2, \cdots, \text{'}, S\_1; \\ G & j = i + Q, & i = 0, 1, 2, \cdots, \text{'}, \text{'}; \\ \mathbf{0} & \text{otherwise} \end{array} \end{array} \end{array}$$

The first equation of the above framework yields the following set of equations:

$$\begin{array}{lll}\chi^{(i+1)}E\_{i+1} + \chi^{(i)}E\_{i} &=& \mathbf{0}, i = 0, 1, \cdots, Q-1\\\chi^{(i+1)}E\_{i+1} + \chi^{(i)}E\_{i} + \chi^{(i-Q)}G &=& \mathbf{0}, i = Q, Q+1, \cdots, S\_{1}-1\\\chi^{(i)}E\_{i} + \chi^{(i-Q)}G &=& \mathbf{0}, i = S\_{1}\end{array}$$

Solving the system of equations, we will get the stated result.

**Theorem 2.** *The stability condition of the system at the truncation point K is given by*

$$r\_1 q \lambda < r\_2 K \lambda\_I$$
  $where \ r\_1 = \sum\_{\wp\_2=0}^{S\_1} \sum\_{\wp\_3=1}^{S\_2} \chi^{(\wp\_2, \wp\_3, N)} \text{ and } r\_2 = \sum\_{\wp\_2=0}^{S\_1} \sum\_{\wp\_3=1}^{S\_2} \sum\_{\wp\_4=1}^{N} \chi^{(\wp\_2, \wp\_3, \wp\_3)} \text{.}$ 

**Proof.** From the well known result of Neuts on the positive recurrence of B*<sup>K</sup>* we have

$$
\chi^{(K)}\mathbb{B}\_{01}\mathbf{e} < \chi^{(K)}\mathbb{B}\_{K0}\mathbf{e},
$$

and by exploiting the structure of the matrices B<sup>01</sup> and B*K*0, we get, for ℘<sup>2</sup> = 0, 1, 2, ··· , *<sup>S</sup>*1, ℘<sup>3</sup> = 1, 2, ··· , *S*<sup>2</sup> and ℘<sup>4</sup> = 0, 1, 2, ··· , *N*.

$$
\chi^{(\wp\_{2\text{-}\ell\triangleright 3\text{-}\ell\text{-}4})} \mathbb{B}\_{01} \mathbf{e} < \chi^{(\wp\_{2\text{-}\ell\triangleright 3\text{-}\ell\text{-}4)}} \mathbb{B}\_{K0} \mathbf{e}.
$$

First, [*χ*(0), *<sup>χ</sup>*(1), ··· , *<sup>χ</sup>*(*S*1)]B01**<sup>e</sup>** < [*χ*(0), *<sup>χ</sup>*(1), ··· , *<sup>χ</sup>*(*S*1)]B*K*0**<sup>e</sup>** Due to the structure of *B*0, L.H.S becomes

$$
\chi^{(\wp\_2)} B\_0 = \chi^{(\wp\_2, \wp\_3, N)} q \lambda\_\dots
$$

On the other hand, due to structure of *B*<sup>1</sup> R.H.S becomes

$$\chi^{(\wp\_2)}B\_1 = [\chi^{(\wp\_{2\prime\wp\_3\prime 1})}, \chi^{(\wp\_{2\prime\wp\_3 2})}, \dots, \chi^{(\wp\_{2\prime\wp\_3\prime 3})}] \mathbb{K}\lambda\_r$$

Therefore, the last inequality becomes

$$\sum\_{\wp\_2=0}^{S\_1} \sum\_{\wp\_3=1}^{S\_2} \chi^{(\wp\_2, \wp\_3, N)} q \lambda < \sum\_{\wp\_2=0}^{S\_1} \sum\_{\wp\_3=1}^{S\_2} \sum\_{\wp\_4=1}^{N} \chi^{(\wp\_2, \wp\_3, \wp\_4)} K \lambda\_{\varGamma}$$

Hence,

$$r\_1 q \lambda < r\_2 K \lambda\_r$$

where *<sup>r</sup>*<sup>1</sup> <sup>=</sup> *<sup>S</sup>*<sup>1</sup> ∑ ℘2=0 *S*2 ∑ ℘3=1 *<sup>χ</sup>*(℘2,℘3,*N*) and *<sup>r</sup>*<sup>2</sup> <sup>=</sup> *<sup>S</sup>*<sup>1</sup> ∑ ℘2=0 *S*2 ∑ ℘3=1 *N* ∑ ℘4=1 *χ*(℘2,℘3,℘4) as desired.

**Remark 1.** *Using r*1*qλ* < *r*2*Kλr, we get,*

$$\frac{r\_1 q \lambda}{r\_2 K \lambda\_r} < 1$$

 $\text{If } \frac{r\_1}{r\_2 K} = m(say), \text{ then }$  
$$\frac{\lambda}{\lambda\_r} < \frac{1}{m}.$$

#### *3.3. Limiting Probability Distribution*

It can be seen from the structure of the rate matrix *B* and from the Theorem 3, that the Markov process {*X*1(*t*), *X*2(*t*), *X*3(*t*), *W*(*t*), *t* ≥ 0} with the state space *E* is regular. Henceforth, the limiting probability distribution

*Υ*(℘1,℘2,℘3,℘4) = lim *<sup>t</sup>*→<sup>∞</sup> *Pr*[*X*1(*t*) = <sup>℘</sup>1, *<sup>X</sup>*2(*t*) = <sup>℘</sup>2, *<sup>X</sup>*3(*t*) = <sup>℘</sup>3, *<sup>W</sup>*(*t*) = <sup>℘</sup>4|*X*1(0), *<sup>X</sup>*2(0), *<sup>X</sup>*3(0), *<sup>W</sup>*(0)],

exists and is independent of the initial state.

Let *Υ* = *<sup>Υ</sup>*(0),*Υ*(1),..., satisfies

*ΥB* = **0**, *Υ***e** = 1.

We can partition the vector *Υ*(℘1), as

$$\mathcal{Y}^{(\wp\_1)} = \left( \mathcal{Y}^{(\wp\_1,0)}, \mathcal{Y}^{(\wp\_1,1)}, \dots, \mathcal{Y}^{(\wp\_1,S\_1)} \right), \wp\_1 \ge 0$$

and

$$Y^{(\wp\_1,\wp\_2)} = \left( Y^{(\wp\_1,\wp\_2,1)}, Y^{(\wp\_1,\wp\_2,2)}, \dots, Y^{(\wp\_1,\wp\_2,S\_2)} \right), \wp\_1 \ge 0, 0 \le \wp\_2 \le S\_1$$

$$Y^{\left(\wp\_1,\wp\_2,\wp\_3\right)} = \left(Y^{\left(\wp\_1,\wp\_2,\wp\_3,0\right)}, Y^{\left(\wp\_1,\wp\_2,\wp\_3,1\right)},\ldots,Y^{\left(\wp\_1,\wp\_2,\wp\_3,N\right)}\right), \wp\_1 \ge 0, 0 \le \wp\_2 \le S\_1, 0 \le \wp\_3 \le N$$

**Theorem 3.** *Utilizing the vector Υ* = (*Υ*(0),*Υ*(1), ... ,) *and the specific structure of B, R can be determined by*

$$R^2 \mathbb{B}\_{K0} + R\mathbb{B}\_{K1} + \mathbb{B}\_{01} = \mathbf{0}$$

*where R is the minimal non-negative of the matrix quadratic equation.*

**Proof.** Assume

$$R = \begin{pmatrix} D\_{00} & D\_{01} & \cdots & D\_{0S\_1} \\ D\_{10} & D\_{11} & \cdots & D\_{1S\_1} \\ D\_{20} & D\_{21} & \cdots & D\_{2S\_1} \\ \vdots & \vdots & \ddots & \vdots \\ D\_{S\_10} & D\_{S\_11} & \cdots & D\_{S\_1S\_1} \end{pmatrix}$$

where

$$D\_{ij} = \begin{pmatrix} R\_{11} & R\_{12} & \cdots & R\_{1S\_2} \\ R\_{21} & R\_{22} & \cdots & R\_{2S\_2} \\ R\_{31} & R\_{32} & \cdots & R\_{3S\_2} \\ \vdots & \vdots & \ddots & \vdots \\ R\_{S\_21} & R\_{S\_22} & \cdots & R\_{S\_2S\_2} \end{pmatrix}, i, j \in 1, 2, \cdots, S\_2$$

Since the block matrix B<sup>01</sup> has (*S*<sup>1</sup> + <sup>1</sup>)(*S*2) number of nonzero rows, the assumed *<sup>R</sup>* matrix also has the same number of nonzero rows. Now, due to the specified structure of B01, the structure of the block matrix *Rmn* is of the form

$$R\_{mn} = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & \cdots & 0 \\ 0 & 0 & 0 & 0 & 0 & \cdots & 0 \\ 0 & 0 & 0 & 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ l\_{mn}^0 & l\_{mn}^1 & l\_{mn}^2 & l\_{mn}^3 & l\_{mn}^4 & \cdots & l\_{mn}^N \end{pmatrix}, m, n \in 1, 2, \cdots, S\_2$$

is also square matrix of dimension *N* + 1.

Now, exploiting the coefficient matrices B*K*0,B*K*1,B<sup>01</sup> with *<sup>R</sup>*<sup>2</sup> and *<sup>R</sup>* equating with **<sup>0</sup>**, we obtain a system of (*N* + 1)-dimensional vector as follows:

For *<sup>i</sup>* <sup>=</sup> 0, 1, ··· , *<sup>S</sup>*1, *<sup>j</sup>* <sup>=</sup> 0, *<sup>m</sup>* <sup>=</sup> 1, 2, ··· , *<sup>S</sup>*2, *<sup>n</sup>* <sup>=</sup> 1, 2, ··· , *<sup>S</sup>*<sup>2</sup> <sup>−</sup> 1, *<sup>C</sup>*(*j*) *<sup>h</sup>* be the diagonal elements of the corresponding matrix *Ej* where *j* = 0, 1, 2, ··· , *S*<sup>1</sup> and *h* = 1, 2, ··· , *S*2(*N* + 1).

(*l* 0 *mnC*(*j*) (*n*−1)(*N*+1)+<sup>1</sup> + *l* 0 *m*(*n*+1) (*n* + 1)*γ*<sup>2</sup> + *l* 0 *mn*(*j* + 1)*γ*<sup>1</sup> + *l* 1 *m*(*n*+1) 1*μ*, *l* 0 *mnλ* + *l* 1 *mnC*(*j*) (*n*−1)(*N*+1)+2+ *l* 1 *m*(*n*+1) (*n* + 1)*γ*<sup>2</sup> + *l* 1 *mn*(*j* + 1)*γ*<sup>1</sup> + *l* 2 *m*(*n*+1) 2*μ*, ··· , *l <sup>N</sup>*−<sup>2</sup> *mn λ* + *l <sup>N</sup>*−<sup>1</sup> *mn <sup>C</sup>*(*j*) (*n*−1)(*N*+1)+*<sup>N</sup>* + *l N*−1 *m*(*n*+1) (*n*+ 1)*γ*<sup>2</sup> + *l <sup>N</sup>*−<sup>1</sup> *mn* (*j* + 1)*γ*<sup>1</sup> + *l N <sup>m</sup>*(*n*+1)*Nμ*, *l <sup>N</sup>*−<sup>1</sup> *mn λ* + *l N mnC*(*j*) *<sup>n</sup>*(*N*+1) + *l N m*(*n*+1) (*n* + 1)*γ*<sup>2</sup> + *l N mn*(*j* + 1)*γ*1+ *δi*0*qλ*) = **0**

$$\text{For } j = 0 \text{, } m = 1 \text{, } 2 \text{, } \cdot \cdot \text{, } \text{S}\_2 \text{, } n = \text{S}\_2$$

(*l* 0 *mnC*(*n*−1)(*N*+1)+<sup>1</sup> + *l* 0 *<sup>m</sup>*11*γ*<sup>2</sup> + *l* 0 *mn*(*j* + 1)*γ*<sup>1</sup> + *l* 1 *<sup>m</sup>*11*μ*, *l* 0 *mnλ* + *l* 1 *mnC*(*n*−1)(*N*+1)+<sup>2</sup> + *l* 1 *<sup>m</sup>*11*γ*<sup>2</sup> + *l* 1 *mn*(*j*+ 1)*γ*<sup>1</sup> + *l* 2 *<sup>m</sup>*12*μ*, ··· , *l <sup>N</sup>*−<sup>2</sup> *mn λ* + *l <sup>N</sup>*−<sup>1</sup> *mn <sup>C</sup>*(*n*−1)(*N*+1)+*<sup>N</sup>* + *<sup>l</sup> <sup>N</sup>*−<sup>1</sup> *<sup>m</sup>*<sup>1</sup> <sup>1</sup>*γ*<sup>2</sup> <sup>+</sup> *<sup>l</sup> <sup>N</sup>*−<sup>1</sup> *mn* (*j* + 1)*γ*<sup>1</sup> + *l N <sup>m</sup>*1*Nμ*, *l <sup>N</sup>*−<sup>1</sup> *mn λ*+ *l N mnCn*(*N*+1) + *l N <sup>m</sup>*11*γ*<sup>2</sup> + *l N mn*(*j* + 1)*γ*<sup>1</sup> + *δi*0*qλ*) = **0**

For *j* = 1, 2, ··· , *S*<sup>1</sup> − 1, *m* = 1, 2, ··· , *S*2, *n* = 1, 2, ··· , *S*<sup>2</sup> − 1 (*l* 0 *mnC*(*j*) (*n*−1)(*N*+1)+<sup>1</sup> + *l* 0 *m*(*n*+1) (*n* + 1)*γ*<sup>2</sup> + *l* 0 *mn*(*j* + 1)*γ*<sup>1</sup> + *l* 1 *m*(*n*+1) 1*μ* + *H*(*s* − *j*)*l* 0 *mnβ S*2 ∑ *d*=1 *l* 0 *mdl* 0 *dnNλ<sup>r</sup>* + *l* 0 *mnλ* + *l* 1 *mnC*(*j*) (*n*−1)(*N*+1)+<sup>2</sup> + *l* 1 *m*(*n*+1) (*n* + 1)*γ*<sup>2</sup> + *l* 1 *mn*(*j* + 1)*γ*<sup>1</sup> + *l* 2 *m*(*n*+1) 2*μ*+ *H*(*s* − *j*)*l* 1 *mnβ*, ··· , *S*2 ∑ *d*=1 *l N*−2 *md ldn<sup>N</sup>*−2*Nλ<sup>r</sup>* <sup>+</sup> *<sup>l</sup> <sup>N</sup>*−<sup>2</sup> *mn λ* + *l <sup>N</sup>*−<sup>1</sup> *mn <sup>C</sup>*(*j*) (*n*−1)(*N*+1)+*<sup>N</sup>* + *l N*−1 *m*(*n*+1) (*n* + 1)*γ*2+ *l <sup>N</sup>*−<sup>1</sup> *mn* (*j* + 1)*γ*<sup>1</sup> + *l N <sup>m</sup>*(*n*+1)*Nμ*+ *H*(*s* − *j*)*l <sup>N</sup>*−<sup>1</sup> *mn β*, *S*2 ∑ *d*=1 *l N*−1 *md l N*−1 *dn Nλ<sup>r</sup>* + *l <sup>N</sup>*−<sup>1</sup> *mn λ* + *l N mnC*(*j*) *<sup>n</sup>*(*N*+1) + *l N m*(*n*+1) (*n* + 1)*γ*<sup>2</sup> + *l N mn*(*j* + 1)*γ*<sup>1</sup> +*δijqλ* + *H*(*s* − *j*)*l N mnβ*) = **0** For *j* = 1, 2, ··· , *S*<sup>1</sup> − 1, *m* = 1, 2, ··· , *S*2, *n* = *S*<sup>2</sup>

$$\begin{aligned} & (l\_{mn}^{0}\mathbb{C}\_{(n-1)(N+1)+1} + l\_{m1}^{0}\mathbbm{1}\gamma\_{2} + l\_{mn}^{0}(j+1)\gamma\_{1} + l\_{m1}^{1}\mathbbm{1}\mu + H(s-j)l\_{mn}^{0}\mathscr{S}\_{r}\sum\_{d=1}^{S\_{2}}l\_{md}^{0}M\lambda\_{r} + l\_{mn}^{0}\lambda\_{r} + l\_{m1}^{0}\mathscr{S}\_{r}\mathscr{S}\_{d}\mathscr{S}\_{r} \\ & l\_{mn}^{1}\mathscr{C}\_{(n-1)(N+1)+2}^{(j)} + l\_{m1}^{1}\mathbbm{1}\gamma\_{2} + l\_{mn}^{1}(j+1)\gamma\_{1} + l\_{m1}^{2}\mathscr{A}\mu + H(s-j)l\_{mn}^{1}\mathscr{S}\_{r}\mathscr{S}\_{r} \cdot \sum\_{d=1}^{S\_{2}}l\_{md}^{N-2}l\_{dn}^{N-2}N\lambda\_{r} + l\_{m1}^{0}\mathscr{S}\_{r}\mathscr{S}\_{d}\mathscr{S}\_{r} \\ & l\_{mn}^{N-2}\lambda + l\_{mn}^{N-1}\mathscr{C}\_{(n-1)(N+1)+N} + l\_{m1}^{N-1}\mathscr{I}\gamma\_{2} + l\_{mn}^{N-1}(j+1)\gamma\_{1} + l\_{m1}^{N}N\mu + \end{aligned}$$

$$\begin{cases} H(s-j)l\_{mn}^{N-1}\mathcal{B}\_{\prime}\sum\_{d=1}^{S\_2}l\_{md}^{N-1}l\_{dn}^{N-1} \\ N\lambda\_r + l\_{mn}^{N-1}\lambda + l\_{mn}^N\mathbb{C}\_{n(N+1)} + l\_{m1}^N1\gamma\_2 + l\_{mn}^N(j+1)\gamma\_1 + \delta\_{lj}q\lambda + H(s-j)l\_{mn}^N\mathcal{B}) = \mathbf{0} \end{cases}$$
 
$$\text{For } j = S\_1, m = 1, 2, \dots, \dots, S\_2, n = 1, 2, \dots, S\_2 - 1$$

$$\begin{split} & (l\_{mn}^{0}\mathbf{C}\_{(n-1)(N+1)+1}^{(j)} + l\_{m(n+1)}^{0}(n+1)\gamma\_{2} + H(s-j)l\_{m\eta}^{0}\mathbf{\hat{\beta}}\_{\cdot} \sum\_{d=1}^{S\_{2}} l\_{md}^{0}l\_{dn}^{0}N\lambda\_{r} + \\ & l\_{mn}^{0}\lambda\_{r} + l\_{mn}^{1}\mathbf{C}\_{(n-1)(N+1)+2}^{(j)} + l\_{m(n+1)}^{1}(n+1)\gamma\_{2} + H(s-j)l\_{m\eta}^{1}\mathbf{\hat{\beta}}\_{\cdot} \dots \sum\_{d=1}^{S\_{2}} l\_{md}^{N-2}l\_{dn}^{N-2}N\lambda\_{r} + l\_{mn}^{N-2}\lambda + l\_{mn}^{N} \\ & l\_{mn}^{N-1}\mathbf{C}\_{(n-1)(N+1)+N}^{(j)} + \\ & l\_{m(n+1)}^{N-1}(n+1)\gamma\_{2} + H(s-j)l\_{mn}^{N-1}\beta\_{2} \sum\_{d=1}^{S\_{2}} l\_{md}^{N-1}l\_{dn}^{N-1}N\lambda\_{r} + l\_{mn}^{N-1}\lambda + l\_{mn}^{N}\mathbf{C}\_{n(N+1)}^{(j)} + l\_{m(n+1)}^{N}(n+1)\gamma\_{2} + M(s-j)l\_{m\eta}^{N} \\ & 1(\gamma\_{2} + \delta\_{ij}\eta\lambda + H(s-j)l\_{mn}^{N}\beta) = \mathbf{0} \\ & \text{For } j = S\_{1}, m = 1, 2, \cdots, S\_{2}, n = S\_{2} \end{split}$$

$$\begin{split} & (l\_{mn}^{0}\mathbb{C}\_{(n-1)(N+1)+1} + l\_{m1}^{0}\mathbbm{1}\gamma\_{2} + H(s-j)l\_{mn}^{0}\mathscr{S}\_{\prime}\sum\_{d=1}^{S\_{2}}l\_{md}^{0}N\lambda\_{n}\lambda\_{\prime} + l\_{mn}^{0}\lambda\_{\prime} + l\_{mn}^{1}\mathscr{C}\_{(n-1)(N+1)+2}^{(j)} + \\ & l\_{m1}^{1}\mathbbm{1}\gamma\_{2} + H(s-j)l\_{mn}^{1}\mathscr{S}\_{\prime}\cdot\cdots\sum\_{d=1}^{S\_{2}}l\_{md}^{N-2}N\lambda\_{\prime} + l\_{mn}^{N-2}\lambda + l\_{mn}^{N-1}\mathscr{C}\_{(n-1)(N+1)+N} + l\_{m1}^{N-1}\mathbbm{1}\gamma\_{2} + \\ & l(\textbf{s}-j)l\_{mn}^{N-1}\mathscr{S}\_{\prime}\cdot\sum\_{d=1}^{S\_{2}}l\_{md}^{N-1}N\lambda\_{\prime} + l\_{mn}^{N-1}\lambda + l\_{mn}^{N}\mathscr{C}\_{n(N+1)} + l\_{m1}^{N}\mathscr{I}\gamma\_{2} + \delta\_{ij}\eta\lambda + H(\textbf{s}-j)l\_{mn}^{N}\mathscr{B}\_{\prime}\textbf{0} = \textbf{0} \end{split}$$

Equating the (*N* + 1)-dimensional vector to zero vector we obtain a set of equations, after solving such equations, one can obtain the elements of *R* matrix.

**Theorem 4.** *Due to the specific structure of B the vector Υ can be determined by*

$$Y^{(i+K-1)} = Y^{(K-1)} \mathbb{R}^i; i \ge 0 \tag{2}$$

*where R is the solution of*

$$R^2 \mathbb{B}\_{K0} + R\mathbb{B}\_{K1} + \mathbb{B}\_{01} = \mathbf{0} \tag{3}$$

*and the vector Υ*(*i*), *i* ≥ 0

$$\mathcal{Y}^{(i)} = \begin{cases} \sigma X^{(0)} \prod\_{j=i+1}^{K} \mathbb{B}\_{j0}(-\mathbb{B}\_{j-1}), 0 \le i \le K-1 \\\ \sigma X^{(0)} R^{(i-K)}, i \ge K \end{cases} \tag{4}$$

*where*

$$\sigma = \left[1 + X^{(0)} \sum\_{i=0}^{K-1} \prod\_{j=i+1}^{K} \mathbb{B}\_{j0} (-\mathbb{B}\_{j-1}) \mathbf{e} \right]^{-1}. \tag{5}$$

**Proof.** The sub vector (*Υ*(0),*Υ*(1),...,*Υ*(*K*−1)) and the block partitioned matrix of *B* satisfies the following set of equations

$$\begin{array}{rcl} \mathcal{Y}^{(0)}\mathbb{B}\_{00} + \mathcal{Y}^{(1)}\mathbb{B}\_{10} &=& \mathbf{0} \\ \mathcal{Y}^{(i-1)}\mathbb{B}\_{01} + \mathcal{Y}^{(i)}\mathbb{B}\_{i1} + \mathcal{Y}^{(i+1)}\mathbb{B}\_{(i+1)0} &=& \mathbf{0}; 1 \le i \le K-1 \\ \mathcal{Y}^{(K-2)}\mathbb{B}\_{01} + \mathcal{Y}^{(K-1)}(\mathbb{B}\_{(K-1)1} + R\mathbb{B}\_2) &=& \mathbf{0}; \end{array} \tag{6}$$

using Equation (6), we get,

$$Y^{(i)} = \,^t Y^{(i+1)} \mathbb{B}\_{(i+1)0} (-B\_i)^{-1}, 0 \le i \le K - 1 \tag{7}$$

where

$$B\_i = \begin{cases} \mathbb{B}\_{i0}, & i = 0\\ \left(\mathbb{B}\_{i1} - \mathbb{B}\_{i0}(-B\_{i-1})^{-1}\mathbb{B}\_{01}\right), & 1 \le i \le K \end{cases}$$

Then

$$(\chi^{(K)}, \chi^{(K+1)}, \chi^{(K+2)} \dots) \begin{pmatrix} \mathbb{B}\_K & \mathbb{B}\_{01} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \dots \\ \mathbb{B}\_{K0} & \mathbb{B}\_{K1} & \mathbb{B}\_{01} & \mathbf{0} & \mathbf{0} & \dots \\ \mathbf{0} & \mathbb{B}\_{K0} & \mathbb{B}\_{K1} & \mathbb{B}\_{01} & \mathbf{0} & \dots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \end{pmatrix} = \mathbf{0} \tag{8}$$

Assume,

$$\sigma\_{\parallel} = \sum\_{i=K}^{\infty} \mathcal{Y}^{(i)} \mathbf{e} \tag{9}$$

$$X^{(i)} = \sigma^{-1} \Upsilon^{(K+i)}, i \ge 0 \tag{10}$$

where *X* is also a 4-dimensional continuous time Markov chain such as *Υ*. The similar partitions are also applicable to *X*. From (8) we get

$$\begin{array}{rcl} \mathcal{Y}^{(K)}B\_N + \mathcal{Y}^{(K+1)}\mathbb{B}\_{K0} &=& \mathbf{0} \\ \mathcal{Y}^{(K+i)} &=& \mathcal{Y}^{(K+i-1)}\mathbb{R}, i \ge 1 \end{array}$$

This can be written as

$$\begin{array}{rcl}X^{(0)}B\_N + X^{(1)}\mathbb{B}\_{01} &=& \mathbf{0} \\ X^{(i)} &=& X^{(i-1)}\mathbb{R}, i \ge 1 \end{array}$$

Since <sup>∞</sup> ∑ *i*=0 *X*(*i*)**e** = 1, then

*X*(0) (*I* − *R*)−1**e** = 1

Hence,

$$\mathcal{Y}^{(i)} = \sigma X^{(0)} R^{(i-K)}\text{,} i \ge \mathcal{K} \tag{11}$$

Again by (7) and (10),

$$\sigma^{(i)} = \sigma X^{(0)} \prod\_{j=i+1}^{K} \mathbb{B}\_{j0}(-B\_{j-1}), 0 \le i \le K - 1 \tag{12}$$

Therefore, combining (11) and (12), we get (4) and *X*(0) is the unique solution of the system

$$X^{(0)}(B\_K + R\mathbb{B}\_{K0}) = \mathbf{0} \tag{13}$$

$$X^{(0)} (I - \mathbb{R})^{-1} \mathbf{e} = 1 \tag{14}$$

Since <sup>∞</sup> ∑ *i*=0 *Υ*(*i*)**e** = 1 and (13),

$$\sigma X^{(0)} \sum\_{i=0}^{K-1} \prod\_{j=i+1}^{K} \mathbb{B}\_{j0} (-B\_{j-1}) \mathbf{e} + \sigma X^{(0)} \sum\_{K}^{\infty} R^{(i-K)} \mathbf{e} = 1$$

which gives (5).

#### **4. Waiting Time Analysis**

The waiting time of a customer is defined as the time interval between the customer enters into the waiting hall and leaves the system after the service completion. Using the Laplace–Stieltjes transform (LST), we look at the waiting time of demand in the queue as well as the orbit independently for each node. Since the state space of the proposed model is infinite, the analytical work on finding the waiting time distribution is a difficult task to do so naturally. Thus, we restrict the orbit size to be finite say *L*(*L* > *K*) to find the waiting time of a customer in the waiting hall and orbit. We denote *Wp* and *Wo* are continuous random variables to represent the waiting time distribution of a customer in the queue and orbit, respectively. The objective is to describe the probability that a customer has to wait, the distribution of the waiting time and nth order moments.

*Waiting Time of a Demand in Queue*

To enable waiting time distribution of *Wp*, we shall define some complementary variables. Suppose that the QIS is at state (℘1, ℘2, ℘3, ℘4), ℘<sup>4</sup> > 0 at an arbitrary time *t*:


$$4. \quad ^\*W\_p(\wp\_1, \wp\_2, \bar{\wp\_3}, \wp\_4)(y) = E\left[e^{yW\_p(\wp\_1, \wp\_2, \wp\_3, \wp\_4)}\right] \\ \text{LST of conditional waiting time (CWT)}.$$

**Theorem 5.** *The expected waiting time of a demand in the queue is defined by*

$$E[\mathcal{W}\_p] = \sum\_{\substack{\wp 1 = 0 \ \wp 2 = 0}}^L \sum\_{\wp 2 = 0}^{S\_1} \sum\_{\wp 3 = 0}^{S\_2} \sum\_{\wp 4 = 0}^{N-1} \mathcal{Y}^{(\wp\_1, \wp\_2, \wp\_3, \wp\_4)} E[\mathcal{W}\_p(\wp\_1, \wp\_{2\cdot}, \wp\_3, \wp\_4 + 1)] \tag{15}$$

**Proof.** To obtain the CWT, we apply first step analysis as follows:

For 0 ≤ ℘<sup>1</sup> ≤ *L*, ℘<sup>2</sup> = 0, 1 ≤ ℘<sup>3</sup> ≤ *S*2, 1 ≤ ℘<sup>4</sup> ≤ *N*

$$\begin{aligned} \,^\*W\_p(\wp\_1, 0, \wp\_3, \wp\_4)(y) &= \frac{\bar{\delta}\_{N\wp 4} \Lambda}{a} \,^\*W\_p(\wp\_1, 0, \wp\_3, \wp\_4 + 1)(y) + \\ \frac{\frac{\wp\_3 \gamma\_2}{a} \,^\*W\_p(\wp\_1, 0, \wp\_3 - 1, \wp\_4)(y) + \frac{\beta}{a} \,^\*W\_p(\wp\_1, 0, \wp\_3, \wp\_4)(y) + \\ &\qquad \frac{\bar{\delta}\_{N\wp 4} \wp\_1 \lambda\_r}{a} \,^\*W\_p(\wp\_1 - 1, 0, \wp\_3, \wp\_4 + 1)(y) + \\ &\qquad \frac{\delta\_{N\wp 4} \bar{\delta}\_{L\wp 1} q \lambda}{a} \,^\*W\_p(\wp\_1 + 1, 0, \wp\_3, \wp\_4)(y), \end{aligned} \tag{16}$$

$$\begin{array}{c} \text{where } a = (y + \bar{\delta}\_{N \wp\_4} \lambda + \wp\_3 \gamma\_2 + \beta + \bar{\delta}\_{N \wp\_4} \wp\_1 \lambda\_r + \delta\_{N \wp\_4} \bar{\delta}\_{L \wp\_1} q \lambda). \\ \text{For } 0 \le \wp\_1 \le L, 1 \le \wp\_2 \le S\_1, 1 \le \wp\_3 \le S\_2, 1 \le \wp\_4 \le N\_r \end{array}$$

$$\begin{aligned} \,^\*W\_{\mathbb{P}}(\wp\_1,\wp\_2,\wp\_3,\wp\_4)(y) &= \frac{\delta\_{N\wp\_4}\lambda}{b}^\*W\_{\mathbb{P}}(\wp\_1,\wp\_2,\wp\_3,\wp\_4+1)(y) + \\ &\quad \frac{H(s-\wp\_4)\beta}{b}^\*W\_{\mathbb{P}}(\wp\_1,\wp\_2+Q,\wp\_3,\wp\_4)(y) + \\ &\quad \frac{\not p\_2\gamma\_1\*}{b}^\*W\_p(\wp\_1,\wp\_2-1,\wp\_3,\wp\_4)(y) + \frac{w\gamma\_2\*}{b}^\*W\_p(\wp\_1,\wp\_2,\wp\_3-1,\wp\_4)(y) + \\ &\quad \frac{(\wp\_4-1)\mu^\*}{b}^\*W\_p(\wp\_1,\wp\_2-1,\wp\_3-1,\wp\_4-1)(y) + \\ &\quad \frac{\delta\_{N\wp\_4}\wp\_1\lambda\_r}{b}^\*W\_p(\wp\_1-1,\wp\_2,\wp\_3,\wp\_4+1)(y) + \frac{\delta\_{N\wp\_1}\delta\_{L\wp\_1}\eta\lambda}{b}^\*W\_p(\wp\_1+1,\wp\_2,\wp\_3,\wp\_4)(y) + \frac{\mu}{b} \end{aligned} \tag{17}$$

where, *b* = (*y* + ¯ *<sup>δ</sup>N*℘4*<sup>λ</sup>* + *<sup>H</sup>*(*<sup>s</sup>* − ℘2)*<sup>β</sup>* + ℘2*γ*<sup>1</sup> + ℘3*γ*<sup>2</sup> + ℘4*<sup>μ</sup>* + ¯ *δN*℘4℘1*λ<sup>r</sup>* + *δN*℘<sup>4</sup> ¯ *δL*℘<sup>1</sup> *qλ*). Now, we differentiate the Equations (16) and (17) for (*n* + 1) times and computing at

 $y = 0$ , we have,

$$\begin{aligned} \text{For } 0 \le \wp\_1 \le L, \wp\_2 = 0, 1 \le \wp\_3 \le S\_2, 1 \le \wp\_4 \le N, \\ E[W\_p^{n+1}(\wp\_1, 0, \wp\_3, \wp\_4)] &= \frac{\bar{\delta}\_{N\wp}\mu}{a} E[W\_p^{n+1}(\wp\_1, 0, \wp\_3, \wp\_4 + 1)] + \\ \frac{w\nu\_{\varproj}}{b} E[W\_p^{n+1}(\wp\_1, \wp\_2, \wp\_3 - 1, \wp\_4)] + \frac{\beta}{a} E[W\_p^{n+1}(\wp\_1, 0, \wp\_3, \wp\_4)] + \\ &\qquad \frac{\bar{\delta}\_{N\wp\_4}\wp\_1 \lambda\_r}{a} E[W\_p^{n+1}(\wp\_1 - 1, 0, \wp\_3, \wp\_4 + 1)] + \\ \frac{\delta\_{N\wp\_4}\bar{\delta}\_{L\wp}q\lambda}{a} E[W\_p^{n+1}(\wp\_1 + 1, 0, \wp\_3, \wp\_4)] + (n+1)E[W\_p^{n+1}(\wp\_1, \wp\_2, \wp\_3, \wp\_4)], \end{aligned} \tag{18}$$

where*a* = (*y* + ¯ *<sup>δ</sup>N*℘4*<sup>λ</sup>* + ℘3*γ*<sup>2</sup> + *<sup>β</sup>* + ¯ *δN*℘4℘1*λ<sup>r</sup>* + *δN*℘<sup>4</sup> ¯ *δL*℘<sup>1</sup> *qλ*). For 0 ≤ ℘<sup>1</sup> ≤ *L*, 1 ≤ ℘<sup>2</sup> ≤ *S*1, 1 ≤ ℘<sup>3</sup> ≤ *S*2, 1 ≤ ℘<sup>4</sup> ≤ *N*, *<sup>E</sup>*[*Wn*+<sup>1</sup> *<sup>p</sup>* (℘1, <sup>℘</sup>2, <sup>℘</sup>3, <sup>℘</sup>4)] = ¯ *δN*℘4*λ <sup>b</sup> <sup>E</sup>*[*Wn*+<sup>1</sup> *<sup>p</sup>* (℘1, <sup>℘</sup>2, <sup>℘</sup>3, <sup>℘</sup><sup>4</sup> <sup>+</sup> <sup>1</sup>) + *H*(*s* − *v*)*β <sup>b</sup> <sup>E</sup>*[*Wn*+<sup>1</sup> *<sup>p</sup>* (*u*, *<sup>v</sup>* <sup>+</sup> *<sup>Q</sup>*, *<sup>w</sup>*, *<sup>x</sup>*)] + *<sup>v</sup>γ*<sup>1</sup> *<sup>b</sup> <sup>E</sup>*[*Wn*+<sup>1</sup> *<sup>p</sup>* (℘1, <sup>℘</sup><sup>2</sup> <sup>−</sup> 1, <sup>℘</sup>3)] + *wγ*<sup>2</sup> *<sup>b</sup> <sup>E</sup>*[*Wn*+<sup>1</sup> *<sup>p</sup>* (℘1, <sup>℘</sup>2, <sup>℘</sup><sup>3</sup> <sup>−</sup> <sup>1</sup>)] + (℘<sup>4</sup> <sup>−</sup> <sup>1</sup>)*<sup>μ</sup> <sup>b</sup> <sup>E</sup>*[*Wn*+<sup>1</sup> *<sup>p</sup>* (℘1, <sup>℘</sup><sup>2</sup> <sup>−</sup> 1, <sup>℘</sup><sup>3</sup> <sup>−</sup> 1, <sup>℘</sup><sup>4</sup> <sup>−</sup> <sup>1</sup>)] + (19) ¯ *δN*℘4℘1*λ<sup>r</sup> <sup>b</sup> <sup>E</sup>*[*Wn*+<sup>1</sup> *<sup>p</sup>* (℘<sup>1</sup> <sup>−</sup> 1, <sup>℘</sup>2, <sup>℘</sup>3, <sup>℘</sup><sup>4</sup> <sup>+</sup> <sup>1</sup>)] + *δN*℘<sup>4</sup> ¯ *δL*℘<sup>1</sup> *qλ <sup>b</sup> <sup>E</sup>*[*Wn*+<sup>1</sup> *<sup>p</sup>* (℘<sup>1</sup> <sup>+</sup> 1, <sup>℘</sup>2, <sup>℘</sup>3, <sup>℘</sup>4)+(*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)*E*[*Wn*+<sup>1</sup> *<sup>p</sup>* (℘1, <sup>℘</sup>2, <sup>℘</sup>3, <sup>℘</sup>4)]

where, *b* = (*y* + ¯ *<sup>δ</sup>N*℘1*<sup>λ</sup>* + *<sup>H</sup>*(*<sup>s</sup>* − ℘2)*<sup>β</sup>* + ℘2*γ*<sup>1</sup> + ℘3*γ*<sup>2</sup> + ℘4*<sup>μ</sup>* + ¯ *δN*℘4℘1*λ<sup>r</sup>* + *δN*℘<sup>4</sup> ¯ *δL*℘<sup>1</sup> *qλ*). The LST of UWT of a demand in the queue is given by

$${}^{\*}W\_{p}(y) = 1 - \sum\_{\wp\_{1}=0}^{L} \sum\_{\wp\_{2}=0}^{S\_{1}} \sum\_{\wp\_{3}=1}^{S\_{2}} \sum\_{\wp\_{4}=0}^{N-1} {}^{Y}\_{\wp\_{1}\wp\_{2},\wp\_{3}\wp\_{4}}(+ \\ \tag{20}$$

$$\sum\_{\wp\_{1}=0}^{L} \sum\_{\wp\_{2}=0}^{S\_{1}} \sum\_{\wp\_{3}=1}^{S\_{2}} \sum\_{\wp\_{4}=0}^{N-1} {}^{Y}\_{\wp\_{4}=0} {}^{Y}\_{\wp\_{1}\wp\_{2},\wp\_{3}\wp\_{4}}(+ \\ \tag{21}$$

The *n*th moments of UWT, using (20), is given by

$$E[\mathcal{W}\_p^n] = \delta\_{0n} + (1 - \delta\_{0n}) \sum\_{\wp\_1=0}^L \sum\_{\wp\_2=0}^{S\_1} \sum\_{\wp\_3=1}^{S\_2} \sum\_{\wp\_4=0}^{N-1} \mathcal{Y}^{(\wp\_1, \wp\_2, \wp\_3, \wp\_4)} E[\mathcal{W}\_p^n(\wp\_1, \wp\_2, \wp\_3, \wp\_4 + 1)] \tag{21}$$

Using Equation (21) and substitute *n* = 1, we get the desired result as in (15).

**Corollary 1.** *The expected waiting time of a orbital demand is defined by*

$$E[\mathcal{W}\_o] = \sum\_{\substack{\wp\_1=0 \ \wp\_2=0}}^{L-1} \sum\_{\wp\_3=1}^{S\_2} \sum\_{\wp\_4=0}^{N} \mathcal{Y}^{(\wp\_1,\wp\_2,\wp\_3,\wp\_4)} E[\mathcal{W}\_o(\wp\_1+1,\wp\_2,\wp\_3,\wp\_4)] \tag{22}$$

#### **5. System Characteristics**

This section explores the necessary and sufficient system performance of the proposed model.

#### *5.1. Mean Inventory Level for First Commodity*

The mean inventory of the first commodity is defined as

$$\Delta\_1 = \sum\_{\wp\_1=0}^{\infty} \sum\_{\wp\_2=1}^{S\_1} \sum\_{\wp\_3=1}^{S\_2} \sum\_{\wp\_4=0}^{N} \wp\_2 Y^{(\wp\_1, \wp\_2, \wp\_3, \wp\_4)}$$

#### *5.2. Mean Inventory Level for Second Commodity*

The mean inventory of the second commodity is defined as

$$\Delta\_2 = \sum\_{\wp\_1=0}^{\infty} \sum\_{\wp\_2=0}^{S\_1} \sum\_{\wp\_3=1}^{S\_2} \sum\_{\wp\_4=0}^{N} \wp\_3 \, \mathcal{Y}^{(\wp\_1, \wp\_2, \wp\_3, \wp\_4)} \, ^\*$$

*5.3. Mean Reorder Rate for First Commodity*

The mean reorder rate of first commodity is defined as

$$\Delta\_3 = \sum\_{\wp\_1=0}^{\infty} \sum\_{\wp\_3=1}^{S\_2} \sum\_{\wp\_4=1}^{N} \wp\_4 \mu \mathcal{Y}^{(\wp\_1,s+1,\wp\_3,\wp\_4)} + \sum\_{\wp\_1=0}^{\infty} \sum\_{\wp\_3=1}^{S\_2} \sum\_{\wp\_4=0}^{N} (s+1) \gamma\_1 \mathcal{Y}^{(\wp\_1s+1,\wp\_3,\wp\_4)}$$

#### *5.4. Mean Reorder Rate for Second Commodity*

The mean reorder rate of second commodity is defined as

$$\Delta\_4 = \sum\_{\wp\_1=0}^{\infty} \sum\_{\wp\_2=1}^{S\_1} \sum\_{\wp\_4=1}^{N} \wp\_4 \mu Y^{(\wp\_1,\wp\_2,1,\wp\_4)} + \sum\_{\wp\_1=0}^{\infty} \sum\_{\wp\_2=0}^{S\_1} \sum\_{\wp\_4=0}^{N} \gamma\_2 Y^{(\wp\_1,\wp\_2,1,\wp\_4)}$$

#### *5.5. Mean Perishable Rate for First Commodity*

The mean perishable rate of first commodity is given by

$$\Delta\_5 = \sum\_{\wp\_1=0}^{\infty} \sum\_{\wp\_2=1}^{S\_1} \sum\_{\wp\_3=1}^{S\_2} \sum\_{\wp\_4=0}^{N} \wp\_2 \gamma\_1 \mathcal{Y}^{(\wp\_1, \wp\_2, \wp\_3, \wp\_4)}$$

#### *5.6. Mean Perishable Rate for Second Commodity*

The mean perishable rate of second commodity is given by

$$\Delta\_6 = \sum\_{\wp\_1=0}^{\infty} \sum\_{\wp\_2=0}^{S\_1} \sum\_{\wp\_3=1}^{S\_2} \sum\_{\wp\_4=0}^{N} w \gamma\_2 Y^{(\wp\_1, \wp\_2, \wp\_3, \wp\_4)}$$

#### *5.7. Mean Number of Primary Customer Lost*

The mean number of primary customer lost in the system is defined as

$$\Delta\_{11} = \sum\_{\wp\_1=0}^{\infty} \sum\_{\wp\_2=0}^{S\_1} \sum\_{\wp\_3=1}^{S\_2} (1-q)\lambda \mathcal{Y}^{(\wp\_1,\wp\_2,\wp\_3,N)}$$

*5.8. Overall Rate of Retrial*

The expected overall rate of retrial of the orbit customer in the system is given by

$$\Delta\_{12} = \sum\_{\wp\_1=1}^{\infty} \sum\_{\wp\_2=0}^{S\_1} \sum\_{\wp\_3=1}^{S\_2} \sum\_{\wp\_4=0}^{N} \wp\_1 \lambda\_r \mathcal{Y}^{(\wp\_1, \wp\_2, \wp\_3, \wp\_4)}$$

#### *5.9. Successful Rate of Retrial*

The expected successful rate of retrial of the orbit customer in the system is given by

$$\Delta\_{13} = \sum\_{\wp\_1=1}^{\infty} \sum\_{\wp\_2=0}^{S\_1} \sum\_{\wp\_3=1}^{S\_2} \sum\_{\wp\_4=0}^{N-1} \wp\_1 \lambda\_I \Upsilon^{(\wp\_1, \wp\_2, \wp\_3, \wp\_4)}$$

#### *5.10. Fraction of Successful Rate of Retrial*

The expected fraction of successful rate of retrial of the orbit customer in the system is given by

$$
\Delta\_{14} = \frac{\Delta\_{13}}{\Delta\_{12}}
$$

#### **6. Cost Analysis and Numerical Illustrations**

The mean total cost (*MTC*) is given by

*MTC* = (*Ch*<sup>1</sup> ∗ *Δ*1)+(*Ch*<sup>2</sup> ∗ *Δ*2)+(*Cs*<sup>1</sup> ∗ *Δ*3)+(*Cs*<sup>2</sup> ∗ *Δ*4)+(*Cp*<sup>1</sup> ∗ *Δ*5)+(*Cp*<sup>2</sup> ∗ *Δ*6) + (*Cw*<sup>1</sup> ∗ *E*[*Wp*]) + (*Cw*<sup>2</sup> ∗ *E*[*Wo*]) + (*Cb* ∗ *Δ*11).

To compute the *MTC* per unit time, the following costs are considered


*Cb* = Cost due to loss of customers per unit time.

#### *Numerical Illustrations*

Numerical analysis is an applied mathematics technique that allows staggeringly large amount of data to be processed and analyzed for trends, thereby aiding in forming conclusions. This is done nowadays in a computer environment, providing massive increases in speed and usefulness of calculations. This numerical work is carried out by fixing the parameters and cost values which are involved in the proposed mathematical model. Due to the nature of the proposed model, the determination of the parameter values assumed randomly under the satisfaction of the stability condition. The numerical value of stationary probability vector whose sum gives one. After such verification's (stability condition and normalising condition), we proceed the numerical works. To perform these numerical illustrations, we fix the parameter values as *Ch*<sup>1</sup> = 0.001, *Ch*<sup>2</sup> = 0.01,*s* = 4, *Cs*<sup>1</sup> = 4.1, *Cs*<sup>2</sup> = 3.5, *Cp*<sup>1</sup> = 2, *Cp*<sup>2</sup> = 1.8, *Cw*<sup>1</sup> = 1.6, *Cw*<sup>2</sup> = 1.1, *Cb* = 3.6, *N* = 5, *λ* = 3.6, *β* = 0.93, *γ*<sup>1</sup> = 0.21, *γ*<sup>2</sup> = 0.2, *λ<sup>r</sup>* = 4.46, *q* = 0.5, *μ* = 25.5. The proposed study is only valid if the assumed parameters must satisfy the stability condition and normalising condition.

**Example 1.** *From Table 1, for varying S*1*, the economical capacity S*1(*S*<sup>∗</sup> <sup>1</sup> ) *is determined at each S*<sup>2</sup> *and is identified by the expected total cost which is column minimum with underline. Similarly, the economical capacity S*2(*S*<sup>∗</sup> <sup>2</sup> ) *is determined at each S*<sup>1</sup> *and is identified by the expected total cost which is row minimum and is given in bold form. When we increase the number of first commodity, the MTC holds the both decreasing and increasing property. It shows that the MTC is not monotonic regarding the change in S*1*. This will produce the minimum MTC at some S*1*. Similarly, the same monotonic property holds for the second commodity. Therefore, we also obtain the minimum MTC for S*2*. Since both S*<sup>1</sup> *and S*<sup>2</sup> *holds the minimum MTC, the proposed model produces the convex point at some* (*S*1, *S*2) *locally. Finally the local minimum expected total cost is determined at S*<sup>∗</sup> 1 *and S*<sup>∗</sup> <sup>2</sup> *where the row minimum and column minimum matches, so that the cost will be identified with both underlined and in bold script. Using Table 1, S*<sup>∗</sup> <sup>1</sup> <sup>=</sup> 25, *<sup>S</sup>*<sup>∗</sup> <sup>2</sup> <sup>=</sup> 7, *MTC*<sup>∗</sup> <sup>=</sup> 8.991429*.*

*The graphical illustration of the optimum MTC is shown in Figure 1. If the values of S*<sup>1</sup> *increases, then the value of MTC have both increasing and increasing quality. It means that S*<sup>1</sup> *will give the optimum MTC. When we work on the S*<sup>2</sup> *with the similar manner as S*1*, this also gives the minimum MTC. Finally, both S*<sup>1</sup> *and S*<sup>2</sup> *are increased together, we obtain the convex at some point of this proposed mathematical model. This convexity also helps us to determine the fixation of parameters and cost values of the model. Under the obtained convex result, we perform the further numerical work.*

**Figure 1.** *S*<sup>1</sup> vs. *S*<sup>2</sup> on *MTC*.

**Table 1.** Mean total cost rate as a function of *S*<sup>1</sup> and *S*2.


**Example 2.** *From the exploration of Figures 2–6, the responses of Δ*3,*Δ*9,*Δ*11,*Δ*<sup>13</sup> *and MTC under the subject of queue capacity and primary arrival rate are explained below.*


**Figure 2.** *λ* vs. *N* on *Δ*3.

**Figure 3.** *λ* vs. *N* on *Δ*9.

**Figure 4.** *λ* vs. *N* on *Δ*11.

**Figure 5.** *λ* vs. *N* on *Δ*13.

**Figure 6.** *λ* vs. *N* on *MTC*.

**Example 3.** *Table 2 shows the different characteristics of MTC and Δ*<sup>11</sup> *with varying combinations Q*,*s*, *S*<sup>2</sup> *and q with set up cost dependent lead time for commodity 1.*



**Table 2.** Mean total cost vs. Customer lost.

**Example 4.** *From Table 3, Under the given cost structure and some parameters associated with service and any kind of arrivals, the comparisons of results of queue dependent with non-queue dependent service policies, the merits of the proposed model are studied and stated below:*



**Table 3.** Queue dependent service time Vs Non queue dependent service time.


**Table 3.** *Cont.*

#### **7. Conclusions**

This paper analyses a single server two commodity inventory system with queuedependent services for finite queue and an optional retrial facility. We applied a Neuts matrix geometric approach to resolve the infinitesimal generator matrix and further we derived a stability condition and the steady state probability vector of the system. Upon computing the necessary system characteristics of the system, we obtained the mean total cost of the considered model. This model explored the queue dependent and non-queue dependent service polices in a two commodity retrial inventory system. The minimal optimum mean total cost is obtained for the queue-dependent service policy. Indeed, this policy also helped to reduce the rate of customer lost and his/her expected waiting time. Apart from that the retrial customers successful retrial rate also increased as we expected. As for the consequences, the queue dependent service is more profitable when the major product attached with its complimentary. As a result of the numerical survey, maximum stock levels of both products are economically controlled. Furthermore, the fixed stocking capacities of inventory and the preference level of customer entering into the orbit, the system may control the size of queue with the ordering quantity which corresponds to lead time and set up cost. However, the comparison of queue dependent and non-queue dependent service polices will give a great impact on the inventory management. The readers can easily understand these different service polices which one is adoptable to bring out the profitable business. This model will be extended into a multi server service system in the future. In such an environment, one can analyse the impact of customer lost and their waiting time with the single server model. This analysis also helps to the business people to develop their business.

**Author Contributions:** Conceptualization, K.J.; data curation, K.J.; formal analysis, K.J., M.A.R. and S.A.; funding acquisition, G.P.J. and C.S.; investigation, N.A. and D.J.; methodology, M.A.R., S.S. and S.A.; project administration, G.P.J. and C.S.; software, M.A.R.; supervision, C.S.; validation, G.P.J. and D.J.; visualization, S.S. and N.A.; writing—review and editing, G.P.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Institute of Information and Communications Technology Planning and Evaluation (IITP) Grant by the Korean Government through MSIT (Development of User identity certification and management technology for self-sovereign identity applications.) under Grant 2021-0-00565.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** Anbazhagan and Amutha would like to thank RUSA Phase 2.0 (F 24-51/2014-U), DST-FIST ( SR/FIST/MS-I/2018/17), DST-PURSE 2nd Phase programme (SR/PURSE Phase 2/38), Govt. of India.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Notations**


#### **References**


**Yu-Kai Weng 1, Cheng-Han Li 1, Chia-Chun Lai <sup>1</sup> and Ching-Wei Cheng 2,\***


**Abstract:** In the egg industry, it is necessary to estimate the egg volume accurately when estimating egg quality or freshness in a non-destructive method. Egg volume and weight could obtain egg density and could be used to determine egg freshness. Therefore, the egg geometric must be obtained first to establish a volume equation with a geometric shape. This research proposes an innovative idea to derive the mathematical model and volume equation of egg shape, calculate its volume, and verify the accuracy of the mathematical equation proposed using the volume displacement method. Using the proposed equation, the minimum error between the calculated egg volume) and actual egg volume is 0.01%. The maximum volume error does not exceed 2%. The egg shape equation can accurately draw the outer contour curve of the egg by the half-length of the maximum long axis and maximum breadth of the short axis, and the distance from the center point of the egg to the maximum breadth (*xm*).

**Keywords:** egg shape equation; displacement of volume method; egg volume

**MSC:** 14-11

#### **1. Introduction**

Egg geometry is often used in food research, agricultural engineering, biological sciences, mechanical engineering, architecture, and so on. This involves research on the classification and ecological morphology of poultry populations [1,2], predicting the weight relationship of chicks after egg hatching [3], and egg hatchability [4,5]. An egg can withstand external forces exerted during grading and transportation. Therefore, the mechanical properties of simulated eggshells are usually analyzed [6], and the best protection method is proposed. Eggshell shapes have also been used to design underwater installations, containers [7], and buildings [8]. Therefore, mathematical equations for egg shape, egg volume, and related parameters are required to perform related research and in applications.

Egg shape can be divided into oval, pyriform, circular, and elliptical forms [1,9]. The shape classification of eggs is mainly based on the ratio of the maximum breadth (*B*) to the maximum length (*L*) of the egg multiplied by 100, which is called the shape index (SI) [10,11]. The SI is used to judge whether the egg is approximately round or oval (i.e., degree of shape). Mathematically speaking, eggs are prolate spheroids [5], which approximate the volume of ellipsoids. Therefore, for egg volume calculation, the deformation equation with the volume of ellipsoids as the base is often used for calculation, as shown in Equation (1) [5,10].

$$V = k\_v \frac{4}{3} \pi \left(\frac{L}{2}\right) \left(\frac{B}{2}\right)^2 \approx k\_v L B^2. \tag{1}$$

**Citation:** Weng, Y.-K.; Li, C.-H.; Lai, C.-C.; Cheng, C.-W. Equation for Egg Volume Calculation Based on Smart's Model. *Mathematics* **2022**, *10*, 1661. https://doi.org/10.3390/ math10101661

Academic Editors: Wen Zhang, Xiaofeng Xu, Jun Wu and Kaijian He

Received: 14 April 2022 Accepted: 10 May 2022 Published: 12 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

*kv* is the compensating coefficient, *B* is the maximum width equatorial, and *L* is the maximum egg length. Eggs shape of various breeds is different, and the *kv* value of the volume calculation Equation (1) will also change accordingly [12]. The commonly used *kv* = 0.5236 is a standard ellipsoid. To increase the accuracy of the egg volume calculation, Narushin (1997) set *kv* = 0.496 [13]. A follow-up study by Narushin (2005) proposed a more accurate correction equation for egg volume, *kv* = (0.6057 − 0.0018B), and the coefficient of determination reaching 0.958, as shown in Equation (2) [14]. It was subsequently proved that *kv* is not a constant but a function of the linear parameters of the egg shape (that is, the maximum length *L* and maximum breadth *B* of the egg) [14–16]. Subsequently, Narushin (2021) again proposed that *kv* = 0.5136 and Equation (2) can be more accurate and closer to the actual egg volume when calculations are performed [15].

$$V = \left(0.5202LB^2 - 0.4065\right). \tag{2}$$

In 1948, the German mathematician Fritz Hügelschäffer used two non-concentric circles to construct a deformed ellipse to form an egg-shaped curve [8,17]. Preston (1953) multiplied the basic ellipse equation by a cubic polynomial so that the deformed ellipse equation can describe a variety of egg shapes [18]. Smart (1969) found the tangent angle (taper angle) of the edge of the egg-shaped contour and introduced the short axis of the ellipse equation to describe the outer shape of the egg [19]. Narushin et al. (2020) applied Hügelschäffer's egg-shaped contour model to calculate the volume of an egg [20].

This study aims to derive a simplified equation that can quickly, directly, and accurately calculate the egg volume (V) based on the egg-shaped contour equation proposed by Smart (1969) and use the displacement of volume method to calculate volume. The egg volume calculated by the equation established in this paper is compared to the actual measured volume. This helps compare the theoretical equation and actual measurement error value as well as the accuracy of the theoretical equation. Further verifications were made by comparing the present results to the previously reported egg volume equations. The egg volume equation established in this study can be used as one of the reference methods for quick and accurate egg volume calculation. In addition, the egg-shaped contour equation was drawn and compared with the actual shooting shape, which further confirmed the reliability of the egg-shaped contour equation drawn in this research.

#### **2. Materials and Methods**

#### *2.1. Egg Sample*

The theoretical equation for egg volume proposed in this paper uses 47 eggs and measures the respective volume of each egg by the displacement of volume method. It is used to verify the difference between the egg volumes found using the theoretical equation and the actual experimental measurement.

#### *2.2. Egg Dimension Parameter and Volume Measurement*

For egg dimensions, a digital caliper is used to measure the maximum long (*L*) and maximum short (*B*) axes of eggs and the distance from the long axis half-length to the plane of the maximum breadth of eggs (*xm*). The resolution of the digital caliper can reach 0.01 mm.

The actual volume of the egg is measured according to the volume displacement method, with the following practical guideline steps, as shown in Figure 1. Place the egg to be measured in a 500 mL graduated empty cylinder, and add water up to the 400 mL mark. Remove the egg from the measuring cylinder, reset the precision electronic balance (HG-2000, Shinko Denshi Co. Ltd., Tokyo, Japan) to zero, and then add the water in the measuring cylinder to 400 mL. As the density of water approaches 1 g/cm3, the weight of the added water corresponds to the volume of the eggs. Removing the eggs from the graduated cylinder requires preventing the water in the graduated cylinder from overflowing. Therefore, let the water droplets adhering to the egg's surface remain in the measuring cylinder before removing the egg for higher precision. The water droplet adhering to the surface of the eggshell has a significant influence on the measurement of egg volume.

**Figure 1.** Flow chart of egg volume measurement.

#### *2.3. Egg Volume Equation*

Preston (1953) multiplied the basic ellipse equation by a polynomial to modify the ellipse shape to achieve a variety of egg shape equations for describing the geometric shapes of various bird eggs [18]. The follow-up work by Smart (1969) also used the ellipse equation as the basis to measure the tangent angle of the outer edge of the egg and introduced the ellipse equation to describe the shape of the egg [19], as follows:

$$\frac{\alpha^2}{a^2} + \frac{y^2}{\left(b + x \tan \theta\right)^2} = 1.\tag{3}$$

Köller [21] multiplied the *y*<sup>2</sup> term by *t*(*x*) based on the ellipse Equation:

$$\frac{x^2}{a^2} + \frac{y^2}{b^2} \cdot t(x) = 1,\tag{4}$$

where *a* and *b* are the long and short axes of the ellipse, respectively, as shown in Figure 2. *t*(*x*) = 1 + *kx* is a simple function, which converts the ellipse into an egg-shaped curve. When the value of *k* is large, the *y*-axis of the ellipse is not symmetric. When the value of *k* is small, it will be closer to the standard ellipse shape, as shown in Figure 2b. *k* = 0 corresponds to the standard elliptic curve.

**Figure 2.** (**a**) Smart (1969) egg-shaped curve parameter location diagram. (**b**) Köller's egg-shaped curve diagram of the change of the function *k* value in Equation (4).

We have found that when (*b* + *x* tan *θ*) <sup>−</sup><sup>2</sup> in Smart (1969) equation is converted by power series, (1 + 2*kx*) can be obtained, and the conversion process is as follows:

Equation (5) can be obtained by expanding Equation (3) to power series as

$$\frac{\alpha^2}{a^2} + \frac{y^2}{b^2} \left[ \left( 1 - 2 \frac{\tan \theta}{b} x \right) + \frac{-2(-3)}{2!} \left( \frac{\tan \theta}{b} x \right)^2 + \dotsb \right] = 1. \tag{5}$$

When *x* = *xm*, Equation (6) can be obtained by neglecting the higher-order terms in Equation (5) as

$$\left(\frac{x^2}{a^2} + \frac{y^2}{b^2}\left(1 - 2\frac{\tan\theta}{b}x\right) = 1\right) \tag{6}$$

and *<sup>x</sup>*<sup>2</sup>

$$\frac{\chi^2}{a^2} + \frac{y^2}{b^2}(1+kx) = 1,\tag{7}$$

where *k* = −2 tan *<sup>θ</sup> <sup>b</sup>* .

From the above process, it can be seen that the egg-shaped contour equations of Köller (2000) and Smart (1969) are consistent. To calculate the egg volume, we use the Smart (1969) egg contour equation and find the egg volume from the equation by rotating the integral around the long axis. The process is as follows:

$$V = \int\_{-a}^{a} \pi y^2 dx \tag{8}$$

$$=\int\_{-a}^{a} \pi \left(1 - \frac{\mathbf{x}^2}{a^2}\right) (b + \mathbf{x}\tan\theta)^2 d\mathbf{x} \tag{9}$$

$$\hat{\theta} = \frac{\pi}{a^2} \left[ \frac{4}{3} a^3 b^2 + \frac{4}{15} a^5 \tan^2 \theta \right]. \tag{10}$$

Finally, the egg volume equation is obtained:

$$V = \pi (\frac{4}{3}ab^2 + \frac{4}{15}a^3 \tan^2 \theta),\tag{11}$$

where *a* is the half-length of the long axis of the egg, *b* is the half-length of the short-axis breadth of point *O* in the center of the egg, and *θ* is the egg taper angle, as shown in Figure 2.

The short-axis breadth of the egg center point of the long axis (2*b*) and edge taper angle (*θ*) are not easy to measure. In this study, the half-length of the long axis of the egg (*a*), maximum breadth of the short axis of the egg (*bm*), and distance from the long axis half-length to the plane of maximum breadth (*xm*) are known conditions. The correct egg shape can be obtained by calculating the theoretical values *b* and tan *θ* of the egg shape parameters.

Equation (12) can be obtained by sorting out Equation (6) as follows:

$$y^2 = \left(1 - \frac{\mathbf{x}^2}{a^2}\right) \frac{b^2}{\left(1 - \frac{2}{b} \cdot \mathbf{x} \cdot \tan\theta\right)}.\tag{12}$$

Equation (12) is differentiated such that:

$$y\frac{dy}{dx} = \frac{-\mathbf{x}}{a^2} \frac{b^2}{\left(1 - \frac{2}{b} \cdot \mathbf{x} \cdot \tan\theta\right)} + \frac{1}{b} \tan\theta \cdot \left(1 - \frac{\mathbf{x}^2}{a^2}\right) \cdot \frac{b^2}{\left(1 - \frac{2}{b} \cdot \mathbf{x} \cdot \tan\theta\right)^2}.\tag{13}$$

When *dy dx* = 0 is the horizontal tangent of the maximum width of the egg, *x* = *xm*, the result is as follows:

$$1 - \chi\_m + \frac{a^2}{b} \tan\left(1 - \left(\frac{\chi\_m}{a}\right)^2\right) \cdot \frac{1}{1 - \frac{\chi\_m}{b} \cdot \tan\theta} = 0,\tag{14}$$

According to the actual measurements in Table 1, putting the *xm* and egg length of halflength ( *<sup>L</sup>* <sup>2</sup> ) average into 1 <sup>−</sup> *xm a* <sup>2</sup> = 0.99 can be obtained. The short-axis breadth half-length (*b*) and egg maximum breadth half-length ( *<sup>B</sup>* <sup>2</sup> ) are close, so we use the maximum breadth of half-length directly for calculation; alos obtaining 1 − *xm <sup>b</sup>* · tan *θ* = 0.99.

**Table 1.** Egg shape geometry parameters.


Therefore assume when, 1 − *xm a* <sup>2</sup> ∼= 1 and 1 − *xm <sup>b</sup>* · tan *<sup>θ</sup>* <sup>∼</sup><sup>=</sup> 1, then:

$$\propto\_m \cong \frac{a^2}{b} \tan \theta. \tag{15}$$

From Figure 2a, the following can be obtained:

$$\frac{b\_m - b}{\varkappa\_m} \cong \tan \theta. \tag{16}$$

From Equations (15) and (16), the short-axis breadth half-length (*b*) of point *O* in the egg can be obtained as follows (see Figure 2a):

$$b = \frac{a^2 b\_m}{(\varkappa\_m^2 + a^2)}.\tag{17}$$

By substituting Equations (15) and (17) into Equation (11), the final egg volume can be obtained as:

$$V = \pi a^3 \left(\frac{b\_m}{x\_m^2 + a^2}\right)^2 \cdot \left(\frac{4}{3}a^2 + \frac{4}{15}x\_m^2\right) \,,\tag{18}$$

where *bm* = *<sup>B</sup>* <sup>2</sup> and *<sup>a</sup>* <sup>=</sup> *<sup>L</sup>* 2 .

#### **3. Results and Discussion**

This paper uses 47 eggs to verify the egg volume equation (Equations (11) and (18)) established based on Smart's equation (Equation (3)). The actual measured parameter results of 47 eggs are shown in Table 1. The most significant variation of egg shape parameters is the actual volume of the egg and the coefficient of variation of *xm*, which are 23.53% and 8.11%, respectively. Egg maximum breadth has the least variability.

The egg volume Equation (11) is introduced by Equations (15) and (17) to obtain the final egg theoretical volume Equation (18). The angle of taper (*θ*) of the egg and half the breadth of the center point (*b*) of the egg can be described by Equations (15)–(17). To confirm the accuracy of the calculated taper angle (*θ*) and the short-axis breadth (*b*) of the egg center point, we selected the egg SIs as large: 81.08, medium: 76.11, small: 70.97. The profile of actual eggs of different specifications was compared with the calculated curve profile, as shown in Figure 3 The volume percent errors of the three eggs were 0.32%, −0.58%, and 0.33%, respectively, within ±1%. MATLAB was used to draw the

egg-shaped contour curve of Equation (3). *b* and tan *θ* in Equation (3) were calculated using Equations (15)–(17). Actual pictures of the eggs were added. This helps to check whether the theoretical equation egg contour, taper angle (*θ*) and half of the short axis width of the egg center point (*b*) are consistent with the actual egg contour. The results are shown in Figure 3. The calculated egg taper angle (*θ*) and half of the short-axis width of the center point (*b*) are almost consistent with the actual egg contour.

**Figure 3.** The actual egg shape and egg-shaped curve drawn by Equations (3) and (15)–(17). The SI (**a**) is 81.08; (**b**) 76.11; (**c**) 70.97.

The egg-shaped contour curve is converted into descriptions by *a*, *bm*, and *xm* from *b* and tan *θ* in Equation (3). Therefore, the egg-shaped contour generated by the mathematical equation can approximate the real egg-shaped contour curve. However, according to the research and investigation in this study, there are still errors in the contour curves drawn by Equation (3) and Equations (15)–(17). In addition, because eggs are biological samples, they may be extruded deformed during production, resulting in an uneven shape.

For egg volume calculation, Equation (1) was used in the past. The *kv* value was mainly set to 0.5236 for the elliptical volume equation. To make it more accurate, Narushin (2021) found that for *kv* = 0.5163, calculation error variation can be reduced to a range of 0–3.7%, with an average of 1.1% [15]. Additionally, from the volume of the Hügelschäffer egg-shaped body and ellipsoids volume Equation (1) *kv* = 0.5236, another Equation (2) for calculating the egg-shaped volume is obtained [15]. Equations (2) and (19) calculate the egg-shaped volume, with no significant difference between the two values obtained [15].

$$V = 0.5163LB^2.\tag{19}$$

In this study, the actual volume of eggs was measured by the volume displacement method. The shape parameters of 47 eggs were substituted into the theoretical egg volume Equation (18), ellipse Equation (1), and egg volume Equations (2) and (19) proposed by Narushin (2021). Errors between the egg volume calculated by the four equations and the actual measured egg volume were compared. The root mean square (rms) error values obtained were 0.47, 0.96, 0.43, and 0.63, respectively. This indicates that our theoretical equation and that proposed by Narushin (2021) (*kv* = 0.5163) can calculate relatively accurate egg volumes. With *kv* = 0.5236 ellipse volume equation, the rms error of egg volume calculated is only 0.96. However, because eggs are not in a standard ellipse shape, the error between the egg volume calculated by the ellipse volume Equation (1) (*kv* = 0.5236) and the actual measured volume will be more significant.

We use the actual measured volume of the egg and the calculated volume of the equation to express the percentage error, which is a standard used to measure the accuracy of the measurement results. It is calculated as the estimate minus the actual value divided by the actual value, multiplied by 100%. It is used to express the accuracy of the volume converted by its equation.

Figure 4a shows that from Narushin (2021) Equation (19), and theoretical calculation Equation (18) in this study, the calculated volume percentage errors are mainly distributed within ±2%. Therefore, from the volume obtained by the theoretical volume equation in this study and that measured by the volume of displacement method, the volume percentage errors of ±0.5%, ±1%, ±1.5%, and ±2% correspond to 22 eggs, 16 eggs, 7 eggs, and 2 eggs, respectively. Therefore, the standard deviation is 0.47%. The volume percentage errors calculated by Equation (19) of Narushin (2021) are ±0.5%, ±1%, ±1.5%, and 2%, corresponding to 22 eggs, 17 eggs, 7 eggs, and 1 egg, respectively. The standard deviation is 0.4%. For the theoretical volume equation (i.e., Equation (18)) in this study, the maximum volume error of the 47 eggs is 1.91%, minimum volume error is 0.01%. Narushin's equation (i.e., Equation (19)) maximum error is 1.59%, and the minimum volume error is 0.05%. The resulting data point distribution and volume calculation accuracy are roughly the same for our theoretical volume Equation (18) and Narushin's Equation (19). However, the maximum volume percentage error of Equation (18) proposed in this paper is 1.91%, which may be caused by the fact that Equation (5) is expanded by the power series method, and the higher-order terms of Equation (6) are omitted for the convenience of calculation. Another reason may be that tan *θ* and *b* are approximate values obtained by converting Equations (15)–(17) after measuring *a*, *bm*, and *xm*, resulting in a slight error in using the theoretical volume Equation (18) in this study. The percentage of error in the volume calculation of the ellipse volume equation is mostly greater than 1.5%, the maximum volume error can reach 3%, and the standard deviation is 0.74%, while the egg volume calculated by Narushin's Equation (2) is mostly evenly distributed between 0 and 3%, as shown in Figure 4b. Figure 4a shows that the egg volume obtained by this theoretical Equation (18) and Narushin's (2021) Equation (19) is quite accurate. Therefore, both the theoretical equation and Equation (19) can accurately and quickly estimate the egg volume.

**Figure 4.** The volume percentage error corresponds to the actual weight of the egg. (**a**) theoretical volume Equations (18) and (19); (**b**) theoretical volume Equation (18), elliptical volume equation, Equation (2).

The theoretical egg volume Equation (18) established in this study is derived from Smart's egg-shaped curve Equation (3), with the parameters shown in Figure 2a including the taper angle (*θ*) and short-axis width (b) of the egg's geometric center point. The egg volume equation was verified using 47 eggs, and the egg volume could be accurately estimated when the SI of the egg ranged from 83.81 to 70.97. The research equation in this paper can be applied to analytical software such as the finite element method in the future to model eggshell or egg-shaped structures. Furthermore, it helps in improving the accuracy of the analysis in engineering or agriculture for analyzing parameters such as the mechanical strength of egg-shaped structures.

#### **4. Conclusions**

This study used the egg-shaped curve equation proposed by Smart (1969) [19] to extend the calculation equation of egg volume as:

$$V = \pi a^3 \left(\frac{b\_m}{x\_m^2 + a^2}\right)^2 \cdot \left(\frac{4}{3}a^2 + \frac{4}{15}x\_m^2\right). \tag{20}$$

The egg volume equation proposed in this study only requires the measurement of the half maximum length of the long axis (*a*), half maximum breadth of the short axis (*bm*), and distance from the center of the egg to the maximum breadth of the short axis (*xm*). Thus, it can quickly estimate the volume of the egg. It was verified that the egg volume with a SI of 70.97 to 83.81, the error range of the calculated theoretical volume and actual volume is within ±2%.

Another contribution of this study is that we do not need to measure the egg's outer taper angle (*θ*) and breadth (*b*) of the center point. Instead, it can be converted from *xm*, *bm*, and *a* using Equations (15)–(17) established in this work. By obtaining the length distance between *θ* and *b*, Smart's (1969) egg-shaped curve Equation (3) can be directly used to draw the egg-shaped curve that is almost the same as the actual egg shape. In addition, the equation proposed in this paper can be applied to the fields of agriculture or food to analyze and simulate the curve modeling of eggs, which will be of great help and can effectively improve the simulation accuracy.

**Author Contributions:** Y.-K.W., C.-H.L., C.-C.L. and C.-W.C. conceived and designed the experiments; Y.-K.W., C.-H.L., C.-C.L. and C.-W.C. performed the experiments; Y.-K.W., C.-H.L., C.-C.L. and C.-W.C. analyzed the data; Y.-K.W., C.-H.L., C.-C.L. and C.-W.C. wrote the paper; Y.-K.W., C.-H.L., C.-C.L. and C.-W.C. reviewed and edited the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

