Abstract

The stochastic uncapacitated lot-sizing problems with incremental quantity discount have been studied in this paper. First, a multistage stochastic mixed integer model is established by the scenario analysis approach and an equivalent reformulation is obtained through proper relaxation under the decreasing unit order price assumption. The proposed reformulation allows us to extend the production-path property to this framework, and furthermore we provide a more accurate characterization of the optimal solution. Then, a backward dynamic programming algorithm is developed to obtain the optimal solution and considering its exponential computation complexity in term of time stages, we design a new rolling horizon heuristic based on the proposed property. Comparisons with the commercial solver CPLEX and other heuristics indicate better performance of our proposed algorithms in both quality of solution and run time.

1. Introduction

The lot-sizing problems have been the subject of intensive research in the last few decades. The basic definition of single-item lot-sizing problems can be stated as follows: the order (production), inventory and backlog quantities in each period should be determined to meet the deterministic but dynamic demand over a finite time horizon. The objective is to minimize the total costs, which consist of the fixed setup cost, order cost and inventory cost. Different quantity discount policies, such as all-units quantity discount and incremental quantity discount, have been widely executed in practice and thus have also been introduced into the lot-sizing problems.

Although the deterministic planning and scheduling models have been intensively studied, in practice there are many different sources of uncertainties, such as customer demand, production lead-time, and product price, which make information that will be needed in subsequent decision stages unavailable to the decision maker. In such cases, the solution provided by a deterministic model may be of little value in terms of applicability of the model’s recommendations, see Beraldi et al. [1]. Thus the stochastic version of lot-sizing problems has been studied recently and with the advent of stochastic programming, the classical deterministic lot-sizing models have been extended to scenario-based multistage stochastic mixed integer programming.

Wagner and Whitin [2] first introduced the single-item dynamic lot-sizing problems without backlogging. They proposed a dynamic programming approach based on the Wagner-Whitin (W-W) property, that is, no production is undertaken if inventory is available. Because W-W property holds for the nondiscount problem, many heuristics have been developed based on this property; however, in the all-units quantity discount case, the property does not necessarily hold due to its discontinuous piecewise linear cost structure by Hu and Munson [3]. For the modified all-units discount problem, Chan et al. [4] demonstrated that a zero-inventory-ordering policy based on the W-W property exists, whose cost is no more than times the optimal cost. Federgruen and Lee [5] characterized structural properties of optimal solutions for both all-units and incremental quantity discount, and proposed dynamic programming algorithms with complexity and , respectively, and with being the number of periods in the planning horizon. For deterministic capacitated lot-sizing problem with general piecewise linear cost function, Shaw and Wagelmans [6] presented a dynamic programming procedure with complexity , where is the number of periods, is the average demand and is the average number of pieces required to represent the production cost function. Note that this pseudopolynomial time algorithm is based on the assumption that the demand is an integral value.

For the stochastic version problem, although Ahmed et al. [7] showed that the W-W property does not hold for the stochastic lot-sizing problems several modified W-W properties have been presented for different versions of the stochastic lot-sizing problems: Guan and Miller [8] studied the stochastic uncapacitated lot-sizing problems with uncertain parameters and introduced the production path property for the optimal solution. Further, the production path property was extended to the stochastic capacitated lot-sizing problems with backlogging in Guan and Miller [8] and the dynamic programming algorithms based on this property have been presented. For the stochastic lot-sizing problems with random lead times, Huang and Küçükyavuz [9] presented the Semi-Wagner-Whitin property under assumption that no order crosses in time.

Besides the algorithms based on the extended W-W properties, Lulli and Sen [10] and Ahmed et al. [7] have presented branch-and-price and branch-and-bound algorithms, respectively to solve the proposed multistage integer programming. Although such branch-and-bound-based (B&B-based) methods have broad application prospects for general integer programming, special structure properties of the stochastic lot-sizing problems have not been explored in order to design customized algorithms, and only computational results for small-size stochastic batch-sizing problems have been reported in Lulli and Sen [10]. Other heuristic methods, such as the fix and relax, have also been redesigned to solve particular stochastic lot-sizing problems, see [1, 11]. We also refer the reader to the recent literature review of [1214].

