# A Stochastic Model for an Input Control Problem in a Two-Level Supply Chain for Production-Time-Dependent Products with Random Demands

Authors: Woo-Sung Kim; Jun-Hee Han
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Inspired by the poultry farming process, we studied an input control problem in a two-level supply chain for production-time-dependent products with random demands. Poultry farms deliver chicks in batches, and the raising process begins at timed intervals. Chicks become broilers after a predetermined raising time. The broilers in the process are shipped to manufacturing plants to satisfy the demand. The remaining chicks grow to the next product size to satisfy the demands for larger chickens. This procedure is repeated until the chicks are fully grown. After the chicks are grown to satisfy the demand for the largest size, the remaining chicks are discarded. In this paper, a stochastic model is presented to study an input control problem in the poultry farming process. Because of the production density, feed, and temperature control, one important issue in the operation of a poultry company is the determination of the raising interval and quantity of input (chicks). While existing mathematical models can provide effective information on the production-planning problem of systems, research has not been conducted on cases of random demand. Identifying a recursive structure and Markovian property for the number of raw materials (chicks) and the unfulfilled demand for each product type in the system, we demonstrate that embedded Markov chain models can be obtained. The equilibrium probabilities of the models can be calculated using matrix analytic methods or probability generating functions. Various numerical experiments are conducted to analyze how performance measures such as amount of disposal, unsatisfied demands, and total cost (considering disposal cost and opportunity cost) change with system parameters.


---

## Body

## 1. Introduction

As the global production and demand for poultry increases, it is necessary to operate the poultry supply chain effectively [1]. According to a report published by U.S. Department of Agriculture (USDA), the U.S. livestock industry has continued to grow since 2001 (Figure 1). The report states that both production and demand for beef, pork, and poultry are likely to continue increasing until 2031 assuming stable live-animal prices. Owing to this increase in demand, the wholesale price of a broiler is likely to be 98.8 cents per pound in 2023 and is projected to increase to 106.5 cents per pound by 2031. In addition, U.S. poultry exports are predicted to grow. Also, as shown in Table 1, according to a report of Food and Agriculture Organization (FAO) [2], the production of poultry meats has increased by 1.39% for all continents of the world even in the pandemic of coronavirus disease 2019 (COVID-19). The scenario in the poultry industry in South Korea is similar. According to U.S. Department of Agriculture [3], Korea’s poultry production and consumption are also projected to increase in 2022. Both import and export of chicken would increase by approximately 3% and 57%, respectively. The large increase in exports is because the exports to Vietnam would resume as the Korean poultry industry recovers from the damage caused by highly pathogenic avian influenza (HPAI) in 2021. Another reason is the growing demand for chicken. The number of parent stock broilers in Korea increased by 6.3% in the first half of 2021. This implies that chicken production would increase over the next year. Meanwhile, competition among chicken companies is intensifying in Korea. It is stated in [3] that vertically integrated chicken companies continued to compete to increase market share. It is essential to operate supply chains effectively to survive in a scenario of severe market competition among major chicken companies.Figure 1
U.S. production of animal product production, 2001–31 (source: Department of Agriculture, USDA agricultural projections from 2001 to 2031).Table 1
Production of poultry meat (thousand tons) in the world continents in 2019 and 2020 (source: FAO, 2021).
Continents of the worldYearIncreasing rate (%)20192020Europe22,08922,2890.91Asia49,66950,3671.41Africa6,5686,7582.89South Africa1,8161,9658.20Central America and the Caribbean5,2795,3992.27South America22,03022,2631.06North America24,36124,5920.95Oceania1,5651,5972.04Total133,377135,2301.39Several studies have been conducted on the effective operation of poultry production. Manning and Baines [4] attempted to identify the key factors that have caused the globalization of the poultry supply chain and their impacts on the chain performance. Shamsuddoha [5] studied the poultry industry in Bangladesh to explore a reverse supply chain management model. Bukhori et al. [6] employed the supply chain operations reference (SCOR) and analytic hierarchy process (AHP) methods to evaluate the performance measures of a poultry slaughterhouse. Taube-Netto [7] employed a mathematical model to increase the operational efficiency of Sadia Concórdia, the largest poultry producer in Brazil. According to this study, the company saved over USD 50 million over three years by improving its operation. Oladokun and Johnson [8] employed a linear programming approach to develop an optimization feed formulation model using local feed ingredients for the Nigerian poultry industry. Further details on the poultry supply chain and its optimization problems can be found in [9, 10]. Refer to [11] for further details on the research into quantitative models for evaluating the performance of a general supply chain.Recently, there are some studies to determine operation strategies of several types of supply chain models and applications. These previous studies introduce new considerations such as ordering automation [12] and green cost [13], but the concept of production-time-dependent products is considered in [14] for the first time. It reflects the product characteristics of the poultry production process. Chicks become chickens of different sizes (weights) depending on the time for which these are raised in the poultry farm. Thus, a product is called production-time-dependent product if its raw material transforms into different types of the product according to the production time (time spent in the process). Even before this concept was proposed, numerous studies had been conducted on production systems with different production times according to the product type. In previous studies, it was assumed that the production time differs depending on the type of job or product [15–18]. The authors in [19, 20] considered a system with one facility. It was assumed that the production time of the product varies depending on the resources performing the work. For systems with stochastic production time, refer to [21–25].The production-time-dependent products proposed in [14] differ from those in previous studies in that even when only one type of raw material is used in the process, different types of products are obtained depending on the production time. With regard to systems that consider production-time-dependent products, Han and Kim [14] solved the network design problem of a two-level supply chain. They developed a Lagrangian-based heuristic algorithm for the deterministic information. Han et al. [26] studied production-planning problems involving outsourcing and considering production-time-dependent products. In [27], a mixed-integer programming model and a heuristic algorithm were developed for the production-planning problem considering multiple production lines for each supplier.While the papers [14, 26, 27] that studied the supply chain for production-time-dependent products consider the characteristics of the product well in mathematical models, they have limitations in that they assume dynamic but deterministic demand. In the actual supply chain, it is almost impossible to predict the exact demand for each product since it is generated from the market, and it is difficult for the company to respond to the demand immediately due to the production time. In other words, it is difficult for poultry production companies to supply the products instantly after demand arises because it takes time for each product to be completed. For example, if the company starts raising chicks after the demand for large-sized chickens arises, it is difficult to satisfy the demand because it takes a considerable amount of time for the chicks to grow. Since poultry production companies produce products in response to market demand that is not determined in advance, it is appropriate to assume that demand occurs randomly. Although the models in the existing research can provide useful information on the production-planning problem of systems, research has not been conducted on cases of random demand. Our research is motivated by the deficiency of results. On the other hand, there are several studies that analyzed the supply chain considering random demand, but these studies did not consider the characteristics of production-time-dependent products [28–32]. In this paper, we analyze the supply chain considering both random demand and production-time-dependent product characteristics, and this is the difference between our study and previous studies.The problem of our interest is the input control problem that determines the amount of raw material input and the input interval for a poultry farm. In poultry farms, chicks are supplied in batches, and the raising process begins at timed intervals. Once the farms begin to raise a set of chicks, a few chicks in the process are shipped to manufacturing plants after these attain the appropriate size (weight) to satisfy the demand. The remaining chicks grow to the next product size to satisfy the demand for that size (production-time-dependent products). This procedure is repeated until the chicks grow fully. After the chicks are delivered to satisfy the demand for the final product type (the largest size), the remaining chicks are discarded. Once the raising process of a group of chicks starts, the farm cannot start to raise the next group of chicks due to the production density, feed, and temperature control. In other words, a new group of chicks cannot enter the system until the current chicks are fully grown and leave the system. Since only one group of chicks can be serviced at a time, the farm must determine the appropriate number of chicks to be fed to meet the demand for a variety of products over a period. If too many chicks are put in or if chicks are entered frequently, there is a high probability that more items will be produced than demanded, and waste can occur. Conversely, if a small number of chicks enter the system, the demand is likely to outweigh the supply and a loss in sales may occur. Thus, considering these two costs (disposal costs and costs for unfilled demand), the appropriate batch size and the interval at which the raising of batches of chicks should be started need to be determined. Our objective is to develop a stochastic model for the input control problem in poultry farms to seek a balance between these costs. Although there are several prior studies that have studied the operational problems of poultry farms, to our best knowledge, no stochastic models have been proposed to deal with the abovementioned problem. Embedded Markov chain models are employed for the analysis. The following are the results of our study:(i)
A recursive structure and a Markovian property for the number of raw materials (chicks) and unfulfilled demand for each type of product in the system are investigated.(ii)
Embedded Markov chain models are used for the number of raw materials (chicks) and unsatisfied demand (for each type) when each type of product demand occurs.(iii)
We demonstrate that the equilibrium probabilities for the Markov chain models could be obtained numerically using matrix analytic methods.(iv)
A probability-generating function can be obtained for the convolution distribution of unsatisfied demands.Through the proposed model, distributions regarding the number of discarded raw materials and the unmet demand of each type of product can be derived, which provides useful information on the relevant decision-making problems. In addition, our analytical model provides information on how distributions and related system metrics (such as mean disposal costs or mean opportunity costs) change according to system parameters. This helps to understand the system behavior theoretically.The remainder of this paper is organized as follows. In Section2, the input control problems in a two-level supply chain for production-time-dependent products are described. In addition, key notations are defined. Section 3 presents the results. We show that the system can be modeled as a Markov chain and that the equilibrium probabilities of the system can be derived. In Section 4, numerical experiments are presented. Finally, the concluding remarks are presented in Section 5.


