Product decomposition strategy for optimization of supply chain planning

Abstract

Optimization of large-scale supply chain planning models requires the application of decomposition strategies to reduce the computational expense. Two major options are to use either spatial or temporal Lagrangean decomposition. In this paper, to further reduce the computational expense a novel decomposition scheme by products is presented. The decomposition is based on a reformulation of knapsack constraints in the problem. The new approach allows for simultaneous decomposition by products and by time periods, enabling the generation of a large number of subproblems, that can be solved by using parallel computing. The case study shows that the proposed product decomposition exhibits similar performance as the temporal decomposition, and that selecting different orders of products and aggregating the linking constraints can improve the efficiency of the algorithm.

Keywords

supply chain planning / Lagrangean decomposition / mixed-integer programming

Introduction

The supply chains of the 21st century are highly complex and spread all across the world (^{Kearney, 2004}). Globalization has brought both access to new markets and increased competition. In this scenario, making the best possible decisions is critical to remain competitive. Making optimal tactical plans involves balancing several trade-offs including inventory, transportation, and production costs. The main decisions are the assignment of production volume to manufacturing facilities, including third parties, inventory levels, and shipments between facilities. The planning horizons addressed range from 6 to 24 months, divided into monthly or even weekly time periods. At the same time, a typical supply chain handles hundreds or thousands of products. The resulting optimization models can be very large, and often require the use of decomposition techniques to obtain good solutions in reasonable time.

We propose a novel Lagrangean decomposition scheme based on decomposing the problem by products, and also a decomposition scheme combining temporal and product decomposition. Lagrangean decomposition (^{Guignard and Kim, 1987}) has been effectively applied to many large-scale optimization problems (^{Wang, 2003}). ^{Graves (1982)} used these techniques for hierarchical planning, decomposing a large problem into two subproblems. ^{Gupta and Maranas (1999)} analyzed the importance of choosing the appropriate relaxation. Along the same lines ^{Jackson and Grossmann (2003)} used temporal decomposition for solving a multi-site, multi-period planning problem. ^{Terrazas-Moreno et al. (2011)} compared temporal and spatial decomposition for supply chain planning. They concluded that temporal decomposition was superior because the relaxation of inventory continuity, generated in temporal decomposition, is closer to the optimal solution compared to the relaxation of material balances obtained when using spatial decomposition. ^{Rodriguez et al. (2014)} used temporal decomposition to solve a planning problem under uncertainty. ^{Calfa et al. (2013)} combined temporal Lagrangean decomposition with bilevel decomposition for the integration of planning and scheduling. The opportunity of decomposing the problem by subproblems is also identified by ^{van Elzakker et al. (2014)}. However, they propose a heuristic procedure, while we generalize these ideas using Lagrangean decomposition.

The update of the Lagrange multipliers is the most critical step of the algorithm, hence, it has been a matter of intensive research. The update methods are based on three basic approaches: Subgradient (^{Fisher, 1985}), cutting planes (^{Kelley, 1960}), and bundle methods that combine the first two (^{Lemarechal et al., 1981}). The cutting plane methods return many unbounded solutions, especially at the beginning of the first iterations. The bundle methods handle these issues adding stabilization terms. Within the cutting plane methods, ^{Mouret et al. (2011)} allow the multipliers to take values within a band of the value given by the subgradient, while ^{Oliveira et al. (2013)} allow the multipliers to take values in a band defined by the direction and the opposite direction given by the subgradient method. Both report convergence improvements compare to plain cutting plane methods. The improvement of the subgradient method has also been a matter of research. ^{Erlenkotter (1978)} proposed a multiplier adjustment method for the Lagrangean relaxation of the uncapacitated location problem, while ^{Wu and Ierapetritou (2006)} use the Nelder-Mead method to identify search directions for adjusting the Lagrangean multipliers. We propose a probing method to improve the step size when updating the multipliers using the subgradient method.

The paper is organized as follows: Section 2 presents an overview of Lagrangean decomposition. In Section 3, a representative instance of the addressed problem is described. In Section 4, temporal decomposition, product decomposition, and simultaneous product and temporal decomposition are derived. The probing subgradient method is described in Section 5. A case study is presented in Section 6, and the main results of the paper are summarized in Section 7.

Lagrangean decomposition overview

^{Held and Karp (1970)} originally proposed the Lagrangean Relaxation for the traveling salesman problem in which the subtour elimination constraints are dualized. The basic idea is to exploit the structure of many optimization problems that involves complicating constraints as shown in Fig. 1. These constraints are added to the objective function with a penalty term proportional to the amount of constraint violation. The resulting relaxation reduces the computational complexity of the solution procedure, because the Lagrangean problem is easier to solve than the original problem, and each subproblem can be solved independently.

Consider the following optimization problem (^{Guignard, 2003}):

(1)P=\underset{x}{\max }=\{z=c^{\mathrm{T}}x|Ax\le b, Dx\le e, x\in X\},

where, X={ℝ}_{+}^{n-p}x{\{\mathrm{0}, \mathrm{1}\}}^pand where Dx\le e are considered to be the complicating constraints. Let l be a vector of nonnegative Lagrangean multipliers. The Lagrangean relaxation of (P) is the following:

