NEW ALGORITHMS FOR (Q, r) SYSTEMS WITH COMPLETE BACKORDERING USING A FILL-RATE CRITERION Candace Arai Yano Department of Industrial & Operations Engineering The University of Michigan Ann Arbor, Michigan 48109, USA January 1985 Technical Report 85-1 REVISED TECHNICAL REPORT 84-17

NEW ALGORITHMS FOR (Q, r) SYSTEMS WITH COMPLETE BACKORDERING USING A FILL-RATE CRITERION Candace Arai Yano University of Michigan Department of Industrial & Operations Engineering Ann Arbor, Michigan

New Algorithms for (Q, r) Systems With Complete Backordering Using a Fill-Rate Criterion ABSTRACT We present an algorithm which determines optimal parameter values for order quantity-reorder point systems with complete backordering. The service level is measured as fraction of demand satisfied directly from shelf, also known as "fill-rate." This algorithm differs from existing algorithms because an exact cost function is used rather than an approximation. We also present a new heuristic algorithm which is more efficient computationally than the optimal procedure and provides excellent results. Results of extensive computational experience also are reported.

New Algorithms for (Q, r) Systems with Complete Backordering Using a Fill-Rate Criterion 1.0 INTRODUCTION Continuous review order quantity-reorder point inventory models, typically called (Q, r) models, have been studied extensively. These models differ in the assumptions about backordering, and in the service level criterion used. In many cases, closed form solutions cannot be obtained or iterative algorithms must be used to converge on optimal solutions. This difficulty of implementation for the practitioner has led to proposals for and evaluation of simple heuristic or approximate solutions. Although the cost of reviewing inventory continuously may be expensive, there are many situations in which (1) computer systems already exist to provide continuous or nearly continuous review of inventory levels; (2) the order quantity itself is important because of quantity discounts, package sizes, minimum order quantities, or administrative considerations; or (3) leadtimes are long and/or highly uncertain, making even approximate solutions for periodic review policies difficult to compute using dynamic programming techniques. We note that it is possible to model a periodic review system as a continuous review system in which the random leadtime includes the period between the time the order point is reached and the time an order is placed. Continuous review systems have the desirable characteristics of being fairly robust to errors in cost estimation and relatively easily incorporating variable leadtimes. Generally operating parameters for continuous review systems are much easier to compute than for periodic review systems, even if an iterative routine must be used. Therefore, continuous review systems will 1

continue to be used and may become more prevalent as the cost of review decreases with advancing technology. Hadley and Whitin [3] analyze the case of complete backordering with a penalty for each unit backlogged. They use an estimate of on-hand inventory which does not incorporate the effect of backorders which must be filled by incoming orders. Their formulation also uses a very simple approximation for average on-hand inventory during the leadtime. Hadley and Whitin suggest an iterative algorithm which derives from two non-linear equations in two unknowns. Holt, Modigliani, Muth and Simon [4] study two situations with complete backlogging. The first is the same as the Hadley-Whitin problem. For this problem they use an estimate of on-hand inventory which distinguishes between cycles in which stockouts do occur and those in which they do not occur. The second formulation charges a shortage cost per unit per unit time using the same approach to estimate inventory levels. For both problems, they derive the optimal reorder point as a function of the order quantity. Parker [6] studies single and multi-item continuous review systems in which (1) a penalty per unit short is assessed and (2) a constraint is placed on fraction of demand backordered. Using an approximate cost function, he derives formulas for the order quantity and reorder point(s) for the two objectives. Presutti and Trepp [7] develop explicit formulas for order quantities and reorder points using models with two different aggregate multi-item service constraints: (1) total number of units backordered per year, and (2) average number of units in a backorder position. The importance of each item is reflected as a weight in the aggregate constraints. Formulas are developed for models in which holding costs are assessed on either on hand inventory or the inventory position. Schrady and Choe [9] develop techniques for multi-item inventory systems in which the objective is to minimize time-weighted shortages. The problem is 2