To the best of our knowledge, little research has been reported on the stochastic lot-sizing problems with incremental quantity discount (SLP-IQD). However as it is reported in the survey by Munson and Rosenblatt [15], of the buyers received quantity discounts for most of the items they purchased and of orders involved either the offer or receipt of some incremental quantity discounts; thus, the study on quantity discount is of great importance in practice. In this paper, a multistage stochastic mixed-integer programming model is established and under decreasing unit order price assumption, an equivalent relaxed formulation is obtained. The reformulation provides the possibility of extending the production path property for optimal solution of SLP-IQD. The extended production path property is not only a direct extension to the case with incremental quantity discount, but also provides a more accurate characterization for the optimal solution. Then, a backward dynamic programming algorithm has been developed. Although it can obtain optimal solutions in polynomial time in terms of the number of nodes, it has an exponential computational complexity in terms of the number of stages. Thus, a new rolling horizon heuristic which makes use of the extended production path property and has flexible parameters settings is presented in order to balance the desired solution quality and run time. Numerical experiments have been implemented to explore the proper parameters settings and validate the effectiveness of the proposed algorithms by comparison with the CPLEX 11.1 solver and other heuristics.

The remainder of the paper is organized as follows. In Section 2, we first introduce the deterministic formulation and then formulate the general multistage stochastic mixed integer model for the stochastic uncapacitated lot-sizing problems with incremental quantity discount (SULP-IQD). An equivalent reformulation is proposed under the decreasing unit order price assumption. In Section 3, the extended production path property is proven and a backward dynamic programming algorithm and a rolling horizon heuristic with flexible parameters settings are developed. Computational results are reported and discussed in Section 4. Section 5 presents conclusions.

2. Mathematical Model

2.1. Deterministic Lot-Sizing Problems with Incremental Quantity Discount

First, we will establish an mathematical model for the deterministic uncapacitated lot-sizing problems with incremental quantity discount (DULS-IQD). Considering a planning horizon of time periods (stages), at each period , the nonnegative demand , variable inventory holding cost , the variable setup cost , and piecewise order cost function at period are given. Variable denotes the inventory quantity at the end of period , and variable and Boolean variable denote the order quantity and fixed charge variable at period , respectively.

The incremental quantity discount cost structure is given as follows: where denotes the number of price intervals. Suppose that and . The decreasing unit order price assumption, that is, with always naturally holds in practice. The piecewise order cost function is given as

Thus, the DULS-IQD can be formulated as subject to

The objective function (2.3) minimizes the sum of inventory cost, setup cost and piecewise order cost. Constraints (2.4) guarantee that the dynamic demands in every period are met. Constraints (2.5) are the capacity constraints, and here we assume that is a sufficiently large upper bound on . Constraints (2.6) ensure that there is no backlogging and order quantity variables are nonnegative. Without loss of generality, we suppose that there is no initial inventory. Setup variable is defined as a binary variable in constraints (2.7).

2.2. Stochastic Lot-Sizing Problems with Incremental Quantity Discount

In this subsection, a stochastic model is established by scenario analysis approach. We assume that the problem parameters follow discrete time stochastic processes with a finite probability space and evolve as a multistage stochastic scenario tree. Let node at stage represent all the realization of the system leading up to and including stage , and set denotes the set of all nodes. Let set denote the set of leaf nodes. Thus a scenario represents a path from the root node (indexed as ) to a leaf node. The related notation is given as follows:: parent of node ,: time stage of node , ,: the set of all descendants of node and ,: the set of nodes on a path from node to node , where , and we assume that includes and ,: children nodes of node , that is, ,: the likelihood (probability) assigned to node , .

Here using the above notation, the deterministic model in the above subsection can be extended to the stochastic environments by replacing the stage index with node index . The multistage stochastic problem formulation, denoted by SP, is given as follows: subject to

An equivalent mixed integer programming formulation can be easily obtained by introducing auxiliary order quantity variables and corresponding Boolean variables. A group of variables are assigned for each order quantity interval as follows:

For model brevity, we introduce the following notations: , , and . Note that the setup cost also increases correspondingly in index , that is, for each node . An equivalent MIP formulation, ESP1, is given as follows: subject to

Proposition 2.1. Formulation ESP1 is equivalent to the original formulation SP.

