2.1. Notation

Symbols that are used throughout the sections are summarized here, whereas those that are used briefly will be defined locally.

2.2. A Periodic Review Inventory System

We consider an infinite horizon, periodic-review inventory system that faces discrete-valued independent and identically distributed (i.i.d.) demands. Demand that is not met immediately is fully backlogged and will be fulfilled on a first-come–first-served basis. The system replenishes inventory according to an (s,S) policy with a constant lead time of L periods. (The random lead time will be introduced after the basic model has been formulated.) The sequence of events in each period is as follows: The inventory is reviewed, an order that was placed L periods ago is received, an ordering decision is made, and the demand in the period is then subtracted out of the inventory. These events are supposed to happen at the beginning of each period.

It is well known (e.g., Veinott and Wagner [18], and ZF [20]) that

For ease of analysis, we impose a minor condition on G(y), which is also assumed in most of the previous work:

Assumption 1: lim_y→±∞ G(y) > K + G(y₀).

It should be noted that when L is a random variable, G(y) should be understood as the expected loss function over demand and the lead time

Under the conditions that are explicated by Zipkin [21], the distribution of demand over the deterministic lead time L is just replaced by the distribution of total demand over the stochastic lead time. The conditions are met when the system faces an exogenous, sequential supply source in which the lead times are i.i.d. over time and no order-crossing occurs.

We also assume that demand takes a discrete value, but that continuous distributions can be discretized with an arbitrary precision.

2.3. An Example of Non-quasiconvex Loss Function

In a recent article, Rosling [11] considered loss functions with nonlinear shortage costs. His basic model is a periodic-review model that is almost identical to ours; in his model, two nonlinear shortage (or backorder) costs are added to the traditional linear shortage cost. Here, we only give their mathematical expressions and refer the reader to Rosling [11] for a further explanation.

Suppose that the inventory position after possible ordering is y; then

where

, which is the shortage cost that accumulates at a constant rate proportional to the number of backorders on the books;

, which is the cost that occurs only for the incremental backorders that arise in a period (Hadley and Whitin [8], Chen and Zheng [3], and Bashyam and Fu [2]); B_b(y) = 1 − Φ^L+1(ξ), which is the cost that is incurred whenever there are any backorders in a period (Silver, Pyke, and Peterson [15]); and H(y) is the expected holding cost (per period). (For brevity, we simply denote G(y|L) by G(y).) Setting π = b = 0, we have a loss function with a traditional linear backorder cost. Generally, b and π are viewed as Lagrangian multipliers, or shadow costs, of the ready and fill rates, respectively.

Assuming that the holding cost incurs at a constant rate and that the lead time for replenishment orders is constant, Rosling [11] showed that when p,π,b ≥ 0, G(y) is quasiconvex if Φ(ξ) is log-concave. This condition is met by many demand distributions, including the Normal, Poisson, Negative Binomial, gamma family, Weibull, and beta with parameters (r,q) such that both r and q ≥ 1. However, there are significant exceptions. Take an example that is given in Eppen and Martin [5]: The demand of each time period is normally distributed with mean 20 units and variance 10 (units)². The lead time (in periods) can take on the values 2 and 4 with probabilities of 0.8 and 0.2, respectively. In such a case, the (expected) loss function fails to be quasiconvex when the backlogging cost is nonlinear. Figure 1 illustrates a case in which p = $2/unit/period, π = 0, b = $200/period, and the holding cost rate is $1/unit/period. (The line that is associated with C_min will be explained later.)

Although the loss function can be non-quasiconvex in the inventory position, G(y) generally increases monotonically after y becomes sufficiently large, and it is nonincreasing when y < 0. For example, a test can be carried to see whether G(y) will become monotonic. Suppose that H(y) is convex. If H(y₁ + 1) − H(y₁) ≥ pB_p(y₁) + πB_π(y₁) + bB_b(y₁) (see (1)), then G(y) is increasing for all y > y₁.

3.1. Auxiliary Function Method and Preliminary Results

To elaborate our computational procedure, we introduce a dummy cost γ, which is assumed to be in the range γ ∈ A = {γ : G(y₀) ≤ γ ≤ K + G(y₀)}. As will be shown, the assumption about γ can be trivially met (see Section 3.4). We then define an auxiliary function that relates to this γ as follows:

In fractional programming, the problem of minimizing l_γ(s,S) over s < S is called a parametric problem. The use of such an approach in inventory theory was first proposed by Feng and Sun [6].