(2)LR(\lambda )=\underset{x}{\max }=\{z=c^{\mathrm{T}}x+{\lambda }^{\mathrm{T}}(e-Dx)|Ax\le b, x\in X\}.

Since the feasible region of the relaxation is larger than the original one, the objective is overestimated, and hence its optimal value provides an upper bound of the original maximization problem (P). The best bound that can be obtained by solving the Lagrangean Dual is:

(3)LD=\underset{\lambda }{\mathrm{m}\mathrm{i}\mathrm{n}}LR\left(\lambda \right)=\underset{\lambda }{\mathrm{m}\mathrm{i}\mathrm{n}}\left\{\underset{x}{\mathrm{m}\mathrm{a}\mathrm{x}}=\left\{z=c^{\mathrm{T}}x+{\lambda }^{\mathrm{T}}\left(e-Dx\right)|Ax\le b, x\in X\right\}\right\}.

Furthermore, ^{Geoffrion (1974)} demonstrated the equivalence between the Lagreangean dual (LD) and primal relaxation of (P) denoted by (P^{*}).

(4)P^{\ast }=\underset{x}{\mathrm{m}\mathrm{a}\mathrm{x}}=\left\{z=c^{\mathrm{T}}x\left|Dx\le e, x\in Conv\left(X\right)\right|Ax\le b, x\in X\right\},

where Conv(X) denotes the convex hull of the set X, and v{·} is the optimal solution of {·}. Then,

(5)v\left\{LP\right\}\ge v\left\{LR\right\}\ge LD=v\left\{P^{\ast }\right\}\ge v\left\{P\right\},

where LP is the LP relaxation of (P). This relationship is explained by the modification of the original structure of the problem in Lagrangean Relaxation, since the complicating constraints are removed from the constraint set and added to the objective function as a penalty term. To overcome this drawback, ^{Guignard and Kim (1987)} proposed Lagrangean Decomposition, which is a special case of the Lagrangean Relaxation. In Lagrangean Decomposition, instead of relaxing the complicating constraints, the set of variables that connects the set of complicating constraint to the other set of constraints is duplicated, generating new equality constraints. Then, the Lagrangean Relaxation is applied to the new problem by dualizing the new set of complicating variables, namely the new equality constraints. Now every constraint in the original problem appears in one of the sub-problems. Figure 2 shows the process to obtain the Lagrangean Decomposition for an instance of problem (P), where the set X is the positive integers. It can also be proved that the upper bound generated from the Lagrangean Decomposition is always at least as tight as the one from the Lagrangean relaxation (^{Guignard and Kim, 1987}). The bound of Lagrangean Decomposition can be further tightened by adding surrogate constraints, resulting in subproblems with overlapping constraints.

Another way of looking at Lagrangean decomposition is as a way to reach consensus in a disagreement between one or more parts. The subproblems share a connection expressed in terms of the linking constraints. When dualizing these constraints, the subproblems do not necessarily agree on the value of the shared variables. To resolve this disagreement they ask a third party to put a price on their disagreement. This price is the Lagrangean multiplier. The iterative procedure proposes new prices at each iteration until the parts agree as much as possible. This way of solving the disagreements makes Lagrangean decomposition the algorithm of choice for decentralized decision-making (^{Kelly and Zyngier, 2008}; ^{Boyd et al., 2010}; ^{Granada et al., 2012}).

Problem description

To clarify the presentation of the algorithms described in the following sections, it is convenient to define a motivating problem. The problem is later used as a case study in Section 6. Consider a set of production facilities, i∈ I, serving a number of customers, j∈ J as shown in Fig. 3. A maximum demand forecast, D_jpt, is given at each customer for every product, p∈ P, at every time period, t∈ T. The length of each time period, not necessarily uniform, is L_t. The objective is to determine the production and shipping plan to maximize the profit. The production plan is described by the production time q_ipt and the production amount x_ipt. Both variables are linked by a fixed production rate R_ip. The unit production cost per amount of product is γ_ip. When a product is produced in a given time period t, a setup time σ_ip with corresponding setup cost α_p is incurred. The binary variable indicating the occurrence of a production run is y_ipt. After production, the amount shipped to each customer, f_ijpt, is decided, with the corresponding transportation cost C_ij. The product that is not sent to customers is kept as inventory, s_ipt, with unit holding cost H_p. All the product arriving at a customer z_jpt is sold at price b_p.

The problem is formulated as a mixed-integer linear programming model (MILP), defined by Eqs. (6)–(13).

(6)-{\displaystyle \sum _i{\displaystyle \sum _j{\displaystyle \sum _p{\displaystyle \sum _t{{\displaystyle C}}_{ij}{{\displaystyle f}}_{ijpt}}}}},-{\displaystyle \sum _i{\displaystyle \sum _p{\displaystyle \sum _t({{\displaystyle \alpha }}_p{{\displaystyle y}}_{ipt}+{{\displaystyle \gamma }}_{ipt}{{\displaystyle x}}_{ipt}+{{\displaystyle H}}_p{{\displaystyle s}}_{ipt})}}}\mathrm{M}:\mathrm{M}\mathrm{a}\mathrm{x}{\displaystyle \sum _j^{}{\displaystyle \sum _p^{}{\displaystyle \sum _t^{}{{\displaystyle \beta }}_p{{\displaystyle z}}_{jpt}}}}