Proof. Suppose that is an optimal solution of ESP1. For each , there exists at most one such that and by constraints (2.13), (2.14), and the optimality. Note that if , then . Thus, by checking constraints in SP, we construct a corresponding feasible solution of SP as . By the definition of , and , the constructed feasible solution of SP has the same objective cost to .
For each optimal solution of SP, , we define the corresponding solution for ESP1 as follows: for given , if , then and if , then . By checking its feasibility and the objective cost in ESP1, we conclude that the constructed solution is also feasible in ESP1 and has same objective cost to .

Further, by relaxing the constraints on Boolean variables and order quantity variables in (2.13)-(2.14), the following formulation, denoted by ESP2, and proposition can be obtained. subject to

Proposition 2.2. Under the decreasing unit order price assumption, formulation ESP2 is equivalent to ESP1 and SP.

Proof. Because the feasible set of the relaxed formulation ESP2 contains the feasible set of ESP1, the optimal solutions of ESP2 must be proven feasible for ESP1. Let be an optimal solution of ESP2. First note that if , we must have otherwise it contradicts the optimality since .
Second, it is asserted that the relaxed constraints always hold for every optimal solution of ESP2. Suppose that ; thus, and . Without loss of generality, suppose that and , then where in the third inequality.
From the above analysis, a better solution can always be obtained by setting , where , and , which contradicts the optimality.
Third, we prove that the relaxed constraints also hold for the optimal solution. Assume that and for all , but the inequality constraint does not hold. Without loss of generality, suppose that where , then Therefore, we reach a contradiction because we obtain a better solution by setting , and for . Since the optimal solution of ESP2 satisfies all the constraints of ESP1, and ESP2 is obtained from ESP1 by relaxing its constraints, this implies that ESP2 is equivalent to ESP1 and SP.

3. Optimality Condition and Algorithms

In this section, we explore the property of the SULS-IQD and design-efficient algorithms. It is necessary to highlight the differences between the deterministic problems and the stochastic problems. First, Ahmed et al. [7] and Huang and Küçükyavuz [9] have shown that W-W property does not hold for the stochastic version problems. The reason is that one production or positive order, which is made to exactly cover the demand along certain path, will inevitably influence all its descendant nodes. The violation happens when this positive order could not cover some nodes and their inventory level is nonzero. Thus, we develop an extended production-path property for the stochastic problem and design a dynamic programming algorithm. Second, since the number of nodes grows exponentially as the number of time stages increases for stochastic version problems, the traditional dynamic programming methods fail to run efficiently. So, we further design an approximate heuristic based on the proposed property to obtain good enough solutions efficiently in this section.

3.1. Extended Production-Path Property

Formulation ESP2 presented in Section 2 is essential to handle the piecewise linear objective function and extend the production path property to the SULS-IQD. The following proposition is not only a direct extension of Guan and Miller [8] when incremental quantity discount is provided, but also a more accurate characterization of the optimal solution. It shows that there always exists an optimal solution such that if a positive order is made at node , the order quantity exactly covers the demand along the path from node to a certain descendant node ; moreover, no positive order will be made for these nodes along the path from node to the parent node of .

Proposition 3.1 (extended production-path property). For any instance of SULS-IDQ, there exists an optimal solution , such that for each node , if , then there exist and such that(1); (2); (3)

Note that under assumption that all lead time is equal to 1 and by similar arguments, the second optimal condition can be regarded as an extension of the Semi-Wagner-Whitin Property in Huang and Küçükyavuz [9]; however, the third optimal condition provides new restrictions which narrow the scope of the expected optimal solutions. In the next section, an improved backward dynamic programming algorithm and a rolling horizon heuristic based on the extended production-path property will be presented.

For any optimal solution of ESP2, , if , we introduce the following definition: where set contains the first node with positive order quantity after node , and set contains these leaf nodes from which to node no positive order has been placed except node . Note that if a certain positive order quantity is properly transferred between node and all the nodes in set , the feasibility can be kept and only the cost for nodes in the following set is changed:

