# Inventory Model with Partial Backordering When Backordered Customers Delay Purchase after Stockout-Restoration

Authors: Ren-Qian Zhang; Yan-Liang Wu; Wei-Guo Fang; Wen-Hui Zhou
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Many inventory models with partial backordering assume that the backordered demand must be filled instantly after stockout restoration. In practice, however, the backordered customers may successively revisit the store because of the purchase delay behavior, producing a limited backorder demand rate and resulting in an extra inventory holding cost. Hence, in this paper we formulate the inventory model with partial backordering considering the purchase delay of the backordered customers and assuming that the backorder demand rate is proportional to the remaining backordered demand. Particularly, we model the problem by introducing a new inventory cost component of holding the backordered items, which has not been considered in the existing models. We propose an algorithm with a two-layer structure based on Lipschitz Optimization (LO) to minimize the total inventory cost. Numerical experiments show that the proposed algorithm outperforms two benchmarks in both optimality and efficiency. We also observe that the earlier the backordered customer revisits the store, the smaller the inventory cost and the fill rate are, but the longer the order cycle is. In addition, if the backordered customers revisit the store without too much delay, the basic EOQ with partial backordering approximates our model very well.


---

## Body

## 1. Introduction

It is observed usually in daily life that customers who have incurred a stockout might not be avid to pay and fetch the backordered item definitely at the next replenishment point, when the stocked out item is available again. Consider that a customer visits an electronics store for a computer but incurs a stockout. If the customer is willing to backorder the demand, the manager will inform her/him of the date when the computer is available again. In an ideal situation, as assumed in the existing model of economic order quantity with partial backordering (EOQ-PBO), the customer will instantly revisit the store and fetch the backordered computer exactly at the next replenishment point [1, 2]. It implies that the backordered customers generate an unlimited demand rate and the store need not keep any inventory for the backordered demand.However, in some practical situations, it need take time to meet the backordered customer and the customer who is not in urgency is more possible to delay the purchase or decrease the order frequency to the store after incurring a stockout [3–5]. Thus, among the customers who have backordered the computer because of the stockout, some of them may repurchase the item immediately as soon as the inventory is replenished again, but some of them will revisit the store relatively later after the stockout-restoration, although they have incurred a stockout. Similarly, some other backordered customers may revisit the store more lately and so on.Because of the purchase delay behavior, the backordered customers will repurchase the item generating a limited, not an unlimited, demand rate. Consequently, the vendor has to keep inventory for the backordered customers, resulting in a new cost component in the inventory model, which therefore leads to a different inventory policy for the vendor as a result of the different cost structure. This paper will reinvestigate the retailer-consumer inventory problem, considering that the backordered customers do not repurchase the item exactly at the next replenishment point but rather revisit the store successively during the next in-stock period because of the purchase delay behavior.


## 2. Literature Review

Assuming that a fraction of unsatisfied customers incurring stockouts is willing to backorder their demand and be satisfied by the next replenishment led to the problem of EOQ with partial backordering (EOQ-PBO) [1]. The basic EOQ-PBO model assumes that a constant percentage β of the demand during the stockout period will be backordered, and the remaining fraction 1-β of the unsatisfied demand is lost [6–10]. Among different variants of the model [2], one of the most frequent extensions is to introduce a changing pattern for the backordering rate rather than assuming that it is a constant percentage [11–15]. These models have considered the psychology regarding how the backordered customers will wait the stockout-restoration before the next replenishment point, which is often varying with the waiting time.Some researches considered new demand patterns rather than using the constant demand rate. The often adopted pattern is that demand is a time-varying function [16–22]. Another demand pattern adopted in EOQ-PBO is inventory-level-dependent demand, which is also used in inventory model of perishable items [23–27]. Besides, some researches assumed that the demand rate was dependent on price and time [28, 29]. We refer the reader to the survey of Pentico and Drake [2] for a more extensive survey. Recently, the problem was extended to consider prepayment, delayed payment, linear and fixed backordering costs, discounts, and so forth [30–33].The above literature shows that most variants of the EOQ-PBO model are focused on different patterns of demand rate or backordering rate.The backordering rate describes how the unsatisfied customers will wait for the next replenishment. However,it does not consider the “sequel” after the customer completes the wait, that is, how the backordered customers revisit the store and fetch their backordered items after the next replenishment arrives and the item is available again.Duran et al. [34] pointed out that unsatisfied customers may delay their purchase in the next replenishment period because they wish to obtain greater customer utility (either increased service or decreased price). More generally, Greenleaf and Lehmann [35] and G.-R. Chen and M.-L. Chen [36] reported that customers may delay their purchase just because they do not have enough time, feel unpleasant for shopping the product, experience perceived risk, rely on advice from others, wait for a lower price, and so forth. It is also possible that one may delay occasionally a purchase because of unexpected events [37]; for example, unexpected business meeting obstructs the backordered customer revisiting the store in time, or even the backordered customer delays the purchase due to avoidance and indecision [38]. Moore and Fitzsimons [39] further showed that consumers may respond negatively to the vendor after the stockout-restoration. For online purchase situations, although retailers can promptly inform consumers of product availability, consumers do not immediately recover the negative effects of the stockout [40]. We can reasonably speculate that all these factors will considerably delay the repurchasing action of the backordered customers.On the other hand, in some practical situations, the retailer may also delay delivering goods to the customers who have lower priorities [41]. Even if the backordered customers require that the vendor sends out the item as soon as it is available, they need still wait a number of days until the product is delivered [42]. In case of cash on delivery, it means that the repurchase is delayed and the retailer has to bear the inventory cost for keeping the backordered item on passage.In sum, despite the possible delay of vendor or customer, an extra inventory of backordered demand should be held with longer time compared with the rigid assumption in most of the existing EOQ-PBO models that the backordered demand must be fulfilled completely as soon as the new order arrives. We have checked the above-mentioned literatures that are related to the customer’s purchase delay behavior, but none of them had introduced the factor of purchase delay into inventory model [3–5, 34–41].This paper extends the EOQ-PBO problem considering the purchase delay of the backordered customer; that is,the backordered demand will be fulfilled gradually during the subsequent replenishment period rather than being met outright at the very beginning of the contiguous inventory period. The key of the problem is to model the endogenous demand rate derived from the backordered customers and to formulate the inventory cost function considering new cost structure.The rest of the paper is organized as follows. In Section3, we present the problem definition and further list the new contributions of our study. Section 4 formulates the model by deriving the backorder demand rate and the total inventory cost function. Section 5 first briefly introduces Lipschitz Optimization (LO) used in the proposed algorithm. Subsequently, the internal-layer algorithm for obtaining the optimal order cycle with given fill rate is presented, based on which the external-layer algorithm for finding the globally optimal solution is proposed. Section 6 evaluates the performance of the proposed algorithm and presents managerial insights by numerical computations. Conclusions are given in Section 7.


## 3. Problem Definition

The traditional EOQ-PBO model assumes one constant demand rate, but the vendor in fact must handle two kinds of demand during each inventory period: (1) a stable demand for ordinary customers who have not incurred any stockouts (we call the corresponding demand per unit timeroutine demand rate); (2) an endogenous demand derived from the previous-period backordered customers, which is characterized by thebackorder demand rate.Specifically, we assume that the endogenous demand is dependent on the amount of the backordered customers who have not been satisfied, which thereby generates a decreasing and convex backorder demand rate. Different from the existing inventory models, the inventory holding cost in this case will include two parts: the inventory holding cost of the item for meeting the routine demand and the inventory holding cost of the item for satisfying theendogenous backordered demand. Moreover, considering ordering cost and the opportunity cost of lost sales, we will determine the optimal fill rate and the optimal replenishment interval so as to minimize the total inventory cost.Our study is based on the literature regarding two streams. One is focused on the optimal inventory lot-sizing determined by minimizing the total inventory cost (Stream I in Table1). The other is focused on customer purchasing behavior, referring mainly to the reason-explanation for delaying purchase and the related statistical analysis (Stream II in Table 1). Because of the endogenous backordered demand, the components of inventory cost are different from those of the existing inventory models. We compare our model with the existing models in Table 1, for clearly highlighting the new contribution of this paper.Table 1
Cost components in related literature and our model.
Inventory cost components/factors
Papers reviewed
Our model
Stream I
Stream II
Ordering cost
√
×
√
Inventory holding cost of in-stock items
√
×
√
Backorder cost of stockouts
√
×
√
Opportunity cost of lost sales
√
×
√
Repurchase delay of backordered customers
×
√
√
Describe/analyze the behavior of purchase delay
×
√
√
Endogenous demand rate of backordered customers
×
×
√
Inventory holding cost of backordered demand
×
×
√
Stream I
Montgomery et al.  [1], Rosenberg [6], Park  [7, 8], Wee  [9], Padmanabhan and Vrat  [23], Wee  [16], Abad  [28, 29], Abboud and Sfairy  [11], Skouri and Papachristos  [17], Zhou et al.  [18], Ghosh and Chaudhuri  [19], San José et al.  [12–14], Hou  [24], Jain et al.  [25], Uthayakumar and Geetha  [26], Pentico and Drake  [10], Shah  [20], Yang et al.  [27], Dye and Hsieh  [15], Taleizadeh et al.  [30–32], and Wee et al.  [33]
Stream II
Ajzen [37], Greenleaf and Lehmann  [35], Darpy  [38], Zinn and Liu [3], Duran et al.  [34], Liberopoulos and Tsikis  [4], Li and Chen  [41], G.-R. Chen and M.-L. Chen [36], Rao et al.  [5], Pizzi and Scarpi [40], and Moore and Fitzsimons [39]To our knowledge, this is the first time that the backordered customer’s purchase delay is introduced into the inventory model to minimize the total inventory cost. Model with this assumption is more complicate on account that we must handle with a new inventory cost component, that is, the inventory holding cost of the item kept for the backordered customers. What is more, a two-layer global optimization algorithm is needed for the solution because of the mathematical intractability.


## 4. Formulation of the Model

### 4.1. Notations and Assumptions