constrained by limits on the total average inventory investment and number of orders placed per unit time. Gross and Ince [2] evaluate the algorithms proposed by Hadley and Whitin, and by Wagner et al. [12], along with a single-pass (non-iterative) algorithm which uses the standard Wilson lot-size formula. They assume a Poisson demand process which, because of its memory-less property, satisfies the assumption of (Q, r) models that demand over the leadtime is independent of when the reorder point is reached. Their results show that the "single-pass" approach performs poorly. The Hadley-Whitin and Wagner, et al. approaches perform well, generally within 5% of optimal. The Hadley-Whitin algorithm performed slightly better than the Wagner algorithm in many cases, because of compensating errors in estimating on-hand inventory. Silver and Wilson [10] develop approaches for determining the percentage cost penalty of using a "single-pass" solution rather than iterative or simultaneous solutions. This is done for four service criteria: (1) fixed cost per stockout occasion, (2) fixed cost per unit short, (3) probability of stockout per replenishment cycle, and (4) fraction of demand satisfied directly from shelf. Nahmias [5] demonstrates the equivalence of three approximate (Q, r) models. One model assesses a penalty per unit short, the second model imposes a constraint on the fill-rate, and the third model specifies a constraint on the probability of stockout during the leadtime. He also provides streamlined versions of related algorithms proposed by Brown [1], Hadley and Whitin [3], and Wagner [11]. Section 2 describes the formulation of the exact cost model unde the assumptions of complete backordering and a service level criterion of fraction of demand satisfied directly from shelf, or "fill-rate." Section 3 details the 3

development of a new heuristic algorithm. Section 4 compares the new heuristic algorithm with the Silver-Wilson heuristic from an analytical viewpoint. Section 5 presents computational experience with both algorithms, and Section 6 concludes with a summary and discussion. 2.0 DEVELOPMENT OF THE EXACT COST MODEL AND ALGORITHM The problem objective involves minimizing the average annual sum of setup and inventory holding costs subject to a constraint on the level of service. Service is measured here as the fraction of demand satisfied directly from shelf, also known as "fill-rate." We assume that the leadtime is known and constant. We also assume that demand is stationary and that demand during the leadtime can be described by some probability distribution. Let Q = order quantity (a decision variable) r = reorder point (a decision variable) D = annual demand in units S = setup costs T = leadtime in years f (.) = density of demand during the leadtime F (.) = cumulative distribution of leadtime demand b (r) = expected number of backorders per replenishment cycle, given reorder point r o = standard deviation of demand during the leadtime 1- a= service level required The optimization problem using an exact cost representation (see [3] for exact representations of average inventory and expected backorders) can be stated as follows: 4

Minimize SD/Q + h{Q/2 + r-p + [8(r)-6 (r+Q)]/Q} r+Q s.t. [l1-F(x)] dx/Q O a r where 6(v) = (t-v) [1-F(t)] dt The cost function is obtained using steady-state analysis of continuous time Markov chains, as done by Hadley and Whitin ([3], pp. 181-188). Minimization of costs will drive the constraint to an equality. Simplifying the objective function and restating the problem using a Lagrangian yields: r+Q Minimize L = SD/Q + h {Q+r - / (x-r)F(x)dx/Q r oo r+Q I [l-F(x)]dx} + x [ [l-F(x)]dx/Q - a] (1)! r+Q J r The first order necessary conditions are: *r+Q / r+Q?L/DQ = {-SD + h] (x-r)F(x)dx + X [ F(x)dx - QF(r+Q)]}/Q2=O (2) r+Q DL/Dr = {h F(x)dx - X [F(r+Q) - F(r)]}/Q = 0 (3) r fr+Q aL/DX = [l-F(x)]dx /Q = a It can be shown that the original problem involves minimization of a convex function on a convex set. Convexity of B(Q,r) = [S(r) - 3(r + Q)]/Q has been established by Zipkin [13], over the range in which safety stock is positive and convexity of the remainder of the objective function is well-known. Convexity of the feasible set can be demonstrated graphically, as shown in Figure 1 for x - N(p,o2). Under these conditions, first order conditions are sufficient for optimality. 5