(7)s.t.R_{ip}{\theta }_{ipt}=x_{ipt}\forall i, p, t,

(8){\displaystyle \sum _p({\theta }_{ipt}+{\sigma }_{ip}{y}_{ipt})}\le L_t\forall i, t,

(9){\theta }_{ipt}\le y_{ipt}{L}_t\forall i, p, t,

(10)s_{ipt}=s_{ipt-\mathrm{1}}+x_{ipt}-{\displaystyle \sum _j{f}_{ijpt}}\forall i, p, t,

(11){\displaystyle \sum _i{f}_{ijpt}}=z_{jpt}\forall j, p, t,

(12)z_{jpt}\le D_{jpt}\forall j, p, t,

(13)x_{ipt}, {\theta }_{ipt}, s_{ipt}, f_{ijpt}, z_{jpt}\ge \mathrm{0}, y_{ipt}\in \left\{\mathrm{0}, \mathrm{1}\right\}.

The profit to be maximized is given by Eq. (6), and it is equal to the difference between the income and the costs, which includes setup costs, production costs, inventory holding costs, and transportation costs. The production time and production amount are related by Eq. (7).Constraint (8) ensures the production time limit is not exceeded, while constraint (9) ensures that the binary variable indicating when an item is produced is set to one if there is production time consumed. Equation (10) represents the system inventory balance.

Finally, Eqs. (11) and (12) define the amount sold and an upper bound for sales, respectively.

Lagrangean decomposition schemes

Temporal decomposition

Temporal decomposition is a well-established algorithm to solve large instances of problem (M). The goal of temporal decomposition is to break the problem into individual subproblems of one or more time periods. For example, a model with a planning horizon of 6 months can be decomposed into 6 one-month subproblems, or three subproblems of two months each. In the problem, time periods are linked by inventory balance constraints (Eq. (10)), which relate the inventory at the end of period t with the inventory at the end of period t−1. These constraints can be dualized to obtain the Lagrangean relaxation (^{Geoffrion, 1974}) of the problem, i.e. transferring the constraint to the objective function with a penalty term referred as Lagrange multiplier. Figure 4 is a graphical representation of the inventory balance constraint.

To obtain the Lagrangean decomposition (^{Guignard and Kim, 1987}) of the problem, the model is reformulated to include auxiliary variables si_ipt for the initial inventory at period t, and sf_ipt for the inventory at the end of the period. In this case, the linking constraint is the inventory continuity equation (Eq. (15)) enforcing that the initial inventory at period t must be equal to the inventory at the end of period t−1. The problem is decomposed by dualizing the continuity constraints. Figure 5 illustrates the formulation for this decomposition.

(14)s{f}_{ipt}=s{i}_{ipt}+x_{ipt}-{\displaystyle \sum _j{f}_{ijpt}}\forall i, p, t,

(15)s{f}_{ipt-\mathrm{1}}=s{i}_{ipt}\forall t>\mathrm{1.}

Applying this reformulation it is possible to obtain the temporal Lagrange dual (TLD(t)) of problem (M).

TLD\left(t\right)：\mathrm{M}\mathrm{a}\mathrm{x}{\displaystyle \sum _j^{}{\displaystyle \sum _p^{}{\displaystyle {{\displaystyle \beta }}_p{{\displaystyle z}}_{jpt}}}}

{\displaystyle {\displaystyle {\displaystyle -{\displaystyle \sum _i{\displaystyle \sum _p{\displaystyle ({{\displaystyle \alpha }}_p{{\displaystyle y}}_{ipt}+{{\displaystyle \gamma }}_{ipt}{{\displaystyle x}}_{ipt}+{{\displaystyle H}}_p)}}}}}}

(16)-{\displaystyle \sum _i{\displaystyle \sum _j{\displaystyle \sum _p{C}_{ij}{f}_{ijpt}}}+{\displaystyle \sum _i{\displaystyle \sum _p({\lambda }_{t}s{f}_{ipt}-{\lambda }_{t-\mathrm{1}}s{i}_{ipt})}}},

(17)\begin{array}{ll}s.t.\hfill & R_{ip}{\theta }_{ipt}=x_{ipt}\hfill \end{array}\forall i, p,

(18){\displaystyle \sum _p({\theta }_{ipt}+{\sigma }_{ip}{y}_{ipt})}\le L_t\forall i,

(19){\theta }_{ipt}\le y_{ipt}{L}_t\forall i, p,

(20)s{f}_{ipt}=s{i}_{ipt}+x_{ipt}-{\displaystyle \sum _j{f}_{ijpt}}\forall i, p,

(21){\displaystyle \sum _i{f}_{ijpt}}=z_{jpt}\forall j, p,

(22)z_{jpt}\le D_{jpt}\forall j, p,

(23)x_{ipt}, {\theta }_{ipt}, f_{ijpt}, z_{jpt}\ge \mathrm{0}, y_{ipt}\in \left\{\mathrm{0}, \mathrm{1}\right\}.