We keep all assumptions of the basic EOQ-PBO model [2, 10] except assuming that the backordered demand is filled instantly at the beginning of the very next order cycle with an infinitely large demand rate. If the backordered customers delay their repurchasing and therefore do not fetch the backordered items immediately at the beginning of the next order cycle, the corresponding demand rate is limited. Therefore, the store will incur an inventory holding cost of the backordered demand because of keeping inventory for the backordered customers, as in right side of Figure 1.Figure 1
Two different cases of satisfying the backordered demand.The notations in Notations section are used to model the problem.Obviously, the total order quantity is equal toDFT+B-. It can be virtually divided into two parts: the first part DFT is used for satisfying theroutine demand of the first time visiting (first-visiting) customers by the shelf stock during in-stock period. The corresponding demand rate is D. The second part, B-=βD(1-F)T, is ordered for satisfying the demand that is backordered during the previous stockout period and will be consumed according to thebackorder demand rateDb(t). We can see from the left side of Figure 1 that if the backordered demand is satisfied instantly, then the inventory of the order quantity B- for the backordered customers becomes 0 immediately as soon as the inventory is replenished, without causing any inventory holding cost; otherwise, the corresponding inventory level will gradually decrease to 0 and the store must bear an extra inventory holding cost (see the right side of Figure 1).Different repurchasing behavior of backordered customers generates different backorder demand rateDb(t). Because of the complexity of the purchase delay factors, we consider an “average effect” of the purchase delay for facilitating modeling the inventory problem. That is, all the backordered customers are regarded as homogeneous so that at any time during the next in-stock period, each backordered customer has the same probability of revisiting the store and buying the backordered item. It means that during any unit time, the amount of the backordered customers who implement the repurchasing action should be proportional to the remaining backordered customers, which is equally expressed in Assumption 1, as follows.Assumption 1.
Thebackorder demand rate is proportional to the remained unsatisfied backordered customers; that is, Dbt=αB(t) and α is a positive real number.Moreover, for the convenience in management, a vendor often commits one-time backorder for a backordered customer. It is possible that a so-called “backordered” customer might miss the entire very next in-stock period due to excessive purchase delay (the delay time is longer thanFT). In this case, the customer will suffer a second stockout and it is highly possible that the corresponding demand will be lost forever due to customer’s impatience [43, 44]. But in this case, we should not fault the vendor again for the “second-time” stockout, because it is completely caused by the customer’s excessive purchase delay. Thus, this kind of “backordered” customers does not lead to any inventory cost and therefore does not impact the inventory model. In fact, if a “backordered” customer misses the next replenishment again and fades forever, the store should not replenish for the corresponding demand. This is equal to the case that the customer does not backorder the item at all when they incur the first time stockout. Thus, we can calibrate the backordering rate β to eliminate the “backordered” customer who incurs a second stockout because of the excessive purchase delay time longer than FT. Hence, to facilitate modeling the inventory problem, we have Assumption 2.Assumption 2.
All the backordered customers must be completely satisfied during all the next in-stock period of[0,FT].It is also possible that all backordered customers are satisfied during a period of time that is shorter thanFT. If it happens, our formulation will not precisely model the inventory problem. But we can approximate this case by adopting a very large attenuation coefficient, which ensures that almost all the backordered customers are satisfied before FT. From this point, we approximately assume, in Assumption 2, that the backordered customers are satisfied during all the in-stock period of [0,FT].


### 4.2. The Backorder Demand Rate

Based on the backorder demand rate assumption (Assumption1), we have a differential equation as(1)Dbt=-dBtdt=αBt,α>0.The solution to the differential equation is Bt=B0e-αt, where B0 is a constant that is dependent on boundary conditions. Thus, the backorder demand rate at time t is(2)Dbt=-dBtdt=αB0e-αt.The backorder demand rate Dbt is an exponentially decreasing function in time t, where α in (2) is theattenuation coefficient of the backorder demand rate. αB0 represents the initial demand rate at time t=0, which is dependent on how many demand has been backordered. We can see that a larger α corresponds to a larger backorder demand rate at the beginning of order cycle, but with a faster attenuation.Since all the backordered customers are satisfied during the next in-stock period (Assumption2), the total backordered demand B- meets(3)B-=∫t=0FTDbtdt=∫t=0FTαB0e-αtdt=B01-e-αFT.The total quantity of backordered demand during each order cycle is B-=βS=βD(1-F)T [10]. Therefore, B0 is determined by the following equation:(4)B0=βD1-FT1-e-αFT.


### 4.3. The Inventory Cost Function

#### 4.3.1. The Inventory Holding Cost of the Backordered Demand

Figure2 presents the inventory level curves, where the upper part of the figure indicates how the backordered demands are satisfied in the coming order cycle and the lower part presents how the routine demand is satisfied.Figure 2
Inventory levels over time.The backorder demand rateDb(t) exponentially attenuates and gradually reaches 0, so that the backordered customers successively revisit the store during the in-stock period [0,FT]. Obviously, the inventory level corresponding to the backorder demand changes as(5)Ibt=B--∫τ=0tDbτdτ=B--∫τ=0tαB0e-ατdτ=B0e-αt-e-αFT.Regarding the area under the inventory level curve of the backorder demand, we have (6)Sb=∫t=0FTIbtdt=B01α1-e-αFT-e-αFTFT.Considering that B0=βD(1-F)T/1-e-αFT (see (4)), the inventory holding cost per unit time of the backordered demand is(7)ChSbT=ChTB01α1-e-αFT-e-αFTFT=βDCh1-Fα-βDCh1-FFTeαFT-1.


#### 4.3.2. The Total Inventory Cost Function

From the existing models, it is easy to know that the ordering cost of one unit time isA/T, the backordering cost per unit time is βDCb(1-F)2T/2, the opportunity cost of lost sales is CoD1-β1-F, and the holding cost of the inventory for satisfying the first-visiting customers is DChF2T/2. Summing up all the cost terms yields the total inventory cost per unit time as in the following:(8)ΓT,F=AT+DChF2+βDCb1-F22T-βDCh1-FααFTeαFT-1+βDCh1-Fα+CoD1-β1-F.Introduce auxiliary functions as follows:(9)uF=DChF2+βDCb1-F22,vF=βDCh1-Fα,wF=vF+CoD1-β1-F,θx=xex-1-1.Then, the total cost function is recast as(10)ΓT,F=AT+uFT-vFθx+wF,where x=αFT. For the convenience, we denote f′x=dfx/dx in this paper.Compared with the existing EOQ-PBO model, the inventory total cost function of (10) has a new term of βDCh1-F/α-βDCh1-F/ααFT/eαFT-1. It leads to the intractability in finding the closed-form solution, which forces us to develop an algorithm based on Lipschitz Optimization (LO).Note that the backordered customer may prepay a deposit for the backordered item. In this case, the store will incur less financial cost for the capital of inventory, where the unit inventory holding cost per unit time of the backordered item should be smaller thanCh, say ch<Ch. To adapt the model for this case, we need only to substitute ch for Ch in the auxiliary function vF in the model. To be concise, we do not differentiate ch and Ch in this paper.


## 4.1. Notations and Assumptions

We keep all assumptions of the basic EOQ-PBO model [2, 10] except assuming that the backordered demand is filled instantly at the beginning of the very next order cycle with an infinitely large demand rate. If the backordered customers delay their repurchasing and therefore do not fetch the backordered items immediately at the beginning of the next order cycle, the corresponding demand rate is limited. Therefore, the store will incur an inventory holding cost of the backordered demand because of keeping inventory for the backordered customers, as in right side of Figure 1.Figure 1
Two different cases of satisfying the backordered demand.The notations in Notations section are used to model the problem.Obviously, the total order quantity is equal toDFT+B-. It can be virtually divided into two parts: the first part DFT is used for satisfying theroutine demand of the first time visiting (first-visiting) customers by the shelf stock during in-stock period. The corresponding demand rate is D. The second part, B-=βD(1-F)T, is ordered for satisfying the demand that is backordered during the previous stockout period and will be consumed according to thebackorder demand rateDb(t). We can see from the left side of Figure 1 that if the backordered demand is satisfied instantly, then the inventory of the order quantity B- for the backordered customers becomes 0 immediately as soon as the inventory is replenished, without causing any inventory holding cost; otherwise, the corresponding inventory level will gradually decrease to 0 and the store must bear an extra inventory holding cost (see the right side of Figure 1).Different repurchasing behavior of backordered customers generates different backorder demand rateDb(t). Because of the complexity of the purchase delay factors, we consider an “average effect” of the purchase delay for facilitating modeling the inventory problem. That is, all the backordered customers are regarded as homogeneous so that at any time during the next in-stock period, each backordered customer has the same probability of revisiting the store and buying the backordered item. It means that during any unit time, the amount of the backordered customers who implement the repurchasing action should be proportional to the remaining backordered customers, which is equally expressed in Assumption 1, as follows.Assumption 1.
Thebackorder demand rate is proportional to the remained unsatisfied backordered customers; that is, Dbt=αB(t) and α is a positive real number.Moreover, for the convenience in management, a vendor often commits one-time backorder for a backordered customer. It is possible that a so-called “backordered” customer might miss the entire very next in-stock period due to excessive purchase delay (the delay time is longer thanFT). In this case, the customer will suffer a second stockout and it is highly possible that the corresponding demand will be lost forever due to customer’s impatience [43, 44]. But in this case, we should not fault the vendor again for the “second-time” stockout, because it is completely caused by the customer’s excessive purchase delay. Thus, this kind of “backordered” customers does not lead to any inventory cost and therefore does not impact the inventory model. In fact, if a “backordered” customer misses the next replenishment again and fades forever, the store should not replenish for the corresponding demand. This is equal to the case that the customer does not backorder the item at all when they incur the first time stockout. Thus, we can calibrate the backordering rate β to eliminate the “backordered” customer who incurs a second stockout because of the excessive purchase delay time longer than FT. Hence, to facilitate modeling the inventory problem, we have Assumption 2.Assumption 2.
All the backordered customers must be completely satisfied during all the next in-stock period of[0,FT].It is also possible that all backordered customers are satisfied during a period of time that is shorter thanFT. If it happens, our formulation will not precisely model the inventory problem. But we can approximate this case by adopting a very large attenuation coefficient, which ensures that almost all the backordered customers are satisfied before FT. From this point, we approximately assume, in Assumption 2, that the backordered customers are satisfied during all the in-stock period of [0,FT].


## 4.2. The Backorder Demand Rate

Based on the backorder demand rate assumption (Assumption1), we have a differential equation as(1)Dbt=-dBtdt=αBt,α>0.The solution to the differential equation is Bt=B0e-αt, where B0 is a constant that is dependent on boundary conditions. Thus, the backorder demand rate at time t is(2)Dbt=-dBtdt=αB0e-αt.The backorder demand rate Dbt is an exponentially decreasing function in time t, where α in (2) is theattenuation coefficient of the backorder demand rate. αB0 represents the initial demand rate at time t=0, which is dependent on how many demand has been backordered. We can see that a larger α corresponds to a larger backorder demand rate at the beginning of order cycle, but with a faster attenuation.Since all the backordered customers are satisfied during the next in-stock period (Assumption2), the total backordered demand B- meets(3)B-=∫t=0FTDbtdt=∫t=0FTαB0e-αtdt=B01-e-αFT.The total quantity of backordered demand during each order cycle is B-=βS=βD(1-F)T [10]. Therefore, B0 is determined by the following equation:(4)B0=βD1-FT1-e-αFT.


## 4.3. The Inventory Cost Function

### 4.3.1. The Inventory Holding Cost of the Backordered Demand