FIGURE 1 We observe that solving for X in equation (3) and substituting for r+Q F (x) by rearranging equation (4) permits us to eliminate X and to simplify J r equation (1) to r+Q (x-r)F(x)dx = r SD/h + Q2 (1-a) [F(r+Q) + a-l]/[F(r+Q)-F(r)] (5) Given any value of r, we can find Q using a one-dimensional search. Similarly, given any value of Q, we can find r satisfying the constraint using equation (4) or its simplified equivalent: rr+Q F(x)dx = (1-a) Q (6) I r The generalized reduced gradient technique is ideally suited to solve this problem. Starting with any feasible value of r, one first solves for a new Q using equation (5), then a new value of r using equation (6). Iteration continues until the two equations are satisfied simultaneously. Solving equation (5) is equivalent to a steepest descent step while equation (6) insures that the constraint is satisfied. 2.1 GLOBAL CONVERGENCE OF THE ALGORITHM If the algorithm described above converges, it converges to the global minimum since the first order conditions are sufficient. We next demonstrate that the algorithm indeed converges. Let g(Q,r) = {(Q,r) satisfying (5)} and 6

h(Q,r) = {(Q,r) satisfying (6)}. It can be shown (see Appendix) that aQ/Dr (g(Q,r)) < aQ//r(h(Q,r)) | (7) Therefore, a situation depicted in Figure 2 in which a "hog cycle" (divergent sequence of points) would occur is guaranteed not to arise. Further, since (7) holds as a strict inequality, convergence is guaranteed, as depicted in Figure 3. FIGURES 2 AND 3 3.0 A NEW HEURISTIC ALGORITHM Much of the computational effort required for the exact procedure described in Section 2 involves computing numerically integrals for which tabled values are not available, even for standardized distributions. Using an approximate but simpler representation of average inventory significantly reduces the computational effort. We therefore chose to approximate average cycle stock as (Q-b(r))2/2Q This approximation is based on the assumptions that demand occurs at a constant rate and that b(r) backorders are outstanding, on average, prior to the arrival of an order. An approximate cost model can be formulated as follows: Minimize SD/Q + h[(Q-b(r))2/2Q + (r-x)f(x)dx] (8) -00 subject to 5(r)/Q < a (9) 00 where b(r) = (x-r) f(x) dx r Minimization of costs will drive the constraint to an equality. Therefore, without loss of generality, (2) can be replaced by 7