One might think of the dummy cost γ as a budget level per period. Because

is the expected total cost per ordering cycle, the difference

gives the cumulative deficit in operating the system per renewal cycle; that is, if the budgeted cost exceeds the expected operating cost, then there will be a surplus, and if the expected operating cost exceeds the budgeted cost, then there will be a fiscal shortfall. The objective is therefore converted to one that minimizes l_γ(s,S).

Our interest in the auxiliary function is due to the fundamental properties of fractional programming regarding the value γ* and the function w_γ, which are defined respectively as γ* = {γ ∈ A : min_s≤S l_γ(s,S) = 0} and w_γ = min_s<S l_γ(s,S). In what follows, we denote by s_γ and S_γ the pair that minimize l_γ(s,S) under a given γ. We also term (s_γ,S_γ) optimal with respect to the auxiliary function; that is,

Lemma 1:

In view of Lemma 1, the minimization of C(s,S) is equivalent to finding a γ ≥ G(y₀), such that min_s<S l_γ(s,S) = 0. Specifically, we begin with an initial dummy cost γ₀ ≥ C_min and then iteratively lower it until the condition min_s<S l_γ(s,S) = 0 is satisfied.

For later use, we rearrange the auxiliary function

As y₀ is a global minimizer of G(y) given a γ ≥ G(y₀), it follows that

Let P(γ) = {y : G(y + 1) ≤ γ < G(y)}. In Figure 2 (where the curve is smoothed and γ = γ₀), there are two elements in P(γ₀): s_γ0¹ and s_γ0². Lemma 2 shows that P(γ) is the set that contains all potential candidates for s_γ. If we denote these candidates as s_γ¹,s_γ²,…, in ascending order of value, then P(γ) = {s_γ¹,s_γ²,…}.

An illustration of candidates for sγ0 and ranges for Sγ0.

Lemma 2:

Proof: Part (a) follows from the fact that

For part (b), we assume to the contrary that G(S_γ) > γ. Then, following (3),

which implies that

If

, then we have 0 > 0, which is an impossible event. Thus,

, leading to

which is a contradiction of (4). Thus, the lemma follows. █

Lemma 2 implies that s_γ ∈ P(γ). Because of Assumption 1 and the assumption that G(y₀) ≤ γ ≤ K + G(y₀), the set P(γ) is clearly nonempty and finite. It can easily be seen that the lower bound of s* is s_γ0¹. Furthermore, not every value of y ≤ S^γ can be an effective order-up-to level. For example, in Figure 2, it can be seen that for an (s_γⁱ,S) to be optimal, y₀ ≤ S ≤ d₁ or S > s_γ², where d₁ = min{y > s_γ¹ : G(y) ≤ γ} and i = 1 or 2. There are more unattractive candidates for S_γ that can be eliminated. However, for notational simplicity, we will compromise some efficiency by ignoring these properties in the sequel algorithm design.

Lemma 2 represents the major extension of the work of FX. In the rest of the section, we discuss a method for the upgrading of order-up-to levels and when to end the search. Interestingly, the properties that are revealed in these discussions are direct extensions of those of FX and in fact these discussions do not engage the loss function.

We close this subsection with a brief comparison of the average-cost approach of ZF and our auxiliary-function approach. By Eq. (6) in ZF and by the proof of Lemma 1 in Zheng [19], one can derive a lower bound for s* and an upper bound for S*, respectively. If their algorithm starts with an initial policy that gives γ₀ and ends with the lowest average cost C_min being identified, then these two bounds coincide with our s_γ0¹ and S^C_min. However, a full enumeration within these bounds is required in their approach, whereas our approach allows a more efficient search. The difference arises from the fact that they directly use the cost function C(s,S) to design their algorithm, whereas we use the auxiliary function l_γ(s,S). Although both functions vary with the policy (s,S), l_γ(s,S) has a simpler form.

3.2. Updating Order-up-to Levels