The inventory available at the beginning of the planning horizon is the initial inventory. In other words, for t = 1, si_ip1 = InitInv_ip.

Product decomposition

Temporal decomposition exploits the natural dynamic structure of time periods, one time period is linked to the next by the inventory balance. On the other hand, products do not exhibit the same structure, products usually share capacity of a limited resource, e.g.: Production, transportation, and inventory. The underlying dynamic structure is not clear at first sight, but it can be exposed through reformulation.

In general, products are linked by a knapsack constraint of the following form:

(24){\displaystyle \sum _p{\eta }_p{x}_p}\le C,

where x_p is the amount of product p, and h_p is the unit consumption factor of x_p with respect to resource C. In model (M), Eq. (8) has this structure, in which the production time is limited by the period length. To reformulate this constraint, it is convenient to think of the planning process as making optimal decisions about one product at a time. For the first product, the full capacity is available. After the first product has been optimized, it is required to know the remaining capacity for the next product. If we denote as c_p the available capacity for product p, then c_p is given by Eq. (25).

(25)c_p=c_{p-\mathrm{1}}-u_{p-\mathrm{1}}\forall p>\mathrm{1},

(26)u_p={\eta }_p{x}_p\forall p.

To make this reformulation possible, we have also assigned an arbitrary order to products. Later in the paper, we explore the impact of the order of the products in the computational efficiency.

As with temporal decomposition, it is convenient to define auxiliary variables for the initial and ending capacity for a given product, ci_p and cf_p, the reformulation obtained is shown in Fig. 6, and the constraints are shown in Eqs. (27) and (28).

(27)c{f}_p=c{i}_p-u_p\forall p,

(28)c{f}_{p-\mathrm{1}}=c{i}_p\forall p>\mathrm{1},

(29)c{i}_{\mathrm{1}}=C.

In problem (M), Eq. (8) ensures that the sum of the production time of all products does not exceed the available production time. Even though it is expressed in time units, the available production time limits the production capacity. It is possible to decompose the problem using the same procedure described in the previous section for temporal decomposition. To accomplish this, auxiliary variables hi_ip and hf_ip are defined to represent the number of production hours available at plant i before and after optimizing a given product. The reformulated constraints are Eqs. (30)–(33).

(30){\theta }_{ipt}+{\sigma }_{ip}{y}_{ipt}\le h{i}_{ipt}\forall ipt,

(31)h{f}_{ipt}=h{i}_{ipt}-({\theta }_{ipt}+{\sigma }_{ip}{y}_{ipt})\forall ipt,

(32)h{f}_{i\left(p-\mathrm{1}\right)t}=h{i}_{ipt}\forall i, p>\mathrm{1},

(33)h{i}_{i\mathrm{1}}=L_{\mathrm{1}}.

The resulting product Lagrange dual (PLD(p)) is given by Eqs. (34)–(43).

PLD\left(p\right):\mathrm{M}\mathrm{a}\mathrm{x}{\displaystyle \sum _j{\displaystyle \sum _t{\beta }_p{z}_{jpt}}}

(34)-{\displaystyle \sum _i{\displaystyle \sum _j{\displaystyle \sum _t{C}_{ij}{f}_{ijpt}}}}+{\displaystyle \sum _i{\displaystyle \sum _t({\mu }_{ipt}h{f}_{ipt}-{\mu }_{i\left(p-\mathrm{1}\right)t}h{i}_{ipt})}},-{\displaystyle \sum _i{\displaystyle \sum _t({\alpha }_p{y}_{ipt}+{\gamma }_{ipt}{x}_{ipt}+H_p{s}_{ipt})}}

(35)\begin{array}{rr}\hfill s.t.& \hfill R_{ip}{\theta }_{ipt}=x_{ipt}\end{array}\forall i, t,

(36){\theta }_{ipt}\le y_{ipt} L_t\forall i, t,

(37){\theta }_{ipt}+{\sigma }_{ip}{y}_{ipt}\le h{i}_{ipt}\forall it,

(38)h{f}_{ipt}=h{i}_{ipt}-({\theta }_{ipt}+{\sigma }_{ip}{y}_{ipt})\forall it,

(39)s_{ipt}=s_{ipt-\mathrm{1}}+x_{ipt}-{\displaystyle \sum _j{f}_{ijpt}}\forall it,

(40){\displaystyle \sum _i{f}_{ijpt}}=z_{jpt}\forall jt,

(41)z_{jpt}\le D_{jpt}\forall jt,

(42)x_{ipt}, {\theta }_{ipt}, s_{ipt}, f_{ijpt}, z_{jpt}\ge \mathrm{0}, y_{ipt}\in \left\{\mathrm{0}, \mathrm{1}\right\}.

The derivation from this section has been motivated by the large number of products considered in supply chain planning. However, the reformulation can be applied to any set sharing a resource in a knapsack constraint. For example, in the optimization of HVAC systems (Heating, Ventilation and AirConditioning), the capacity of the cooling unit is shared among all rooms (^{Alvarez et al., 2013}).

Simultaneous product and temporal decomposition