Proof. First, it is asserted that there exists at least an optimal solution by Weierstrass’ theorem since the feasible set is compact (note that can be bounded by ) and the objective function is continuous. Then for any given optimal solution of ESP2, , if , the first part holds from Proposition 2.2. Next, we show that an optimal solution holding the above property can always be constructed from any given optimal solution of ESP2 by adjusting the variable’s value.
We scan the optimal solution from stage to stage . Assume that there exists a node with not holding the property, but for nodes at stage the conclusion holds, the following approach adjusts variables’ value assigned to nodes in to satisfy the property without changing variables’ value before stage . For any feasible solution of ESP2, we define the object cost function as
The objective cost for is where is the unique positive order quantity at node , and the node set is divided into four subsets: , , , and . For any node , no order is placed.
Since for (by definition) and for (otherwise it contradicts the assumption), we choose small positive scalar satisfying such that the following solutions are feasible for ESP2 and keep the value of other variables unchanged. For given satisfying constraints (3.5)–(3.7), the objective costs of and are Note that the first equality comes from the definition of , , and function. The second one is obtained by comparing with the value of and the last one is obtained by rearranging the terms in the former, where .
Since is an optimal solution of ESP2, thus we must have and , are also optimal solution for ESP2. Now we increase so that at least one of the following cases occur.
Case  1. If equality in (3.5) holds, then we have eliminated the undesired positive order at node and will be used to construct an eventual solution by a similar method.
Case  2. If equality in (3.6) for a certain holds, then there exists such that and . Next, will be used to construct an eventual solution by a similar method.
Case  3. If equality in (3.7) for a certain holds, then we apply the above analysis to . Since there are finite nodes in , eventually case 1 or case 2 will occur (in the worst case after finite steps).
Thus, the optimal solution holding the proposed property can always be constructed after finite steps.

3.2. Dynamic Programming Algorithm

To recursively calculate the optimal solution, the following functions are introduced as in Guan and Miller [8]:: value function at node when its initial inventory is ,  that is,   subject to constraints (2.18) and :: production value function at when its initial inventory is and :: nonproduction value function at when its initial inventory is and .

From Proposition 3.1, if a positive order quantity is made at node when the initial inventory is , then there exists a node , such that . Therefore, the following equations hold:

To obtain the exact optimal solution of SP, it is not necessary to completely characterize the value function . We only need to calculate its values at possible positive discontinuous points. That is, with zero initial inventory assumption for node , we only need to evaluate for since in (3.11). Thus, we give the following dynamic programming for SULS-IQD:

Without loss of generality, we assume for all in the following analysis. Dynamic programming demonstrates a straight method to obtain the optimal solution. Although Guan and Miller [8] showed that general SULS without incremental quantity discount can be solved in , the above algorithm runs in exponential time in terms of since . The exponential computational time encourages us to develop a more effective algorithm for the problems with large numbers of stages (see Algorithm 1).

For each stage to
 For each node at stage
  For each possible initial inventory ,  
  (if , set )
   step 1: Calculate by (3.11)
   step 2: Calculate by (3.12) if ,
    otherwise
   step 3: Calculate by (3.13)
  End For Iteration ( )
 End For Iteration ( )
End For Iteration ( )

3.3. Rolling Horizon Heuristic

In dynamic lot-sizing and scheduling problems with a large planning horizon, rolling horizon heuristics have been developed to decompose the original large-scale problem into a set of smaller subproblems. See, for example, [11, 16]. In contrast with the existing heuristics, the proposed RHH based on the extended production-path property stems from the proposed dynamic programming algorithm. In the heuristic, nonproduction strategy is developed to take advantage of the accurate characterization of the optimal solution, and then look-forward and look-backward heuristic strategies are developed to reduce the computation complexity of the complete enumeration in evaluation of the value function. The computation complexity is also analyzed to demonstrate its advantage for the problems with a large planning horizon.

Nonproduction Strategy
At each iteration for in DP algorithm for given , if , step 1 can be skipped and set since it violates part 3 of Proposition 3.1. In such case, even if , and there exist optimal solutions with , we will not lose the optimality because it is guaranteed that there exists an optimal solution satisfying the extended production path property in Proposition 3.1. The nonproduction strategy can be stated as

Look-Forward Strategy
At step 1 for given node at stage and initial inventory , (3.11) requires calculating all the positive orders that cover demand from node to a certain node . These complete enumeration calculations are very time-consuming, thus we define the following look-forward strategy with parameter FT (forward time stage):
where and denotes the set of node ’s descendants within FT generations. The problem is how to select proper FT. It is obvious that increasing FT will not only improve the quality of solution, but will also increase computation time. Thus, the quality of solution and run time can be balanced by properly selecting FT. The performance of the proposed RHH with different FT settings will be tested by numerical experiments in the next section.