In this subsection, we discuss a method for updating the optimal order-up-to level in conjunction with the upgrading of the dummy cost. We use Figure 2 to demonstrate how the dummy cost will be upgraded in our algorithm. Under γ = γ₀, suppose that for all S < S₁, l_γ(s_γⁱ,S) ≥ 0, i = 1 or 2, but l_γ(s_γ¹,S₁) < 0 (or l_γ(s_γ²,S₁) < 0). We then immediately upgrade the dummy cost from γ to γ₁ = C(s_γ¹,S₁) (or C(s_γ²,S₁)). In the sequel, when we say that the algorithm updates the dummy cost γ at pair (s₁,S₁), we mean that under γ, for all the potential (s,S) policies with s < S < S₁, l_γ(s,S) ≥ 0, but by raising S to S₁, a pair (s₁,S₁) emerges such that l_γ(s₁,S₁) < 0, where s₁ < S₁,s₁ ∈ P(γ) and S₁ ≤ S^γ. This is how and when the dummy cost γ is upgraded in the proposed algorithm: As soon as l_γ(s,S) turns negative, the dummy cost is lowered from γ to γ₁. The following discussion is based on such an upgrading scheme.

Now suppose that the dummy cost γ is upgraded to γ₁ by a given policy (s₁,S₁), where s₁ ∈ P(γ), and a new set P(γ₁) is formed as s_γ1¹ < s_γ1² < ···. Do we need to evaluate all (s_γ1ⁱ,S₁ − j) policies, where j = 1,2,…? The answer is no, which is deduced from the following lemma.

Lemma 3: Suppose that γ is upgraded to γ₁ < γ by a policy (s₁,S₁), where s₁ ∈ P(γ), S₁ < S^γ, and

Then for any s_γ1ⁱ ∈ P(γ₁) and y < S₁, we have

Proof: Apparently, given s₁ < S₁, l_γ(s₁,S₁), as a function of γ > 0, strictly decreases in γ. As a result, for s_γ1ⁱ ∈ P(γ₁) and y < S₁,

The last inequality results from the updating scheme of the dummy cost that was described earlier. █

3.3. Identifying C_min and an Optimal (s,S) Policy

Once γ is upgraded to γ* = C_min, the algorithm needs to be terminated and an optimal (s,S) policy identified. This subsection discusses when the algorithm should terminate.

Lemma 4: Suppose that the algorithm upgrades the dummy cost γ at (s₁,S₁) to γ₁ = C(s₁,S₁) where s₁ ∈ P(γ). In addition, suppose that for all s ∈ P(γ₁) and S₁ ≤ y ≤ S^γ₁,

Then, γ* = C_min = γ₁, s* = s_γ* = s₁, and S* = S_γ* = S₁.

Proof: In view of Lemma 3, for a case in which s ∈ P(γ₁) and y < S₁, l_γ1(s,y) ≥ 0. With these assumptions, we have

and, thus, γ₁ ≤ C_min. However, as γ₁ = C(s₁,S₁) ≥ C_min = γ*, we have γ₁ = γ*. Furthermore, as l_γ*(s₁,S₁) = 0, we can conclude that s* = s_γ* = s₁ and S* = S_γ* = S₁. █

This lemma indicates that in the final γ upgrading attempt, the search will continue until S^γ*, which is the largest order-up-to level to be evaluated, is reached.

3.4. The Algorithm

We now present the descent-search algorithm, which starts with an initial dummy cost γ₀ = C(s,S) ≤ K + G(y₀). This kind of dummy cost can be easily obtained. For instance, when s = y₀ − 1 and S = y₀,

can be selected as γ₀. It is generally easy to locate s_γ0¹; for example, we can start with y = 0. If G(0) > γ₀, then s_γ0¹ ≥ 0, and if G(0) ≤ γ₀, then s_γ0¹ < 0. This is because G(y) is decreasing when y < 0. We assume that a number y is also available so that for y ≥ y (see the remark at the end of Section 2.3), G(y) is increasing. In finding y, a global minimizer of G(y), y₀, also becomes known. In the algorithm, we do not describe how S^γ is computed, because such a procedure is straightforward. The algorithm will be called the descent-search algorithm, because the dummy cost is lowered after each iteration.

Descent-search algorithm

The algorithm works as follows. It starts with γ₀. The set of potential reorder levels associated with the dummy cost γ₀ consists of s_γ0ⁱ, i = 1,2,…,τ. The set I contains the levels in P(γ₀) that are currently used for evaluation. At the very beginning, the set I contains just one reorder level s = s_γ0¹, and the algorithm begins with the computation of l_γ0(s_γ0¹,s_γ0¹ + 1).

