Lagrangian based approach to solve a two level capacitated lot sizing problem

Two-level, multi-item, multi-period-capacitated dynamic lot-sizing problem with inclusions of backorders and setup times, TL_CLSP_BS, a well-known NP-hard problem, is solved using a novel procedure. Lagrangian relaxation of the material balance constraint reduces TL_CLSP_BS to a single-constraint continuous knapsack problem. Reduced problem is solved using bounded variable linear programs (BVLPs). We obtain promising bounds, which provides a better start to the branch-and-bound procedure. Limited empirical investigations are carried out on four problem sizes. In terms of computational time, the developed procedure is efficient than the CPLEX solver of GAMS. Further, while GAMS could not solve the largest problem size considered here, our procedure could solve the same in around one second time. This clearly highlights the efficacy of the developed procedure. The solution technique is applicable to any problem structure, which is reducible to the application of BVLPs. Subjects: Operations Research; Production Systems; Manufacturing & Processing; Supply Chain Management


PUBLIC INTEREST STATEMENT
A two-level-, multi-item-, multi-period-capacitated dynamic lot-sizing problem with inclusions of backorders and setup times-TL_CLSP_BS-a well-known NP-hard problem, is solved using a novel procedure. Lagrangian relaxation of the material balance constraint reduces TL_CLSP_BS to a single-constraint continuous knapsack problem. Reduced problem is solved using bounded variable linear programs. We obtain promising bounds, which provides a better start to the branch-andbound procedure. Limited empirical investigations are carried out on four problem sizes. In terms of computational time, the developed procedure is efficient than the CPLEX solver of GAMS. Further, while GAMS could not solve the largest problem size considered here, our procedure could solve the same in around one second time. This clearly highlights the efficacy of the developed procedure. The solution technique is applicable to any problem structure, which is reducible to the application of bounded variable linear programs.