Figure2 presents the inventory level curves, where the upper part of the figure indicates how the backordered demands are satisfied in the coming order cycle and the lower part presents how the routine demand is satisfied.Figure 2
Inventory levels over time.The backorder demand rateDb(t) exponentially attenuates and gradually reaches 0, so that the backordered customers successively revisit the store during the in-stock period [0,FT]. Obviously, the inventory level corresponding to the backorder demand changes as(5)Ibt=B--∫τ=0tDbτdτ=B--∫τ=0tαB0e-ατdτ=B0e-αt-e-αFT.Regarding the area under the inventory level curve of the backorder demand, we have (6)Sb=∫t=0FTIbtdt=B01α1-e-αFT-e-αFTFT.Considering that B0=βD(1-F)T/1-e-αFT (see (4)), the inventory holding cost per unit time of the backordered demand is(7)ChSbT=ChTB01α1-e-αFT-e-αFTFT=βDCh1-Fα-βDCh1-FFTeαFT-1.


### 4.3.2. The Total Inventory Cost Function

From the existing models, it is easy to know that the ordering cost of one unit time isA/T, the backordering cost per unit time is βDCb(1-F)2T/2, the opportunity cost of lost sales is CoD1-β1-F, and the holding cost of the inventory for satisfying the first-visiting customers is DChF2T/2. Summing up all the cost terms yields the total inventory cost per unit time as in the following:(8)ΓT,F=AT+DChF2+βDCb1-F22T-βDCh1-FααFTeαFT-1+βDCh1-Fα+CoD1-β1-F.Introduce auxiliary functions as follows:(9)uF=DChF2+βDCb1-F22,vF=βDCh1-Fα,wF=vF+CoD1-β1-F,θx=xex-1-1.Then, the total cost function is recast as(10)ΓT,F=AT+uFT-vFθx+wF,where x=αFT. For the convenience, we denote f′x=dfx/dx in this paper.Compared with the existing EOQ-PBO model, the inventory total cost function of (10) has a new term of βDCh1-F/α-βDCh1-F/ααFT/eαFT-1. It leads to the intractability in finding the closed-form solution, which forces us to develop an algorithm based on Lipschitz Optimization (LO).Note that the backordered customer may prepay a deposit for the backordered item. In this case, the store will incur less financial cost for the capital of inventory, where the unit inventory holding cost per unit time of the backordered item should be smaller thanCh, say ch<Ch. To adapt the model for this case, we need only to substitute ch for Ch in the auxiliary function vF in the model. To be concise, we do not differentiate ch and Ch in this paper.


## 4.3.1. The Inventory Holding Cost of the Backordered Demand

Figure2 presents the inventory level curves, where the upper part of the figure indicates how the backordered demands are satisfied in the coming order cycle and the lower part presents how the routine demand is satisfied.Figure 2
Inventory levels over time.The backorder demand rateDb(t) exponentially attenuates and gradually reaches 0, so that the backordered customers successively revisit the store during the in-stock period [0,FT]. Obviously, the inventory level corresponding to the backorder demand changes as(5)Ibt=B--∫τ=0tDbτdτ=B--∫τ=0tαB0e-ατdτ=B0e-αt-e-αFT.Regarding the area under the inventory level curve of the backorder demand, we have (6)Sb=∫t=0FTIbtdt=B01α1-e-αFT-e-αFTFT.Considering that B0=βD(1-F)T/1-e-αFT (see (4)), the inventory holding cost per unit time of the backordered demand is(7)ChSbT=ChTB01α1-e-αFT-e-αFTFT=βDCh1-Fα-βDCh1-FFTeαFT-1.


## 4.3.2. The Total Inventory Cost Function

From the existing models, it is easy to know that the ordering cost of one unit time isA/T, the backordering cost per unit time is βDCb(1-F)2T/2, the opportunity cost of lost sales is CoD1-β1-F, and the holding cost of the inventory for satisfying the first-visiting customers is DChF2T/2. Summing up all the cost terms yields the total inventory cost per unit time as in the following:(8)ΓT,F=AT+DChF2+βDCb1-F22T-βDCh1-FααFTeαFT-1+βDCh1-Fα+CoD1-β1-F.Introduce auxiliary functions as follows:(9)uF=DChF2+βDCb1-F22,vF=βDCh1-Fα,wF=vF+CoD1-β1-F,θx=xex-1-1.Then, the total cost function is recast as(10)ΓT,F=AT+uFT-vFθx+wF,where x=αFT. For the convenience, we denote f′x=dfx/dx in this paper.Compared with the existing EOQ-PBO model, the inventory total cost function of (10) has a new term of βDCh1-F/α-βDCh1-F/ααFT/eαFT-1. It leads to the intractability in finding the closed-form solution, which forces us to develop an algorithm based on Lipschitz Optimization (LO).Note that the backordered customer may prepay a deposit for the backordered item. In this case, the store will incur less financial cost for the capital of inventory, where the unit inventory holding cost per unit time of the backordered item should be smaller thanCh, say ch<Ch. To adapt the model for this case, we need only to substitute ch for Ch in the auxiliary function vF in the model. To be concise, we do not differentiate ch and Ch in this paper.


## 5. Solution

### 5.1. Solution Scheme and Lipschitz Optimization (LO)

#### 5.1.1. Solution Scheme

The convexity of the cost functionΓT,F cannot be obtained and to derive a closed-form solution for the problem of min⁡ΓT,F is very difficult, so we consider a two-layer strategy to find the optimal solution. Specifically, in the internal layer we minimize function ΓT,F supposing that F is given, say F=F~. This yields the solution to the problem of minT⁡Γ~(T,F~) for any given F=F~∈[0,1]. Define a new function(11)Γ^F~=minT⁡Γ~T,F~,representing the minimal value of Γ(T,F) for F=F~.In the external layer of the algorithm, we search the optimal value ofF~ that minimizes function Γ^(F~), which finally yields the globally optimal solution to the problem of min⁡ΓT,F. However, it is still very difficult to derive a closed-form solution even optimizing the single variable problem minT⁡Γ~(T,F~) for F=F~, so we use Lipschitz Optimization (LO) to solve the problem of minT⁡Γ~(T,F~) and minF~Γ^(F~), respectively.Functionf(x) is Lipschitzian on interval [a,b] if, for any two points x,y∈[a,b], the absolute difference of the function is bounded by a constant multiplying the distance between x and y. That is, |fx-fy|≤L|x-y|, where L is theLipschitz constant. We have found that it is easy to determine theLipschitz constant for function Γ~(T,F~) given F=F~, and therefore the LO algorithm with known Lipschitz constant can be applied to optimizing function Γ~(T,F~). However, the Lipschitz constant of function Γ^(F~) is very difficult to be obtained, so we have to optimize function Γ^(F~) without knowing Lipschitz constant. Below, we briefly introduce the basic idea of LO algorithms and the reader is referred to the related literature for the details [45, 46].


#### 5.1.2. Shubert’s Algorithm with a Known Lipschitz Constant

There are several LO algorithms for optimizing a Lipschitz function if the Lipschitz constantL is known. Thepassive algorithm is the simplest one that evaluates the function on points a+ε/L,a+3ε/L,a+5ε/L,…∈[a,b], where the function value on the best point does not differ more than ε from the global optimal value.A more efficient LO method isShubert’s algorithm. Shubert’s algorithm uses piecewise linear functions with slopes of L and -L to approximate the objective function from below on interval [a,b] (see Figure 3).Figure 3
The basic idea of Shubert’s algorithm.Based on the piecewise linear function approximation, the algorithm continuously selects the interval with the minimal lower bound and partitions it into subintervals until optima are obtained under a given preciseness. We refer the reader to Shubert [45] for the details. In this paper, we will adopt Shubert’s algorithm to minimize Γ~(T,F~), obtaining the optimal value of T with given F.


#### 5.1.3. The DIRECT Algorithm without Lipschitz Constant

Jones et al. [46] proposed a new variant of LO algorithm, named DIRECT, which does not need to know the Lipschitz constant. It is different from Shubert’s algorithm that the DIRECT algorithm approximates the objective function from the middle point of the interval and partitions the selected interval (the interval with minimal lower bound) into three equal parts, as shown in Figure 4.Figure 4
Interval partition of the DIRECT algorithm.Figure4 shows that we first partition the interval into three equal parts whose centers are c1,c2, and c3, respectively. The two lines crossing the middle points c1,c2, and c3 with slopes ±K provide lower bounds for the objective function on intervals a,a2,a2,b2, and b2,b, respectively. Then select the interval with minimal lower bound, that is, [a,a2], for the further trisection while updating K with K′.Generally, letci=(ai+bi)/2 denote the center of interval [ai,bi]. The two lines crossing the middle point ci with slopes ±K provide lower bound for the objective function on interval [ai,bi]. Denote c-i=bi-ai, the width of interval [ai,bi]. It is noticeable that there are multiple intervals corresponding to each value of c-i. Each pair of c-i and f(ci) corresponds to a point with coordinates (c-i,f(ci)) on a two-dimensional plane. As long as DIRECT selects the interval for the further trisection, the value of K can be updated for the new iteration by finding the lower right of the convex hull of the set of points (c-i,f(ci)), which can be realized by theconhull algorithm proposed in Björkman and Holmström [47]. The above procedure is repeated until the convergence basin of the optimum is found. Obviously, the DIRECT algorithm continuously renews K to determine the slope of the crossing lines for approximating the objective function, so the Lipschitz constant is not needed. To save space, we refer the reader to Jones et al. [46] for the details of the DIRECT algorithm.


### 5.2. The Optimal Value ofT When F Is Given

#### 5.2.1. The First-Order Condition of Optimality

IfF is given, say F=F~, we introduce two functions rT and qT as follows:(12)rT=AT+uF~T,qT=vF~θαF~T. Then, the cost function with given F is recast as(13)Γ~T=Γ~T,F~=rT-qT+wF~.The optimal value of Γ~T,F~ defined as Γ^(F~)=minT⁡Γ~(T,F~) in (11) is(14)Γ^F~=minT⁡Γ~T,F~=minT⁡Γ~T.Proposition3 gives the first-order necessary condition for the minima of function Γ~T.Proposition 3.
The minima of functionΓ~T solve the following equation:(15)Γ~′T=r′T-q′T=0,where(16)r′T=-AT2+uF~,q′T=αF~·vF~·θ′αF~T.​Proof.
WhenT→0, r(T)→∞ and qT>0, so Γ~(T)→∞. On the other hand, if T→∞, we also have Γ~(T)→∞. Moreover, Γ~(T) is continuous and differentiable, and for a given T,  Γ~(T) has a finite value. Hence, the minima of function Γ~(T) meet the first-order condition. This completes the proof.


#### 5.2.2. Bounds of the Optimal Value ofT with Given F