## 2. Problem Statement

### 2.1. Input Control Problems in a Poultry Company

In this section, we describe our problems, which are based on a poultry company in Korea. Many companies in the poultry industry operate manufacturing plants and poultry farms. In manufacturing plants, demand for a variety of products is realized, and finished products are produced to meet this demand. Each product is manufactured from chickens of different sizes. This implies that different products require different semi-finished products (raw materials). After the plants review the number of chickens (of each size) in their stock, they order semi-finished products from poultry farms.Poultry farms supply semi-finished products to manufacturing plants. Thus, from the perspective of the supply chain, poultry farms play the role of suppliers that produce various types (weights) of chickens (semi-finished products). The weight (type) of a chicken is proportional to the length of time it has been raised on the farm (production time). Chickens raised on farms are larger in size and weight. Thus, after chicks are supplied to the farms, these transform into a variety of products depending on the length of time over which these are raised there. As mentioned in the previous section, products with such characteristics are termed as being production-time dependent. The semi-finished products are delivered to manufacturing plants to address the orders that the plants have. Because the farm and production plant are operated by the same company, the production plan of the farms is linked to that of the manufacturer. Thus, we can assume that the relevant information is shared fully and that a business organization controls all the units in the supply chain. This enables us to assume one aggregated poultry farm. Since one business organization has all the information about the farms and controls them, it estimates the demand for each product first and then decides how many of each product should be produced over a period of time. After the aggregate plan is finalized, the amount each farm must produce is determined based on the plan. Our study focuses on how the business organization determines its aggregate plan when the demand for each product is random. Therefore, in this study, one aggregated poultry farm is assumed. Note that our model can be applied to individual farms as well, given information about the demand occurring on each farm. However, in this case, other costs such as logistics costs and setup costs should be considered depending on the distance from the manufacturing plant in general. Thus, determining the production by varying the operation of each farm under the assumption of random demand is also an interesting topic, and it is beyond the scope of this study. This topic will be the future direction of this study.An important issue in the operation of a poultry company is the determination of the raising interval and quantity of input (chicks). Because the production density, feed, and temperature need to be controlled during the raising process, the process of raising a new set cannot be started until the production of the previous set is complete. A cleaning process is carried out after the process for a set is complete to prevent infection. New chicks can then enter this process. Because it is difficult to process multiple sets simultaneously, to respond to random demand, it is necessary to determine the appropriate raising interval and the number of chicks whose rearing needs to be started. When an excessively large number of chicks need to be reared in a set, disposal costs may be incurred because of an increase in the number of discarded chicks. Meanwhile, when an excessively small number of chicks are raised, the demand for a specific product may not be satisfied (opportunity cost is incurred). The opportunity cost of not satisfying demand in a timely manner differs for each product type (chicken size). The stochastic model proposed in this study aims to achieve a trade-off between the two costs.The following assumptions are made in this study and are based on a real system:(1)
The supply chain consists of a manufacturing plant and a poultry farm. A company operates multiple poultry farms in an integrated manner, and the farms share their information completely. Therefore, it can be assumed to be an aggregated poultry farm.(2)
A manufacturing plant produces multiple types of products. Each product requires different semi-finished products (chickens of different sizes). The output of each product is determined by the demand, which is random. Thus, the plant orders a random number of semi-finished products from the farm.(3)
The demand for each type of product is random. It is assumed that the demand that occurs for a specific type of product during a unit of time follows a general discrete distribution (independent and identically distributed).(4)
A farm that is a supplier produces multiple types of products to satisfy demand.(5)
The semi-finished products produced on the farm have properties of production-time-dependent products. Once a raw material is input to the process, it transforms into a different type of product depending on the production time. In the farm, after the chicks are received and rearing begins, these become chickens of different sizes depending on the length of time that they stay in the system.(6)
On the farm, the production of a batch of chicks begins simultaneously. Note that the products in a batch may result in different types of products depending on the time for which these remain in the process. In addition, a new set cannot enter the process until the production of the previous batch is complete.(7)
We assume that the farm should determine the batch size and time interval to start batch production. These correspond to the decision variables for the problem. A farm may decide to produce a batch size ofQ at time intervals of T.(8)
When a certain type of semi-finished product is produced, it is used immediately to satisfy the unmet demand that has accumulated until that point. The demand that is not satisfied by the produced item is carried over if the quantity produced is less than the accumulated demand (demand is not lost). If leftovers that satisfy the demand are available, these continue to grow and become the next type of product. After the demand for the final type of product is satisfied, the remaining product is discarded. Figure2 shows the behavior of the system.(9)
We consider two types of cost. The first type of cost is incurred when demand is not satisfied. It is assumed that unsatisfied demand carries over but incurs costs. The second type of cost is incurred when the product is discarded. In the case of excessive product input, products may remain after the demand for the final type of product is satisfied. Such products incur a disposal cost.Figure 2
Example of system behavior: case 1 (a); case 2 (b).
(a)(b)As mentioned in the previous section, the assumption of production-time-dependent product is new that cannot be found in previous studies, and other assumptions are often used in other inventory models. Note that we assume that demand can carry over. Figure2 depicts an example of system behavior over time. Consider the case in which a farm produces three types of products and inputs a batch consisting of Q chicks. After a certain amount of time μ1 has passed since the beginning of production, chicks become the first type of product and are used to satisfy the demand d1 that has accumulated until this point. The remaining products Q−d1 continue to grow and become the next type of products. In the example, after the batch is placed, time μi is required for the chicks to become units of type-i product. The demand for the type-i product is represented by di (in the example, i = 1, 2, 3). Note that the product is used immediately to satisfy the unsatisfied demand that has accumulated until that point as soon as a certain type of semi-finished product is completed. After the demands for all the product types are satisfied, the remaining products are discarded (case 1). If the demand is not satisfied by the current input, the unsatisfied demand is carried over and is addressed by the next input (case 2).


### 2.2. Notation and System Parameters