After decomposing the problem (M) by time periods the resulting subproblems have multiple products. Similarly, when decomposing problem (M) by products, the resulting subproblems are multiperiod. Applying the same concepts reviewed in the previous sections, it is possible to decompose problem (M) by products and time periods simultaneously. The network structure of the resulting subproblems is shown in Fig. 7.

In Fig. 7, node n₁ corresponds to the subproblem for the first product in the first period. In the structure, some nodes such as node n₅ have multiple inputs and multiple outputs. The problem is decomposed dualizing all the linking constraints, which are represented by the arcs in the diagram.

The resulting Lagrange dual is given by problem (PTLD(p,t)).

(43)PTLD(p,t):\mathrm{M}\mathrm{a}\mathrm{x}{\displaystyle \sum _j{\beta }_p{z}_{jpt}}-{\displaystyle \sum _i({\alpha }_p{y}_{ipt}+{\gamma }_{ipt}{x}_{ipt}+H_p{s}_{ipt})}-{\displaystyle \sum _i{\displaystyle \sum _j{C}_{ij}{f}_{ijpt}}}+{\displaystyle \sum _i({\mu }_{ipt}h{f}_{ipt}-{\mu }_{i(p-\mathrm{1})t}h{i}_{ipt})}+{\displaystyle \sum _i({\lambda }_{ipt}s{f}_{ipt}-{\lambda }_{ipt-\mathrm{1}}s{i}_{ipt})},

(44)s.t.R_{ip}{\theta }_{ipt}=x_{ipt}\forall i,

(45){\theta }_{ipt}\le y_{ipt}{L}_t\forall i,

(46){\theta }_{ipt}+{\sigma }_{ip}{y}_{ipt}\le h{i}_{ipt}\forall i,

(47)h{f}_{ipt}=h{i}_{ipt}-({\theta }_{ipt}+{\sigma }_{ip}{y}_{ipt})\forall i,

(48)s{f}_{ipt}=s{i}_{ipt}+x_{ipt}-{\displaystyle \sum _j{f}_{ijpt}}\forall i,

(49){\displaystyle \sum _i{f}_{ijpt}}=z_{jpt}\forall j,

(50)z_{jpt}\le D_{jpt}\forall j,

(51)x_{ipt}, {\theta }_{ipt}, s{i}_{ipt}, s{f}_{ipt}, h{i}_{ipt}, h{f}_{ipt}, f_{ijpt}, z_{jpt}\ge \mathrm{0}, y_{ipt}\in \{\mathrm{0}, \mathrm{1}\}.

The ability to simultaneously decompose a problem into both time periods and products gives the flexibility to decide on the number of subproblems. For example, if a planning problem has only 12 time periods, the maximum number of subproblems obtained when applying temporal decomposition is 12. However, if the problem has also 50 products, the number of subproblems can be up to 600. This can be helpful for very large scale problems.

An improved subgradient method

The subgradient method has been the most popular choice to update the multipliers in Lagrangean decomposition. The update formula proposed by ^{Fisher (1985)} has been used in a large number of applications because of its simplicity and ease of computation.

The main drawback of the formula is that its application does not guarantee a descent in the Lagrange dual (in case of maximization). Even though the subgradient gives a descent direction, the step size used might be too large. This could lead a decomposition algorithm diverging from the optimal Lagrange dual for some iterations. Several more iterations are required to bring the bound value back to the initial values. In the traditional subgradient, the step is taken regardless of the reduction in the Lagrange dual, and if the function fails to decrease in a specified number of iterations, the step size is decreased. Usually, the initial step is chosen to be between 0 and 2, and it is halved when the Lagrange dual fails to decrease in 3 consecutive iterations.

We propose a simple algorithm to avoid taking ascending steps in the search process. The idea is to probe the value of the candidate step and then decide the length of the step based on the difference between the current value of the Lagrange dual (z_k) and its value at the candidate step (z_c). We can distinguish the following cases:

1) The candidate step improves the bound (z_k-z_c>\delta). Take full step α.

2) The candidate step is similar to the current solution\left(\left|z_k-z_c\right|\le \delta \right) . Take half step, {\displaystyle \frac{\alpha }{\mathrm{2}}} .

3) The candidate step worsens the bound\left(z_c-z_k>\delta \right) . Take a quarter step, {\displaystyle \frac{\alpha }{\mathrm{4}}} .

The decisions on the step size rely on the piecewise-convex nature of the Lagrange dual function. Clearly, if the value of the function at the candidate step is similar to the current value, there is a point between the current and the candidate solution that minimizes the Lagrange dual function in the current search direction. Following the same idea, if the bound worsens at the candidate step the point that minimizes the Lagrange dual function in the current search direction is located closer to the current point. Figure 8 summarizes the cases. The method is computationally more intensive than the traditional subgradient because more Lagrange dual evaluations are required. Each evaluation is conducted by solving all the subproblems for a given value of the Lagrange multipliers. To partially alleviate this issue, when the case 1 is found, i.e. the candidate step improves the bound, the value of z_k is immediately updated to save one evaluation. The extra effort in the computation is compensated by having optimization trajectories with less divergence, which can be especially useful when the time budget for the iterations is low, for example, in problems whose subproblems are hard to solve.