DenoteT^F the optimal value of T with given F. But even if F is known, the convexity of function Γ~T cannot be obtained, and it is also difficult to solve the first-order necessary condition to obtain T^F for minimizing Γ~T. To use LO to find the optimal solution for a given F, we first determine the bounds of T^F.Proposition 4.
LetT-F denote the upper bound of T^F and define a sequence {T-Fk} as(17)T-Fk=AuF~-αF~·vF~·θ′x-k-1,where x-(k-1)=αF~T-Fk-1, (k=1,…,∞), and  T-F(0)=A/u(F~). Then, T^F≤T-F(k) such that, for any k, T-Fk provides an upper bound for T^F and T-F=limk→∞⁡T-F(k) is the minimal upper bound of T^F.Proof.
See AppendixB.Proposition 5.
LetT_F denote the lower bound of T^F and define a sequence {T_Fk} as(18)T_Fk=AuF~-αF~·vF~·θ′x_k-1,where x_(k-1)=αF~T_Fk-1, (k=1,…,∞), and T_F0=A/uF~+αF~·vF~/2. Then, T^F≥T_Fk such that, for any k, T_Fk provides a lower bound for T^F and T_F=limk→∞⁡T_F(k) is the maximal lower bound of T^F.Proof.
See AppendixB.Propositions4 and 5 give the bounds of T^F. For the limits of upper and lower bounds when k→∞, we have Proposition 6.Proposition 6.
The limits of upper and lower bounds both solve the first-order necessary condition given by (15).Proof​.
See AppendixB.Corollary 7.
Ifα→∞, that is, the backorder demand rate is infinitely large, the upper and lower bounds are identical and both solve the first-order condition of the basic EOQ-PBO model.​Proof.
It is easy to verify thatlimx→∞⁡θ′x=0 (see Lemma A.2 in Appendix A). Thus, when α→∞, we have x=αFT→∞ and so limk→∞T-F(k)=limk→∞T_Fk=2A/u(F)=2A/D[ChF2+βCb(1-F)2]. This is equation (9) of Pentico and Drake [10], which solves the first-order condition of the basic EOQ-PBO model. The proof is completed.Corollary7 shows that the basic EOQ-PBO is a special case of our model. If α is sufficiently large, we can use the basic EOQ-PBO model as an approximation. The numerical experiment in Section 6 suggests that if α>30, the basic EOQ-PBO can approximate our model very well.Corollary 8.
Forα<∞, if the limits of upper and lower bounds are identical, that is, T-F=T_F, then they are the unique optimal solution to the problem of minTΓ~T; that is, T^F=T-F=T_F.Proof.
It is straightforward from Propositions3 and 6.Corollary8 allows us to improve the efficiency of finding the optimal solution T^F. For a given F, we first check if the limits of the upper and lower bounds are identical. Only when the two limits are not equal, should we further use some optimization algorithms, such as the LO algorithm, to minimize function Γ~T; otherwise, the optimal value T^F=T-F=T_F is attained. Proposition 9 presents a sufficient condition under which T-F and T_F are equal to each other.Proposition 9.
The optimal value of order cycleT with given F can be determined by T^F=T-F=T_F if the following inequality holds:(19)α2F~2vF~A<12uF~-αF~vF~θ′αF~T-F3/2.Proof.
See AppendixB.


#### 5.2.3. Lipschitz Optimization onT with Given F

Based on the bounds ofT determined by Propositions 4 and 5, we can use LO to minimize Γ~(T) when F is given, for which Proposition 10 presents the Lipschitz constant.Proposition 10.
ForT1,T2∈[T_F,T-F], one has |Γ~T1-Γ~(T2)|≤LF~·|T1-T2|, where LF~ is the Lipschitz constant that is determined by(20)LF~=max⁡uF~+αF~·vF~2-AT-F2,uF~-AT_F2.Proof.
See AppendixB.Additionally, Proposition9 shows that if the limits of the upper and lower bounds are identical, then the limit directly gives the optimal value of T as T^F=T-F=T_F. In this case, we do not need to use the LO method. The procedure of the algorithm optimizing T with given F is detailed as follows. We name this algorithm Opt_Local for the convenience.Algorithm Opt_LocalStep 1.
For givenF=F~, calculate the limit of the upper bound T-F from (17), by iterations.Step 2.
Check whether (19) holds. If yes, the optimal value of T with given F is obtained by T^F=T-F and go to Step 6; otherwise, go to Step 3.Step 3.
Calculate the limit of the lower boundT_F from (18) by iterations.Step 4.
IfT_F=T-F, then the optimal value of T with given F is obtained by T^F=T_F=T-F and go to Step 6; otherwise, go to Step 5.Step 5.
Determine theLipschitz constant from (20). Apply Lipschitz Optimization (Shubert’s algorithm) to the problem of minT⁡Γ~T and the optimal value T^F is obtained.Step 6.
Output the optimal solutionT^F, and Γ^F~=minT⁡Γ~(T)=Γ~(T^F).In Step5, we set the stop criterion based on relative improvements. Specifically, if the relative improvement on the objective function Γ~T between two successive iterations keeps smaller than a very small positive number ε for p consecutive iterations, then stop. In this paper, we set ε=10-6 and p=10. Thus, if the relative improvement keeps smaller than 0.0001/% for 10 consecutive iterations, stop.


### 5.3. Algorithm for Finding the Globally Optimal Solution

To obtain the global minima, we need minimize functionΓ^(F~)=minT⁡Γ~(T,F~) to find the optimal fill rate F. The convexity of Γ^(F~) cannot be obtained, so the LO algorithm is employed again. Unfortunately, it is very difficult to derive closed-form expression for the Lipschitz constant of Γ^(F~), so we adopt the DIRECT algorithm that does not need the Lipschitz constant. The condition of using DIRECT is that the objective function should at least be continuous within the neighborhood of global optima [46].Consider two values ofF, say F=F~ and F=F~′. As mentioned above, with given F, the inventory cost function ΓT,F is optimized by algorithm Opt_Local, yielding Γ^F~=ΓT^F~,F~ and Γ^F~′=ΓT^F~′,F~′, respectively. Since function ΓT,F is differentiable and continuous in T,F, if the difference between F~′ and F~ is small enough, Γ^F~′ can be sufficiently close to Γ^F~; that is, limF~′→F~Γ^F~′=Γ^F~. Therefore, Γ^(F~) is continuous in F~.The one-dimensional DIRECT algorithm optimizing functionΓ^(F~), named AlgorithmOpt_Global, is depicted as follows.Algorithm Opt_GlobalStep 1.
Divide the selected intervals within[0,1] by the DIRECT algorithm into three equal sections and yield candidate points of F~∈[0,1].Step 2.
Determine the value ofΓ^(F~) on each candidate point of F~ by using AlgorithmOpt_Local and obtain the current best solution.Step 3.
Check whether the stop criterion is met. If yes, go to Step5; otherwise, go to Step 4.Step 4.
UpdateK, the slope of the crossing linear functions, by calling the algorithm ofconhull [47], and approximate the lower bound of function Γ^(F~) for each interval by using the crossing linear function with slope K. Intervals with the smallest lower bound of function Γ^(F~) are selected for the further trisection. Go to Step 1.Step 5.
Output the globally optimal solution.Note that, in the first iteration, after partitioning the original interval[0,1] into three equal subintervals, we select the subinterval for the next trisection directly based on the function value of the middle point of each subinterval, without approximating the objective function by the crossing lines with slope K. K will be obtained and updated starting from the second iteration. Moreover, if the relative improvement on the objective function between two successive iterations keeps smaller than ε=10-6 for p=10 consecutive iterations, then stop. In addition, as Zhang [48] showed that we must further compare Γ∗, the optimal value of ΓT,F, with CoD, the opportunity cost of not stocking at all, to determine the optimal inventory policy. The smaller one is the optimal decision.


## 5.1. Solution Scheme and Lipschitz Optimization (LO)

### 5.1.1. Solution Scheme

The convexity of the cost functionΓT,F cannot be obtained and to derive a closed-form solution for the problem of min⁡ΓT,F is very difficult, so we consider a two-layer strategy to find the optimal solution. Specifically, in the internal layer we minimize function ΓT,F supposing that F is given, say F=F~. This yields the solution to the problem of minT⁡Γ~(T,F~) for any given F=F~∈[0,1]. Define a new function(11)Γ^F~=minT⁡Γ~T,F~,representing the minimal value of Γ(T,F) for F=F~.In the external layer of the algorithm, we search the optimal value ofF~ that minimizes function Γ^(F~), which finally yields the globally optimal solution to the problem of min⁡ΓT,F. However, it is still very difficult to derive a closed-form solution even optimizing the single variable problem minT⁡Γ~(T,F~) for F=F~, so we use Lipschitz Optimization (LO) to solve the problem of minT⁡Γ~(T,F~) and minF~Γ^(F~), respectively.Functionf(x) is Lipschitzian on interval [a,b] if, for any two points x,y∈[a,b], the absolute difference of the function is bounded by a constant multiplying the distance between x and y. That is, |fx-fy|≤L|x-y|, where L is theLipschitz constant. We have found that it is easy to determine theLipschitz constant for function Γ~(T,F~) given F=F~, and therefore the LO algorithm with known Lipschitz constant can be applied to optimizing function Γ~(T,F~). However, the Lipschitz constant of function Γ^(F~) is very difficult to be obtained, so we have to optimize function Γ^(F~) without knowing Lipschitz constant. Below, we briefly introduce the basic idea of LO algorithms and the reader is referred to the related literature for the details [45, 46].


### 5.1.2. Shubert’s Algorithm with a Known Lipschitz Constant

There are several LO algorithms for optimizing a Lipschitz function if the Lipschitz constantL is known. Thepassive algorithm is the simplest one that evaluates the function on points a+ε/L,a+3ε/L,a+5ε/L,…∈[a,b], where the function value on the best point does not differ more than ε from the global optimal value.A more efficient LO method isShubert’s algorithm. Shubert’s algorithm uses piecewise linear functions with slopes of L and -L to approximate the objective function from below on interval [a,b] (see Figure 3).Figure 3
The basic idea of Shubert’s algorithm.Based on the piecewise linear function approximation, the algorithm continuously selects the interval with the minimal lower bound and partitions it into subintervals until optima are obtained under a given preciseness. We refer the reader to Shubert [45] for the details. In this paper, we will adopt Shubert’s algorithm to minimize Γ~(T,F~), obtaining the optimal value of T with given F.


### 5.1.3. The DIRECT Algorithm without Lipschitz Constant

Jones et al. [46] proposed a new variant of LO algorithm, named DIRECT, which does not need to know the Lipschitz constant. It is different from Shubert’s algorithm that the DIRECT algorithm approximates the objective function from the middle point of the interval and partitions the selected interval (the interval with minimal lower bound) into three equal parts, as shown in Figure 4.Figure 4
Interval partition of the DIRECT algorithm.Figure4 shows that we first partition the interval into three equal parts whose centers are c1,c2, and c3, respectively. The two lines crossing the middle points c1,c2, and c3 with slopes ±K provide lower bounds for the objective function on intervals a,a2,a2,b2, and b2,b, respectively. Then select the interval with minimal lower bound, that is, [a,a2], for the further trisection while updating K with K′.Generally, letci=(ai+bi)/2 denote the center of interval [ai,bi]. The two lines crossing the middle point ci with slopes ±K provide lower bound for the objective function on interval [ai,bi]. Denote c-i=bi-ai, the width of interval [ai,bi]. It is noticeable that there are multiple intervals corresponding to each value of c-i. Each pair of c-i and f(ci) corresponds to a point with coordinates (c-i,f(ci)) on a two-dimensional plane. As long as DIRECT selects the interval for the further trisection, the value of K can be updated for the new iteration by finding the lower right of the convex hull of the set of points (c-i,f(ci)), which can be realized by theconhull algorithm proposed in Björkman and Holmström [47]. The above procedure is repeated until the convergence basin of the optimum is found. Obviously, the DIRECT algorithm continuously renews K to determine the slope of the crossing lines for approximating the objective function, so the Lipschitz constant is not needed. To save space, we refer the reader to Jones et al. [46] for the details of the DIRECT algorithm.