Step 1.1 evaluates all l_γ1(s,S), with S being fixed but s being taken from I. The minimum of l_γ1(s,S) and l_γ1(s′,S) is then obtained. Once this minimum becomes negative, the algorithm immediately upgrades γ₁ in Step 2 by the (s′,S) policy.

If the minimum of l_γ1() is nonnegative, then the algorithm proceeds to Step 1.2. The current order-up-to level S is replaced by S + 1; meanwhile, a test is carried out to see whether S > S^γ. If this is the case, then the algorithm jumps to Step 3 in terms of Lemma 2b. Otherwise, another test is undertaken to see whether G(S − 1) > γ ≥ G(S) and S − 1 can be added to I (by Lemma 2a). Afterward, the algorithm goes back to Step 1.1.

In Step 2.1, the dummy cost γ₁ is upgraded to a “new” γ₁ by (s′,S) according to Lemma 1b. The “default” optimal policy (s*,S*) and dummy cost γ* are also updated at this step. In view of Lemma 3, the search for a “new” order-up-to level begins with S* in Step 2.2. With the new dummy cost, the algorithm regenerates the first candidate for the optimal reorder level and adds it into set I.

Finally, the computation comes to an end in Step 3, in which the optimal (s,S) policy and C_min are identified. The reason why it is possible to have C_min in Step 3 is given by Lemma 4.

We use Figures 2 and 3 to illustrate how the algorithm is executed. At γ₀, the algorithm first searches for the smallest candidate for s_γ0 (i.e., s_γ0¹). It then sets S₀ = s_γ0¹ + 1 to start the process that will upgrade γ₀. As S₀ is increased to S₀′, s_γ0²(= S₀′ − 1) is encountered, and both (s_γ0¹,S₀′) and (s_γ0²,S₀′) are evaluated. As S is further increased to S₁, l_γ0(s_γ0ⁱ,S₁) turns to a negative value (i = 1 or 2) and γ₀ is lowered to γ₁. Under γ₁, two candidates of s_γ1 are identified: s_γ1¹ and s_γ1². The next around of dummy-cost upgrading then begins with the evaluation of both (s_γ1ⁱ,S₁), i = 1 and 2. (For the reason for skipping the order-up-to level below S₁, refer to the beginning of Section 3.2.) This process continues until γ₁ is upgraded, and eventually the optimal policy (s*,S*) is confirmed after S^γ* has been reached.

A typical execution of the algorithm.

It is well recognized that an evaluation of C(s_γ1¹,S), and hence l_γ1(s_γ1¹,S), is a major operation. We thus call each computation of l_γ1(s_γ1¹,S) a critical evaluation. Each time S advances by a unit, l_γ1(s_γ1¹,S) has to be “completely” recalculated, with virtually nothing from the previous evaluation of l_γ1(s_γ1¹,S − j), j ≥ 1, being usable. Therefore, the complexity of the algorithm is mainly determined by the number of critical evaluations (NCE) of l_γ1(s_γ1¹,S) pertaining to different S. Because the NCE is determined by the cost structure and initial dummy cost γ₀, the complexity of the algorithm is dependent on the holding-shortage cost and the initial dummy cost. Given an initial dummy cost of γ₀, the NCE can be greater than (S^γ* − s_γ0¹), because if the dummy cost is upgraded at (s′,S) from γ to γ₁, then S will probably be evaluated twice. Nevertheless, the NCE is bounded by 2(S^γ* − s_γ0¹).

Although the main computational effort relates to the NCE, after each updating of γ, a set of new candidates for the optimal reorder level must be relocated. For this reason, we introduce another criterion: the number of dummy profit cost upgrades (NDC), which measures the number of times that γ is upgraded. Both the NCE and NDC will be obtained later in numerical examples.

We conclude this section with a brief discussion of the complexity of the algorithm. Note that the algorithm ends in a finite number of iterations, because each γ upgrading will push s and/or S forward by at least one. Also note that when the loss function is quasiconvex, the algorithm essentially collapses to that proposed by FX. Note that FX's algorithm has the same complexity as ZF's—about twice the effort that is required to evaluate a specific (s,S) policy. When the loss function is non-quasiconvex, additional effort is required to evaluate perhaps several (s_γⁱ,S) policies simultaneously for each S. However, according to ZF, in calculating C(s_γ¹,S), the C(s_γⁱ,S) are readily available, where i = 2,… and s_γⁱ > s_γ¹. Thus, our algorithm has roughly the same complexity as the two other algorithms (see Chen and Feng [4] for more details).