Note that the decision variables in this study are the input quantity and input interval. LetQ and T denote the input quantity (batch size) and interval, respectively. Thus, we place a batch of Q products into the process after every T time units. Our objective is to determine the optimal values of Q and T considering discard cost and opportunity cost.We assume that a supplier (farm) producesM types of products and supplies these to the manufacturing plant. After the batch is input, the type-i product is completed after the passage of time μi (production time). Thus, if the input time is t, the type-i product is completed after t+μi. For convenience, it is assumed that the first input occurs at t=0. Subsequently, the input of the k-th batch occurs at t = (k − 1) T. The demand for the type-i product during k−2T+μi,k−1T+μi is denoted as Dik. It is independent and identically distributed (i.i.d.) and follows the discrete probability distribution Di. Note that the system strives to satisfy the demand for a particular type of product immediately upon the completion of production. The mean of the random variable Di is denoted as λi. We assume a Q value that satisfies Q>∑l=1Mλl owing to the stability condition. Table 2 presents the variables used in this study.Table 2
Notation.
QInput quantity (batch size), decision variableTInput interval, decision variableMNumber of product typesμiTime required to produce type-i product (production time)DikDemands for type-i product during k−2T+μi,k−1T+μiDiDik–Di for all k. Distribution of demand for type-i productλiEDiaijPrDi=jbijPr∑l=1iDl=jXikNumber of raw materials (chicks) remaining in thek-th input immediately after time t=k−1T+μixijLimiting probability ofXik.xij≔limk⟶∞PXik=jj=0,…,QYikRemaining demand for type-i product immediately after t=k−1T+μi (demand not satisfied even with the k-th input product)yijLimiting probability ofYik.yij≔limk⟶∞PYik=jUikThe convolution ofY1k,Y2k,…,YMk. that is, Uik=∑j=1iYjkuijLimiting probability ofUik.uij≔limk⟶∞PUik=jciCost incurred when the demand for a unit of type-i product is not satisfiedcdisCost of disposing one unit of raw material


## 2.1. Input Control Problems in a Poultry Company

In this section, we describe our problems, which are based on a poultry company in Korea. Many companies in the poultry industry operate manufacturing plants and poultry farms. In manufacturing plants, demand for a variety of products is realized, and finished products are produced to meet this demand. Each product is manufactured from chickens of different sizes. This implies that different products require different semi-finished products (raw materials). After the plants review the number of chickens (of each size) in their stock, they order semi-finished products from poultry farms.Poultry farms supply semi-finished products to manufacturing plants. Thus, from the perspective of the supply chain, poultry farms play the role of suppliers that produce various types (weights) of chickens (semi-finished products). The weight (type) of a chicken is proportional to the length of time it has been raised on the farm (production time). Chickens raised on farms are larger in size and weight. Thus, after chicks are supplied to the farms, these transform into a variety of products depending on the length of time over which these are raised there. As mentioned in the previous section, products with such characteristics are termed as being production-time dependent. The semi-finished products are delivered to manufacturing plants to address the orders that the plants have. Because the farm and production plant are operated by the same company, the production plan of the farms is linked to that of the manufacturer. Thus, we can assume that the relevant information is shared fully and that a business organization controls all the units in the supply chain. This enables us to assume one aggregated poultry farm. Since one business organization has all the information about the farms and controls them, it estimates the demand for each product first and then decides how many of each product should be produced over a period of time. After the aggregate plan is finalized, the amount each farm must produce is determined based on the plan. Our study focuses on how the business organization determines its aggregate plan when the demand for each product is random. Therefore, in this study, one aggregated poultry farm is assumed. Note that our model can be applied to individual farms as well, given information about the demand occurring on each farm. However, in this case, other costs such as logistics costs and setup costs should be considered depending on the distance from the manufacturing plant in general. Thus, determining the production by varying the operation of each farm under the assumption of random demand is also an interesting topic, and it is beyond the scope of this study. This topic will be the future direction of this study.An important issue in the operation of a poultry company is the determination of the raising interval and quantity of input (chicks). Because the production density, feed, and temperature need to be controlled during the raising process, the process of raising a new set cannot be started until the production of the previous set is complete. A cleaning process is carried out after the process for a set is complete to prevent infection. New chicks can then enter this process. Because it is difficult to process multiple sets simultaneously, to respond to random demand, it is necessary to determine the appropriate raising interval and the number of chicks whose rearing needs to be started. When an excessively large number of chicks need to be reared in a set, disposal costs may be incurred because of an increase in the number of discarded chicks. Meanwhile, when an excessively small number of chicks are raised, the demand for a specific product may not be satisfied (opportunity cost is incurred). The opportunity cost of not satisfying demand in a timely manner differs for each product type (chicken size). The stochastic model proposed in this study aims to achieve a trade-off between the two costs.The following assumptions are made in this study and are based on a real system:(1)
The supply chain consists of a manufacturing plant and a poultry farm. A company operates multiple poultry farms in an integrated manner, and the farms share their information completely. Therefore, it can be assumed to be an aggregated poultry farm.(2)
A manufacturing plant produces multiple types of products. Each product requires different semi-finished products (chickens of different sizes). The output of each product is determined by the demand, which is random. Thus, the plant orders a random number of semi-finished products from the farm.(3)
The demand for each type of product is random. It is assumed that the demand that occurs for a specific type of product during a unit of time follows a general discrete distribution (independent and identically distributed).(4)
A farm that is a supplier produces multiple types of products to satisfy demand.(5)
The semi-finished products produced on the farm have properties of production-time-dependent products. Once a raw material is input to the process, it transforms into a different type of product depending on the production time. In the farm, after the chicks are received and rearing begins, these become chickens of different sizes depending on the length of time that they stay in the system.(6)
On the farm, the production of a batch of chicks begins simultaneously. Note that the products in a batch may result in different types of products depending on the time for which these remain in the process. In addition, a new set cannot enter the process until the production of the previous batch is complete.(7)
We assume that the farm should determine the batch size and time interval to start batch production. These correspond to the decision variables for the problem. A farm may decide to produce a batch size ofQ at time intervals of T.(8)
When a certain type of semi-finished product is produced, it is used immediately to satisfy the unmet demand that has accumulated until that point. The demand that is not satisfied by the produced item is carried over if the quantity produced is less than the accumulated demand (demand is not lost). If leftovers that satisfy the demand are available, these continue to grow and become the next type of product. After the demand for the final type of product is satisfied, the remaining product is discarded. Figure2 shows the behavior of the system.(9)
We consider two types of cost. The first type of cost is incurred when demand is not satisfied. It is assumed that unsatisfied demand carries over but incurs costs. The second type of cost is incurred when the product is discarded. In the case of excessive product input, products may remain after the demand for the final type of product is satisfied. Such products incur a disposal cost.Figure 2
Example of system behavior: case 1 (a); case 2 (b).
(a)(b)As mentioned in the previous section, the assumption of production-time-dependent product is new that cannot be found in previous studies, and other assumptions are often used in other inventory models. Note that we assume that demand can carry over. Figure2 depicts an example of system behavior over time. Consider the case in which a farm produces three types of products and inputs a batch consisting of Q chicks. After a certain amount of time μ1 has passed since the beginning of production, chicks become the first type of product and are used to satisfy the demand d1 that has accumulated until this point. The remaining products Q−d1 continue to grow and become the next type of products. In the example, after the batch is placed, time μi is required for the chicks to become units of type-i product. The demand for the type-i product is represented by di (in the example, i = 1, 2, 3). Note that the product is used immediately to satisfy the unsatisfied demand that has accumulated until that point as soon as a certain type of semi-finished product is completed. After the demands for all the product types are satisfied, the remaining products are discarded (case 1). If the demand is not satisfied by the current input, the unsatisfied demand is carried over and is addressed by the next input (case 2).


## 2.2. Notation and System Parameters