## 5.1.1. Solution Scheme

The convexity of the cost functionΓT,F cannot be obtained and to derive a closed-form solution for the problem of min⁡ΓT,F is very difficult, so we consider a two-layer strategy to find the optimal solution. Specifically, in the internal layer we minimize function ΓT,F supposing that F is given, say F=F~. This yields the solution to the problem of minT⁡Γ~(T,F~) for any given F=F~∈[0,1]. Define a new function(11)Γ^F~=minT⁡Γ~T,F~,representing the minimal value of Γ(T,F) for F=F~.In the external layer of the algorithm, we search the optimal value ofF~ that minimizes function Γ^(F~), which finally yields the globally optimal solution to the problem of min⁡ΓT,F. However, it is still very difficult to derive a closed-form solution even optimizing the single variable problem minT⁡Γ~(T,F~) for F=F~, so we use Lipschitz Optimization (LO) to solve the problem of minT⁡Γ~(T,F~) and minF~Γ^(F~), respectively.Functionf(x) is Lipschitzian on interval [a,b] if, for any two points x,y∈[a,b], the absolute difference of the function is bounded by a constant multiplying the distance between x and y. That is, |fx-fy|≤L|x-y|, where L is theLipschitz constant. We have found that it is easy to determine theLipschitz constant for function Γ~(T,F~) given F=F~, and therefore the LO algorithm with known Lipschitz constant can be applied to optimizing function Γ~(T,F~). However, the Lipschitz constant of function Γ^(F~) is very difficult to be obtained, so we have to optimize function Γ^(F~) without knowing Lipschitz constant. Below, we briefly introduce the basic idea of LO algorithms and the reader is referred to the related literature for the details [45, 46].


## 5.1.2. Shubert’s Algorithm with a Known Lipschitz Constant

There are several LO algorithms for optimizing a Lipschitz function if the Lipschitz constantL is known. Thepassive algorithm is the simplest one that evaluates the function on points a+ε/L,a+3ε/L,a+5ε/L,…∈[a,b], where the function value on the best point does not differ more than ε from the global optimal value.A more efficient LO method isShubert’s algorithm. Shubert’s algorithm uses piecewise linear functions with slopes of L and -L to approximate the objective function from below on interval [a,b] (see Figure 3).Figure 3
The basic idea of Shubert’s algorithm.Based on the piecewise linear function approximation, the algorithm continuously selects the interval with the minimal lower bound and partitions it into subintervals until optima are obtained under a given preciseness. We refer the reader to Shubert [45] for the details. In this paper, we will adopt Shubert’s algorithm to minimize Γ~(T,F~), obtaining the optimal value of T with given F.


## 5.1.3. The DIRECT Algorithm without Lipschitz Constant

Jones et al. [46] proposed a new variant of LO algorithm, named DIRECT, which does not need to know the Lipschitz constant. It is different from Shubert’s algorithm that the DIRECT algorithm approximates the objective function from the middle point of the interval and partitions the selected interval (the interval with minimal lower bound) into three equal parts, as shown in Figure 4.Figure 4
Interval partition of the DIRECT algorithm.Figure4 shows that we first partition the interval into three equal parts whose centers are c1,c2, and c3, respectively. The two lines crossing the middle points c1,c2, and c3 with slopes ±K provide lower bounds for the objective function on intervals a,a2,a2,b2, and b2,b, respectively. Then select the interval with minimal lower bound, that is, [a,a2], for the further trisection while updating K with K′.Generally, letci=(ai+bi)/2 denote the center of interval [ai,bi]. The two lines crossing the middle point ci with slopes ±K provide lower bound for the objective function on interval [ai,bi]. Denote c-i=bi-ai, the width of interval [ai,bi]. It is noticeable that there are multiple intervals corresponding to each value of c-i. Each pair of c-i and f(ci) corresponds to a point with coordinates (c-i,f(ci)) on a two-dimensional plane. As long as DIRECT selects the interval for the further trisection, the value of K can be updated for the new iteration by finding the lower right of the convex hull of the set of points (c-i,f(ci)), which can be realized by theconhull algorithm proposed in Björkman and Holmström [47]. The above procedure is repeated until the convergence basin of the optimum is found. Obviously, the DIRECT algorithm continuously renews K to determine the slope of the crossing lines for approximating the objective function, so the Lipschitz constant is not needed. To save space, we refer the reader to Jones et al. [46] for the details of the DIRECT algorithm.


## 5.2. The Optimal Value ofT When F Is Given

### 5.2.1. The First-Order Condition of Optimality

IfF is given, say F=F~, we introduce two functions rT and qT as follows:(12)rT=AT+uF~T,qT=vF~θαF~T. Then, the cost function with given F is recast as(13)Γ~T=Γ~T,F~=rT-qT+wF~.The optimal value of Γ~T,F~ defined as Γ^(F~)=minT⁡Γ~(T,F~) in (11) is(14)Γ^F~=minT⁡Γ~T,F~=minT⁡Γ~T.Proposition3 gives the first-order necessary condition for the minima of function Γ~T.Proposition 3.
The minima of functionΓ~T solve the following equation:(15)Γ~′T=r′T-q′T=0,where(16)r′T=-AT2+uF~,q′T=αF~·vF~·θ′αF~T.​Proof.
WhenT→0, r(T)→∞ and qT>0, so Γ~(T)→∞. On the other hand, if T→∞, we also have Γ~(T)→∞. Moreover, Γ~(T) is continuous and differentiable, and for a given T,  Γ~(T) has a finite value. Hence, the minima of function Γ~(T) meet the first-order condition. This completes the proof.


### 5.2.2. Bounds of the Optimal Value ofT with Given F

DenoteT^F the optimal value of T with given F. But even if F is known, the convexity of function Γ~T cannot be obtained, and it is also difficult to solve the first-order necessary condition to obtain T^F for minimizing Γ~T. To use LO to find the optimal solution for a given F, we first determine the bounds of T^F.Proposition 4.
LetT-F denote the upper bound of T^F and define a sequence {T-Fk} as(17)T-Fk=AuF~-αF~·vF~·θ′x-k-1,where x-(k-1)=αF~T-Fk-1, (k=1,…,∞), and  T-F(0)=A/u(F~). Then, T^F≤T-F(k) such that, for any k, T-Fk provides an upper bound for T^F and T-F=limk→∞⁡T-F(k) is the minimal upper bound of T^F.Proof.
See AppendixB.Proposition 5.
LetT_F denote the lower bound of T^F and define a sequence {T_Fk} as(18)T_Fk=AuF~-αF~·vF~·θ′x_k-1,where x_(k-1)=αF~T_Fk-1, (k=1,…,∞), and T_F0=A/uF~+αF~·vF~/2. Then, T^F≥T_Fk such that, for any k, T_Fk provides a lower bound for T^F and T_F=limk→∞⁡T_F(k) is the maximal lower bound of T^F.Proof.
See AppendixB.Propositions4 and 5 give the bounds of T^F. For the limits of upper and lower bounds when k→∞, we have Proposition 6.Proposition 6.
The limits of upper and lower bounds both solve the first-order necessary condition given by (15).Proof​.
See AppendixB.Corollary 7.
Ifα→∞, that is, the backorder demand rate is infinitely large, the upper and lower bounds are identical and both solve the first-order condition of the basic EOQ-PBO model.​Proof.
It is easy to verify thatlimx→∞⁡θ′x=0 (see Lemma A.2 in Appendix A). Thus, when α→∞, we have x=αFT→∞ and so limk→∞T-F(k)=limk→∞T_Fk=2A/u(F)=2A/D[ChF2+βCb(1-F)2]. This is equation (9) of Pentico and Drake [10], which solves the first-order condition of the basic EOQ-PBO model. The proof is completed.Corollary7 shows that the basic EOQ-PBO is a special case of our model. If α is sufficiently large, we can use the basic EOQ-PBO model as an approximation. The numerical experiment in Section 6 suggests that if α>30, the basic EOQ-PBO can approximate our model very well.Corollary 8.
Forα<∞, if the limits of upper and lower bounds are identical, that is, T-F=T_F, then they are the unique optimal solution to the problem of minTΓ~T; that is, T^F=T-F=T_F.Proof.
It is straightforward from Propositions3 and 6.Corollary8 allows us to improve the efficiency of finding the optimal solution T^F. For a given F, we first check if the limits of the upper and lower bounds are identical. Only when the two limits are not equal, should we further use some optimization algorithms, such as the LO algorithm, to minimize function Γ~T; otherwise, the optimal value T^F=T-F=T_F is attained. Proposition 9 presents a sufficient condition under which T-F and T_F are equal to each other.Proposition 9.
The optimal value of order cycleT with given F can be determined by T^F=T-F=T_F if the following inequality holds:(19)α2F~2vF~A<12uF~-αF~vF~θ′αF~T-F3/2.Proof.
See AppendixB.


### 5.2.3. Lipschitz Optimization onT with Given F

Based on the bounds ofT determined by Propositions 4 and 5, we can use LO to minimize Γ~(T) when F is given, for which Proposition 10 presents the Lipschitz constant.Proposition 10.
ForT1,T2∈[T_F,T-F], one has |Γ~T1-Γ~(T2)|≤LF~·|T1-T2|, where LF~ is the Lipschitz constant that is determined by(20)LF~=max⁡uF~+αF~·vF~2-AT-F2,uF~-AT_F2.Proof.
See AppendixB.Additionally, Proposition9 shows that if the limits of the upper and lower bounds are identical, then the limit directly gives the optimal value of T as T^F=T-F=T_F. In this case, we do not need to use the LO method. The procedure of the algorithm optimizing T with given F is detailed as follows. We name this algorithm Opt_Local for the convenience.Algorithm Opt_LocalStep 1.
For givenF=F~, calculate the limit of the upper bound T-F from (17), by iterations.Step 2.
Check whether (19) holds. If yes, the optimal value of T with given F is obtained by T^F=T-F and go to Step 6; otherwise, go to Step 3.Step 3.
Calculate the limit of the lower boundT_F from (18) by iterations.Step 4.
IfT_F=T-F, then the optimal value of T with given F is obtained by T^F=T_F=T-F and go to Step 6; otherwise, go to Step 5.Step 5.
Determine theLipschitz constant from (20). Apply Lipschitz Optimization (Shubert’s algorithm) to the problem of minT⁡Γ~T and the optimal value T^F is obtained.Step 6.
Output the optimal solutionT^F, and Γ^F~=minT⁡Γ~(T)=Γ~(T^F).In Step5, we set the stop criterion based on relative improvements. Specifically, if the relative improvement on the objective function Γ~T between two successive iterations keeps smaller than a very small positive number ε for p consecutive iterations, then stop. In this paper, we set ε=10-6 and p=10. Thus, if the relative improvement keeps smaller than 0.0001/% for 10 consecutive iterations, stop.