Case study

The computational performance of the Lagrangean decomposition algorithms proposed is assessed by applying them to the case study presented in Section 3. The considered problem involves 6 manufacturing sites and 10 markets, while the number of products and time periods was varied from 6 to 20. The number of time periods considered is equal to the number of products to ensure a fair comparison between temporal and product decomposition. In the entire section, this number will be referred as the size of the problem. Each time period represents 720 h. The MILP formulations are modeled in JuMP (^{Dunning et al., 2017}) and solved with Gurobi 7.5 (^{Gurobi Optimization, 2015}) in an Intel i7 machine with 16 Gb of RAM. Model sizes are shown in Table 1. Since the intention is to evaluate the decomposition schemes, the optimal solution was provided as a lower bound instead of using a Lagrangean heuristic. The method used to update the Lagrangean multipliers is the standard subgradient with the initial step size of 2, decreasing the step to half of the current step when no improvement is obtained in 3 consecutive iterations.

Decomposition schemes comparison

The different ways of decomposing the problem, by product (P), time periods (T), and both, products and time periods (PT), were compared applying the algorithms to the case study. The effect of the initial multipliers in the formulation was also studied. Two possible values for the initial multipliers were considered: 1) zero (l₀ = 0), and 2) the multipliers obtained from solving the LP relaxation (l₀ from LP). The results for solution time and optimality gap after 90 s are shown in Table 2.

In the case study, the temporal decomposition converged to better bounds than the ones obtained with product decomposition. Even though the number of products and time periods was set to be the same, the product subproblems took longer to solve than the time period subproblems, which can be seen in the larger average iteration time. Allowing more iterations, showed that the minimum value for product decomposition in a problem of size 6 was 0.23%, which is still larger than the gap obtained with temporal decomposition. In the case study, the termination criterion was chosen to be a time limit, because when optimizing supply chain planning there is usually a limited time to obtain a solution, in order to allow the analysis of different scenarios. The temporal decomposition being more efficient than the product decomposition for the case study cannot be generalized. The value of the Lagrangean relaxation for products decomposition strongly depends on the capacity. In the limiting case, when there is unlimited production capacity the problem for each product can be solved independently. Figure 9 shows the sensitivity analysis of the gap with respect to capacity for product decomposition of size 20. For large values of capacities, the gap decreases to reach very low gaps. The values do not reach zero gap because of the relative gap setting for solving the subproblems was set at 0.5%.

The sensitivity analysis explicitly shows that the gap obtained using product decomposition depends on the parameters of the model. For large capacity values, the gaps obtained are even lower than the gap obtained with temporal decomposition (0.78%). For small values of capacity, the gap is also small, because the problem becomes too constrained and finding the optimal solution becomes easier.

Another noteworthy aspect is that the number of products was set to be the same as the number of time periods for comparison purposes. However, the number of products in supply chain planning models usually exceeds the number of time periods. In this situation, product decomposition leads to a larger number of smaller subproblems than temporal decomposition, making it better suited to take advantage of parallel computing capabilities. The simultaneous product and temporal decomposition for 6 products and time periods is an example of such situation; the number of subproblems obtained is much larger (O (n²) where n is the problem size), but with a much lower average iteration time. The tradeoff is that a larger number of subproblems comes with a larger number of Lagrangean multipliers that need to be optimized, which could slow down the overall solution time. For the problem with 20 products and time periods, the average iteration time for simultaneous product and temporal decomposition becomes similar to the average iteration time for temporal decomposition. A balance is reached in which solving 400 small subproblems takes almost the same time as solving 20 larger subproblems. This indicates that the number of subproblems must be carefully chosen to balance the average iteration time and the number of iterations required to converge to the optimal solution.

The success of Lagrangean decomposition is based on the exponential increase of the solution time with problem size. For example, if a problem is decomposed into two subproblems, the solution time required to solve both subproblems is much lower than the time required to solve the full problem. However, several iterations are required and subproblems are solved one time per iteration. If the problem is decomposed into more subproblems, more iterations are required to converge. The number of subproblems that minimizes the solution time is problem-dependent. The benefit of the decomposition schemes proposed is that it provides more options to the modeler with the goal of finding the most efficient algorithm to solve a given problem.

The choice of the initial value of the Lagrangean multipliers is also an important factor impacting the efficiency of the algorithm. In all the cases, but the product decomposition of size 6, choosing the initial multipliers as the multipliers obtained from the LP relaxation allowed to reach a smaller gap. Because the value of the initial Lagrange dual is smaller for the initial multipliers from the LP relaxation, the gap that needs to be closed is also smaller. Figure 10 compares the Lagrange dual (ZLD) with respect to time for the temporal decomposition of size 20.

Aggregation of linking constraints

When solving a Lagrangean decomposition problem, it is necessary to find a bound, upper bound in the case of maximization and lower bound in the case of minimization that is as tight as the LP relaxation of the problem. Thus, if the decomposition is applied to a relaxation of the original problem, the result of the decomposition algorithm also leads to a valid bound for the original problem. Clearly, the quality of the bound obtained depends on how close/near is the solution of the relaxation chosen to the optimal solution. As mentioned before, a larger number of multipliers that need to be optimized translates into a larger number of iterations to converge to the optimal Lagrange dual. A relaxation of the problem that can help to reduce the number of multipliers can be obtained by creating surrogate constraints from the linking constraints, i.e. aggregating the linking constraints.