Note that the decision variables in this study are the input quantity and input interval. LetQ and T denote the input quantity (batch size) and interval, respectively. Thus, we place a batch of Q products into the process after every T time units. Our objective is to determine the optimal values of Q and T considering discard cost and opportunity cost.We assume that a supplier (farm) producesM types of products and supplies these to the manufacturing plant. After the batch is input, the type-i product is completed after the passage of time μi (production time). Thus, if the input time is t, the type-i product is completed after t+μi. For convenience, it is assumed that the first input occurs at t=0. Subsequently, the input of the k-th batch occurs at t = (k − 1) T. The demand for the type-i product during k−2T+μi,k−1T+μi is denoted as Dik. It is independent and identically distributed (i.i.d.) and follows the discrete probability distribution Di. Note that the system strives to satisfy the demand for a particular type of product immediately upon the completion of production. The mean of the random variable Di is denoted as λi. We assume a Q value that satisfies Q>∑l=1Mλl owing to the stability condition. Table 2 presents the variables used in this study.Table 2
Notation.
QInput quantity (batch size), decision variableTInput interval, decision variableMNumber of product typesμiTime required to produce type-i product (production time)DikDemands for type-i product during k−2T+μi,k−1T+μiDiDik–Di for all k. Distribution of demand for type-i productλiEDiaijPrDi=jbijPr∑l=1iDl=jXikNumber of raw materials (chicks) remaining in thek-th input immediately after time t=k−1T+μixijLimiting probability ofXik.xij≔limk⟶∞PXik=jj=0,…,QYikRemaining demand for type-i product immediately after t=k−1T+μi (demand not satisfied even with the k-th input product)yijLimiting probability ofYik.yij≔limk⟶∞PYik=jUikThe convolution ofY1k,Y2k,…,YMk. that is, Uik=∑j=1iYjkuijLimiting probability ofUik.uij≔limk⟶∞PUik=jciCost incurred when the demand for a unit of type-i product is not satisfiedcdisCost of disposing one unit of raw material


## 3. Analysis

In this section, the recursive structure of the system to obtain the embedded Markov chain model is introduced. After the model is obtained, several methods are used to derive the equilibrium probability distribution of the models. Performance metrics that can provide effective information while designing a system are also presented.


### 3.1. Recursive Structure

Recall that thek-th batch is input to the system at t = (k − 1) T. In addition, the production of the type-i product is completed in μi time units after the batch is input. It is used to satisfy the demand accumulated until then. We define the following notation:(i)
Xik: amount of raw material (number of chicks) remaining in the system immediately after time t=k−1T+μi.(ii)
Yik: remaining demand for type-i products immediately after t=k−1T+μi (demand not satisfied even with the k-th input product).For the first type of product,(1)X1k+1=max0,Q−Y1k+D1k+1,(2)Y1k+1=max0,Y1k+D1k+1−Q,for k≥1. After the (k + 1)-th batch (of Q units) is input, the units in the batch become units of the type-1 product at t=kT+μ1. Equation (1) indicates that the finished units of type-1 product are then used to address the demand that was not satisfied by the previous batch Y1k and the demand that occurs during k−1T+μ1,kT+μ1D1k+1.X1k+1=0 if the unsatisfied demand at time t=kT+μ1 exceeds QQ<Y1k+D1k+1. This is because all the products in the batch are used to address the demand. In this case, Y1k+1>0 because unsatisfied demand for the type-1 product exists immediately after t=kT+μ1. Note that the unsatisfied demand is carried over and is addressed with subsequent batches.Similarly, for the type-i product i>1,(3)Xik+1=max0,Xi−1k+1−Yik+Dik+1,(4)Yik+1=max0,Yik+Dik+1−Xi−1k+1,for k≥1. After the products in the (k + 1)-th batch are used to satisfy the demands for all the previous types, the remaining chicks Xi−1k+1 become type-i product to satisfy the demand for this product type Yik+Dik+1. When the number of chicks remaining is insufficient, only part of the demand for the type-i product is satisfied, and Xik+1=0 (equation (3)). In this case, the demand Yik+1 that is not satisfied by the (k + 1)-th batch is carried over, and the amount is determined by equation (4).A performance measure that is of interest is the amount of input discarded. After the demand for the final type of product is satisfied, the remaining input is discarded. This incurs disposal costs. Note thatXMk represents the amount of inventory remaining after the k-th input satisfies the demand. To analyze the limiting distribution of XMk, we define Uik=∑j=1iYjk. Then, the following can be derived based on equations (3) and (4):(5)Uik+1=max0,Uik+∑j=1iDjk+1−Q,(6)Xik+1=max0,Q−Uik+∑j=1iDjk+1.The proof for (5) and (6) is provided in Appendix A.Note thatUik+1 and Xik+1 are determined by Uik. Furthermore, the convolution of D1,…,andDi.Uik+1 obeys a Markovian property and has discrete values. Thus, an embedded Markov chain model can be constructed. Its embedding points are t=mT+μim=1,2,… immediately after the demand for type-i products occurs. The Markov chain is denoted by Uik,k≥1. Its steady-state probability is defined as(7)uij=limk⟶∞PUik=j,i=1,…,M.Then, the average amount of input discarded can be expressed as follows:(8)limk⟶∞EXMk=∑i+j≤QQ−i−juMibMj,where bij represents the convolution of probability distributions of D1,D2,…,Di. That is, bij=P∑j=1iDk=j. In addition, the limiting distribution of Xik can be derived as follows using the Markov chain model:(9)xi0=∑m+n≥Quimbin,(10)xiQ−l=∑m+n=l<Quimbin,where xij≔limk⟶∞PXik=jj=0,…,Q.Another performance measure of interest is the distribution of unsatisfied demand for each type of product. It can be derived from equation (4):(11)Yik+1=max0,Yik+Dik+1−Xi−1k+1.The unsatisfied demand is affected by two factors: the remaining number of (k + 1)-th input Xi−1k+1 immediately after time t=kT+μi−1 and the demands that occurred during k−1T+μi,kT+μi. Assuming that the systems are in a steady state, the above recursive structure enables us to derive the limiting distribution of unsatisfied demands.


### 3.2. Derivation of Equilibrium Probabilities