Look-Backward Strategy
At step 2 for given node at stage and initial inventory , if stage is smaller than , it is likely and in such cases step 2 can be skipped. For the same reason, if is larger than a certain value, we can expect and this iteration for can be skipped too. Look-backward strategy with parameter BF (backward time stage) can be stated as follows:
where BF depends on the distribution of demand. We will set proper based on distribution of demand by numerical experiments in the next section.
In the above strategies, the major modification comes from the iteration for and calculation of at step 1 for given node . In order to give a brief characterization of the above heuristic rules and analyze its computational complexity, we give the following lemma.

Lemma 3.2. For each node , the rolling horizon heuristic with and only needs to evaluate the following values: where .

Proof. This lemma is proven by induction from nodes at stage to leaf nodes at stage . Since positive order must be made at node due to the zero initial inventory assumption, thus for node at stage by (3.15) possible initial inventory at node could only be a certain value in set . Suppose the lemma holds for node at stage , where , that is, the set of initial inventory for node is given by .
Now consider the possible initial inventory at node .
Case 1. If a positive order has been made at node , by (3.15) in the look-forward strategy, the set of initial inventory at nodes is given by .
Case 2. if no order has been made at node , the set of possible initial inventory at node is given by . However, note that by look-backward strategy (3.16) we skip the calculation of for those such that ; thus, in this case we only need to consider initial inventory value in set   =  ,   at node .
Combine the above two cases, the conclusion holds for node .

Figure 1 gives an intuitive example of a balanced scenario tree with at each node. For the rolling horizon heuristic with and , given for node , consider the initial inventory set for node . If a positive order has been made at node , will be introduced by (3.15). By the look-backward strategy (3.16), we do not need to calculate , where . Other values will be inherited if no order has been made at node . Note that is no more than where we assume that for all .

By summarizing the above analysis, the rolling horizon heuristic is given as in Algorithm 2.

For each stage to
 For each node at stage
  For each initial inventory
  Set ; (if ,  set )
   step 1: If , go  to step 3; otherwise calculate
    by (3.15)
   step 2: If or , goto step 3; otherwise
    calculate by (3.12)
   step 3: Calculate by (3.13)
  End For Iteration ( )
 End For Iteration ( )
End For Iteration ( )

Next the computation complexity of RHH is analyzed. For given node at stage , there exist at most possible initial inventories . For each given , it takes time to complete calculation of at step 1. Step 2 and step 3 will be completed in time. Thus the total run time can be estimated as (assume that for all )

Proposition 3.3. For any instance of SULS-IDQ, the rolling horizon heuristic with parameter and runs in time.

The above analysis can be applied to the dynamic programming algorithm, thus the total run time for DP is given by By comparison of computational complexity in (3.18) and (3.19), we observe that RHH works more efficiently for large . Effective parameters settings of RHH will be further explored by numerical experiments in the next section.

4. Computational Experiment

In this section, the computational results on both DP and RHH are reported. In the computational analysis, we first concentrate on identifying proper settings of parameters FT and BT for RHH by comparison of its relative error gap and run time with DP’s. Then, DP and RHH are compared with the CPLEX solver and other heuristics for the lot-sizing problems with fixed charge.

4.1. Problem Instance

In order to evaluate performance of the proposed DP and RHH, and explore the proper parameters settings of and for RHH, 8 groups of problem instances are generated by varying the number of stages from 5 to 12 and the number of branches from 2 to 4. The number of incremental quantity discount intervals is fixed to 3. Table 1 gives the number of branches , the number of stages , the number of nodes , and the number of integer and continuous variables in formulation ESP2 for each group of problem instances. Other parameters are generated randomly and we assume that they are a sequence of i.i.d. random variables which are subject to truncated normal distribution or uniform distribution. We report in Table 2 part of the problem parameters and the variance of will be further defined by (coefficient of variation) in the following subsection. Due to decreasing unit order price assumption, the unit order price and setup cost at node is generated by and , where is the discount factor and we set . Our experiments are conducted on a workstation clocked at 2.33 GHz and equipped with 11.9 GB of core memory.