## 5.2.1. The First-Order Condition of Optimality

IfF is given, say F=F~, we introduce two functions rT and qT as follows:(12)rT=AT+uF~T,qT=vF~θαF~T. Then, the cost function with given F is recast as(13)Γ~T=Γ~T,F~=rT-qT+wF~.The optimal value of Γ~T,F~ defined as Γ^(F~)=minT⁡Γ~(T,F~) in (11) is(14)Γ^F~=minT⁡Γ~T,F~=minT⁡Γ~T.Proposition3 gives the first-order necessary condition for the minima of function Γ~T.Proposition 3.
The minima of functionΓ~T solve the following equation:(15)Γ~′T=r′T-q′T=0,where(16)r′T=-AT2+uF~,q′T=αF~·vF~·θ′αF~T.​Proof.
WhenT→0, r(T)→∞ and qT>0, so Γ~(T)→∞. On the other hand, if T→∞, we also have Γ~(T)→∞. Moreover, Γ~(T) is continuous and differentiable, and for a given T,  Γ~(T) has a finite value. Hence, the minima of function Γ~(T) meet the first-order condition. This completes the proof.


## 5.2.2. Bounds of the Optimal Value ofT with Given F

DenoteT^F the optimal value of T with given F. But even if F is known, the convexity of function Γ~T cannot be obtained, and it is also difficult to solve the first-order necessary condition to obtain T^F for minimizing Γ~T. To use LO to find the optimal solution for a given F, we first determine the bounds of T^F.Proposition 4.
LetT-F denote the upper bound of T^F and define a sequence {T-Fk} as(17)T-Fk=AuF~-αF~·vF~·θ′x-k-1,where x-(k-1)=αF~T-Fk-1, (k=1,…,∞), and  T-F(0)=A/u(F~). Then, T^F≤T-F(k) such that, for any k, T-Fk provides an upper bound for T^F and T-F=limk→∞⁡T-F(k) is the minimal upper bound of T^F.Proof.
See AppendixB.Proposition 5.
LetT_F denote the lower bound of T^F and define a sequence {T_Fk} as(18)T_Fk=AuF~-αF~·vF~·θ′x_k-1,where x_(k-1)=αF~T_Fk-1, (k=1,…,∞), and T_F0=A/uF~+αF~·vF~/2. Then, T^F≥T_Fk such that, for any k, T_Fk provides a lower bound for T^F and T_F=limk→∞⁡T_F(k) is the maximal lower bound of T^F.Proof.
See AppendixB.Propositions4 and 5 give the bounds of T^F. For the limits of upper and lower bounds when k→∞, we have Proposition 6.Proposition 6.
The limits of upper and lower bounds both solve the first-order necessary condition given by (15).Proof​.
See AppendixB.Corollary 7.
Ifα→∞, that is, the backorder demand rate is infinitely large, the upper and lower bounds are identical and both solve the first-order condition of the basic EOQ-PBO model.​Proof.
It is easy to verify thatlimx→∞⁡θ′x=0 (see Lemma A.2 in Appendix A). Thus, when α→∞, we have x=αFT→∞ and so limk→∞T-F(k)=limk→∞T_Fk=2A/u(F)=2A/D[ChF2+βCb(1-F)2]. This is equation (9) of Pentico and Drake [10], which solves the first-order condition of the basic EOQ-PBO model. The proof is completed.Corollary7 shows that the basic EOQ-PBO is a special case of our model. If α is sufficiently large, we can use the basic EOQ-PBO model as an approximation. The numerical experiment in Section 6 suggests that if α>30, the basic EOQ-PBO can approximate our model very well.Corollary 8.
Forα<∞, if the limits of upper and lower bounds are identical, that is, T-F=T_F, then they are the unique optimal solution to the problem of minTΓ~T; that is, T^F=T-F=T_F.Proof.
It is straightforward from Propositions3 and 6.Corollary8 allows us to improve the efficiency of finding the optimal solution T^F. For a given F, we first check if the limits of the upper and lower bounds are identical. Only when the two limits are not equal, should we further use some optimization algorithms, such as the LO algorithm, to minimize function Γ~T; otherwise, the optimal value T^F=T-F=T_F is attained. Proposition 9 presents a sufficient condition under which T-F and T_F are equal to each other.Proposition 9.
The optimal value of order cycleT with given F can be determined by T^F=T-F=T_F if the following inequality holds:(19)α2F~2vF~A<12uF~-αF~vF~θ′αF~T-F3/2.Proof.
See AppendixB.


## 5.2.3. Lipschitz Optimization onT with Given F

Based on the bounds ofT determined by Propositions 4 and 5, we can use LO to minimize Γ~(T) when F is given, for which Proposition 10 presents the Lipschitz constant.Proposition 10.
ForT1,T2∈[T_F,T-F], one has |Γ~T1-Γ~(T2)|≤LF~·|T1-T2|, where LF~ is the Lipschitz constant that is determined by(20)LF~=max⁡uF~+αF~·vF~2-AT-F2,uF~-AT_F2.Proof.
See AppendixB.Additionally, Proposition9 shows that if the limits of the upper and lower bounds are identical, then the limit directly gives the optimal value of T as T^F=T-F=T_F. In this case, we do not need to use the LO method. The procedure of the algorithm optimizing T with given F is detailed as follows. We name this algorithm Opt_Local for the convenience.Algorithm Opt_LocalStep 1.
For givenF=F~, calculate the limit of the upper bound T-F from (17), by iterations.Step 2.
Check whether (19) holds. If yes, the optimal value of T with given F is obtained by T^F=T-F and go to Step 6; otherwise, go to Step 3.Step 3.
Calculate the limit of the lower boundT_F from (18) by iterations.Step 4.
IfT_F=T-F, then the optimal value of T with given F is obtained by T^F=T_F=T-F and go to Step 6; otherwise, go to Step 5.Step 5.
Determine theLipschitz constant from (20). Apply Lipschitz Optimization (Shubert’s algorithm) to the problem of minT⁡Γ~T and the optimal value T^F is obtained.Step 6.
Output the optimal solutionT^F, and Γ^F~=minT⁡Γ~(T)=Γ~(T^F).In Step5, we set the stop criterion based on relative improvements. Specifically, if the relative improvement on the objective function Γ~T between two successive iterations keeps smaller than a very small positive number ε for p consecutive iterations, then stop. In this paper, we set ε=10-6 and p=10. Thus, if the relative improvement keeps smaller than 0.0001/% for 10 consecutive iterations, stop.


## 5.3. Algorithm for Finding the Globally Optimal Solution

To obtain the global minima, we need minimize functionΓ^(F~)=minT⁡Γ~(T,F~) to find the optimal fill rate F. The convexity of Γ^(F~) cannot be obtained, so the LO algorithm is employed again. Unfortunately, it is very difficult to derive closed-form expression for the Lipschitz constant of Γ^(F~), so we adopt the DIRECT algorithm that does not need the Lipschitz constant. The condition of using DIRECT is that the objective function should at least be continuous within the neighborhood of global optima [46].Consider two values ofF, say F=F~ and F=F~′. As mentioned above, with given F, the inventory cost function ΓT,F is optimized by algorithm Opt_Local, yielding Γ^F~=ΓT^F~,F~ and Γ^F~′=ΓT^F~′,F~′, respectively. Since function ΓT,F is differentiable and continuous in T,F, if the difference between F~′ and F~ is small enough, Γ^F~′ can be sufficiently close to Γ^F~; that is, limF~′→F~Γ^F~′=Γ^F~. Therefore, Γ^(F~) is continuous in F~.The one-dimensional DIRECT algorithm optimizing functionΓ^(F~), named AlgorithmOpt_Global, is depicted as follows.Algorithm Opt_GlobalStep 1.
Divide the selected intervals within[0,1] by the DIRECT algorithm into three equal sections and yield candidate points of F~∈[0,1].Step 2.
Determine the value ofΓ^(F~) on each candidate point of F~ by using AlgorithmOpt_Local and obtain the current best solution.Step 3.
Check whether the stop criterion is met. If yes, go to Step5; otherwise, go to Step 4.Step 4.
UpdateK, the slope of the crossing linear functions, by calling the algorithm ofconhull [47], and approximate the lower bound of function Γ^(F~) for each interval by using the crossing linear function with slope K. Intervals with the smallest lower bound of function Γ^(F~) are selected for the further trisection. Go to Step 1.Step 5.
Output the globally optimal solution.Note that, in the first iteration, after partitioning the original interval[0,1] into three equal subintervals, we select the subinterval for the next trisection directly based on the function value of the middle point of each subinterval, without approximating the objective function by the crossing lines with slope K. K will be obtained and updated starting from the second iteration. Moreover, if the relative improvement on the objective function between two successive iterations keeps smaller than ε=10-6 for p=10 consecutive iterations, then stop. In addition, as Zhang [48] showed that we must further compare Γ∗, the optimal value of ΓT,F, with CoD, the opportunity cost of not stocking at all, to determine the optimal inventory policy. The smaller one is the optimal decision.


## 6. Computational Study

### 6.1. Benchmark Algorithms

To evaluate the performance ofAlgorithmOpt_Global, we use two algorithms as benchmarks. The first one is a grid search (GS) algorithm for evaluating the optimality of AlgorithmOpt_Global, and the second one is a two-dimensional DIRECT algorithm [46] for evaluating the computational efficiency.


#### 6.1.1. Grid Search onF~

We use grid search (GS) to approximately optimize function Γ^(F~) and compare the solution with that obtained by AlgorithmOpt_Global. For implementing the GS algorithm, we divide the interval F~∈[0,1] by steps of 10-4. For each problem instance, we need to optimize function Γ~T=Γ~T,F~ for ten thousand times by calling AlgorithmOpt_Local and keep the value of F~ with minimal inventory cost as the optimal solution.


#### 6.1.2. Two-Dimensional DIRECT Algorithm

The problem can also be solved by using the two-dimensional DIRECT algorithm (TD) proposed by Jones et al. [46]. However, the two-dimensional split of the space can be very time-consuming if requiring a high precision. We here compare the performance of the Opt_Global algorithm with that of the TD algorithm to show how much Opt_Global will reduce the computational effort. The stop condition of the TD algorithm is the same as that of AlgorithmOpt_Global; that is, the relative improvement keeps smaller than 10-6 for 10 consecutive iterations.