This section presents the procedure for deriving the equilibrium probabilities of the embedded Markov chain. The probabilities are used to calculate the performance measures of the system.First, we analyzeUik,k≥1. The state space of the Markov chain is 0,1,…. Its one-step transition probability matrix is given by(12)piU=BA2A3A4⋯A0A1A2A3⋯OQA0A1A2⋯OQOQA0A1⋯OQOQOQA0⋯⋮⋮⋮OQ⋯⋮⋮⋮⋮⋯,whereB=∑j=0QbijbiQ+1biQ+2biQ+3⋯bi2Q−1∑j=0Q−1bijbiQbiQ+1biQ+2⋯bi2Q−2∑j=0Q−2bijbiQ−1biQbiQ+1⋯bi2Q−3∑j=0Q−3bijbiQ−2biQ−1biQ⋯bi2Q−4⋮⋮⋮⋮⋯⋮⋮⋮⋮⋮⋯⋮∑j=01bijbi2bi3bi4⋯biQ,A0=bi0bi1bi2bi3⋯biQ−10bi0bi1bi2⋯biQ−200bi0bi1⋯biQ−3000bi0⋯biQ−40000⋯biQ−5⋮⋮⋮⋮⋯⋮0000⋯bi0,Am=bimQbimQ+1bimQ+2bimQ+3⋯bim+1Q−1bimQ−1bimQbimQ+1bimQ+2⋯bim+1Q−2bimQ−2bimQ−1bimQbimQ+1⋯bim+1Q−3bimQ−3bimQ−2bimQ−1bimQ⋯bim+1Q−4bimQ−4bimQ−3bimQ−2bimQ−1⋯bim+1Q−5⋮⋮⋮⋮⋯⋮bim−1Q+1bim−1Q+2bim−1Q+3bim−1Q+4⋯bimQform=1,2,….OQ is a Q×Q matrix with all elements equal to zero. The submatrices of PiU are Q×Q matrices. The matrix can be obtained from recursive equation (5). For example, the first row of the matrix gives the probabilities of transitioning from state 0 to another in a single time unit. When Uik=0, if ∑j=1iDjk+1≤Q (with probability ∑j=0Qbij, Uik+1 becomes 0 (the first element of the first row). With probability biQ+1,∑j=1iDjk+1=Q+1; then, Uik+1=1 (the second element of the first row). The other elements of the matrix (transition probabilities) can also be obtained from equation (5) in a similar way. The structure of the matrix is called the upper Hessenberg matrix and is similar to that of an M/G/1-type Markov chain in [33]. Hence, two well-known approaches can be employed to calculate the equilibrium probabilities: the matrix analytic methods in [34] and the root-finding methods in [35].Equilibrium probabilities can be calculated numerically by employing matrix analytic methods. LetR be the smallest nonnegative solution of the nonlinear matrix equation:(13)R=∑i=0∞AiRi.Then, the equilibrium probability vector of the Markov chain,π, obeys Ramaswami’s formula:(14)πi=π0B¯i+∑j=1i−1πiA¯i+1−jI−A¯1−1,for i>0, with A¯i=∑j≥iAjRj−i and B¯i=∑j≥iAj+1Rj−i. Various algorithms have been developed to calculate R and the initial condition vector. For the algorithms, refer to [34] (functional iterations methods), [36] (Newton iterations), [37] (cyclic reduction), and [38] (invariant subspace). These algorithms can be used to calculate the probabilities.In addition, root-finding methods can be employed. By derivingπ=πP, we observe that the balance equations have the following recursive structure:(15)pi1=pi0IQ−BA0−1,(16)pik+1=pikIQ−A1−∑m=1kpik−mA1+mA0−1,where pik is a row vector comprising uikQ,uikQ+1,…,uik+1Q−1. That is, uik≔uikQ,uikQ+1,…,uik+1Q−1. IQ is a Q×Q identity matrix. Note that the equilibrium probabilities pik+1 in (16) are expressed in terms of pim for m≤k. Thus, the equilibrium probabilities can be calculated recursively if the values pi0 are known.The initial conditionpi0 can be derived from the probability generating function (PGF) in certain cases. The balance equations are as follows:(17)ui0=∑m=0Quim∑n=0Q−mbin,uij=∑m=0Q+juimbiQ+j−m.We obtain the following results by deriving the PGF of the Markov chainUi∗z:(18)Ui∗z=zQ∑m=0Q−1uim∑n=0Q−mbin−∑m=0Q−1uim∑n=0Q−mzm+nbinzQ−Biz,where Biz=∑n=0∞binzn. There are Q unknown probabilities in the numerator, ui0,ui1,…,uiQ−1. These are elements of pi0. The unknowns need to be evaluated to uniquely determine the PGF.We employ Rouche’s theorem in [39] to determine the unknowns. Let fz and giz denote zQ and Biz, respectively. In addition, we set C=z:z=1+ε. Then,(19)fz=zQ=1+εQ=1+Qε+oε=1+Qε+oε,giz=∑n=0∞binzn≤∑n=0∞binzn=∑n=0∞bin1+εn=∑n=0∞bin1+kε+oε=1+ε∑l=1iλl+oε.Note thatfz>giz for all i because of the stability condition Q>∑l=1Mλl. Because the Q roots of fz are zero (all are inside C), by Rouche’s theorem, all the Q zeros of zQ−Biz are inside C. When ε⟶0, zQ−Biz has Q roots on and inside the unit circle. Because the PGF converges in the unit circle, the solutions inside the unit circle that render the denominator equal to zero should be the roots of the numerator. Hence, the Q linear equations can be obtained by substituting the roots of the denominator into the numerator, and the initial conditions can be obtained by solving these equations. For the equations to have a unique solution, the Q roots must be different. Appendix B shows that the roots of the denominator of the PGF are different when ∑k=1iDk follows Poisson, binomial, and geometric distributions. For the other distributions, we can determine ui0,ui1,…,uiQ−1 when zQ−Biz=0 has Q distinct roots. Once the initial condition is determined, we can calculate the equilibrium probabilities by employing equations (13) and (14).Next, we analyzeYik. Note that Y1k = U1k. The limiting distribution of Y1k can be derived using this method. Here, we consider the limiting distribution of Yik with i>1. Rewriting equation (4),(20)Yik+1=max0,Yik+Dik+1−Xi−1k+1.We can obtain the following by deriving the probability mass function ofDik+1−Xi−1k+1:(21)gij≔limk⟶∞PDik+1−Xi−1k+1=j,j=−Q,−Q+1,…=∑l=0Q+jailxi−1l−j,if−Q≤j<0,=∑l=0Qaij+lxi−1l,ifj≥0.Subsequently, the state space of the Markov chain (Yik,k≥1) is 0,1,…. Furthermore, its one-step transition probability matrix is(22)piY=HG2G3G4⋯G0G1G2G3⋯OQG0G1G2⋯OQOQG0G1⋯OQOQOQG0⋯⋮⋮⋮OQ⋯⋮⋮⋮⋮⋯,H=∑j=−Q0gijgi1gi2gi3⋯giQ−1∑j=−Q−1gijgi0gi1gi2⋯giQ−2∑j=−Q−2gijgi−1gi0gi1⋯giQ−3∑j=−Q−3gijgi−2gi−1gi0⋯giQ−4⋮⋮⋮⋮⋯⋮⋮⋮⋮⋮⋯⋮∑j=−Q−Q+1gijgi−Q+2gi−Q+3gi−Q+4⋯gi0,G0=gi−Qgi−Q+1gi−Q+2gi−Q+3⋯gi−10gi−Qgi−Q+1gi−Q+2⋯gi−200gi−Qgi−Q+1⋯gi−3000gi−Q⋯gi−40000⋯gi−5⋮⋮⋮⋮⋯⋮0000⋯gi−Q,Gm=gim−1Qgim−1Q+1gim−2Q+2gim−2Q+3⋯gimQ−1gim−1Q−1gim−1Qgim−2Q+1gim−2Q+2⋯gimQ−2gim−1Q−2gim−1Q−1gim−1Qgim−2Q+1⋯gimQ−3gim−1Q−3gim−1Q−2gim−1Q−1gim−1Q⋯gimQ−4gim−1Q−4gim−1Q−3gim−1Q−2gim−1Q−1⋯gimQ−5⋮⋮⋮⋮⋯⋮gim−2Q+1gim−2Q+2gim−2Q+3gim−2Q+4⋯gim−1Qform=1,….Note that the structure of the matrix is identical to that ofPiU. Of the two well-known approaches for calculating equilibrium probabilities, it is difficult to employ the root-finding methods in [35] to Yik,k≥1. This is because unlike Uik,k≥1, we cannot ensure that the initial probabilities can be obtained explicitly. It is difficult to obtain an explicit expression for the limiting probabilities of Dik+1−Xi−1k+1. Hence, Rouche’s theorem is inapplicable. Nonetheless, the equilibrium probabilities of Yik,k≥1 can be calculated numerically using the matrix analysis method. Equations (13) and (14) can be modified conveniently for PiY. The procedure for calculating the steady-state probabilities of PiY is identical to that for calculating the probability of PiU.


## 3.1. Recursive Structure

Recall that thek-th batch is input to the system at t = (k − 1) T. In addition, the production of the type-i product is completed in μi time units after the batch is input. It is used to satisfy the demand accumulated until then. We define the following notation:(i)
Xik: amount of raw material (number of chicks) remaining in the system immediately after time t=k−1T+μi.(ii)
Yik: remaining demand for type-i products immediately after t=k−1T+μi (demand not satisfied even with the k-th input product).For the first type of product,(1)X1k+1=max0,Q−Y1k+D1k+1,(2)Y1k+1=max0,Y1k+D1k+1−Q,for k≥1. After the (k + 1)-th batch (of Q units) is input, the units in the batch become units of the type-1 product at t=kT+μ1. Equation (1) indicates that the finished units of type-1 product are then used to address the demand that was not satisfied by the previous batch Y1k and the demand that occurs during k−1T+μ1,kT+μ1D1k+1.X1k+1=0 if the unsatisfied demand at time t=kT+μ1 exceeds QQ<Y1k+D1k+1. This is because all the products in the batch are used to address the demand. In this case, Y1k+1>0 because unsatisfied demand for the type-1 product exists immediately after t=kT+μ1. Note that the unsatisfied demand is carried over and is addressed with subsequent batches.Similarly, for the type-i product i>1,(3)Xik+1=max0,Xi−1k+1−Yik+Dik+1,(4)Yik+1=max0,Yik+Dik+1−Xi−1k+1,for k≥1. After the products in the (k + 1)-th batch are used to satisfy the demands for all the previous types, the remaining chicks Xi−1k+1 become type-i product to satisfy the demand for this product type Yik+Dik+1. When the number of chicks remaining is insufficient, only part of the demand for the type-i product is satisfied, and Xik+1=0 (equation (3)). In this case, the demand Yik+1 that is not satisfied by the (k + 1)-th batch is carried over, and the amount is determined by equation (4).A performance measure that is of interest is the amount of input discarded. After the demand for the final type of product is satisfied, the remaining input is discarded. This incurs disposal costs. Note thatXMk represents the amount of inventory remaining after the k-th input satisfies the demand. To analyze the limiting distribution of XMk, we define Uik=∑j=1iYjk. Then, the following can be derived based on equations (3) and (4):(5)Uik+1=max0,Uik+∑j=1iDjk+1−Q,(6)Xik+1=max0,Q−Uik+∑j=1iDjk+1.The proof for (5) and (6) is provided in Appendix A.Note thatUik+1 and Xik+1 are determined by Uik. Furthermore, the convolution of D1,…,andDi.Uik+1 obeys a Markovian property and has discrete values. Thus, an embedded Markov chain model can be constructed. Its embedding points are t=mT+μim=1,2,… immediately after the demand for type-i products occurs. The Markov chain is denoted by Uik,k≥1. Its steady-state probability is defined as(7)uij=limk⟶∞PUik=j,i=1,…,M.Then, the average amount of input discarded can be expressed as follows:(8)limk⟶∞EXMk=∑i+j≤QQ−i−juMibMj,where bij represents the convolution of probability distributions of D1,D2,…,Di. That is, bij=P∑j=1iDk=j. In addition, the limiting distribution of Xik can be derived as follows using the Markov chain model:(9)xi0=∑m+n≥Quimbin,(10)xiQ−l=∑m+n=l<Quimbin,where xij≔limk⟶∞PXik=jj=0,…,Q.Another performance measure of interest is the distribution of unsatisfied demand for each type of product. It can be derived from equation (4):(11)Yik+1=max0,Yik+Dik+1−Xi−1k+1.The unsatisfied demand is affected by two factors: the remaining number of (k + 1)-th input Xi−1k+1 immediately after time t=kT+μi−1 and the demands that occurred during k−1T+μi,kT+μi. Assuming that the systems are in a steady state, the above recursive structure enables us to derive the limiting distribution of unsatisfied demands.


## 3.2. Derivation of Equilibrium Probabilities

This section presents the procedure for deriving the equilibrium probabilities of the embedded Markov chain. The probabilities are used to calculate the performance measures of the system.First, we analyzeUik,k≥1. The state space of the Markov chain is 0,1,…. Its one-step transition probability matrix is given by(12)piU=BA2A3A4⋯A0A1A2A3⋯OQA0A1A2⋯OQOQA0A1⋯OQOQOQA0⋯⋮⋮⋮OQ⋯⋮⋮⋮⋮⋯,whereB=∑j=0QbijbiQ+1biQ+2biQ+3⋯bi2Q−1∑j=0Q−1bijbiQbiQ+1biQ+2⋯bi2Q−2∑j=0Q−2bijbiQ−1biQbiQ+1⋯bi2Q−3∑j=0Q−3bijbiQ−2biQ−1biQ⋯bi2Q−4⋮⋮⋮⋮⋯⋮⋮⋮⋮⋮⋯⋮∑j=01bijbi2bi3bi4⋯biQ,A0=bi0bi1bi2bi3⋯biQ−10bi0bi1bi2⋯biQ−200bi0bi1⋯biQ−3000bi0⋯biQ−40000⋯biQ−5⋮⋮⋮⋮⋯⋮0000⋯bi0,Am=bimQbimQ+1bimQ+2bimQ+3⋯bim+1Q−1bimQ−1bimQbimQ+1bimQ+2⋯bim+1Q−2bimQ−2bimQ−1bimQbimQ+1⋯bim+1Q−3bimQ−3bimQ−2bimQ−1bimQ⋯bim+1Q−4bimQ−4bimQ−3bimQ−2bimQ−1⋯bim+1Q−5⋮⋮⋮⋮⋯⋮bim−1Q+1bim−1Q+2bim−1Q+3bim−1Q+4⋯bimQform=1,2,….OQ is a Q×Q matrix with all elements equal to zero. The submatrices of PiU are Q×Q matrices. The matrix can be obtained from recursive equation (5). For example, the first row of the matrix gives the probabilities of transitioning from state 0 to another in a single time unit. When Uik=0, if ∑j=1iDjk+1≤Q (with probability ∑j=0Qbij, Uik+1 becomes 0 (the first element of the first row). With probability biQ+1,∑j=1iDjk+1=Q+1; then, Uik+1=1 (the second element of the first row). The other elements of the matrix (transition probabilities) can also be obtained from equation (5) in a similar way. The structure of the matrix is called the upper Hessenberg matrix and is similar to that of an M/G/1-type Markov chain in [33]. Hence, two well-known approaches can be employed to calculate the equilibrium probabilities: the matrix analytic methods in [34] and the root-finding methods in [35].Equilibrium probabilities can be calculated numerically by employing matrix analytic methods. LetR be the smallest nonnegative solution of the nonlinear matrix equation:(13)R=∑i=0∞AiRi.Then, the equilibrium probability vector of the Markov chain,π, obeys Ramaswami’s formula:(14)πi=π0B¯i+∑j=1i−1πiA¯i+1−jI−A¯1−1,for i>0, with A¯i=∑j≥iAjRj−i and B¯i=∑j≥iAj+1Rj−i. Various algorithms have been developed to calculate R and the initial condition vector. For the algorithms, refer to [34] (functional iterations methods), [36] (Newton iterations), [37] (cyclic reduction), and [38] (invariant subspace). These algorithms can be used to calculate the probabilities.In addition, root-finding methods can be employed. By derivingπ=πP, we observe that the balance equations have the following recursive structure:(15)pi1=pi0IQ−BA0−1,(16)pik+1=pikIQ−A1−∑m=1kpik−mA1+mA0−1,where pik is a row vector comprising uikQ,uikQ+1,…,uik+1Q−1. That is, uik≔uikQ,uikQ+1,…,uik+1Q−1. IQ is a Q×Q identity matrix. Note that the equilibrium probabilities pik+1 in (16) are expressed in terms of pim for m≤k. Thus, the equilibrium probabilities can be calculated recursively if the values pi0 are known.The initial conditionpi0 can be derived from the probability generating function (PGF) in certain cases. The balance equations are as follows:(17)ui0=∑m=0Quim∑n=0Q−mbin,uij=∑m=0Q+juimbiQ+j−m.We obtain the following results by deriving the PGF of the Markov chainUi∗z:(18)Ui∗z=zQ∑m=0Q−1uim∑n=0Q−mbin−∑m=0Q−1uim∑n=0Q−mzm+nbinzQ−Biz,where Biz=∑n=0∞binzn. There are Q unknown probabilities in the numerator, ui0,ui1,…,uiQ−1. These are elements of pi0. The unknowns need to be evaluated to uniquely determine the PGF.We employ Rouche’s theorem in [39] to determine the unknowns. Let fz and giz denote zQ and Biz, respectively. In addition, we set C=z:z=1+ε. Then,(19)fz=zQ=1+εQ=1+Qε+oε=1+Qε+oε,giz=∑n=0∞binzn≤∑n=0∞binzn=∑n=0∞bin1+εn=∑n=0∞bin1+kε+oε=1+ε∑l=1iλl+oε.Note thatfz>giz for all i because of the stability condition Q>∑l=1Mλl. Because the Q roots of fz are zero (all are inside C), by Rouche’s theorem, all the Q zeros of zQ−Biz are inside C. When ε⟶0, zQ−Biz has Q roots on and inside the unit circle. Because the PGF converges in the unit circle, the solutions inside the unit circle that render the denominator equal to zero should be the roots of the numerator. Hence, the Q linear equations can be obtained by substituting the roots of the denominator into the numerator, and the initial conditions can be obtained by solving these equations. For the equations to have a unique solution, the Q roots must be different. Appendix B shows that the roots of the denominator of the PGF are different when ∑k=1iDk follows Poisson, binomial, and geometric distributions. For the other distributions, we can determine ui0,ui1,…,uiQ−1 when zQ−Biz=0 has Q distinct roots. Once the initial condition is determined, we can calculate the equilibrium probabilities by employing equations (13) and (14).Next, we analyzeYik. Note that Y1k = U1k. The limiting distribution of Y1k can be derived using this method. Here, we consider the limiting distribution of Yik with i>1. Rewriting equation (4),(20)Yik+1=max0,Yik+Dik+1−Xi−1k+1.We can obtain the following by deriving the probability mass function ofDik+1−Xi−1k+1:(21)gij≔limk⟶∞PDik+1−Xi−1k+1=j,j=−Q,−Q+1,…=∑l=0Q+jailxi−1l−j,if−Q≤j<0,=∑l=0Qaij+lxi−1l,ifj≥0.Subsequently, the state space of the Markov chain (Yik,k≥1) is 0,1,…. Furthermore, its one-step transition probability matrix is(22)piY=HG2G3G4⋯G0G1G2G3⋯OQG0G1G2⋯OQOQG0G1⋯OQOQOQG0⋯⋮⋮⋮OQ⋯⋮⋮⋮⋮⋯,H=∑j=−Q0gijgi1gi2gi3⋯giQ−1∑j=−Q−1gijgi0gi1gi2⋯giQ−2∑j=−Q−2gijgi−1gi0gi1⋯giQ−3∑j=−Q−3gijgi−2gi−1gi0⋯giQ−4⋮⋮⋮⋮⋯⋮⋮⋮⋮⋮⋯⋮∑j=−Q−Q+1gijgi−Q+2gi−Q+3gi−Q+4⋯gi0,G0=gi−Qgi−Q+1gi−Q+2gi−Q+3⋯gi−10gi−Qgi−Q+1gi−Q+2⋯gi−200gi−Qgi−Q+1⋯gi−3000gi−Q⋯gi−40000⋯gi−5⋮⋮⋮⋮⋯⋮0000⋯gi−Q,Gm=gim−1Qgim−1Q+1gim−2Q+2gim−2Q+3⋯gimQ−1gim−1Q−1gim−1Qgim−2Q+1gim−2Q+2⋯gimQ−2gim−1Q−2gim−1Q−1gim−1Qgim−2Q+1⋯gimQ−3gim−1Q−3gim−1Q−2gim−1Q−1gim−1Q⋯gimQ−4gim−1Q−4gim−1Q−3gim−1Q−2gim−1Q−1⋯gimQ−5⋮⋮⋮⋮⋯⋮gim−2Q+1gim−2Q+2gim−2Q+3gim−2Q+4⋯gim−1Qform=1,….Note that the structure of the matrix is identical to that ofPiU. Of the two well-known approaches for calculating equilibrium probabilities, it is difficult to employ the root-finding methods in [35] to Yik,k≥1. This is because unlike Uik,k≥1, we cannot ensure that the initial probabilities can be obtained explicitly. It is difficult to obtain an explicit expression for the limiting probabilities of Dik+1−Xi−1k+1. Hence, Rouche’s theorem is inapplicable. Nonetheless, the equilibrium probabilities of Yik,k≥1 can be calculated numerically using the matrix analysis method. Equations (13) and (14) can be modified conveniently for PiY. The procedure for calculating the steady-state probabilities of PiY is identical to that for calculating the probability of PiU.


## 4. Numerical Experiments

The numerical experiments are described in this section. Recall that the structure of the matricesPiU and PiU is the same as that of the M/G/1-type Markov chain. The equilibrium probabilities of the M/G/1 type Markov chain can be calculated numerically by employing the matrix analytic method [33]. To implement the matrix analytic method, several algorithms have been introduced to calculate the probabilities in the previous literature: functional iteration methods [34], Newton iterations [36], cyclic reduction [37], and invariant subspace [38]. Methods for implementing each algorithm are presented in [39]. Here, we implemented a cycle reduction method that is presented in [37, 40].First, we investigated the distribution ofXi (k) according to the variation in the loading. Cases 1–4 are to examine the distribution of the remaining amount of input after the demand for each product is satisfied as the total demand rate changes. The loading is defined as the ratio of demand to supply or capacity. Note that here the average of the total demand for each product during T is ∑i=1Mλi, and the quantity supplied during T is Q. Thus, loading is expressed as ∑i=1Mλi/Q. Note that Xi (k) represents the remaining number of raw materials used to satisfy demand and the state space is 0,1,⋯,Q. Figure 3 depicts the distribution of Xi (k) according to loading. Xi in the figure indicates the limiting distribution of Xi (k). The limiting distribution of Xi (k) is a discrete probability distribution. However, to illustrate the shape of the graph and express several graphs in one figure, we connect each probability value and express it smoothly. We assume five types of products (M = 5, T = 1, and Q = 50). Then, the loading can be calculated by ∑i=15λi/Q (total average demand during T divided by the amount of input during that time). In addition, in the example, it is assumed that the demand for all the types of products follows a Poisson distribution with an equal rate in each case (the case with different demands for each type of product is addressed in the following numerical example). For example, in cases 1 and 2, if λ=λ1λ2λ3λ4λ5 is defined, λ=2.52.52.52.52.5 and 55555, respectively. As the demand for each type of product increases, the amount of raw input materials tends to decrease. Thereby, the limiting distribution of Xik shifts to the left as i increases (cases 1 and 2). In addition, the distribution flattens as i increases. This is because the demand arrival process is a Poisson process. Note that in the case of a Poisson distribution, the distribution flattens as the mean increases. As the loading increases, the distributions tend to shift to the left overall, and the probability that the remaining input becomes zero (when the accumulated demand is higher than the input) increases. In cases 3 and 4, as the graph approaches zero, the probability tends to decrease and then, increases from zero. In case 4, the probability that the remaining input would be zero after all the types of demand have been satisfied, x50, is 0.3355. In case 3, x50 is 0.0305. It may appear to be effective to set the input as close to the demand as feasible (considering the cost of disposing the input that remains after the entire demand is satisfied). However, the case in which demand is not satisfied should be considered. In case 4, x40 is 0.016. This indicates the likelihood that the remaining input would be zero after the demand for the type-4 product is satisfied. That is, it is likely that the demand for the type-5 product would not be satisfied. The unsatisfied demand is presented as an extreme distribution of Yi (k). This implies that both the distributions should be considered in combination.Figure 3
Limiting distribution ofXi (k) according to loading. (a) Case 1 (loading = 0.25). (b) Case 2 (loading = 0.5). (c) Case 3 (loading = 0.75). (d) Case 4 (loading = 0.9).
(a)(b)(c)(d)We can also examine the distribution ofXi (k) with variations in Q and T. Note that Q and T are decision variables. Cases 5–8 show that when the average demand and distribution of each product are estimated, the distribution can be derived according to the value of the decision variable. Figure 4 shows four examples. Note that cases 5 and 6 have equal loading values. If the unit of T is day, case 5 is to input 50 chicks into the process each day, and case 6 is to input 100 chicks every two days. As is evident from the graph, the distribution is different even for an equal loading value. When only the cost of disposal is considered, it is more effective to input a small amount more frequently than to input a large amount less frequently. x50 is 0.0305 in case 5, whereas the probability is 0.0034. The values of x50 for cases 7 and 8 are 0.3355 and 0.1953, respectively.Figure 4
Limiting distribution ofXi (k) with different Q and T. (a) Case 5 (loading = 0.75, T = 1, and Q = 50). (b) Case 6 (loading = 0.75, T = 2, and Q = 100). (c) Case 7 (loading = 0.9, T = 1, and Q = 50). (d) Case 8 (loading = 0.9, T = 2, and Q = 100).
(a)(b)(c)(d)Next, we analyze the limiting distribution ofYi (k). Examples of Y distribution are shown in Figure 5. Cases 9, 10, and 11 assume λ=55555, and case 12 assumes λ=1591317. From cases 9, 10, and 11, it is unfavorable to input chicks at small time intervals. The graph for Y5 (k) (marked as y5 in the graph) reveals that unsatisfied demand increases as the input interval decreases. Therefore, in contrast to the observation in the distribution of X, it is desirable to input a large number of chicks at long time intervals considering the opportunity cost incurred when the demand is not satisfied. In addition, a comparison of cases 11 and 12 shows that although the average total demand is identical, the results may involve different performance measures if the demand distribution for each type of product is different. In Case 11, the average unsatisfied demand for the final type of product, limk⟶∞EY5k, is 2.1105. In case 12, the average results are 2.8653. Therefore, even if the sum of the averages of the individual product demand is identical, different performance metric values can be obtained according to the parameters of the individual product demand distribution.Figure 5
Examples of limiting distribution ofYi (k). (a) Case 9 (Q = 50, T = 1). (b) Case 10 (Q = 25, T = 0.5). (c) Case 11 (Q = 10, T = 0.2). (d) Case 12 (Q = 10, T = 0.2) (different demand rates).
(a)(b)(c)(d)Note that the discarded quantity,XMk, is determined by the sum of the averages of the individual demand distributions (equations (8)–(10)). The limiting distribution of XMk is determined by the distribution of aggregate demand rather than by that of the demand for individual products. Meanwhile, Yik is affected by individual distributions. This implies that the distributions of individual product types should be considered when a trade-off is to be achieved between the two costs.We define the total cost function as follows to achieve a trade-off between disposal and opportunity costs (when demand is not satisfied):(23)totalcost=limk⟶∞∑i=1MciEYik+cdisEXMk,where ci is the cost incurred when the demand for a unit of type-i product is not satisfied and cdis is the cost of disposing one unit of raw material. We define Ctotal≔c1c2…cMcdis as the cost. The values of the cost function according to the cost and the values of Q and T are shown in Figures 6 and 7.Figure6 shows the values of the total cost function according to Q for the cases with different cost values. In all the cases, it is assumed that the vector representing demand λ is [5 5 5 5 5] (T = 1). However, the opportunity cost of each type of product and cost of discarding the remaining raw materials vary across cases. The vectors Ctotal in cases 13, 14, 15, and 16 are [5 5 5 5 5 5], [1 3 5 7 9 5], [9 7 5 3 1 5], and [5 5 5 5 5 20], respectively. Comparing cases 13, 14, and 15, case 14 has the highest total cost. This is considered to be owing to the large opportunity cost for the final type of product. Over time, raw materials transform into the next type of product while satisfying the demand of each type. During this time, the quantity of remaining raw materials reduces as the demand for the previous type of product is satisfied. Thereby, the probability of shortage of raw materials at the time of production of the next type of product increases. Therefore, case 14 (where the opportunity cost of the final product, c5, is higher than that for the other cases) incurs a higher cost. As Q increases, the average number of discarded raw materials increases, and the probability of not satisfying demand decreases. Thus, the results provide information regarding the Q value that balances the two costs given a demand, a cost vector, and T. In addition, a comparison of cases 13 and 16 reveals that when the disposal cost is high, the increment in total cost increases as Q increases beyond its optimal value.Figure 6
Total cost function according toci and cdis.BothQ and T should be determined to make an optimal decision. Figure 7 shows the value of the cost function according to the Q value for the cases with different T values. Note that λ=55555 and Ctotal=555555. It can be observed that in the example, the cost is lower than those for the other cases when T = 0.5. Case 17 (T = 0.5) has the lowest cost when Q is 15, and the total cost is 17.3461. The minimum costs for cases 18 and 19 are 24.3885 and 35.1237, respectively. These occur when Q is 28 and 55, respectively. The values of T and Q that provide the minimum cost can be calculated using our model given the demand and cost.Figure 7
Total cost function according toT and Q.


## 5. Concluding Remarks

In this study, we present a stochastic model to analyze the input control problem of the poultry farming process. Poultry farms produce multiple types (weights) of semi-finished products (chicken) and supply these to manufacturing plants. The type of product depends on the length of time it has spent in the process. A product with such a characteristic is called a production-time-dependent product. Because an important issue in the operation of a poultry company is the determination of the input interval and input quantity considering random demands, a stochastic model is proposed for the problem.From a theoretical perspective, our mathematical model is the first model that considers both the characteristics of production-time-dependent products and the random demand. Although there are several previous studies on the supply chain for production-time-dependent products, to our best knowledge, no mathematical models have been proposed to deal with the random demands. Since the demand for each product comes from the market, the assumption of random demand is essential in real poultry farms. We demonstrated that embedded Markov chain models can be obtained by identifying a recursive structure and Markovian property in the system. The recursive structure can be used not only to calculate steady-state probabilities but also for simple simulations to understand the system behavior. Also, it was found that the structure of the Markov chains is identical to that of the M/G/1-type Markov chain. Based on the structure, the steady-state probability can be calculated by two well-known methods: probability generating function method and matrix analysis method. We also proved that the probability generating function can be uniquely determined.From a practical perspective, our model can provide useful information on the decision-making problems in the poultry supply chain. Companies that supply poultry products to the market must control their farms and direct production orders. At this time, the company should determine how many of each product should be produced over a period before the demands are realized. Note that it takes a considerable amount of time as the chicks grow into each type of product. After the aggregate plan is decided, the amount each farm must produce is determined based on the plan. Our model provides information on how the company determines its aggregate plan when the demand for each product is random. Note that the equilibrium probabilities completely characterize the steady-state behavior of the system. Various performance measures such as discarded quantity, unmet demand, and total cost can be calculated using the equilibrium probabilities, and how the values of these measures change according to system parameter values can be analyzed. This could help the business organization issuing production orders to the farms.This study could be extended in several ways. One direction is to validate and improve our model based on actual company data. While the present model suggests methods for obtaining the equilibrium probabilities and the performance measures, the closed-form expression for the performance measures has not been derived. If the distribution of demand can be specified based on the actual data rather than assuming a general distribution, it may be possible to obtain an analytic solution for the steady-state probability or a closed-form expression for the performance measure. Also, the data-driven models to determine theQ and T values can be developed based on the actual demand. Another direction is to incorporate the various features or constraints of the practical system in the model. For example, there are no constraints on T in this study. However, there may be constraints on T because various operations such as cleaning are performed after input. This requires collaboration with a company. Also, chicks in the raising process may die due to disease or stress in practice. This means that a part of the input may be lost (i.e., yield ≠ 100%). Because the event wherein part of the input is lost also occurs randomly, the model can be improved by incorporating these issues. Furthermore, unsatisfied demand is carried over, but the costs associated with unmet demand are not considered in the present model. Taking these costs into account would also be an extension of the model. Finally, multiple farms with different behaviors may be considered. Note that one aggregated poultry farm is assumed in this study. After the aggregate plan is determined, the amount each farm must produce is determined based on the plan. However, the company operates multiple farms, and each farm may by operated with different system parameter values or constraints. In this case, in order to solve the decision-making problem, it is necessary to consider not only the system parameter values but also the transportation cost to the manufacturing plant or the setup cost. Developing a stochastic programming model taking these considerations into account will also be one of the directions to expand the study.



---

*Source: 4033469-2022-12-23.xml*