4.2. Numerical Results

In order to evaluate the performance of the proposed RHH method with different parameters, we define two different implementations of RHH and the problem parameter varies from to for each group of problems. with parameters and has been designed to obtain good enough solution in a short time, while with parameters and has been designed to obtain better solutions, where denotes the smallest integer larger than . Table 3 gives the optimal objective function values and CPU time in seconds for DP, and the percentage relative error () and CPU time for and , where For , the computational results show that is able to solve all the instances in the shortest time with no more than relative error except, for instance, . The solution quality and run time are further improved when the number of stages becomes larger and the worst results of can be interpreted as the look-forward time horizon being too small to obtain near optimal solution. For , the relative error has been further decreased compared with . Although the run time has been increased, still has advantage in run time compared with especially for large problem instances. We also note that as the coefficient of variation increases, the solutions quality and run time have been constantly improved for all problem instances.

Then in the second experiment, we concentrate on the comparison of solution quality and computation time with standard MIP solver CPLEX (version 11.1) and another heuristic dynamic slope scaling procedure (DSSP). DSSP proposed by Kim and Pardalos [17] is a newly developed effective heuristic algorithm that provides good-quality solution to the concave piecewise linear optimization problem and among the heuristics it works well in practice [18] The initial solutions and updating schemes are implemented in accordance with the recommendations in [18] and the stopping criterion is given as follows: if or the iteration reaches , then it stops, where . To compare the solution quality for given time limit, we set the time limit of CPLEX solver to time of the total run time of DP. We report in Table 4 the percentage relative error and CPU time of RHH (its parameters are set according to ), DSSP and CPLEX and the optimal objective function value and CPU time of DP. Table 4 shows that the proposed RHH outperforms the DSSP in both the quality of solution and run time for almost all the test instances, and the CPLEX solver fails to find optimal solution within the given time.

In summary, the proposed DP can solve the SULP-IQD efficiently compared with the standard CPLEX solver and by properly setting the parameters, we obtain effective RHH which outperforms the DSSP heuristic for the tested instances. The computational results also show that RHH performs better for problem instances with a larger number of stages and high coefficient of variation.

5. Conclusion

In this paper, we study the stochastic uncapacitated lot-sizing problems with incremental quantity discount where the uncertain parameters are supposed to evolve as discrete time stochastic processes. First, we establish the original stochastic formulation by scenario analysis approach. Another two formulations are built by introducing auxiliary variables and relaxing certain constraints. Then, it is proven that under the decreasing unit order price assumption, the relaxed formulation is equivalent to the original one. Based on this reformulation, the extended production-path property is presented for the SULP-IQD and it enhances the ability to further refine the desired optimal solution by providing a more accurate characterization.

To obtain the exact optimal solution, a dynamic programming algorithm is developed. Although the dynamic programming algorithm has the polynomial-time computational complexity in terms of the number of nodes, it runs exponentially in terms of the number of stages. Thus, a new rolling horizon heuristic is further designed which contains three types of strategies to reduce the computational time. The nonproduction strategy is designed based on the accurate characterization of the optimal solution, and the look-forward and look-backward strategy is developed to overcome the complete enumeration calculations in the production and nonproduction value function. Numerical experiments are carried out to identify proper parameters settings of the proposed RHH and to evaluate the performance of the proposed algorithms by comparison with the CPLEX solver and DSSP heuristic. The computational results of a large group of problem instances with different parameters setting suggest that DP outperforms the CPLEX solver in run time required for obtaining optimal solution and the proposed RHH demonstrates satisfactory run time and solution quality compared with DSSP heuristic; moreover, as the computational complexity analysis suggests, the performance of RHH is better for problems with a greater number of stages.

Since the main difficulties for the stochastic lot-sizing problems stem from the setup cost and uncertain parameters, it will be an area of future research to analyze the properties of the problem and present effective algorithms for the stochastic lot-sizing problems with complex setup requirements, such as setup carryovers by Buschkühl et al. [14] and sequence-dependent setup costs by Beraldi et al. [1].

Acknowledgments

The paper is supported by the National Nature Science Foundation of China (NSFC no. 60874071 and 60834004) and Research Fund for the Doctoral Program of Higher Education of China (RFDP no. 20090002110035).