### 6.2. Outline of Numerical Examples

We investigate the performance of our model and algorithm under different scenarios by solving numerical examples. The parameter of the numerical examples should cover different cases as many as possible. Thus, we enumerate different values of parameters to construct different examples, referring to the parameter scheme of the major item of the numerical examples in Zhang et al. [49], as follows.(i) Ordering Cost. We consider four values, A = $100/order, $1000/order, $2500/order, and $5000/order, respectively, representing the difference in lower and higher ordering cost.(ii) Inventory Holding Cost and Backordering Cost. Set Ch = $5/unit/year, $10/unit/year, $25/unit/year, and $50/unit/year representing the inventory holding cost changes from low to high. Similarly, we set backorder cost as Cb = $5/unit/year, $10/unit/year, $25/unit/year, or $50/unit/year, respectively.(iii) Opportunity Cost. Set Co =$5/unit, $10/unit, $25/unit, and $50/unit, respectively.(iv) Backordering Rate. Backordering rate β is changed from 0.1 to 0.9 with a step of 0.2.(v) Demand. Set the routine demand rate D = 100 units/year, 1000 units/year, 5000 units/year, and 10000 units/year, respectively.(vi) Attenuation Coefficient of the Backorder Demand Rate. Set α=0.1,0.5,1,5,10,50,100,500, representing the backorder demand rate attenuates from slow to fast.By combination ofA×Ch×Cb×Co×β×D×α=4×4×4×4×5×4×8 we generate 40960 different numerical examples for computational experiments.


### 6.3. Performance of the Algorithm

We apply, respectively, the three algorithms:AlgorithmOpt_Global, grid search (GS), and the two-dimensional DIRECT (TD) algorithm to solve the generated 40960 numerical examples. Let ΓOpt, ΓGS, and ΓTD denote the best value of the objective function obtained by Opt_Global, GS, and TD, respectively. Let ΓGS-ΓOpt/ΓOpt indicate the relative deviation of the solution found by GS from that obtained by the Opt_Global algorithm, and ΓTD-ΓOpt/ΓOpt indicate the relative deviation of TD. We group all the numerical examples by the attenuation coefficient α, and the computational results are presented in Table 2.Table 2
Relative deviations of the two benchmarks.
α
Dev. ofΓGS (%)
Dev. ofΓTD (%)
Min.
Avg.
Max.
Min.
Avg.
Max.
0.1
<10−6
8.03 × 10−6
2.06 × 10−2
<10−6
2.43 × 10−5
6.16 × 10−1
0.5
<10−6
3.57 × 10−4
3.88 × 10−2
<10−6
3.25 × 10−4
2.15 × 10−1
1
<10−6
5.05 × 10−4
3.88 × 10−2
<10−6
2.22 × 10−4
4.21 × 10−3
5
<10−6
7.69 × 10−4
3.88 × 10−2
<10−6
4.46 × 10−4
6.95 × 10−1
10
<10−6
2.81 × 10−4
2.56 × 10−2
<10−6
2.42 × 10−4
4.22 × 10−1
50
<10−6
4.58 × 10−5
2.02 × 10−2
<10−6
3.51 × 10−4
3.85 × 10−1
100
<10−6
4.32 × 10−5
3.90 × 10−2
<10−6
4.13 × 10−4
9.12 × 10−2
500
<10−6
6.35 × 10−5
1.72 × 10−2
<10−6
2.70 × 10−5
5.13 × 10−1It can be observed from Table2 that AlgorithmOpt_Global performs better than grid search (GS) and two-dimensional DIRECT (TD). This verifies the optimality of Opt_Global. Moreover, if we naively apply TD to solve the problem, it yields the worst solution among the three algorithms.We also evaluate the computational efficiency through the computational time (on a PC with 3.4 GHz Intel CPU and 4 GB memory), which is presented in Table3.Table 3
Computational times of the three algorithms.
α
CPU time ofOpt_Global (s)
CPU time ofGS (s)
CPU time ofTD (s)
Min.
Avg.
Max.
Min.
Avg.
Max.
Min.
Avg.
Max.
0.1
0.41
1.90
2.85
41.14
44.41
53.01
2.902
17.92
4769
0.5
0.40
1.68
2.74
42.76
44.26
53.13
3.058
18.56
4957
1
0.41
1.92
3.14
42.38
44.19
53.27
3.125
18.66
5023
5
0.41
1.92
3.75
42.15
44.12
53.42
3.117
19.22
5129
10
0.39
1.87
4.00
42.19
42.50
53.85
3.211
17.84
4883
50
0.39
1.59
4.26
42.70
42.92
47.24
4.259
19.56
5212
100
0.25
1.83
4.19
42.23
42.89
52.30
3.568
18.29
5019
500
0.42
1.95
4.00
42.01
43.29
53.13
3.689
18.65
5217Table3 shows that the Opt_Global algorithm is much more efficient than GS and TD, in terms of computational time of CPU. In particular, TD may consume more than one hour to get the solution for some cases but obtain inferior results compared with that of Opt_Global(Table 2). Thus, it is reasonable to conclude that the two-layer structure in the proposed algorithm significantly improves the computational efficiency.Note that we have tried searching the optimal value ofF by smaller steps of 10-6 to improve the preciseness of the GS algorithm, but it cannot improve the solution significantly. What is more, under the search step of 10-6, it costs more than 4000 seconds for the computer to solve one example by GS and we totally need more than three years to complete the computation of all the numerical examples on our computer. Thus, a search step of 10-4 for F is used, with which all the numerical examples were solved within 20 days by the GS algorithm.


### 6.4. Optimal Solution Changing with Attenuation Coefficientα

We have shown that whenα→∞, the optimal solution of our model converges to the optima of the basic EOQ-PBO model (Corollary 7). Let T∞, F∞, and Γ∞ represent the optimal order cycle, the optimal fill rate, and the corresponding optimal inventory cost when α→∞, which can be obtained by solving the basic EOQ-PBO model [10]. Correspondingly, let Tα, Fα, and Γα indicate the optimal order cycle, the optimal fill rate, and the optimal inventory cost for a limited value of α<∞. Densify α by, respectively, setting α=0.1,0.25,0.5,1,3,5,10,25,50,75,100,300,500. Implementing the scheme in Section 6.2 to generate numerical examples, each value of α corresponds to A×Ch×Cb×Co×β×D=4×4×4×4×5×4=5120 problem instances. Solving all the numerical examples by Opt_Global for each α, the average, minimum, and maximum values of the ratio of Fα/F∞ and Tα/T∞ are presented in Figure 5.Figure 5
The ratio of the optimal solution,Fα/F∞ and Tα/T∞, under different α.It can be observed thatFα is decreasing in α until reaching F∞ when α→∞ (Fα/F∞→1), which has an intuitive interpretation: A smaller α corresponds to a smaller backorder demand rate so that the unsatisfied customer of the previous order cycle will fetch the backordered item later. In this case, the store must keep the backordered item for a longer time and incur a larger inventory holding cost. Thus, if α is smaller, the vendor should backorder less resulting in a larger fill rate, to avoid too much inventory holding cost of the backordered demand. Conversely, if α increases, the optimal fill rate Fα decreases.Moreover, another method for avoiding a high inventory holding cost caused by the backordered demand in case of a smallα is to diminish the maximum backorder level by promising a shorter stockout period. Intuitively, this can be realized by replenishing inventory more frequently based on a shorter order cycle. It follows that, under a smaller α, a shorter order cycle is more preferable. In other words, Tα is increasing in α, which can also be observed in Figure 5.In addition, Figure5 shows that, for some numerical examples, the optimal order cycle and the fill rate are always equal to their counterparts of the basic EOQ-PBO (Tα/T∞ and Fα/F∞ are equal to 1). We have checked these examples and found that if the basic EOB-PBO’s optimal policy is F∞=1, that is, to meet all demands without any lost sales, then Fα is also equal to 1, no matter what value of α is. In fact, since Fα is decreasing in α and F∞≤Fα≤1, it is an inevitable outcome that Fα must be equal to 1 in case of F∞=1.Figure6 presents the relative deviations of the optimal inventory cost Γα from Γ∞, where Γ∞ represents the optimal inventory cost of the basic EOQ-PBO model.Figure 6
The optimal inventory cost decreases inα until reaching that of the basic EOQ-PBO.We can see that the inventory cost is decreasing inα, and the optimal solution converges quickly to the case of α→∞ as α increases. In particular, when ln⁡α>3.4, that is, α>30, the optimal solution will be sufficiently close to that of the basic EOQ-PBO with a relative deviation less than 5%. That means if the backordered customers revisit the store without too much purchase delay, the basic EOQ-PBO can approximate our model very well.


## 6.1. Benchmark Algorithms

To evaluate the performance ofAlgorithmOpt_Global, we use two algorithms as benchmarks. The first one is a grid search (GS) algorithm for evaluating the optimality of AlgorithmOpt_Global, and the second one is a two-dimensional DIRECT algorithm [46] for evaluating the computational efficiency.


### 6.1.1. Grid Search onF~

We use grid search (GS) to approximately optimize function Γ^(F~) and compare the solution with that obtained by AlgorithmOpt_Global. For implementing the GS algorithm, we divide the interval F~∈[0,1] by steps of 10-4. For each problem instance, we need to optimize function Γ~T=Γ~T,F~ for ten thousand times by calling AlgorithmOpt_Local and keep the value of F~ with minimal inventory cost as the optimal solution.


### 6.1.2. Two-Dimensional DIRECT Algorithm

The problem can also be solved by using the two-dimensional DIRECT algorithm (TD) proposed by Jones et al. [46]. However, the two-dimensional split of the space can be very time-consuming if requiring a high precision. We here compare the performance of the Opt_Global algorithm with that of the TD algorithm to show how much Opt_Global will reduce the computational effort. The stop condition of the TD algorithm is the same as that of AlgorithmOpt_Global; that is, the relative improvement keeps smaller than 10-6 for 10 consecutive iterations.


## 6.1.1. Grid Search onF~

We use grid search (GS) to approximately optimize function Γ^(F~) and compare the solution with that obtained by AlgorithmOpt_Global. For implementing the GS algorithm, we divide the interval F~∈[0,1] by steps of 10-4. For each problem instance, we need to optimize function Γ~T=Γ~T,F~ for ten thousand times by calling AlgorithmOpt_Local and keep the value of F~ with minimal inventory cost as the optimal solution.


## 6.1.2. Two-Dimensional DIRECT Algorithm

The problem can also be solved by using the two-dimensional DIRECT algorithm (TD) proposed by Jones et al. [46]. However, the two-dimensional split of the space can be very time-consuming if requiring a high precision. We here compare the performance of the Opt_Global algorithm with that of the TD algorithm to show how much Opt_Global will reduce the computational effort. The stop condition of the TD algorithm is the same as that of AlgorithmOpt_Global; that is, the relative improvement keeps smaller than 10-6 for 10 consecutive iterations.


## 6.2. Outline of Numerical Examples