r00 f (x-r) f(x)dx - Qa = 0 (9') r Note that cycle stock has been represented by (Q - b(r))2/2Q The average inventory level is -b(r) immediately before an order arrives. Therefore, stock available upon arrival of an order is Q -b(r) on average. On average, stock is available approximately a fraction of the cycle equal to (Q-b(r))/Q. Average safety stock is expressed as r (r-x) f(x) dx v -00 which is a commonly used representation of safety stock that is somewhat more accurate than r-u. The problem can be restated as a non-linear optimization problem using a Lagrangian as follows: L = SD/Q + h[(Q-b(r))2/2Q + (r-x)f(x) dx] (10) - _00 + x [(r) - Qa] Taking the partial derivatives of L with respect to Q, r, and X and setting them equal to zero yields: 3L/DQ = -SD/Q2 + h/2 - h[b(r)]2/2Q2 -Xa =0 (11) 9L/ r = -hb(r)[1-F(r)]/Q + h -X[l-F(r)] = 0 (12) aL/aX = b(r) - Qa = 0 (13) Solving for X in (12) and substituting for X and b(r) in (11) yields: 8

Q = {2SD[l-F(r)]/h [(la>2) [1-F(r)]-2a}-5 (14) or Q = {2SD[l-F(r)]/h [(1-a)2 - (l+a2)F(r)]}'5 The form of the equation for Q provides for a very simple iterative approach for simultaneously determining Q and r. The algorithm is initialized by setting r = V = mean demand during the leadtime which is equivalent to setting F(r) =.5 when the density of demand during the leadtime is symmetric. Heuristic Algorithm 1. Set j = O Set ko = 0 Set F(ro) =.5 2. Solve Qo = {2SD [l-F(O)]/h[(l+a2) (l-F(O)) - 2a]}' 5 3. Set j = j+l 4. Find ri such that b(rj)/Q__1 = 5. Qj = {2SD [l-F(rj)]/h[(l+a2) (l-F(r)) - 2a]}-5 6. Determine k = (rj- )/G 7. If IQ, - Qj- or rj - rjl or kj - kj is sufficiently small, stop. Otherwise go to step 3. Note that if leadtime demand is distributed N(p,a2), then step 4 can be simplified to read as follows: Find rj such that E [(rj - p)/a] = Qjlao/ where E(.) is the standard normal loss function. 9

3.1 CONVEXITY AND CONVERGENCE Let L represent the Hessian of the Lagrangian. It can be shown (see Appendix) that L is positive definite for all non-negative Q. Therefore, the solution of the first order necessary conditions in (5), (6) and (7) gives a unique global minimum. The algorithm proposed in this paper simply utilizes these conditions to achieve their simultaneous solution. The algorithm is simply an adaptation of the generalized reduced gradient technique, as is the algorithm for the exact cost model. It can be shown that the heuristic algorithm is guaranteed to converge from a procedural viewpoint. (The proof parallels that of the exact algorithm in Section 2). However, additional conditions must be satisfied to insure that at each iteration Qj is positive and real-valued. In order for convergence to be assured, we must have (l+a2) (l-F(rj)) - 2a > 0 or equivalently F(rj) (1-a)2/(1+a2) which, in turn, insures that Q is positive and real valued. This also implies a necessary relationship between the standard deviation of leadtime demand the standard economic order quantity (EOQ). Consider the case of normally distributed leadtime demand. Let k = (r-D)/a where p =T D = mean leadtime demand and a = standard deviation of leadtime demand. Then F(r) < (1-a)2/(1+a2) is equivalent to -(k) < (1-a)2/(l+a2) Let k' be such that m(k') = (1-a)2/(1x2). We know that Qj = a E(kj)/a for all j by virtue of the service level constraint. But 10

Qj = EOQ [(l-m(kj))/(l-$(kj)-2a)].5 Therefore, E(kj) = a(EOQ) [(l-m(kj))/((l+a2)(1-m(kj)-2a))]'5/]LT Now E(k') = E{~-1 [(1-c)2/(l+a2)]} by definition of k'. We must have E(kj) > E(k') for all j, or EOQ, E {~-1 [C(1-a)2/(l+a2)]} a [(1-D (kj))/((l+a2) (1-(kj)-2a)]'5 For constant a, the right hand side achieves its maximum at the smallest permissible value of k (-cc if unrestricted or 0 if safety stock is constrained to be non-negative). In the unrestricted case the condition reduces to: EOQ > E{d-l[(1-a)2/(l+a2)]} (1-a)/ and to the following when k must be non-negative: EOQ > E{V-1 [(1-a)2/(l+a2)} (a2-4a+ l)/a The right hand sides are plotted as a function of a in Figures 4 and 5 respectively. FIGURES 4 AND 5 4.0 COMPARISON WITH SILVER-WILSON HEURISTIC Silver and Wilson [10] develop a heuristic for the problem posed in this paper, using an estimate of on-hand inventory equal to: Q/2 + (r-p) where P = mean demand during the leadtime. The first term overestimates average on-hand cycle stock while the second term underestimates safety stock. Silver and Wilson solve for Q and r using a system of two non-linear equations in two unknowns. It is clear, however, that the same approach as 11

described in Section 3 can be used to develop an algorithm with their approximate cost function. Such an approach gives Q = {2SD(l-F(r))/[h(l-F(r) - 2a)]}.5 This results in a value of Q larger than in our algorithm, hence a smaller value of r as well. Note that the computations required to execute the new algorithm differs from the computational requirements of the Silver-Wilson heuristic by one addition and three multiplications per iteration. Calculation of the term 1 - F(r) -2a requires one multiplication and two subtractions. The term (l1+a2) (l-F(r)) -2a requires one addition, four multiplications and two subtractions. 5.0 COMPUTATIONAL EXPERIENCE We examined 1440 problems with a wide range of cost parameters, leadtimes, leadtime demand variability, and service levels in order to gain some empirical evidence of the performance of the speed of convergence of the new algorithms. For these problems we assume that leadtime demand is distributed N(i,a2). Table 1 lists the possible values of the parameters. The 1440 problems represent all combinations of these parameter values. The item cost, C, is normalized to 1.0, and i, the annual holding cost rate, takes on values.20,.25,.30,.35. Annual demand D, is set so that the Economic Order Quantity, given S and h, is equal to 1000 or 5000. Hence D = 500,000h/S or 12,500,000 h/S. TABLE 1 12

The leadtime may be.02,.04,.08, or.16 of a year, corresponding to approximately, one, two, four, or eight weeks, respectively. The coefficient of variation of leadtime demand is allowed to take a maximum value of.40 to insure that the distributions represented reflect primarily positive demand. There are several objectives of this computational study. The first, of course, is to compare the solutions from the heuristic with the optimal solution. The second objective is to determine how quickly each algorithm converges. It was necessary to eliminate from consideration many of the problems because either (1) the EOQ/a ratio indicated that the heuristic algorithm would not converge, (2) the safety stock factor would be negative so that the optimal algorithm could not be applied, or (3) the combinations of parameters generated two or more essentially equivalent problems (different ps but equal as). We decided to use relevant problems for S=10 to compare the results of the two algorithms, primarily because for S=100 and S=1000, nearly all the problems would be eliminated for one of the 3 reasons above. There were 94 problems available for this comparison. For convergence results, we chose to use all available problems; there were 803 problems for the heuristic and 94 for the exact cost procedure using S=10. It is important to note that because of the form of equation (5) and availability of standard normal probability tables in increments of only.01, it was possible to determine the value of Q using the exact algorithm only to the nearest.01. As such, we probably have not obtained exactly optimal values of Q and r. However, a search would have been computationally prohibitive, as it would be in real applications. Therefore, the comparisons described here may realistically represent what would be achieved in practice. We calculated the values of: 13

Cost heuristic -Cost "optimal" Cost "optimal" using the exact cost function for both cost computations. We found that all but 16 of these values are within the range of (-.01,.01) and of the 16 remaining values the largest was.02 and 12 were negative (i.e., the heuristic outperformed the "exact" procedure). We also observed that the heuristic achieved the specified fill-rate in all cases. Thus, the heuristic appears to perform extremely well relative to the exact procedure as it would be applied in practice. Convergence results for the heuristic are displayed in Table 2 for 803 problems. The algorithm achieved (and verified) convergence in 5 iterations or less in over 94% of the problems. An examination of the problems in which convergence was achieved slowly revealed a correlation between number of iterations and the proximity of the EOQ/o ratio to its maximum permitted value. TABLE 2 The exact procedure generally converged in 3 or 4 iterations, as indicated in Table 3. Unfortunately, however, this procedure required several times the computation time of the heuristic, despite fewer iterations on average. (The problems were run in BASIC on the IBM PC/XT and therefore detailed CPU times are not available. However, it was observed that the heuristic generally solved the problems in much less than a minute, while the optimal procedure required a few to several minutes). TABLE 3 14

6.0 CONCLUSIONS The heuristic developed in Section 3 performs quite well when compared with approximate solutions to the exact cost model whose optimal solutions would be computationally prohibitive to obtain. It has the further advantage of being applicable where negative safety stock is appropriate. However, it cannot be used when a is large relative to the EOQ, and in such situations, it is advisable to use another heuristic. When a is large, only very approximate values of Q order can be obtained from the exact cost procedure. Further research is needed to develop appropriate procedures to determine values of Q and r when a > EOQ. 15

k 2.0 1.5 \ ^c a =.02, a = 500.5 =.01, a = 100 0 k i i i 5 10 15 20 Q (X 00) Figure 1 Feasible Regions 16

g(Q,r) ^^^ ^^^h(Qr) Q Figure 2 Divergent Sequence 17

| g(Q,r) h(Q,r) Q Figure 3 Convergent Sequence 18

EOQ 1.20 1.00,.80,.60.40.20 0 ____i i i _ 0.10.20.30.40 a Figure 4 Maximum Value of EOQ/a Guaranteeing Convergence (Safety Stock Unconstrained) 19

EOQ 1.0.80.60.40.20 0.10.20.30 a Figure 5 Maximum Value of EOQ/o Guaranteeing Convergence (Non-negative Safety Stock) 20

Table 1 Parameter Values Parameter Values S 10, 25, 100, 500, 1000 h.20,.25,.30,.35 a.02,.05 L.02,.04,.08,.16 D 500,000 h/S, 12,500,000 h/S o/p.10,.25,.40 21

Table 2 Convergence Results for Heuristic Procedure No. of Iterations No. of Problems Cumulative % 1 18 2.2 2 134 18.9 3 508 82.2 4 56 89.2 5 41 94.3 6 10 95.5 7 12 97.0 8 11 98.4 9 1 98.5 10 0 11 3 98.9 12 0 13 7 99.8 14 0 15 0 16 2 100.0 Total 803 22

Table 3 Convergence Results for Exact Cost Procedure No. of Iterations No. of Problems Cumulative % 2 5 5.3 3 60 69.1 4 26 96.8 5 3 100.0 Total 94 23

Appendix We want to show that aQ/ar [g(Q,r)] < | aQ/ar[h(Q,r)] | where g(Q,r) - (Q,r) satisfying r+Q SD/h = Q (x-r) F(x)dx r - Q2(1-a) [F(r+Q) + a-l]/[T(r+Q) - F(r)] (15) ana h(Q,r) - (Q,r) satisfying r+Q F(x)dx/Q = 1 - a r Observe that since we want to find aQ/ar (.), it is desirable to express Q as a function of r. It is then straightforward to show that r+Q aQ/ar [g(Q,r)] = - {QF(r+Q) -J F(x)dx r + Q2[F(r+Q) (f(r) - f(r+Q)) - (l-a)f(r)]}/ [Q{F(r+Q) + 2 (1-a)2 [F(r+Q) - F(r)] - 2(1-a) [F(r+Q) - F(r)] F(r+Q)} - Q2 F(r+Q) f(r+Q)] and aQ/Dr [h(Q,r)] = fr+Q Q{F(r+Q) - F(r)]/[QF(r+Q) - F(x)dx] r Some simplification and utilization of inequalities such as QF(r) <f / F(x)dx < QF(r+Q) will yield the desired result. 24

We also need to demonstrate that the Lagrangian of the approximate cost function is convex. L, the Hessian of the Lagrangian is: {2SD + h [b(r)]2}/Q3 hb(r)[l-F(r)]/Q2 hb(r)[l-F(r)]/Q2 h{b(r)f(r) + [1-F(r)]2}/Q + Xf(r) It can be demonstrated easily that L is positive definite for all non-negative Q. 25

References [1] Brown, R. G., Decision Rules for Inventory Management, New York: Holt, Rinehart and Winston, 1967. [2] Gross, Donald and John F. Ince, "Comparison and Evaluation of Approximate Continuous Review Inventory Models," Int. 3. Prod. Res., Vol. 13, No. 1 (1975), pp. 9-23. [3] Hadley, G. and T. M. Whitin, Analysis of Inventory Systems, Englewood Cliffs, NJ: Prentice-Hall, Inc., 1963. [4] Holt, Charles C., Franco Modigliani, John F. Muth, and Herbert A. Simon, Planning Production, Inventories, and Work Force, Englewood Cliffs, NJ: Prentice-Hall, Inc., 1960. [5] Nahmias, Steven, "On the Equivalence of the Approximate Continuous Review Inventory Models, NRLQ, Vol. 23, No. 1 (1976), pp. 31-36. [6] Parker, L. L., "Economical Order Quantities with Reorder Points with Uncertain Demand," NRLQ, Vol. II, No. 4 (1964), pp. 351-358. [7] Presutti, V. and R. Trepp, "More Ado About Economic Order Quantities," NRLQ, Vol. 17, No. 2 (1970), pp. 243-252. [8] Schneider, H. "Methods for Determining the Reorder Point of an (s,S) Ordering Policy when a Service Level is Specified," JORS, Vol. 29, No. 12 (1978), pp. 1181-1194. [9] Schrady, D. A. and V. C. Choe, "Models for Multi-Item Continuous Review Policies Subject to Constraints," NRLQ, Vol. 18, No. 4 (1971), pp. 451-464. [10] Silver, Edward A. and Terry G. Wilson, "Cost Penalties of Simplified Procedures for Selecting Reorder Points and Order Quantities," American Production and Inventory Control Society Conference Proceedings, 1975, pp. 219-233. [11] Wagner, H. M., Principles of Operations Research, Englewood Cliffs, New Jersey: Prentice-Hall, Inc. 1969. [12] Wagner, H., M. O'Hagen, and B. Lundh, "An Empirical Study of Exactly and Approximately Optimal Inventory Policies," Management Science, Vol. 11, No. 7 (1965), pp. 690-723. [13] Zipkin, Paul, "Inventory Service Level Measures: Convexity and Approximation," Research Working Paper 84-9, Graduate School of Business, Columbia University, 1984. 26