Take, for example, the linking constraint for temporal decomposition (Eq. (15)). The constraint is indexed by the manufacturing site i and the product p. The equation actually represents ip linking constraints between node t and node t + 1. Surrogates can be constructed by adding the constraints by either site, product, or both. Figure 11 provides an example of different aggregation possibilities for temporal decomposition.

The different aggregation strategies were applied to the case study, obtaining the results shown in Table 3. In the table none represents the original formulation without aggregation, i represents the constraints that were added by manufacturing site, t by time period, p by product, it by site and time period, and ip by site and product.

The results indicate that the aggregation of linking constraints does not improve the Lagrange dual gap. In the best case, for temporal decomposition of size 6, the results remain similar while decreasing the average iteration time, as for the case of aggregation by sites. For temporal decomposition of size 20, aggregation of products significantly worsens the Lagrange dual. For product decomposition, aggregations by site produce the largest increases in the Lagrange dual. For simultaneous product and temporal decomposition, the aggregation reduces the number of multipliers, but the Lagrange dual also worsens. The observation that the aggregation of linking constraints can lead to the same bound reducing the number of multipliers is reported by ^{Chen and Guignard (1998)}. Consistently, in the case study, we obtained some scenarios where the bound is maintained and somewhere it is worsened.

Effect of product order

To obtain the product decomposition formulation, the products are assigned an arbitrary order. This gives rise to the question about the impact of the order of the products in the results. To answer this question, product Lagrangean decomposition was applied to 100 random orders of products, for the problem of size 6 after 50 iterations. The histogram of the results is shown in Fig. 12. The obtained results have a low dispersion, and there are not clear outliers. This indicates the order of the products does not have a significant impact on the solution gap.

Probing subgradient performance

The efficiency of the improved subgradient method is assessed by comparison of its application to the traditional subgradient and to the three decomposition schemes analyzed: Temporal decomposition, product decomposition, and simultaneous product and temporal decomposition. The results of the first 10 iterations for problems of size 6 is shown in Fig. 13. A direct comparison can be made looking closely at the results of iteration 2. For the case of temporal decomposition, the candidate step resulted in a large increase in the Lagrange dual. By taking a quarter step, it was possible to find a point where the Lagrange dual decreased. For product decomposition, the candidate step also produces an increase in the Lagrange dual. However, in this case, the quarter step also did not result in a decrease of the Lagrange dual. Nevertheless, the effect of picking an ascending value was mitigated.

Both the temporal decomposition, and the product-temporal decomposition show that the traditional subgradient can recover from diverging, but it can take several iterations to catch up with the bounds obtained using the probing subgradient. In the product decomposition, the difference is even clearer. The Lagrange duals obtained after 10 iterations are completely different. These two results combined indicate that the probing subgradient is a good choice in presence of a low budget of iterations.

Conclusions

Product decomposition is a novel decomposition scheme that provides additional options for modelers to find the most efficient optimization algorithms. The idea of decomposing was motivated by the large number of products in some supply chain planning problems. However, the derivation can be applied to other problems with sets sharing resource in knapsack-like constraints. The results are also useful to design algorithms similar to the rolling horizon, in which the problem is solved one product at a time. The results of temporal decomposition applied to case study resulted better than the results of product decomposition. However, the sensitivity analysis indicates that product decomposition can lead to a smaller gap when the capacity is larger, demonstrating that the results depend on the parameters of the model. Furthermore, both approaches are mathematically equivalent. The only difference is the effect of the relaxation obtained from the Lagrange dual in the objective function. This effect is problem-specific and depends on the modeling parameters. The choice between product or temporal decomposition depends on the number of elements on each set (products and time periods) and the structure of the objective function. Modelers are encouraged to experiment with both approaches when facing a new supply chain planning problem. In most applications, the number of products is larger than the number of time periods, making product decomposition better suited to take advantage of parallelization.

The results of the simultaneous product and temporal decomposition show the tradeoff between number of subproblems (or multipliers), and the difficulty to converge to the optimal Lagrange dual. A larger number of subproblems requires a larger number of multipliers, with a lower average iteration time, and a resulting algorithm with a lower probability to converge to the optimal Lagrange dual. The simultaneous decomposition is recommended only for cases in which decomposing only by time periods or products still results in subproblems that are too computational expensive to solve, with large average iteration times. Aggregation of linking constraints should also be considered in a case by case basis, as it can be beneficial in some applications, but not in others. The results indicate that it can lead to reductions in the number of multipliers while maintaining the quality of the bounds. Together with the decomposition schemes, it represents another lever for modelers to design more efficient algorithms. On the other hand, the order of the products did not have much impact on the results.

Finally, the probing subgradient method resulted in a more stable Lagrange dual method, preventing the Lagrange dual to take large ascending steps to make the algorithm diverge in the first iterations. The method requires the solution of more optimization problems per iteration. Therefore, it should only be used for expensive subproblems, i.e. a large average iteration time.