Introduction
To create a production plan over a finite number of periods, where the demand has to be fulfilled while facing finite capacities is known to be a difficult problem. It is even harder in case of a multilevel product structure, where items are related to each other by successor-predecessor relationship, as per the bill of material. For the dynamic multi-item multi-level-capacitated lot-sizing problem (ML_CLSP_BS) considered here, the following assumptions apply (Maes, McClain, & Van Wassenhove, 1991): (1) A general two-level product structure is present, without loops in the bill of material.
(2) The planning horizon is divided into a number of periods.
(3) Items with dynamic external demands are arranged in the product structure.
(4) Inventory holding costs are computed based on the end-of-period inventory.
(5) Backorders are allowed and have a backorder cost.
(6) Once an item is produced, a fixed setup time and a fixed setup cost occur in addition to variable production cost.
(7) The objective is to minimize the sum of holding costs, setup costs, backorder costs, and the variable production costs.
There had been studies addressing the solution to a multi-level CLSP in the past, but certain realistically relevant variables were excluded to simplify the solution procedure. Billington, McClain, and Thomas (1986) consider a general product structure and Billington, Blackburn, Maes, Millen, and Van Wassenhove (1994) assume a linear product structure with uniform processing times; they however did not consider the setup times. Maes and Van Wassenhove (1986) and Maes et al. (1991) developed heuristics for the case of linear product structure, where capacity is not shared in between different levels. Kuik, Salomon, Van Wassenhove, and Maes (1993) worked on a linear product structure, for which only a capacitated level exists among multiple levels, they too did not consider setup times. Tempelmeier and Helber (1994) then developed a heuristic for the case of a general product structure, where multiple resources share their capacities by items appearing in different levels of the product structure. Tempelmeier and Derstroff (1996) applied Lagrangian relaxation, França, Armentano, Berretta, and Clark (1997) developed a heuristic procedure, and Segerstedt (1996) used dynamic programming for the general product structure with setup times. Recent reviews of modeling and solution approaches of lot-sizing problems are provided by Drexl and Kimms (1997), Karimi, Fatemi Ghomi, and Wilson (2003) and Jans and Degraeve (2007). Other recent works on multi-level lot-sizing problems are Tempelmeier and Buschkühl (2009), Sahling, Buschkuhl, Tempelmeier, and Helber (2009), Han, Tang, Kaku, and Mu (2008, Pitakaso, Almeder, Doerner, and Hartl (2004), and Dellaert and Jeunet (2003). This paper has attempted a two-level-capacitated lot-sizing problem. Some of the recent literatures on two-level lot-sizing problems are Melo and Wolsey (2010) and Toledo, Arantes, and França (2011).
In Section 2, we discuss the formulation of multi-level CLSP_BS. Section 3 explains a two-level product structure and formulates two-level CLSP_BS by modifying the formulation of multi-level CLSP_BS. In Section 4, the Lagrangian relaxation is applied to the formulation of TL_CLSP_BS and the problem reduces to a knapsack problem. Section 5 details the solution procedure for the reduced knapsack problem. Section 6 provides the computational details of the described approach. Finally, the work and its significance is concluded in Section 7.

Decision variables
XP it number of items "i" to be produced during the period "t" XINV it number of items "i" carried as inventory at the end of period "t" XBO it number of items "i" that will be backordered from period "t" YS it binary variable for setup of the resource for item "i" during the period "t" ∈ {0, 1} i.e. 0 (if there is no setup required), 1 otherwise Parameters CP it unit cost of producing item "i" in period "t" CS it unit cost of setup, for item "i" in period "t" CINV it unit cost of holding inventory of item "i" for 1 period CBO it unit cost of Backordering item "i", which was demanded during the period "t" CAP it capacity available to produce item "i" during the period "t" CAPT t capacity available in time units, in a period "t" D it external demand of item "i" during the period "t" TP i time required to process the item "i" TS i time required to setup the production for item "i" S i set of successors of item "i" in the product structure N ij number of items "i" required to produce 1 unit of its successor "j" in the product structure YS it ∈ [0, 1] ∀i, t "i," for the production of next level of item "j." Refer Figure 1 for a better idea about connection between levels of "i" and "j." Also note that non-end items may have an external demand in addition to the demand that is generated by its successor (internal demand); hence the presence of D it is justified for non-end items too. Equation 3 is the time capacity constraint ensuring that the total time utilized in doing production of all the items is always less than or equal to the maximum time available in any period. Equation 4 is the production capacity constraint, which ensures the production quantity to be always less than or equal to the maximum production capacity available for all items and time periods. Equation 5 is the non-negativity restriction. Equations 6 and 7 assumes that there are no initial or final inventory and backorders. Equation 8 enforces binary restriction on the setup variables. The setup variable is equal to one when setup takes place for an item in a period, otherwise zero.

Two-level problem formulation
Figure 2 is a typical example of the product structure with four items in two levels, showing the successor-predecessor relationship. Item 1 is the end item or the final assembly, for which an external demand would be known as a result of the forecast. Items 2, 3, and 4 are the non-end items or the intermediate sub-assemblies, required in 6, 4, and 8 numbers, respectively, to produce one unit of Item 1.
For the two-level product structure, constraint (2) can be written as: Equation 2′ is true for end as well as non-end items in a two-level product structure. So, formulation of TL_CLSP_BS is given by P: Minimize Equation 1, subject to Equations 2′, and 3-8.

Relaxation and reduction
It is evident from the formulation of the TL_CLSP_BS that Equation 2 is the only constraint that links the two levels of the product structure. Equation 2 is relaxed using Lagrangian procedure to solve a two-level version of the problem.

Relaxing material balance constraint
We use an unrestricted Lagrangian multiplier λ2 it to relax (2), an equality constraint. The problem (P 2 ) takes the following shape for a TL_CLSP_BS. P 2 : Subject to: (Equations 3-8).

Reduction to continuous knapsack problem
To solve P 2 , we temporarily omit the term D it λ2 it , λ2 it XP 1t N i1 and also the terms containing inventory and backorder variables. The problem takes the following shape: Subject to: (Equations 3-5) and (Equation 8).
Note that the omission of some terms at this step is just to make the problem structure amenable to application of the bounded variable linear program (BVLP) solution procedure. The terms omitted now will be included back to the solution procedure at an appropriate step, so that the intent of their presence in the problem structure is not lost. Equation 4 gives an upper bound (=CAP it ) on the value of the production variable XP it . If a procedure can take care of this upper bound in its steps, we can eliminate Equation 4 from the problem structure. The problem now takes a special structure of a single-constraint continuous knapsack problem (P k ) for each time period "t", where the variable XP it has an upper bound (=CAP it ).
This single-constraint continuous knapsack problem can be solved efficiently using the concept of BVLP as described in the following section.

Solution procedure
We now devise a very efficient and easy to implement solution scheme using BVLP (Murty, 1976). BVLP is solved inside the iterations of sub-gradient optimization. This is detailed out as follows.

Determining production and setup variables
The procedure starts with sub-gradient optimization by assuming an initial value for the Lagrangian multipliers. We arrange the set (i, t) in the sequence of ascending ratio of the term . In this order of (i, t), we replace XP it by its upper bound CAP it till the right-hand resource is available. All XP it , which remain unallocated, are made zero. For each set of (i, t), we now determine the sense (negative or positive) of the term [(CP it − 2 it )XP it + CS it ]. For those set of (i, t), for which the term is negative, we make YS it = 1; and for all other (i, t), we make YS it = 0. In case of two-level problems, capacities of the second-level items are substantially higher, when compared to its individual demand D it , where i ∈ {2, I}. Because of this, the lower bounds obtained with these values of XP it (=CAP it ) may be higher than the optimal solution. It is hence necessary to

Inventory and backorder variables
We now get back the values of inventory and backorder variables, which were lost due to the procedure of Lagrangian relaxation and BVLPs.
In case both α_BO it < 0 and α_INV it < 0, the term corresponding to the lower of the two is assigned a magnitude equal to Similar is the case if both α_BO it > 0 and α_INV it > 0, we then allocate a zero value to both variables XBO it = 0 and XINV it = 0.
In this way, we now have determined the values of XP it , YS it , λ2 it , XINV it , and XBO it . We start the sub-gradient optimization iterations to update the value of λ2 it . At each iteration of the sub-gradient optimization procedure, we have an updated λ2 it , and hence a correspondingly updated XP it , YS it , XINV it , and XBO it . Continue sub-gradient iterations till the improvement is very less (say < 0.001). At the end of sub-gradient iterations, we have the final updated values of XP it , YS it , λ2 it , XINV it , XBO it that are used to obtain the value of the "lower bound".

Determining upper bound
In Section 5.2, in order to determine the variables using BVLPs, we may have ended up getting an infeasible solution. To obtain an upper bound, we modify the solution to obtain a feasible solution. Constraints (3-8) are all taken care by the procedure and only constraint left out is Equation 2, which was relaxed using a Lagrangian multiplier. Mathematically, backorders are to be zero at the start and end of the planning horizon. But the total quantity of an item produced over the complete planning horizon (all time periods) is not equal to the total quantity of that item in demand over all time periods, due to the used procedure. This causes infeasibility, which needs to be appropriately dealt with.
In order to find a feasible solution (upper bound), we check if increase XP it for the time period for which CP it is minimum; till either XP it = CAP it or satisfied. After doing this adjustment, we move period by period from the first time period, for which D it and XP it are both known. If XP it − D it is positive, it indicates an overproduction and hence an inventory will be carried. Otherwise, the negative value of XP it − D it suggests that a backorder is to be carried. In this way, we obtain the values of XINV i1 and XBO i1 by observing the sense (negative or positive) of the term XP it − D it , and numerically equal to |XP it − D it |.
From the second period to all the subsequent periods, we watch the sense of the term So now, we get the values of production, setup, inventory, and backorder variables, which can be substituted in the problem "P" to obtain its objective value, which is the "feasible solution" or the "upper bound." In later part of this paper, this upper bound obtained is referred to as the "Initial" solution obtained by "Procedure".

Improvement by branch and bound
We apply branch-and-bound procedure to search the better solutions and prune the branches containing inferior solutions. We use the upper bound obtained above (referred as "current best" in the enumeration procedure) and the lower bound obtained at the end of sub-gradient optimization procedure as the root node of the search tree and then do branching on this root node.
At root node, we have obtained the initial best solution keeping all the binary setup variables free to take any value (0 or 1). As we move forward, we consider the binary variables one by one and start to fix their value to 0 or 1. The solution procedure is applied keeping these binary variables fixed at each branch. If the solution obtained at a branch is better than the "current best," it is updated with this solution, else the branch is pruned. We keep moving from one node to the other by branching on the binary variables and keeping the previous binary variable fixed at which we obtained a better "current best" solution. In this way, all binary variables are considered one by one; the procedure is stopped when all branches are explored or pruned and all binary variables are fixed to either 0 or 1. At the end of this enumeration procedure, we obtain the "current best" feasible solution to the problem CLSP_BS, which is either equal to the optimal solution or acts as an upper bound to the optimal solution with very low duality gap.
The efficacy of this procedure is checked by implementing it through the developed computer codes on a variety of randomly generated problem sets of different sizes. In later part of this paper (Section 6), this upper bound obtained is referred to as the "After Enumeration" solution obtained by "Procedure".

Minimum cost network flow problem
Minimum cost network flow (MCF) problem is solved by obtaining information about the value of setup variables from the procedure discussed in previous sub-sections. In this way, MCF provides us with yet another and a very good upper bound (feasible solution). When MCF is solved on the values of setup variables obtained in "Initial" procedure, the solution thus obtained is called "Initial" solution obtained from MCF. Also during enumeration, at each node, as soon as we get the information about the value of setup variables, a MCF problem is thus solved. The solution so obtained at the end of this enumeration is called "After Enumeration" solution of MCF. It is expected that MCF will give very promising duality gaps in shorter computational times.

Computational experiences
To check the performance of the solution procedure described in the previous section, the computational experiments on the randomized data-sets are performed. The procedure was coded in C and run on 2.40 GHz, 3.87 GB RAM, 64-bit operating system, Intel (R) Core (TM) i5 CPU M 450 stand-alone PCs. Initially, the code was developed in MATLAB, but after running those codes, it was realized that MATLAB was only marginally efficient, when compared with the GAMS performance. As this was against our expectations, we proceed to test a single common routine written in MATLAB and also in C. It was observed that to get the same output, the computational time taken by MATLAB is 6-20 times that of C. The reason to this inefficiency of MATLAB is probably due to its large overheads, that C does not possess. Therefore, we moved on to convert the developed MATLAB code into C; the results given in this paper are hence those obtained by the C code only.

Data
The cost parameters used in the formulation is generated by coding a random generator in C language. We choose uniform distribution to generate parameters which is widely used in the literature of CLSP for this purpose. The following Table 1 states the range of values of different parameters that are taken to form the random problems. It is clear from Table 1 that we assume backorder costs to be higher as compared to the inventory carrying costs, as backorder costs also account for the loss of goodwill for the customer, who could not be served instantly. It also implies that our model is comparatively more open to inventory as compared to the backorders, as prompt satisfaction of the demand is primarily important in today's competitive environment. Also, setup costs are considerably larger as compared to the production cost; similarly, setup time is assumed notably larger than the production time.
Range of the randomly generated values of demand and production capacity almost overlap each other, but the upper and lower limits of capacity is assumed slightly higher than that of demand, so that the scope of inventory is ensured in the model while at the same time allowing backorders too when demand is more than the production capacity. The relation between demand and production capacity can also be stated in terms of tightness factor, which can be defined as the ratio of average periodic demand and production capacity. For an uncapacitated problem, the value of tightness factor will be 0; while for the case when required capacity is exactly equal to the available capacity, the tightness factor is 1. As one can observe that for the data considered in this work, the tightness factor is roughly about 0.8. With such a rigid tightness factor, it is possible that some infeasible problems are generated, such infeasible problems however are selectively eliminated from the problem sets.

Time complexity
Complexity of a mixed-integer programming problem is hard to measure mathematically, as it depends on many factors for a given model. Hence, experimental approach to determine the complexity is the most favorable and adopted practical indicator. Though number of variables (integer/ continuous) and constraints may not necessarily correspond to the computing burden, we provide a summary of the four sizes of problems considered to show the relationships. Number of constraints and variables in the problem are a fair indicator of the computational time. While the number of binary variables is of the order (IT), the order of continuous variables is (3IT) and the number of constraints is (2IT + T).
As the problem size grows, number of variables and constraints increase. Note that the sizes (5 × 10 and 10 × 5) are almost the same as far as the number of binary and continuous variables is considered; there is, however, a change in the number of constraints for the two sizes. It would be interesting to note the relative behavior of these sizes with each other, when we compare their average computational times on the same platform and similar problems (Table 2). Table 3, we note the average computational time and the number of nodes taken by GAMS to solve the problems of various sizes, using the method of branch and bound and also branch and cut. One may note that as the problem size increases, the computational time as well as the number of nodes increases. Also, it is noted that branch and cut is much more efficient when compared to branch and bound. All problem sizes considered in this work are solvable by branch and cut within a second; however, computational time using branch and bound is 3-4 s. For small-sized problems, while branch and bound appears to be more efficient than branch and cut, as problem size grows, computational time of branch and bound grows to the order of around 4 s. For further large-sized problems, branch and bound cease to perform at all, as the solver went out of memory after trying to solve the problem for a few minutes; "N/A" is indicated for such problems in Table 3. We emphasize here that this table is provided here just to give a general idea to the reader about the order of computational time taken by the standard CPLEX solver of GAMS to solve a particular size of the problem. In order to actually compare the problem sizes on different aspects, we perform t-test, the result of which is shown in the next tables.

Computational time and gap: procedure
Now in Table 4, we show the efficacy of the Procedure and MCF to solve the same problem sets which are solved in GAMS (Table 3). Average computational time taken and average optimality gap attained by Procedure and MCF to solve different sizes of CLSP_BS is shown here. In Table 4, "Initial" means that the Procedure (or MCF) is only implemented without enumerating the solutions by the use of branch and bound. "After Enumeration" in the table, as the name suggests is the computational time and gap after the implementation of branch and bound to the solution obtained from "Initial" of Procedure (or MCF). The solution time increases as the size of the solved problems increase. Procedure-Initial could solve the problem with an average optimality gap of 38-72% in 0.05-0.1 s, while the same problems on enumeration are solved much closer to optimality, with a gap of 0.3-1.1% in 1-2 s. The performance of MCF is much better when compared to Procedure. MCF-Initial could solve the problem with an average optimality gap of 0.2-0.5% in 0.01-0.03 s, while there is not much improvement noted after enumeration. This provides the user of the code the option to decide if a quick approximate solution with some optimality gap is acceptable (use the solution of "Initial") or a solution very close to optimal is required while compromising on the solution time (use the solution of "After Enumeration"). Also note here that Procedure-Initial performed much better in case of single-level CLSP_BS. For TL_CLSP_BS, the performance of Procedure-Initial has hugely deteriorated and branch-and-bound enumeration plays a major role in improving the solution. Table 5 is provided next to directly compare the average solution time taken by different techniques to solve the same set of problems on a common platform. When one compares the computational time of GAMS branch and bound (tgmsbb) with that of the Procedure-After Enumeration (tprobb)  and MCF-After Enumeration (tmcfbb), the efficacy of the developed algorithms is indubitably illustrated. While for the largest sized problems, GAMS branch and bound went out of memory, "Procedure" could solve the problem in around 1 s and "MCF" solved the problem in around 0.6 s. For problem size 10 × 5, GAMS branch and bound solved the problem in around 3 s, but Procedure and MCF after enumeration could solve the problem closer to optimality in around (or less than) 1 s. The computational time taken by Procedure-Initial (tproin) and MCF-Initial (tmcfin) are shown here to compare with the computational time of GAMS branch and cut (tgmsbc), which is probably one of the fastest solving techniques available commercially. While for the largest sized problems, GAMS branch and cut took around 0.28 s to solve the problem, "Procedure-Initial" and "MCF-Initial" could solve the problem (of course with some optimality gap) in around 50th fraction of a second. This clearly highlights the efficacy of the developed algorithms. Though, however, Procedure-Initial and MCF-Initial used for the computation of the upper bounds were there to provide a starting point to the branch-and-bound procedure, it actually generated reasonably good bounds.

Comparing computational time and optimality gap of different problem sizes
In Table 6, we show the t-test results comparing the computational times and the optimality gaps of the different size problems. This is done to statistically compare the observations of one size with the other. Note that in the following table, against 5 × 10-5 × 5 row, the column "tgmsbb" or "tgmsbc" indicates the "t" calculated for difference between the GAMS branch-and-cut computational times of (5 × 10/5 × 5) and 1; and "tgmsbb" is the "t" calculated for difference between GAMS branch-andbound computational time for 5 × 10/5 × 5 and 1. Similarly, for "tproin," "gproin," "tprobb," "tmcfin," "gmcfin," "tmcfbb," and all other rows of the table. The t values in Table 6 are compared with the t-critical values given in Table 7 to determine the significance. One may note that in general, the computational time increases as the size of problem increases. Similarly the optimality gap; smaller problems are solved much closer to optimality as compared to the large-sized problems. That is, smaller the problem, efficient is the procedure devised in this work to solve that problem. Note that  against some columns, we have few negative values, for example, the row 10 × 5-5 × 5. This means that statistically the optimality gap obtained by Procedure-Initial of size 5 × 5 is more than that of the size 5 × 10. Similar is the explanation for other rows/columns with a negative value. "N/A" in the third column (tgmsbb) indicates that the computational time to solve one or both of these problems is not available as the solver went out of memory.

Comparing computational time of procedure and GAMS
Next in Table 8, we compare the computational time taken by the Procedure and MCF (both after enumeration), with that taken by GAMS branch and bound and branch and cut. Note that in the following table, against 5 × 5 row, the column "tprobb − tgmsbc" indicates the "t" calculated for difference between the computational time of Procedure-After Enumeration/computational time of GAMS branch and cut and 1. Similarly for the other columns "tprobb -tgmsbb," "tmcfbb -tgmsbb," "tmcfbb -tmcfbc," and all other rows of the table. The t values in Table 6 are compared with the t-critical values given in Table 7 to determine the significance. It is noted that branch and bound of Procedure and MCF is inferior when compared to branch and cut of GAMS. However, branch and bound of Procedure and MCF clearly outperforms the branch and bound of GAMS. GAMS branch and bound even cease to perform for the larger sizes (as evident from column 3 and 5), but branch and bound of Procedure and MCF solve the problems of largest size considered here. This clearly highlights the efficacy of the developed procedure. One may note that while Table 6 provides comparison of computational times and gaps for the two problem sizes, Table 8 compares time taken by different methods for the same problem size.

Conclusion
Very few researchers attempted multi-level-capacitated lot-sizing problems in the past. Considering two-level lot-sizing problem to be a part of the multi-level problem, we make an effort to solve this problem through a novel technique. Two-level, multi-item, multi-period-capacitated dynamic lot-sizing problem with inclusions of backorders and setup times (TL_CLSP_BS) is well known to be NP hard.
Lagrangian relaxation of the material balance constraint leaves the problem structure reducible to a continuous knapsack problem. This problem is solved by the use of BVLPs. The advantage of applying BVLPs is that, by the use basic algebraic computations, we are able to solve a well-known complex mixed-integer linear programming problem close to optimality. Sub-gradient optimization procedure is used to solve the Lagrangian relaxation. After obtaining lower bound through this Lagrangian procedure, the upper bound is obtained as a feasible solution. Another upper bound is obtained by solving a minimal cost network flow problem. These bounds are then improved by implicit enumeration procedure of branch and bound.
The procedure is verified through limited empirical investigations on four different problem sizes. Computational time taken by the developed procedure is noted to be less than that of the commercial solver GAMS. Also, for the largest sized problem considered, branch and bound of GAMS could