# Dynamic Programming and Heuristic for Stochastic Uncapacitated Lot-Sizing Problems with Incremental Quantity Discount

Authors: Yuli Zhang; Shiji Song; Cheng Wu; Wenjun Yin
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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.


---

## Body

## 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 4/3 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 O(T3) and O(T2), respectively, and with T 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 O(T2q¯d¯), where T is the number of periods, d¯ is the average demand and q¯ 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 [12–14].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], 83% of the buyers received quantity discounts for most of the items they purchased and 37% 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 Section2, 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 ofT time periods (stages), at each period t, the nonnegative demand dt, variable inventory holding cost ht, the variable setup cost ct, and piecewise order cost function ft(·) at period t∈{1,…,T} are given. Variable st denotes the inventory quantity at the end of period t, and variable xt and Boolean variable yt denote the order quantity and fixed charge variable at period t, respectively.The incremental quantity discount cost structure is given as follows:(2.1)theunitordercost={pt,1forthefirstQt,1units,pt,2forthenextQt,2-Qt,1units,⋮⋮pt,KtforeachunitinexcessofQt,Kt-1units,
where Kt denotes the number of price intervals. Suppose that Qt,0=0 and Qt,Kt=+∞. The decreasing unit order price assumption, that is, pt,k>pt,k+1,k=1,…,Kt-1 with 0=Qt,0<Qt,1<⋯<Qt,k<Qt,k+1<⋯<Qt,Kt=+∞ always naturally holds in practice. The piecewise order cost function is given as
(2.2)ft(xt)={pt,1xt,ifQt,0≤xt<Qt,1,∑j=1k-1{(pt,j-pt,j+1)Qt,j}+pt,kxt,ifQt,k-1≤xt<Qt,k,k=2,…,Kt.Thus, the DULS-IQD can be formulated as(2.3)min:∑t=1T‍{ft(xt)+htst+ctyt}
subject to
(2.4)st-1+xt-dt=st,t=1,…,T,(2.5)xt≤Myt,t=1,…,T,(2.6)xt≥0,st≥0,s0=0t=1,…,T,(2.7)yt∈{0,1},t=1,…,T.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 M is a sufficiently large upper bound on xt. 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 yt 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 noden at stage t(t=1,…,T) represent all the realization of the system leading up to and including stage t, and set N denotes the set of all nodes. Let set L denote the set of leaf nodes. Thus a scenario represents a path from the root node (indexed as n=0) to a leaf node. The related notation is given as follows:a(n): parent of node n,t(n): time stage of node n, t(n)∈{1,…,T},D(n): the set of all descendants of node n and n∈D(n),P(n,m): the set of nodes on a path from node n to node m, where m∈D(n), and we assume that P(n,m) includes n and m,C(n): children nodes of node n, that is, C(n)={m∈N:a(m)=n},θn: the likelihood (probability) assigned to node n, θn>0.Here using the above notation, the deterministic model in the above subsection can be extended to the stochastic environments by replacing the stage indext∈{1,…,T} with node index n∈N. The multistage stochastic problem formulation, denoted by SP, is given as follows:
(2.8)(SP)min:∑n∈N‍{θn{fn(xn)+hnsn+cnyn}}
subject to
(2.9)sa(n)+xn-dn=sn,n∈N,xn≤Myn,n∈N,xn≥0,sn≥0,sa(0)=0,n∈N,yn∈{0,1},n∈N.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:(2.10)xn,k≥0,n∈N,k=1,…,Kn,yn,k={0,xn,k=0,1,Qn,k-1<xn,k≤Qn,k.For model brevity, we introduce the following notations:p¯n,k=Δθnpn,k, h¯n=Δθnhn, and c¯n,k=Δθn{cn+∑j=1k-1‍{(pn,j-pn,j+1)Qn,j}}. Note that the setup cost also increases correspondingly in index k, that is, c¯n,k<c¯n,k+1 for each node n∈N. An equivalent MIP formulation, ESP1, is given as follows:
(2.11)(ESP1)min:∑n∈N‍{∑k=1Kn‍{p¯n,kxn,k+c¯n,kyn,k}+h¯nsn}
subject to
(2.12)sa(n)+∑k=1Kn‍xn,k-dn=sn,n∈N,(2.13)Qn,k-1yn,k≤xn,k≤Qn,kyn,k,n∈N,k=1,…,Kn,(2.14)∑k=1Kn‍yn,k≤1,n∈N,(2.15)xn,k≥0,sn≥0,sa(0)=0,n∈N,k=1,…,Kn,(2.16)yn,k∈{0,1},n∈N,k=1,…,Kn.Proposition 2.1.
Formulation ESP1 is equivalent to the original formulation SP.Proof.
Suppose that{xn,k*,yn,k*,sn*∣n∈N,k=1,…,Kn} is an optimal solution of ESP1. For each n∈N, there exists at most one kn∈{1,…,Kn} such that yn,kn*=1 and M≥xn,kn*>0 by constraints (2.13), (2.14), and the optimality. Note that if ∑k′=1Kn‍xn,k′*=xn,kn*=0, then ∑k′=1Kn‍yn,k′*=yn,kn*=0. Thus, by checking constraints in SP, we construct a corresponding feasible solution of SP as {x¯n*=xn,kn*,y¯n*=yn,kn*,s¯n*=sn*∣n∈N}. By the definition of p¯n,k, c¯n,k and h¯n, the constructed feasible solution of SP has the same objective cost to {xn,k*,yn,k*,sn*∣n∈N,k=1,…,Kn}.
For each optimal solution of SP,{xn*,yn*,sn*∣n∈N}, we define the corresponding solution for ESP1 as follows: for given n∈N, if Qn,k-1<xn*≤Qn,k, then {x¯n,k*=xn*,y¯n,k*=1,s¯n*=sn*,x¯n,k′*=y¯n,k′*=0,k′≠k} and if xn*=0, then {x¯n,k*=0,y¯n,k*=0,s¯n*=sn*,k∈{1,…,Kn}}. 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 {xn*,yn*,sn*∣n∈N}.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.
(2.17)(ESP2)min:∑n∈N‍{∑k=1Kn‍{p¯n,kxn,k+c¯n,kyn,k}+h¯nsn}
subject to
(2.18)sa(n)+∑k=1Kn‍xn,k-dn=sn,n∈N,xn,k≤yn,kM,n∈N,k=1,…,Kn,xn,k≥0,sn≥0,sa(0)=0n∈N,k=1,…,Kn,yn,k∈{0,1},n∈N,k=1,…,Kn.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{x~n,k*,y~n,k*,s~n*∣n∈N,k=1,…,Kn} be an optimal solution of ESP2. First note that if y~n,k*=1, we must have x~n,k*>0 otherwise it contradicts the optimality since c¯n,k>0.
Second, it is asserted that the relaxed constraints∑k=1Kn‍yn,k≤1,n∈N always hold for every optimal solution of ESP2. Suppose that y~n,l*=y~n,m*=1,n∈N,1≤l≠m≤Kn; thus, x~n,l*>0 and x~n,m*>0. Without loss of generality, suppose that m>l and Qn,r-1≤x~n,l*+x~n,m*<Qn,r,1≤r≤Kn, then
(2.19)c¯n,r+p¯n,r(x~n,l*+x~n,m*)=θn{cn+∑j=1r-1‍{(pn,j-pn,j+1)Qn,j}+pn,r(x~n,l*+x~n,m*)}=θn{∑j=mr-1‍{(pn,j-pn,j+1)Qn,j}+pn,r(x~n,l*+x~n,m*)-pn,m(x~n,l*+x~n,m*)}+c¯n,m+p¯n,m(x~n,l*+x~n,m*)≤θn{∑j=mr-1‍{(pn,j-pn,j+1)(x~n,l*+x~n,m*)}+pn,r(x~n,l*+x~n,m*)-pn,m(x~n,l*+x~n,m*)}+c¯n,m+p¯n,m(x~n,l*+x~n,m*)=c¯n,m+p¯n,m(x~n,l*+x~n,m*)<c¯n,m+p¯n,mx~n,m*+p¯n,lx~n,l*<c¯n,m+p¯n,mx~n,m*+c¯n,l+p¯n,lx~n,l*,
where ∑j=mr-1‍{(pn,j-pn,j+1)=pn,m-pn,r in the third inequality.
From the above analysis, a better solution can always be obtained by settingx^n,r*=∑j=1Kn‍x~n,j*,y^n,r*=1, where Qn,r-1≤∑j=1Kn‍x~n,j*≤Qn,r, and x^n,k*=0,y^n,k*=0,k≠r, which contradicts the optimality.
Third, we prove that the relaxed constraintsQn,k-1yn,k<xn,k≤Qn,kyn,k,n∈N,k=1,…,Kn also hold for the optimal solution. Assume that x~n,l*>0 and x~n,m*=0 for all m≠l, but the inequality constraint Qn,l-1≤x~n,l*≤Qn,l does not hold. Without loss of generality, suppose that Qn,r-1<x~n,l*≤Qn,r≤Qn,l-1 where r≤l-1, then
(2.20)c¯n,l+p¯n,lx~n,l*=θn{cn+∑j=1l-1‍{(pn,j-pn,j+1)Qn,j}+pn,lx~n,l*}=c¯n,r+p¯n,rx~n,l*+θn{∑j=rl-1‍{(pn,j-pn,j+1)Qn,j}+(pn,l-pn,r)x~n,l*}>c¯n,r+p¯n,rx~n,l*.
Therefore, we reach a contradiction because we obtain a better solution by setting x^n,r*=x~n,l*,y^n,r*=1, and x^n,k*=0,y^n,k*=0 for k≠r. 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.


## 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 ofT time periods (stages), at each period t, the nonnegative demand dt, variable inventory holding cost ht, the variable setup cost ct, and piecewise order cost function ft(·) at period t∈{1,…,T} are given. Variable st denotes the inventory quantity at the end of period t, and variable xt and Boolean variable yt denote the order quantity and fixed charge variable at period t, respectively.The incremental quantity discount cost structure is given as follows:(2.1)theunitordercost={pt,1forthefirstQt,1units,pt,2forthenextQt,2-Qt,1units,⋮⋮pt,KtforeachunitinexcessofQt,Kt-1units,
where Kt denotes the number of price intervals. Suppose that Qt,0=0 and Qt,Kt=+∞. The decreasing unit order price assumption, that is, pt,k>pt,k+1,k=1,…,Kt-1 with 0=Qt,0<Qt,1<⋯<Qt,k<Qt,k+1<⋯<Qt,Kt=+∞ always naturally holds in practice. The piecewise order cost function is given as
(2.2)ft(xt)={pt,1xt,ifQt,0≤xt<Qt,1,∑j=1k-1{(pt,j-pt,j+1)Qt,j}+pt,kxt,ifQt,k-1≤xt<Qt,k,k=2,…,Kt.Thus, the DULS-IQD can be formulated as(2.3)min:∑t=1T‍{ft(xt)+htst+ctyt}
subject to
(2.4)st-1+xt-dt=st,t=1,…,T,(2.5)xt≤Myt,t=1,…,T,(2.6)xt≥0,st≥0,s0=0t=1,…,T,(2.7)yt∈{0,1},t=1,…,T.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 M is a sufficiently large upper bound on xt. 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 yt 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 noden at stage t(t=1,…,T) represent all the realization of the system leading up to and including stage t, and set N denotes the set of all nodes. Let set L denote the set of leaf nodes. Thus a scenario represents a path from the root node (indexed as n=0) to a leaf node. The related notation is given as follows:a(n): parent of node n,t(n): time stage of node n, t(n)∈{1,…,T},D(n): the set of all descendants of node n and n∈D(n),P(n,m): the set of nodes on a path from node n to node m, where m∈D(n), and we assume that P(n,m) includes n and m,C(n): children nodes of node n, that is, C(n)={m∈N:a(m)=n},θn: the likelihood (probability) assigned to node n, θn>0.Here using the above notation, the deterministic model in the above subsection can be extended to the stochastic environments by replacing the stage indext∈{1,…,T} with node index n∈N. The multistage stochastic problem formulation, denoted by SP, is given as follows:
(2.8)(SP)min:∑n∈N‍{θn{fn(xn)+hnsn+cnyn}}
subject to
(2.9)sa(n)+xn-dn=sn,n∈N,xn≤Myn,n∈N,xn≥0,sn≥0,sa(0)=0,n∈N,yn∈{0,1},n∈N.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:(2.10)xn,k≥0,n∈N,k=1,…,Kn,yn,k={0,xn,k=0,1,Qn,k-1<xn,k≤Qn,k.For model brevity, we introduce the following notations:p¯n,k=Δθnpn,k, h¯n=Δθnhn, and c¯n,k=Δθn{cn+∑j=1k-1‍{(pn,j-pn,j+1)Qn,j}}. Note that the setup cost also increases correspondingly in index k, that is, c¯n,k<c¯n,k+1 for each node n∈N. An equivalent MIP formulation, ESP1, is given as follows:
(2.11)(ESP1)min:∑n∈N‍{∑k=1Kn‍{p¯n,kxn,k+c¯n,kyn,k}+h¯nsn}
subject to
(2.12)sa(n)+∑k=1Kn‍xn,k-dn=sn,n∈N,(2.13)Qn,k-1yn,k≤xn,k≤Qn,kyn,k,n∈N,k=1,…,Kn,(2.14)∑k=1Kn‍yn,k≤1,n∈N,(2.15)xn,k≥0,sn≥0,sa(0)=0,n∈N,k=1,…,Kn,(2.16)yn,k∈{0,1},n∈N,k=1,…,Kn.Proposition 2.1.
Formulation ESP1 is equivalent to the original formulation SP.Proof.
Suppose that{xn,k*,yn,k*,sn*∣n∈N,k=1,…,Kn} is an optimal solution of ESP1. For each n∈N, there exists at most one kn∈{1,…,Kn} such that yn,kn*=1 and M≥xn,kn*>0 by constraints (2.13), (2.14), and the optimality. Note that if ∑k′=1Kn‍xn,k′*=xn,kn*=0, then ∑k′=1Kn‍yn,k′*=yn,kn*=0. Thus, by checking constraints in SP, we construct a corresponding feasible solution of SP as {x¯n*=xn,kn*,y¯n*=yn,kn*,s¯n*=sn*∣n∈N}. By the definition of p¯n,k, c¯n,k and h¯n, the constructed feasible solution of SP has the same objective cost to {xn,k*,yn,k*,sn*∣n∈N,k=1,…,Kn}.
For each optimal solution of SP,{xn*,yn*,sn*∣n∈N}, we define the corresponding solution for ESP1 as follows: for given n∈N, if Qn,k-1<xn*≤Qn,k, then {x¯n,k*=xn*,y¯n,k*=1,s¯n*=sn*,x¯n,k′*=y¯n,k′*=0,k′≠k} and if xn*=0, then {x¯n,k*=0,y¯n,k*=0,s¯n*=sn*,k∈{1,…,Kn}}. 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 {xn*,yn*,sn*∣n∈N}.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.
(2.17)(ESP2)min:∑n∈N‍{∑k=1Kn‍{p¯n,kxn,k+c¯n,kyn,k}+h¯nsn}
subject to
(2.18)sa(n)+∑k=1Kn‍xn,k-dn=sn,n∈N,xn,k≤yn,kM,n∈N,k=1,…,Kn,xn,k≥0,sn≥0,sa(0)=0n∈N,k=1,…,Kn,yn,k∈{0,1},n∈N,k=1,…,Kn.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{x~n,k*,y~n,k*,s~n*∣n∈N,k=1,…,Kn} be an optimal solution of ESP2. First note that if y~n,k*=1, we must have x~n,k*>0 otherwise it contradicts the optimality since c¯n,k>0.
Second, it is asserted that the relaxed constraints∑k=1Kn‍yn,k≤1,n∈N always hold for every optimal solution of ESP2. Suppose that y~n,l*=y~n,m*=1,n∈N,1≤l≠m≤Kn; thus, x~n,l*>0 and x~n,m*>0. Without loss of generality, suppose that m>l and Qn,r-1≤x~n,l*+x~n,m*<Qn,r,1≤r≤Kn, then
(2.19)c¯n,r+p¯n,r(x~n,l*+x~n,m*)=θn{cn+∑j=1r-1‍{(pn,j-pn,j+1)Qn,j}+pn,r(x~n,l*+x~n,m*)}=θn{∑j=mr-1‍{(pn,j-pn,j+1)Qn,j}+pn,r(x~n,l*+x~n,m*)-pn,m(x~n,l*+x~n,m*)}+c¯n,m+p¯n,m(x~n,l*+x~n,m*)≤θn{∑j=mr-1‍{(pn,j-pn,j+1)(x~n,l*+x~n,m*)}+pn,r(x~n,l*+x~n,m*)-pn,m(x~n,l*+x~n,m*)}+c¯n,m+p¯n,m(x~n,l*+x~n,m*)=c¯n,m+p¯n,m(x~n,l*+x~n,m*)<c¯n,m+p¯n,mx~n,m*+p¯n,lx~n,l*<c¯n,m+p¯n,mx~n,m*+c¯n,l+p¯n,lx~n,l*,
where ∑j=mr-1‍{(pn,j-pn,j+1)=pn,m-pn,r in the third inequality.
From the above analysis, a better solution can always be obtained by settingx^n,r*=∑j=1Kn‍x~n,j*,y^n,r*=1, where Qn,r-1≤∑j=1Kn‍x~n,j*≤Qn,r, and x^n,k*=0,y^n,k*=0,k≠r, which contradicts the optimality.
Third, we prove that the relaxed constraintsQn,k-1yn,k<xn,k≤Qn,kyn,k,n∈N,k=1,…,Kn also hold for the optimal solution. Assume that x~n,l*>0 and x~n,m*=0 for all m≠l, but the inequality constraint Qn,l-1≤x~n,l*≤Qn,l does not hold. Without loss of generality, suppose that Qn,r-1<x~n,l*≤Qn,r≤Qn,l-1 where r≤l-1, then
(2.20)c¯n,l+p¯n,lx~n,l*=θn{cn+∑j=1l-1‍{(pn,j-pn,j+1)Qn,j}+pn,lx~n,l*}=c¯n,r+p¯n,rx~n,l*+θn{∑j=rl-1‍{(pn,j-pn,j+1)Qn,j}+(pn,l-pn,r)x~n,l*}>c¯n,r+p¯n,rx~n,l*.
Therefore, we reach a contradiction because we obtain a better solution by setting x^n,r*=x~n,l*,y^n,r*=1, and x^n,k*=0,y^n,k*=0 for k≠r. 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 Section2 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 n, the order quantity exactly covers the demand along the path from node n to a certain descendant node m; moreover, no positive order will be made for these nodes along the path from node n to the parent node of m.Proposition 3.1 (extended production-path property).
For any instance of SULS-IDQ, there exists an optimal solution(x*,y*,s*), such that for each node n∈N, if xn,k*>0, then there exist m∈D(n) and Qn,k-1<dn,m-sa(n)*≤Qn,k such that(1)
xn,i*=0,yn,i*=0,foralli=1,…,Kn,i≠k;(2)
xn,k*=dn,m-sa(n)*,yn,k*=1;(3)
xl,i*=0,yl,i*=0,foralll∈P(n,a(m))∖{n},i=1,…,Kl.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,(x,y,s), if xn,k>0, we introduce the following definition:
(3.1)NP(n)={m∣m∈D(n),∑k=1Kl‍xl,k=0,∀l∈P(n,a(m))∖{n},∑k=1Km‍xm,k>0},NNP(n)={m∣m∈D(n)∩L,∑k=1Kl‍xl,k=0,∀l∈P(n,m)∖{n}},
where set NP(n) contains the first node with positive order quantity after node n, and set NNP(n) contains these leaf nodes from which to node n no positive order has been placed except node n. Note that if a certain positive order quantity is properly transferred between node n and all the nodes in set NP(n), the feasibility can be kept and only the cost for nodes in the following set is changed:
(3.2)ϕ(n)={m∣m∈P(n,l),∀l∈NP(n)∪NNP(n)}.Proof.
First, it is asserted that there exists at least an optimal solution by Weierstrass’ theorem since the feasible set is compact (note that∑k=1Kn‍xn,k can be bounded by maxm∈D(n)∪L{dn,m}) and the objective function is continuous. Then for any given optimal solution of ESP2, (x*,y*,s*), if xn,k*>0, 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(x*,y*,s*) from stage 1 to stage T. Assume that there exists a node n with xn,k*>0 not holding the property, but for nodes at stage t<t(n) the conclusion holds, the following approach adjusts variables’ value assigned to nodes in D(n) to satisfy the property without changing variables’ value before stage t(n). For any feasible solution (x,y,s) of ESP2, we define the object cost function as
(3.3)ℱ(x,y,s)=∑n∈N‍{∑k=1Kn‍{p¯n,kxn,k+c¯n,kyn,k}+h¯nsn}.
The objective cost for(x*,y*,s*) is
(3.4)ℱ(x*,y*,s*)=∑m∈N∖ϕ(n)‍{∑j=1Km‍{p¯m,jxm,j*+c¯m,jym,j*}+h¯msm*}+p¯n,kxn,k*+c¯n,kyn,k*+h¯nsn*+∑m∈ϕ(n){{n}∪NP(n)}‍h¯msm*+∑m∈NP(n)‍{p¯m,jmxm,jm*+c¯m,jmym,jm*+h¯msm*},
where jm is the unique positive order quantity at node m,m∈NP(n), and the node set N is divided into four subsets: {n}, NP(n), ϕ(n)∖{{n}∪NP(n)}, and N∖ϕ(n). For any node m∈ϕ(n)∖{{n}∪NP(n)}, no order is placed.
Sincexm,jm*>0 for m∈NP(n) (by definition) and sm*>0 for m∈NNP(n) (otherwise it contradicts the assumption), we choose small positive scalar ɛ satisfying
(3.5)0<ɛ≤xn,k*,(3.6)0<ɛ≤sm*,form∈NNP(n),(3.7)0<ɛ≤xm,jm*,form∈NP(n),
such that the following solutions are feasible for ESP2
(3.8)x¯n,k*=xn,k*-ɛ,s¯m*=sm*-ɛform∈ϕ(n)∖NP(n),x¯m,jm*=xm,jm*+ɛforn∈NP(n),x^n,k*=xn,k*+ɛ,s^m*=sm*+ɛform∈ϕ(n)∖NP(n),(3.9)x^m,jm*=xm,jm*-ɛforn∈NP(n)
and keep the value of other variables unchanged. For given ɛ satisfying constraints (3.5)–(3.7), the objective costs of (x¯*,y¯*,s¯*) and (x^*,y^*,s^*) are
(3.10)ℱ(x¯*,y¯*,s¯*)=∑m∈N∖ϕ(n)‍{∑j=1Km‍{p¯m,jxm,j*+c¯m,jym,j*}+h¯msm*}+p¯n,k(xn,k*-ɛ)+c¯n,kyn,k*+h¯n(sn*-ɛ)+∑m∈ϕ(n)∖{{n}∪NP(n)}‍h¯m(sm*-ɛ)+∑m∈NP(n)‍{p¯m,jm(xm,jm*+ɛ)+c¯m,jmym,jm*+h¯msm*}=ℱ(x*,y*,s*)-p¯n,kɛ-h¯nɛ-∑m∈ϕ(n)∖{{n}∪NP(n)}‍h¯mɛ+∑m∈NP(n)‍p¯m,jmɛ=ℱ(x*,y*,s*)-p¯n,kɛ-∑m∈ϕ(n)∖NP(n)‍h¯mɛ+∑m∈NP(n)‍p¯m,jmɛ,ℱ(x^*,y^*,s^*)=∑m∈N∖ϕ(n)‍{∑j=1Km‍{p¯m,jxm,j*+c¯m,jym,j*}+h¯msm*}+p¯n,k(xn,k*+ɛ)+c¯n,kyn,k*+h¯n(sn*+ɛ)+∑m∈ϕ(n)∖{{n}∪NP(n)}‍h¯m(sm*+ɛ)+∑m∈NP(n)‍{p¯m,jm(xm,jm*-ɛ)+c¯m,jmym,jm*+h¯msm*}=ℱ(x*,y*,s*)+p¯n,kɛ+h¯nɛ+∑m∈ϕ(n)∖{{n}∪NP(n)}‍h¯mɛ-∑m∈NP(n)‍p¯m,jmɛ.=ℱ(x*,y*,s*)+p¯n,kɛ+∑m∈ϕ(n)∖NP(n)‍h¯mɛ-∑m∈NP(n)‍p¯m,jmɛ.
Note that the first equality comes from the definition of (x¯*,y¯*,s¯*), (x^*,y^*,s^*), and ℱ function. The second one is obtained by comparing with the value of ℱ(x*,y*,s*) and the last one is obtained by rearranging the terms in the former, where h¯nɛ+∑m∈ϕ(n)∖{{n}∪NP(n)}‍h¯mɛ=∑m∈ϕ(n)∖NP(n)‍h¯mɛ.
Since(x*,y*,s*) is an optimal solution of ESP2, thus we must have p¯n,k+∑m∈ϕ(n)∖NP(n)‍h¯m-∑m∈NP(n)‍p¯m,jm=0 and (x¯*,y¯*,s¯*), (x^*,y^*,s^*) 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 n and (x¯*,y¯*,s¯*) will be used to construct an eventual solution by a similar method.
Case  2. If equality in (3.6) for a certain m holds, then there exists m∈D(n) such that x¯n,k*=dn,m-sa(n)* and x¯l,k*=0,y¯l,k*=0foralll∈P(n,a(m))∖{n},k=1,…,Kl. Next, (x¯*,y¯*,s¯*) will be used to construct an eventual solution by a similar method.
Case  3. If equality in (3.7) for a certain m holds, then we apply the above analysis to (x^*,y^*,s^*). Since there are finite nodes in D(m), eventually case 1 or case 2 will occur (in the worst case NP(n)=∅ 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]:H(n,s): value function at node n when its initial inventory is sa(n)=s,  that is, H(n,s)=min{∑m∈D(n)‍{∑k=1Km‍{p¯m,kxm,k+c¯m,kym,k}+h¯msm}}  subject to constraints (2.18) and sa(n)=s:HP(n,s): production value function at n when its initial inventory is sa(n)=s and ∑k=1Kn‍xn,k>0:HNP(n,s): nonproduction value function at n when its initial inventory is sa(n)=s and ∑k=1Kn‍xn,k=0.From Proposition3.1, if a positive order quantity is made at node n when the initial inventory is s, then there exists a node j∈D(n), such that xn=dn,j-s. Therefore, the following equations hold:
(3.11)HP(n,s)=minj∈D(n):Qn,k-1<dn,j-s≤Qn,k{∑m∈C(n)‍Hc¯n,k+p¯n,k(dn,j-s)+hn(dn,j-dn)+∑m∈C(n)‍H(m,dn,j-dn)},(3.12)HNP(n,s)=hn(s-dn)+∑m∈C(n)‍H(m,s-dn),(3.13)H(n,s)=min{HP(n,s),HNP(n,s)}.To obtain the exact optimal solution of SP, it is not necessary to completely characterize the value functionH(0,s). We only need to calculate its values at possible positive discontinuous points. That is, with zero initial inventory assumption for node n∈N∖{0}, we only need to evaluate H(n,s) for s=d0,m-d0,a(n)≥0,m∈N since dn,m-dn=d0,m-d0,n in (3.11). Thus, we give the following dynamic programming for SULS-IQD:Without loss of generality, we assume∣C(n)∣=𝒞 for all n∈N 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 O(∣V(0)∣2max{logV(0),𝒞}), the above algorithm runs in exponential time in terms of T since V(0)=(𝒞T-1)/(𝒞-1). The exponential computational time encourages us to develop a more effective algorithm for the problems with large numbers of stages (see Algorithm 1).Algorithm 1:Dynamic programming algorithm.
For each staget=1 to T
For each noden at stage t
For each possible initial inventorysm=d0,m-d0,a(n)≥0,  m∈N
(ift=1, set sm=0)
step 1: CalculateHP(n,sm) by (3.11)
step 2: CalculateHNP(n,sm) by (3.12) if d0,m-d0,n>=0,
otherwiseHNP(n,sm)=+∞
step 3: CalculateH(n) by (3.13)
End For Iteration (sm)
End For Iteration (n)
End For Iteration (t)


### 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.(
1
) Nonproduction Strategy
At each iteration forsm in DP algorithm for given n, if m∈D(n)∖{n}, step 1 can be skipped and set H(n)=HNP(n,sm) since it violates part 3 of Proposition 3.1. In such case, even if H(n)=HP(n,sm)<HNP(n,sm), and there exist optimal solutions with ∑k=1Kn‍xn,k>0, 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
(3.14)ifm∈D(n)∖{n},thenskipstep1andsetH(n)=HNP(n,sm)atstep3.(
2
) Look-Forward Strategy
At step 1 for given noden at stage t and initial inventory sm=d0,m-d0,a(n), (3.11) requires calculating all the positive orders that cover demand from node n to a certain node j∈D(n). These complete enumeration calculations are very time-consuming, thus we define the following look-forward strategy with parameter FT (forward time stage):
(3.15)HP(n,s)=minj∈DFT(n):Qn,k-1≤dn,j-s≤Qn,k{∑m∈C(n)‍c¯n,k+p¯n,k(dn,j-s)+hn(dn,j-dn)+∑m∈C(n)‍H(m,d0,j-d0,n)},
whereDFT(n)={j∣j∈D(n) and t(j)-t(n)≤FT} denotes the set of node n’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.(
3
) Look-Backward Strategy
At step 2 for given noden at stage t and initial inventory sm=d0,m-d0,a(n), if stage t(m) is smaller than t, it is likely sm-dn=d0,m-d0,n<0 and in such cases step 2 can be skipped. For the same reason, if t-t(m) is larger than a certain value, we can expect sm=d0,m-d0,a(n)<0 and this iteration for sm can be skipped too. Look-backward strategy with parameter BF (backward time stage) can be stated as follows:
(3.16)ift-t(m)>BF,thenskipthisiterationforsm;elseift-t(m)=BF,skipstep2andsetH(n)=HP(n,sm)atstep3,
where BF depends on the distribution of demand. We will set properBF based on distribution of demand by numerical experiments in the next section.
In the above strategies, the major modification comes from the iteration forsm and calculation of HP(n) at step 1 for given node n. 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 noden∈N∖{0}, the rolling horizon heuristic with FT and BT only needs to evaluate the following values:
(3.17){H(n,d0,m-d0,a(n)):d0,m-d0,a(n)≥0,m∈RHFT,BT(n)},
where RHFT,BT(n)=∪m∈P(0,a(n)),t(n)-t(m)≤FT+BTDFT(m).Proof.
This lemma is proven by induction from nodes at staget=1 to leaf nodes at stage T. Since positive order must be made at node 0 due to the zero initial inventory assumption, thus for node n at stage t(n)=2 by (3.15) possible initial inventory sa(n) at node n could only be a certain value in set {d0,m-d0,a(n):d0,m-d0,a(n)≥0,m∈RHFT,BT(n)}={d0,m-d0,0:d0,m-d0,0≥0,m∈DFT(0)}. Suppose the lemma holds for node a(n) at stage t, where 2≤t≤T, that is, the set of initial inventory for node a(n) is given by {d0,m-d0,a(a(n)):d0,m-d0,a(a(n))≥0,m∈RHFT,BT(a(n))}.
Now consider the possible initial inventory at noden.
Case1.If a positive order has been made at node a(n), by (3.15) in the look-forward strategy, the set of initial inventory at nodes n is given by {d0,m-d0,a(n):d0,m-d0,a(n)≥0,m∈DFT(a(n))}.
Case2.if no order has been made at node a(n), the set of possible initial inventory at node n is given by {d0,m-d0,a(n):d0,m-d0,a(n)≥0,m∈RHFT,BT(a(n))}. However, note that by look-backward strategy (3.16) we skip the calculation of HNP(a(n),d0,m-d0,a(n)) for those m such that t(a(n))-t(m)=BF; thus, in this case we only need to consider initial inventory value in set {d0,m-d0,a(n)≥0,m∈∪{DFT(m):m∈P(0,a(a(n))),t(a(n))-t(m)≤FT+BT-1}}  =  {d0,m-d0,a(n)≥0,m∈∪{DFT(m):m∈P(0,a(a(n))),  t(n)-t(m)≤FT+BT}} at node n.
Combine the above two cases, the conclusion holds for noden.Figure1 gives an intuitive example of a balanced scenario tree with 𝒞=2 at each node. For the rolling horizon heuristic with FT=2 and BT=1, given RHH2,1(a(n)) for node a(n), consider the initial inventory set for node n. If a positive order has been made at node a(n), {d0,j-d0,a(n):d0,j-d0,a(n)≥0,j∈DFT(a(n))andt(j)=t(n)+FT-1} will be introduced by (3.15). By the look-backward strategy (3.16), we do not need to calculate H(n,sm), where t(m)=t(n)-BT-1. Other values will be inherited if no order has been made at node a(n). Note that ∣RHFT,BT(n)∣ is no more than (FT+BT)𝒞FT where we assume that ∣C(n)∣=𝒞 for all n∈N.Figure 1
Example of RHH withFT=2 and BT=1.By summarizing the above analysis, the rolling horizon heuristic is given as in Algorithm2.Algorithm 2:Rolling horizon heuristic (FT, BT).
For each staget=1 to T
For each noden at stage t
For each initial inventorysm∈{d0,m-d0,a(n)≥0:m∈RHFT,BT(n)}
SetHP(n,sm)=HNP(n,sm)=+∞; (if t=1,  set sm=0)
step 1: Ifm∈D(n)∖{n}, go  to step 3; otherwise calculate HP(n,sm)
by (3.15)
step 2: Ift-t(m)=BF or d0,n-d0,m<0, goto step 3; otherwise
calculateHNP(n,sm) by (3.12)
step 3: CalculateH(n) by (3.13)
End For Iteration (sm)
End For Iteration (n)
End For Iteration (t)Next the computation complexity of RHH is analyzed. For given noden at stage t, there exist at most ∣RHFT,BT(n)∣ possible initial inventories sm. For each given sm, it takes O(Kn∣C(n)∣∣DFT(n)∣) time to complete calculation of HP(n,sm) at step 1. Step 2 and step 3 will be completed in O(C(n)) time. Thus the total run time can be estimated as (assume that Kn=K for all n∈N)
(3.18)∑n∈N‍{K𝒞|RHFT,BT(n)||DFT(n)|}≤K𝒞(FT+BT)𝒞FT∑n∈N‍{|DFT(n)|}=K𝒞(FT+BT)𝒞FT𝒞FT+1-1𝒞-1|D(0)|=K(FT+BT)𝒞FT+1𝒞-1{𝒞FT+1-1}|D(0)|≤K(FT+BT)𝒞2(FT+1)|D(0)|.Proposition 3.3.
For any instance of SULS-IDQ, the rolling horizon heuristic with parameterFT and BT runs in O(K(FT+BT)𝒞2(FT+1)∣D(0)∣) time.The above analysis can be applied to the dynamic programming algorithm, thus the total run time for DP is given by(3.19)∑n∈N‍{K𝒞|D(0)||D(n)|}=K𝒞|D(0)|∑n∈N‍{|D(n)|}=K𝒞|D(0)|∑t=1T‍𝒞t-1{1+𝒞1+⋯+𝒞T-t}=K|D(0)|𝒞𝒞-1{T𝒞T-𝒞T-1𝒞-1}.
By comparison of computational complexity in (3.18) and (3.19), we observe that RHH works more efficiently for large T. Effective parameters settings of RHH will be further explored by numerical experiments in the next section.


## 3.1. Extended Production-Path Property

Formulation ESP2 presented in Section2 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 n, the order quantity exactly covers the demand along the path from node n to a certain descendant node m; moreover, no positive order will be made for these nodes along the path from node n to the parent node of m.Proposition 3.1 (extended production-path property).
For any instance of SULS-IDQ, there exists an optimal solution(x*,y*,s*), such that for each node n∈N, if xn,k*>0, then there exist m∈D(n) and Qn,k-1<dn,m-sa(n)*≤Qn,k such that(1)
xn,i*=0,yn,i*=0,foralli=1,…,Kn,i≠k;(2)
xn,k*=dn,m-sa(n)*,yn,k*=1;(3)
xl,i*=0,yl,i*=0,foralll∈P(n,a(m))∖{n},i=1,…,Kl.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,(x,y,s), if xn,k>0, we introduce the following definition:
(3.1)NP(n)={m∣m∈D(n),∑k=1Kl‍xl,k=0,∀l∈P(n,a(m))∖{n},∑k=1Km‍xm,k>0},NNP(n)={m∣m∈D(n)∩L,∑k=1Kl‍xl,k=0,∀l∈P(n,m)∖{n}},
where set NP(n) contains the first node with positive order quantity after node n, and set NNP(n) contains these leaf nodes from which to node n no positive order has been placed except node n. Note that if a certain positive order quantity is properly transferred between node n and all the nodes in set NP(n), the feasibility can be kept and only the cost for nodes in the following set is changed:
(3.2)ϕ(n)={m∣m∈P(n,l),∀l∈NP(n)∪NNP(n)}.Proof.
First, it is asserted that there exists at least an optimal solution by Weierstrass’ theorem since the feasible set is compact (note that∑k=1Kn‍xn,k can be bounded by maxm∈D(n)∪L{dn,m}) and the objective function is continuous. Then for any given optimal solution of ESP2, (x*,y*,s*), if xn,k*>0, 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(x*,y*,s*) from stage 1 to stage T. Assume that there exists a node n with xn,k*>0 not holding the property, but for nodes at stage t<t(n) the conclusion holds, the following approach adjusts variables’ value assigned to nodes in D(n) to satisfy the property without changing variables’ value before stage t(n). For any feasible solution (x,y,s) of ESP2, we define the object cost function as
(3.3)ℱ(x,y,s)=∑n∈N‍{∑k=1Kn‍{p¯n,kxn,k+c¯n,kyn,k}+h¯nsn}.
The objective cost for(x*,y*,s*) is
(3.4)ℱ(x*,y*,s*)=∑m∈N∖ϕ(n)‍{∑j=1Km‍{p¯m,jxm,j*+c¯m,jym,j*}+h¯msm*}+p¯n,kxn,k*+c¯n,kyn,k*+h¯nsn*+∑m∈ϕ(n){{n}∪NP(n)}‍h¯msm*+∑m∈NP(n)‍{p¯m,jmxm,jm*+c¯m,jmym,jm*+h¯msm*},
where jm is the unique positive order quantity at node m,m∈NP(n), and the node set N is divided into four subsets: {n}, NP(n), ϕ(n)∖{{n}∪NP(n)}, and N∖ϕ(n). For any node m∈ϕ(n)∖{{n}∪NP(n)}, no order is placed.
Sincexm,jm*>0 for m∈NP(n) (by definition) and sm*>0 for m∈NNP(n) (otherwise it contradicts the assumption), we choose small positive scalar ɛ satisfying
(3.5)0<ɛ≤xn,k*,(3.6)0<ɛ≤sm*,form∈NNP(n),(3.7)0<ɛ≤xm,jm*,form∈NP(n),
such that the following solutions are feasible for ESP2
(3.8)x¯n,k*=xn,k*-ɛ,s¯m*=sm*-ɛform∈ϕ(n)∖NP(n),x¯m,jm*=xm,jm*+ɛforn∈NP(n),x^n,k*=xn,k*+ɛ,s^m*=sm*+ɛform∈ϕ(n)∖NP(n),(3.9)x^m,jm*=xm,jm*-ɛforn∈NP(n)
and keep the value of other variables unchanged. For given ɛ satisfying constraints (3.5)–(3.7), the objective costs of (x¯*,y¯*,s¯*) and (x^*,y^*,s^*) are
(3.10)ℱ(x¯*,y¯*,s¯*)=∑m∈N∖ϕ(n)‍{∑j=1Km‍{p¯m,jxm,j*+c¯m,jym,j*}+h¯msm*}+p¯n,k(xn,k*-ɛ)+c¯n,kyn,k*+h¯n(sn*-ɛ)+∑m∈ϕ(n)∖{{n}∪NP(n)}‍h¯m(sm*-ɛ)+∑m∈NP(n)‍{p¯m,jm(xm,jm*+ɛ)+c¯m,jmym,jm*+h¯msm*}=ℱ(x*,y*,s*)-p¯n,kɛ-h¯nɛ-∑m∈ϕ(n)∖{{n}∪NP(n)}‍h¯mɛ+∑m∈NP(n)‍p¯m,jmɛ=ℱ(x*,y*,s*)-p¯n,kɛ-∑m∈ϕ(n)∖NP(n)‍h¯mɛ+∑m∈NP(n)‍p¯m,jmɛ,ℱ(x^*,y^*,s^*)=∑m∈N∖ϕ(n)‍{∑j=1Km‍{p¯m,jxm,j*+c¯m,jym,j*}+h¯msm*}+p¯n,k(xn,k*+ɛ)+c¯n,kyn,k*+h¯n(sn*+ɛ)+∑m∈ϕ(n)∖{{n}∪NP(n)}‍h¯m(sm*+ɛ)+∑m∈NP(n)‍{p¯m,jm(xm,jm*-ɛ)+c¯m,jmym,jm*+h¯msm*}=ℱ(x*,y*,s*)+p¯n,kɛ+h¯nɛ+∑m∈ϕ(n)∖{{n}∪NP(n)}‍h¯mɛ-∑m∈NP(n)‍p¯m,jmɛ.=ℱ(x*,y*,s*)+p¯n,kɛ+∑m∈ϕ(n)∖NP(n)‍h¯mɛ-∑m∈NP(n)‍p¯m,jmɛ.
Note that the first equality comes from the definition of (x¯*,y¯*,s¯*), (x^*,y^*,s^*), and ℱ function. The second one is obtained by comparing with the value of ℱ(x*,y*,s*) and the last one is obtained by rearranging the terms in the former, where h¯nɛ+∑m∈ϕ(n)∖{{n}∪NP(n)}‍h¯mɛ=∑m∈ϕ(n)∖NP(n)‍h¯mɛ.
Since(x*,y*,s*) is an optimal solution of ESP2, thus we must have p¯n,k+∑m∈ϕ(n)∖NP(n)‍h¯m-∑m∈NP(n)‍p¯m,jm=0 and (x¯*,y¯*,s¯*), (x^*,y^*,s^*) 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 n and (x¯*,y¯*,s¯*) will be used to construct an eventual solution by a similar method.
Case  2. If equality in (3.6) for a certain m holds, then there exists m∈D(n) such that x¯n,k*=dn,m-sa(n)* and x¯l,k*=0,y¯l,k*=0foralll∈P(n,a(m))∖{n},k=1,…,Kl. Next, (x¯*,y¯*,s¯*) will be used to construct an eventual solution by a similar method.
Case  3. If equality in (3.7) for a certain m holds, then we apply the above analysis to (x^*,y^*,s^*). Since there are finite nodes in D(m), eventually case 1 or case 2 will occur (in the worst case NP(n)=∅ 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]:H(n,s): value function at node n when its initial inventory is sa(n)=s,  that is, H(n,s)=min{∑m∈D(n)‍{∑k=1Km‍{p¯m,kxm,k+c¯m,kym,k}+h¯msm}}  subject to constraints (2.18) and sa(n)=s:HP(n,s): production value function at n when its initial inventory is sa(n)=s and ∑k=1Kn‍xn,k>0:HNP(n,s): nonproduction value function at n when its initial inventory is sa(n)=s and ∑k=1Kn‍xn,k=0.From Proposition3.1, if a positive order quantity is made at node n when the initial inventory is s, then there exists a node j∈D(n), such that xn=dn,j-s. Therefore, the following equations hold:
(3.11)HP(n,s)=minj∈D(n):Qn,k-1<dn,j-s≤Qn,k{∑m∈C(n)‍Hc¯n,k+p¯n,k(dn,j-s)+hn(dn,j-dn)+∑m∈C(n)‍H(m,dn,j-dn)},(3.12)HNP(n,s)=hn(s-dn)+∑m∈C(n)‍H(m,s-dn),(3.13)H(n,s)=min{HP(n,s),HNP(n,s)}.To obtain the exact optimal solution of SP, it is not necessary to completely characterize the value functionH(0,s). We only need to calculate its values at possible positive discontinuous points. That is, with zero initial inventory assumption for node n∈N∖{0}, we only need to evaluate H(n,s) for s=d0,m-d0,a(n)≥0,m∈N since dn,m-dn=d0,m-d0,n in (3.11). Thus, we give the following dynamic programming for SULS-IQD:Without loss of generality, we assume∣C(n)∣=𝒞 for all n∈N 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 O(∣V(0)∣2max{logV(0),𝒞}), the above algorithm runs in exponential time in terms of T since V(0)=(𝒞T-1)/(𝒞-1). The exponential computational time encourages us to develop a more effective algorithm for the problems with large numbers of stages (see Algorithm 1).Algorithm 1:Dynamic programming algorithm.
For each staget=1 to T
For each noden at stage t
For each possible initial inventorysm=d0,m-d0,a(n)≥0,  m∈N
(ift=1, set sm=0)
step 1: CalculateHP(n,sm) by (3.11)
step 2: CalculateHNP(n,sm) by (3.12) if d0,m-d0,n>=0,
otherwiseHNP(n,sm)=+∞
step 3: CalculateH(n) by (3.13)
End For Iteration (sm)
End For Iteration (n)
End For Iteration (t)


## 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.(
1
) Nonproduction Strategy
At each iteration forsm in DP algorithm for given n, if m∈D(n)∖{n}, step 1 can be skipped and set H(n)=HNP(n,sm) since it violates part 3 of Proposition 3.1. In such case, even if H(n)=HP(n,sm)<HNP(n,sm), and there exist optimal solutions with ∑k=1Kn‍xn,k>0, 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
(3.14)ifm∈D(n)∖{n},thenskipstep1andsetH(n)=HNP(n,sm)atstep3.(
2
) Look-Forward Strategy
At step 1 for given noden at stage t and initial inventory sm=d0,m-d0,a(n), (3.11) requires calculating all the positive orders that cover demand from node n to a certain node j∈D(n). These complete enumeration calculations are very time-consuming, thus we define the following look-forward strategy with parameter FT (forward time stage):
(3.15)HP(n,s)=minj∈DFT(n):Qn,k-1≤dn,j-s≤Qn,k{∑m∈C(n)‍c¯n,k+p¯n,k(dn,j-s)+hn(dn,j-dn)+∑m∈C(n)‍H(m,d0,j-d0,n)},
whereDFT(n)={j∣j∈D(n) and t(j)-t(n)≤FT} denotes the set of node n’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.(
3
) Look-Backward Strategy
At step 2 for given noden at stage t and initial inventory sm=d0,m-d0,a(n), if stage t(m) is smaller than t, it is likely sm-dn=d0,m-d0,n<0 and in such cases step 2 can be skipped. For the same reason, if t-t(m) is larger than a certain value, we can expect sm=d0,m-d0,a(n)<0 and this iteration for sm can be skipped too. Look-backward strategy with parameter BF (backward time stage) can be stated as follows:
(3.16)ift-t(m)>BF,thenskipthisiterationforsm;elseift-t(m)=BF,skipstep2andsetH(n)=HP(n,sm)atstep3,
where BF depends on the distribution of demand. We will set properBF based on distribution of demand by numerical experiments in the next section.
In the above strategies, the major modification comes from the iteration forsm and calculation of HP(n) at step 1 for given node n. 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 noden∈N∖{0}, the rolling horizon heuristic with FT and BT only needs to evaluate the following values:
(3.17){H(n,d0,m-d0,a(n)):d0,m-d0,a(n)≥0,m∈RHFT,BT(n)},
where RHFT,BT(n)=∪m∈P(0,a(n)),t(n)-t(m)≤FT+BTDFT(m).Proof.
This lemma is proven by induction from nodes at staget=1 to leaf nodes at stage T. Since positive order must be made at node 0 due to the zero initial inventory assumption, thus for node n at stage t(n)=2 by (3.15) possible initial inventory sa(n) at node n could only be a certain value in set {d0,m-d0,a(n):d0,m-d0,a(n)≥0,m∈RHFT,BT(n)}={d0,m-d0,0:d0,m-d0,0≥0,m∈DFT(0)}. Suppose the lemma holds for node a(n) at stage t, where 2≤t≤T, that is, the set of initial inventory for node a(n) is given by {d0,m-d0,a(a(n)):d0,m-d0,a(a(n))≥0,m∈RHFT,BT(a(n))}.
Now consider the possible initial inventory at noden.
Case1.If a positive order has been made at node a(n), by (3.15) in the look-forward strategy, the set of initial inventory at nodes n is given by {d0,m-d0,a(n):d0,m-d0,a(n)≥0,m∈DFT(a(n))}.
Case2.if no order has been made at node a(n), the set of possible initial inventory at node n is given by {d0,m-d0,a(n):d0,m-d0,a(n)≥0,m∈RHFT,BT(a(n))}. However, note that by look-backward strategy (3.16) we skip the calculation of HNP(a(n),d0,m-d0,a(n)) for those m such that t(a(n))-t(m)=BF; thus, in this case we only need to consider initial inventory value in set {d0,m-d0,a(n)≥0,m∈∪{DFT(m):m∈P(0,a(a(n))),t(a(n))-t(m)≤FT+BT-1}}  =  {d0,m-d0,a(n)≥0,m∈∪{DFT(m):m∈P(0,a(a(n))),  t(n)-t(m)≤FT+BT}} at node n.
Combine the above two cases, the conclusion holds for noden.Figure1 gives an intuitive example of a balanced scenario tree with 𝒞=2 at each node. For the rolling horizon heuristic with FT=2 and BT=1, given RHH2,1(a(n)) for node a(n), consider the initial inventory set for node n. If a positive order has been made at node a(n), {d0,j-d0,a(n):d0,j-d0,a(n)≥0,j∈DFT(a(n))andt(j)=t(n)+FT-1} will be introduced by (3.15). By the look-backward strategy (3.16), we do not need to calculate H(n,sm), where t(m)=t(n)-BT-1. Other values will be inherited if no order has been made at node a(n). Note that ∣RHFT,BT(n)∣ is no more than (FT+BT)𝒞FT where we assume that ∣C(n)∣=𝒞 for all n∈N.Figure 1
Example of RHH withFT=2 and BT=1.By summarizing the above analysis, the rolling horizon heuristic is given as in Algorithm2.Algorithm 2:Rolling horizon heuristic (FT, BT).
For each staget=1 to T
For each noden at stage t
For each initial inventorysm∈{d0,m-d0,a(n)≥0:m∈RHFT,BT(n)}
SetHP(n,sm)=HNP(n,sm)=+∞; (if t=1,  set sm=0)
step 1: Ifm∈D(n)∖{n}, go  to step 3; otherwise calculate HP(n,sm)
by (3.15)
step 2: Ift-t(m)=BF or d0,n-d0,m<0, goto step 3; otherwise
calculateHNP(n,sm) by (3.12)
step 3: CalculateH(n) by (3.13)
End For Iteration (sm)
End For Iteration (n)
End For Iteration (t)Next the computation complexity of RHH is analyzed. For given noden at stage t, there exist at most ∣RHFT,BT(n)∣ possible initial inventories sm. For each given sm, it takes O(Kn∣C(n)∣∣DFT(n)∣) time to complete calculation of HP(n,sm) at step 1. Step 2 and step 3 will be completed in O(C(n)) time. Thus the total run time can be estimated as (assume that Kn=K for all n∈N)
(3.18)∑n∈N‍{K𝒞|RHFT,BT(n)||DFT(n)|}≤K𝒞(FT+BT)𝒞FT∑n∈N‍{|DFT(n)|}=K𝒞(FT+BT)𝒞FT𝒞FT+1-1𝒞-1|D(0)|=K(FT+BT)𝒞FT+1𝒞-1{𝒞FT+1-1}|D(0)|≤K(FT+BT)𝒞2(FT+1)|D(0)|.Proposition 3.3.
For any instance of SULS-IDQ, the rolling horizon heuristic with parameterFT and BT runs in O(K(FT+BT)𝒞2(FT+1)∣D(0)∣) time.The above analysis can be applied to the dynamic programming algorithm, thus the total run time for DP is given by(3.19)∑n∈N‍{K𝒞|D(0)||D(n)|}=K𝒞|D(0)|∑n∈N‍{|D(n)|}=K𝒞|D(0)|∑t=1T‍𝒞t-1{1+𝒞1+⋯+𝒞T-t}=K|D(0)|𝒞𝒞-1{T𝒞T-𝒞T-1𝒞-1}.
By comparison of computational complexity in (3.18) and (3.19), we observe that RHH works more efficiently for large T. 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 ofFT and BT 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 C, the number of stages T, the number of nodes N, 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 dn,k will be further defined by σ=8*cv (coefficient of variation) in the following subsection. Due to decreasing unit order price assumption, the unit order price and setup cost at node n is generated by p¯n,k=(1-α)p¯n,k-1 and c¯n,k=c¯n,k-1+αp¯n,k-1Q¯n,k-1, where α is the discount factor and we set α=0.05. Our experiments are conducted on a workstation clocked at 2.33 GHz and equipped with 11.9 GB of core memory.Table 1
Problem instances.
Instance
C
T
N
Int.
Cont.
G
1
2
10
1023
3069
7161
G
2
2
11
2047
6141
8188
G
3
2
12
4095
12285
16380
G
4
3
7
1093
3279
4372
G
5
3
8
3280
9840
13120
G
6
3
9
9841
29523
39564
G
7
4
6
1365
4095
5460
G
8
4
7
5461
16383
21844Table 2
Problem parameters.
Parameter
Distribution
h
-
n
,
k
Uniform(1,2)
p
-
n
,
1
Uniform(1,2)
c
-
n
,
1
Uniform(3,5)
Q
-
n
,
k
Uniform(5,15)
d
n
,
k
Truncated normal(8,σ2)


### 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 parametercv varies from 0.05 to 0.25 for each group of problems. RHH1 with parameters FT=⌈T/3⌉ and BT=1 has been designed to obtain good enough solution in a short time, while RHH2 with parameters FT=⌈T/2⌉ and BT=2 has been designed to obtain better solutions, where ⌈x⌉ denotes the smallest integer larger than x. Table 3 gives the optimal objective function values and CPU time in seconds for DP, and the percentage relative error (RE(%)) and CPU time for RHH1 and RHH2, where
(4.1)RE=ValueofRHH-ValueofDPValueofDP×100%.
For RHH1, the computational results show that RHH1 is able to solve all the instances in the shortest time with no more than 3% relative error except, for instance, G7. The solution quality and run time are further improved when the number of stages T becomes larger and the worst results of G7 can be interpreted as the look-forward time horizon being too small to obtain near optimal solution. For RHH2, the relative error has been further decreased compared with RHH1. Although the run time has been increased, RHH2 still has advantage in run time compared with DP 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.Table 3
Performance of RHH with different (FT, BT).
Instance
cv
RH
H
1
RH
H
2
DP
RE (%)
CPU (sec)
RE (%)
CPU (sec)
Value
CPU (sec)
G
1
-
1
0.05
0.289
0.016
0.0
0.094
12968.28
0.100
G
1
-
2
0.15
0.168
0.015
0.0
0.110
12991.93
0.985
G
1
-
3
0.25
0.135
0.015
0.0
0.110
13082.23
0.969
G
2
-
1
0.05
0.039
0.079
0.0
0.219
25694.23
4.672
G
2
-
2
0.15
0.015
0.079
0.0
0.219
26111.22
4.672
G
2
-
3
0.25
0.005
0.078
0.0
0.218
26447.11
4.703
G
3
-
1
0.05
0.009
0.172
0.0
1.109
51470.76
22.047
G
3
-
2
0.15
0.008
0.172
0.0
1.109
51722.50
22.109
G
3
-
3
0.25
0.006
0.156
0.0
1.110
52772.90
22.172
G
4
-
1
0.05
2.431
0.047
0.509
0.156
8322.42
0.656
G
4
-
2
0.15
1.529
0.031
0.195
0.125
8998.02
0.641
G
4
-
3
0.25
1.391
0.047
0.158
0.125
9126.63
0.640
G
5
-
1
0.05
2.877
0.140
0.635
0.516
24898.40
7.344
G
5
-
2
0.15
1.869
0.125
0.409
0.516
26456.16
7.359
G
5
-
3
0.25
1.212
0.141
0.190
0.500
27716.36
7.375
G
6
-
1
0.05
2.950
0.438
0.214
5.937
74120.22
81.484
G
6
-
2
0.15
1.879
0.453
0.064
5.969
79096.31
81.563
G
6
-
3
0.25
1.434
0.453
0.054
6.047
82641.24
81.844
G
7
-
1
0.05
7.894
0.031
1.270
0.093
7814.61
0.828
G
7
-
2
0.15
5.590
0.016
0.783
0.093
8551.74
0.813
G
7
-
3
0.25
4.350
0.015
0.435
0.094
9024.62
0.812
G
8
-
1
0.05
1.697
0.437
0.247
2.219
30614.77
16.969
G
8
-
2
0.15
1.001
0.453
0.075
2.234
33982.89
16.938
G
8
-
3
0.25
0.597
0.437
0.053
2.234
36099.06
17.032Then 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 ∥xk-xk-1∥≤ɛ or the iteration reaches MaxIteration, then it stops, where MaxIteration=20. To compare the solution quality for given time limit, we set the time limit of CPLEX solver to 5 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 RHH2), 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.Table 4
Comparison with DSSP and CPLEX.
Instance
RHH
DSSP
CPLEX
DP
RE (%)
CPU (sec)
RE (%)
CPU (sec)
RE (%)
CPU (sec)
Value
CPU (sec)
G
1
-
1
0.0
0.094
2.225
1.062
11.881
5.110
11452.16
0.985
G
1
-
2
0.0
0.109
2.299
0.907
10.441
5.063
11653.62
0.984
G
1
-
3
0.0
0.094
2.581
1.062
9.555
5.015
11685.286
0.984
G
2
-
1
0.0
0.516
1.986
2.172
0.619
25.750
22602.25
4.719
G
2
-
2
0.0
0.516
2.583
2.438
3.800
24.093
23408.59
4.719
G
2
-
3
0.0
0.500
3.040
2.094
6.103
24.110
23259.73
4.719
G
3
-
1
0.0
1.140
1.947
4.578
0.645
111.891
45820.49
22.187
G
3
-
2
0.0
1.125
2.845
4.250
1.162
114.657
46665.20
22.172
G
3
-
3
0.0
1.125
3.104
4.406
0.858
113.641
47492.82
22.218
G
4
-
1
0.704
0.140
0.721
1.015
3.744
3.391
7285.71
0.656
G
4
-
2
0.110
0.140
0.893
0.828
15.302
3.375
7916.02
0.656
G
4
-
3
0.035
0.140
1.089
0.735
13.440
3.297
8279.95
0.641
G
5
-
1
0.867
0.515
0.484
3.313
3.334
37.516
21619.82
7.453
G
5
-
2
0.174
0.500
0.632
2.796
1.199
38.265
23767.47
7.422
G
5
-
3
0.198
0.516
1.163
3.110
0.814
37.703
24543.10
7.438
G
6
-
1
0.202
6.063
0.401
8.781
0.800
414.046
65133.43
82.640
G
6
-
2
0.068
6.109
0.747
11.453
1.879
412.860
70243.82
82.375
G
6
-
3
0.022
6.110
1.034
11.157
0.843
419.937
74203.32
82.563
G
7
-
1
0.201
0.360
0.716
1.078
4.300
4.250
6681.52
0.828
G
7
-
2
0.0
0.359
0.889
0.859
1.742
4.172
7480.93
0.812
G
7
-
3
0.0
0.344
1.368
1.188
1.188
4.2347
7980.60
0.828
G
8
-
1
0.367
2.204
0.429
4.578
3.918
86.125
26885.65
17.265
G
8
-
2
0.082
2.172
0.579
4.812
0.741
87.125
30213.60
17.234
G
8
-
3
0.0
2.171
1.060
5.594
1.111
89.468
32503.12
17.157In 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.


## 4.1. Problem Instance

In order to evaluate performance of the proposed DP and RHH, and explore the proper parameters settings ofFT and BT 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 C, the number of stages T, the number of nodes N, 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 dn,k will be further defined by σ=8*cv (coefficient of variation) in the following subsection. Due to decreasing unit order price assumption, the unit order price and setup cost at node n is generated by p¯n,k=(1-α)p¯n,k-1 and c¯n,k=c¯n,k-1+αp¯n,k-1Q¯n,k-1, where α is the discount factor and we set α=0.05. Our experiments are conducted on a workstation clocked at 2.33 GHz and equipped with 11.9 GB of core memory.Table 1
Problem instances.
Instance
C
T
N
Int.
Cont.
G
1
2
10
1023
3069
7161
G
2
2
11
2047
6141
8188
G
3
2
12
4095
12285
16380
G
4
3
7
1093
3279
4372
G
5
3
8
3280
9840
13120
G
6
3
9
9841
29523
39564
G
7
4
6
1365
4095
5460
G
8
4
7
5461
16383
21844Table 2
Problem parameters.
Parameter
Distribution
h
-
n
,
k
Uniform(1,2)
p
-
n
,
1
Uniform(1,2)
c
-
n
,
1
Uniform(3,5)
Q
-
n
,
k
Uniform(5,15)
d
n
,
k
Truncated normal(8,σ2)


## 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 parametercv varies from 0.05 to 0.25 for each group of problems. RHH1 with parameters FT=⌈T/3⌉ and BT=1 has been designed to obtain good enough solution in a short time, while RHH2 with parameters FT=⌈T/2⌉ and BT=2 has been designed to obtain better solutions, where ⌈x⌉ denotes the smallest integer larger than x. Table 3 gives the optimal objective function values and CPU time in seconds for DP, and the percentage relative error (RE(%)) and CPU time for RHH1 and RHH2, where
(4.1)RE=ValueofRHH-ValueofDPValueofDP×100%.
For RHH1, the computational results show that RHH1 is able to solve all the instances in the shortest time with no more than 3% relative error except, for instance, G7. The solution quality and run time are further improved when the number of stages T becomes larger and the worst results of G7 can be interpreted as the look-forward time horizon being too small to obtain near optimal solution. For RHH2, the relative error has been further decreased compared with RHH1. Although the run time has been increased, RHH2 still has advantage in run time compared with DP 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.Table 3
Performance of RHH with different (FT, BT).
Instance
cv
RH
H
1
RH
H
2
DP
RE (%)
CPU (sec)
RE (%)
CPU (sec)
Value
CPU (sec)
G
1
-
1
0.05
0.289
0.016
0.0
0.094
12968.28
0.100
G
1
-
2
0.15
0.168
0.015
0.0
0.110
12991.93
0.985
G
1
-
3
0.25
0.135
0.015
0.0
0.110
13082.23
0.969
G
2
-
1
0.05
0.039
0.079
0.0
0.219
25694.23
4.672
G
2
-
2
0.15
0.015
0.079
0.0
0.219
26111.22
4.672
G
2
-
3
0.25
0.005
0.078
0.0
0.218
26447.11
4.703
G
3
-
1
0.05
0.009
0.172
0.0
1.109
51470.76
22.047
G
3
-
2
0.15
0.008
0.172
0.0
1.109
51722.50
22.109
G
3
-
3
0.25
0.006
0.156
0.0
1.110
52772.90
22.172
G
4
-
1
0.05
2.431
0.047
0.509
0.156
8322.42
0.656
G
4
-
2
0.15
1.529
0.031
0.195
0.125
8998.02
0.641
G
4
-
3
0.25
1.391
0.047
0.158
0.125
9126.63
0.640
G
5
-
1
0.05
2.877
0.140
0.635
0.516
24898.40
7.344
G
5
-
2
0.15
1.869
0.125
0.409
0.516
26456.16
7.359
G
5
-
3
0.25
1.212
0.141
0.190
0.500
27716.36
7.375
G
6
-
1
0.05
2.950
0.438
0.214
5.937
74120.22
81.484
G
6
-
2
0.15
1.879
0.453
0.064
5.969
79096.31
81.563
G
6
-
3
0.25
1.434
0.453
0.054
6.047
82641.24
81.844
G
7
-
1
0.05
7.894
0.031
1.270
0.093
7814.61
0.828
G
7
-
2
0.15
5.590
0.016
0.783
0.093
8551.74
0.813
G
7
-
3
0.25
4.350
0.015
0.435
0.094
9024.62
0.812
G
8
-
1
0.05
1.697
0.437
0.247
2.219
30614.77
16.969
G
8
-
2
0.15
1.001
0.453
0.075
2.234
33982.89
16.938
G
8
-
3
0.25
0.597
0.437
0.053
2.234
36099.06
17.032Then 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 ∥xk-xk-1∥≤ɛ or the iteration reaches MaxIteration, then it stops, where MaxIteration=20. To compare the solution quality for given time limit, we set the time limit of CPLEX solver to 5 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 RHH2), 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.Table 4
Comparison with DSSP and CPLEX.
Instance
RHH
DSSP
CPLEX
DP
RE (%)
CPU (sec)
RE (%)
CPU (sec)
RE (%)
CPU (sec)
Value
CPU (sec)
G
1
-
1
0.0
0.094
2.225
1.062
11.881
5.110
11452.16
0.985
G
1
-
2
0.0
0.109
2.299
0.907
10.441
5.063
11653.62
0.984
G
1
-
3
0.0
0.094
2.581
1.062
9.555
5.015
11685.286
0.984
G
2
-
1
0.0
0.516
1.986
2.172
0.619
25.750
22602.25
4.719
G
2
-
2
0.0
0.516
2.583
2.438
3.800
24.093
23408.59
4.719
G
2
-
3
0.0
0.500
3.040
2.094
6.103
24.110
23259.73
4.719
G
3
-
1
0.0
1.140
1.947
4.578
0.645
111.891
45820.49
22.187
G
3
-
2
0.0
1.125
2.845
4.250
1.162
114.657
46665.20
22.172
G
3
-
3
0.0
1.125
3.104
4.406
0.858
113.641
47492.82
22.218
G
4
-
1
0.704
0.140
0.721
1.015
3.744
3.391
7285.71
0.656
G
4
-
2
0.110
0.140
0.893
0.828
15.302
3.375
7916.02
0.656
G
4
-
3
0.035
0.140
1.089
0.735
13.440
3.297
8279.95
0.641
G
5
-
1
0.867
0.515
0.484
3.313
3.334
37.516
21619.82
7.453
G
5
-
2
0.174
0.500
0.632
2.796
1.199
38.265
23767.47
7.422
G
5
-
3
0.198
0.516
1.163
3.110
0.814
37.703
24543.10
7.438
G
6
-
1
0.202
6.063
0.401
8.781
0.800
414.046
65133.43
82.640
G
6
-
2
0.068
6.109
0.747
11.453
1.879
412.860
70243.82
82.375
G
6
-
3
0.022
6.110
1.034
11.157
0.843
419.937
74203.32
82.563
G
7
-
1
0.201
0.360
0.716
1.078
4.300
4.250
6681.52
0.828
G
7
-
2
0.0
0.359
0.889
0.859
1.742
4.172
7480.93
0.812
G
7
-
3
0.0
0.344
1.368
1.188
1.188
4.2347
7980.60
0.828
G
8
-
1
0.367
2.204
0.429
4.578
3.918
86.125
26885.65
17.265
G
8
-
2
0.082
2.172
0.579
4.812
0.741
87.125
30213.60
17.234
G
8
-
3
0.0
2.171
1.060
5.594
1.111
89.468
32503.12
17.157In 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].



---

*Source: 582323-2012-10-03.xml*