In summary, product decomposition is a novel approach that provides modelers with more options to design the most efficient decomposition algorithms and obtain better results.

References

[1] Alvarez J, Redondo J, Camponogara E, Normey-Rico J, Berenguel M, Ortigosa P (2013). Optimizing building comfort temperature regulation via model predictive control. Energy and Building, 57: 361–372

[2] Boyd S, Parikh N, Chu E, Peleato B, Eckstein J (2010). Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends in Machine Learning, 3(1): 1–122

[3] Calfa B A, Agarwal A, Grossmann I E, Wassick J M (2013). Hybrid Bilevel-Lagrangean decomposition scheme for the integration of planning and scheduling of a network of batch plants. Industrial & Engineering Chemistry Research, 52(5): 2152–2167

[4] Chen B, Guignard M (1998). Polyhedral analysis and decompositions for capacitated plant locationtype problems. Discrete Applied Mathematics, 82(1–3): 79–91

[5] Dunning I, Huchette J, Lubin M (2017). JuMP: A modeling language for mathematical optimization. SIAM Review, 59(2): 295–320

[6] Erlenkotter D (1978). A dual-based procedure for uncapacitated facility location. Operations Research, 26(6): 992–1009

[7] Fisher M L (1985). An applications oriented guide to Lagrangian relaxation. Interfaces, 15(2): 10–21

[8] Geoffrion A M (1974). Lagrangean relaxation for integer programming. In: Balinski M, eds. Mathematical Programming Studies. Berlin: Springer

[9] Granada M, Rider M J, Mantovani J, Shahidehpour M (2012). A decentralized approach for optimal reactive power dispatch using a Lagrangian decomposition method. Electric Power Systems Research, 89: 148–156

[10] Graves S C (1982). Using Lagrangean techniques to solve hierarchical production planning problems. Management Science, 28(3): 260–275

[11] Guignard M (2003). Lagrangean relaxation. Top (Madrid), 11(2): 151–200

[12] Guignard M, Kim S (1987). Lagrangean decomposition: A model yielding stronger Lagrangean bounds. Mathematical Programming, 39(2): 215–228

[13] Gupta A, Maranas C D (1999). A hierarchical Lagrangean relaxation procedure for solving midterm planning problems. Industrial & Engineering Chemistry Research, 38(5): 1937–1947

[14] Gurobi Optimization (2015). Gurobi optimizer reference manual.

[15] Held M, Karp R M (1970). The traveling-salesman problem and minimum spanning trees. Operations Research, 18(6): 1138–1162

[16] Jackson J R, Grossmann I E (2003). Temporal decomposition scheme for nonlinear multisite production planning and distribution models. Industrial & Engineering Chemistry Research, 42(13): 3045–3055

[17] Kearney A T (2004). The complexity challenge: A survey on complexity management across the supply chain. AT Kearney Publications,

[18] Kelley J E Jr (1960). The cutting-plane method for solving convex programs. Journal of the Society for Industrial and Applied Mathematics, 8(4): 703–712

[19] Kelly J D, Zyngier D (2008). Hierarchical decomposition heuristic for scheduling: Coordinated reasoning for decentralized and distributed decision-making problems. Computers & Chemical Engineering, 32(11): 2684–2705

[20] Lemarechal C, Strodiot J J, Bihain A (1981). On a bundle algorithm for nonsmooth optimization. In: Managasarian O L, Meyer R R, Robinson S M, eds. Nonlinear programming 4. New York: Academic Press

[21] Mouret S, Grossmann I E, Pestiaux P (2011). A new Lagrangian decomposition approach applied to the integration of refinery planning and crude-oil scheduling. Computers & Chemical Engineering, 35(12): 2750–2766

[22] Oliveira F, Gupta V, Hamacher S, Grossmann I E (2013). A Lagrangean decomposition approach for oil supply chain investment planning under uncertainty with risk considerations. Computers & Chemical Engineering, 50: 184–195

[23] Terrazas-Moreno S, Trotter P A, Grossmann I E (2011). Temporal and spatial Lagrangean decompositions in multi-site, multi-period production planning problems with sequence-dependent changeovers. Computers & Chemical Engineering, 35(12): 2913–2928

[24] van Elzakker M, Zondervan E, Raikar N, Hoogland H, Grossmann I E (2014). An SKU decomposition algorithm for the tactical planning in the FMCG industry. Computers & Chemical Engineering, 62: 80–95

[25] Wang S H (2003). An improved stepsize of the subgradient algorithm for solving the Lagrangian relaxation problem. Computers & Electrical Engineering, 29(1): 245–249

[26] Wu D, Ierapetritou M (2006). Lagrangean decomposition using an improved Nelder–Mead approach for lagrangean multiplier update. Computers & Chemical Engineering, 30(5): 778–789

[27] Rodriguez M A, Harjunkoski I, Grossmann I E (2014). Optimal supply chain design and management over a multi-period horizon under demand uncertainty. Part II: A Lagrangean decomposition algorithm. Computers & Chemical Engineering, 62: 211–224

RIGHTS & PERMISSIONS

The Author(s) 2018. Published by Higher Education Press. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0)