We investigate the performance of our model and algorithm under different scenarios by solving numerical examples. The parameter of the numerical examples should cover different cases as many as possible. Thus, we enumerate different values of parameters to construct different examples, referring to the parameter scheme of the major item of the numerical examples in Zhang et al. [49], as follows.(i) Ordering Cost. We consider four values, A = $100/order, $1000/order, $2500/order, and $5000/order, respectively, representing the difference in lower and higher ordering cost.(ii) Inventory Holding Cost and Backordering Cost. Set Ch = $5/unit/year, $10/unit/year, $25/unit/year, and $50/unit/year representing the inventory holding cost changes from low to high. Similarly, we set backorder cost as Cb = $5/unit/year, $10/unit/year, $25/unit/year, or $50/unit/year, respectively.(iii) Opportunity Cost. Set Co =$5/unit, $10/unit, $25/unit, and $50/unit, respectively.(iv) Backordering Rate. Backordering rate β is changed from 0.1 to 0.9 with a step of 0.2.(v) Demand. Set the routine demand rate D = 100 units/year, 1000 units/year, 5000 units/year, and 10000 units/year, respectively.(vi) Attenuation Coefficient of the Backorder Demand Rate. Set α=0.1,0.5,1,5,10,50,100,500, representing the backorder demand rate attenuates from slow to fast.By combination ofA×Ch×Cb×Co×β×D×α=4×4×4×4×5×4×8 we generate 40960 different numerical examples for computational experiments.


## 6.3. Performance of the Algorithm

We apply, respectively, the three algorithms:AlgorithmOpt_Global, grid search (GS), and the two-dimensional DIRECT (TD) algorithm to solve the generated 40960 numerical examples. Let ΓOpt, ΓGS, and ΓTD denote the best value of the objective function obtained by Opt_Global, GS, and TD, respectively. Let ΓGS-ΓOpt/ΓOpt indicate the relative deviation of the solution found by GS from that obtained by the Opt_Global algorithm, and ΓTD-ΓOpt/ΓOpt indicate the relative deviation of TD. We group all the numerical examples by the attenuation coefficient α, and the computational results are presented in Table 2.Table 2
Relative deviations of the two benchmarks.
α
Dev. ofΓGS (%)
Dev. ofΓTD (%)
Min.
Avg.
Max.
Min.
Avg.
Max.
0.1
<10−6
8.03 × 10−6
2.06 × 10−2
<10−6
2.43 × 10−5
6.16 × 10−1
0.5
<10−6
3.57 × 10−4
3.88 × 10−2
<10−6
3.25 × 10−4
2.15 × 10−1
1
<10−6
5.05 × 10−4
3.88 × 10−2
<10−6
2.22 × 10−4
4.21 × 10−3
5
<10−6
7.69 × 10−4
3.88 × 10−2
<10−6
4.46 × 10−4
6.95 × 10−1
10
<10−6
2.81 × 10−4
2.56 × 10−2
<10−6
2.42 × 10−4
4.22 × 10−1
50
<10−6
4.58 × 10−5
2.02 × 10−2
<10−6
3.51 × 10−4
3.85 × 10−1
100
<10−6
4.32 × 10−5
3.90 × 10−2
<10−6
4.13 × 10−4
9.12 × 10−2
500
<10−6
6.35 × 10−5
1.72 × 10−2
<10−6
2.70 × 10−5
5.13 × 10−1It can be observed from Table2 that AlgorithmOpt_Global performs better than grid search (GS) and two-dimensional DIRECT (TD). This verifies the optimality of Opt_Global. Moreover, if we naively apply TD to solve the problem, it yields the worst solution among the three algorithms.We also evaluate the computational efficiency through the computational time (on a PC with 3.4 GHz Intel CPU and 4 GB memory), which is presented in Table3.Table 3
Computational times of the three algorithms.
α
CPU time ofOpt_Global (s)
CPU time ofGS (s)
CPU time ofTD (s)
Min.
Avg.
Max.
Min.
Avg.
Max.
Min.
Avg.
Max.
0.1
0.41
1.90
2.85
41.14
44.41
53.01
2.902
17.92
4769
0.5
0.40
1.68
2.74
42.76
44.26
53.13
3.058
18.56
4957
1
0.41
1.92
3.14
42.38
44.19
53.27
3.125
18.66
5023
5
0.41
1.92
3.75
42.15
44.12
53.42
3.117
19.22
5129
10
0.39
1.87
4.00
42.19
42.50
53.85
3.211
17.84
4883
50
0.39
1.59
4.26
42.70
42.92
47.24
4.259
19.56
5212
100
0.25
1.83
4.19
42.23
42.89
52.30
3.568
18.29
5019
500
0.42
1.95
4.00
42.01
43.29
53.13
3.689
18.65
5217Table3 shows that the Opt_Global algorithm is much more efficient than GS and TD, in terms of computational time of CPU. In particular, TD may consume more than one hour to get the solution for some cases but obtain inferior results compared with that of Opt_Global(Table 2). Thus, it is reasonable to conclude that the two-layer structure in the proposed algorithm significantly improves the computational efficiency.Note that we have tried searching the optimal value ofF by smaller steps of 10-6 to improve the preciseness of the GS algorithm, but it cannot improve the solution significantly. What is more, under the search step of 10-6, it costs more than 4000 seconds for the computer to solve one example by GS and we totally need more than three years to complete the computation of all the numerical examples on our computer. Thus, a search step of 10-4 for F is used, with which all the numerical examples were solved within 20 days by the GS algorithm.


## 6.4. Optimal Solution Changing with Attenuation Coefficientα

We have shown that whenα→∞, the optimal solution of our model converges to the optima of the basic EOQ-PBO model (Corollary 7). Let T∞, F∞, and Γ∞ represent the optimal order cycle, the optimal fill rate, and the corresponding optimal inventory cost when α→∞, which can be obtained by solving the basic EOQ-PBO model [10]. Correspondingly, let Tα, Fα, and Γα indicate the optimal order cycle, the optimal fill rate, and the optimal inventory cost for a limited value of α<∞. Densify α by, respectively, setting α=0.1,0.25,0.5,1,3,5,10,25,50,75,100,300,500. Implementing the scheme in Section 6.2 to generate numerical examples, each value of α corresponds to A×Ch×Cb×Co×β×D=4×4×4×4×5×4=5120 problem instances. Solving all the numerical examples by Opt_Global for each α, the average, minimum, and maximum values of the ratio of Fα/F∞ and Tα/T∞ are presented in Figure 5.Figure 5
The ratio of the optimal solution,Fα/F∞ and Tα/T∞, under different α.It can be observed thatFα is decreasing in α until reaching F∞ when α→∞ (Fα/F∞→1), which has an intuitive interpretation: A smaller α corresponds to a smaller backorder demand rate so that the unsatisfied customer of the previous order cycle will fetch the backordered item later. In this case, the store must keep the backordered item for a longer time and incur a larger inventory holding cost. Thus, if α is smaller, the vendor should backorder less resulting in a larger fill rate, to avoid too much inventory holding cost of the backordered demand. Conversely, if α increases, the optimal fill rate Fα decreases.Moreover, another method for avoiding a high inventory holding cost caused by the backordered demand in case of a smallα is to diminish the maximum backorder level by promising a shorter stockout period. Intuitively, this can be realized by replenishing inventory more frequently based on a shorter order cycle. It follows that, under a smaller α, a shorter order cycle is more preferable. In other words, Tα is increasing in α, which can also be observed in Figure 5.In addition, Figure5 shows that, for some numerical examples, the optimal order cycle and the fill rate are always equal to their counterparts of the basic EOQ-PBO (Tα/T∞ and Fα/F∞ are equal to 1). We have checked these examples and found that if the basic EOB-PBO’s optimal policy is F∞=1, that is, to meet all demands without any lost sales, then Fα is also equal to 1, no matter what value of α is. In fact, since Fα is decreasing in α and F∞≤Fα≤1, it is an inevitable outcome that Fα must be equal to 1 in case of F∞=1.Figure6 presents the relative deviations of the optimal inventory cost Γα from Γ∞, where Γ∞ represents the optimal inventory cost of the basic EOQ-PBO model.Figure 6
The optimal inventory cost decreases inα until reaching that of the basic EOQ-PBO.We can see that the inventory cost is decreasing inα, and the optimal solution converges quickly to the case of α→∞ as α increases. In particular, when ln⁡α>3.4, that is, α>30, the optimal solution will be sufficiently close to that of the basic EOQ-PBO with a relative deviation less than 5%. That means if the backordered customers revisit the store without too much purchase delay, the basic EOQ-PBO can approximate our model very well.


## 7. Conclusions

It is frequent that if demands are not filled and the unsatisfied customers are willing to accept the backorder, she/he will fetch the demanded item in the next order cycle. The basic EOQ with partial backordering assumes that all backordered customers will fetch their products simultaneously at the beginning of the very next order cycle and the store will not incur inventory holding cost for the backordered demand. However, this assumption may be violated in some cases where the backordered customer may revisit the store successively because of the purchase delay, and therefore the backordered demand is not filled instantly.This paper proposed a new inventory model with partial backordering by assuming that the backordered customers repurchase the backordered products successively and the arriving rate of the backordered customer is proportional to the backordered customers who have not been filled. This leads to an exponentially decreasing demand rate in time regarding the backordered demand, based on which we formulated the total inventory cost function. To solve the model, we developed an efficient algorithm based on Lipschitz Optimization and conduct numerical experiments proving the efficiency and effectiveness of the algorithm.From the theoretical results and numerical computations, we mainly find that (1) the larger the attenuation coefficient of the backorder demand rate, the lower the fill rate with which the store will satisfy the customer to avoid holding backordered items for a longer time resulting in higher inventory holding cost; (2) the more slowly the backordered customer revisits the store, the shorter the order cycle that the store should set is, to diminish the backorder level and therefore decrease the inventory holding cost of the backordered demand; (3) a larger backorder demand attenuation coefficient leads to a smaller total inventory cost, so the store should deliver the item to the backordered customers as soon as possible.Compared with the traditional EOQ-PBO model, the new parameter needed when implementing our model in the real-world is the attenuation coefficientα. α reflects the proportion of the backordered customer who has revisited the store and fetched the backordered item. To estimate this parameter, the vendor need record the number of the backordered customers who revisit the store per unit time, which is compared with the total number of the backordered customers of the previous stockout. Thereby, the attenuation coefficient α can be calibrated.The study can be extended by combining the limited backorder demand rate considered in this paper with other inventory problems, for example, deteriorating items, time-varying demand rate, and the vendor-retailer supply chain, to further enrich the study of EOQ-PBO models. The consumer behavior model describing how the backordered customers repurchase their demanded items also deserves more researches. In addition, this paper mainly concentrates on the theoretical analysis of the new inventory problem. To implement a theoretical inventory model based on practical data is very important for the research on operations management, which will be carried out as the empirical study in the future work.



---

*Source: 6425403-2016-02-10.xml*