Benders’ Decomposition Methods for Structured Optimization, including Stochastic Optimization 1 1 Block Ladder Structure We consider the solution of a linear optimization model of the basic format: cT x + minimizex,y s.t. fT y Ax = Bx + x≥0 b Dy = d y≥0. Here the variables x can be thought of as stage-1 decisions, governed by constraints Ax = b, x ≥ 0 and with cost function cT x. Once the x variables have been chosen, the y variables are chosen next, subject to the constraints Dy = d − Bx, y ≥ 0 and with cost function f T y. We also consider the more complex format associated with two-stage optimization under uncertainty: minimize x, y1 , . . . , yK s.t. cT x + α1 f1T y1 + α2 f2T y2 + · · · Ty + αK fK K Ax B1 x + B2 x = b = d1 D 1 y1 + ... D 2 y2 = d2 .. BK x .. . . + D K yK = dK x, y1 , y2 , . . . , yK ≥ 0 . Actually, each of these formats can be thought of as a special case of the other format. The more complex format above is known as “block-ladder”, and is illustrated in Figure 1. 2 Stage-2 Variables Objectives Stage-1 Variables RHS State 1 State 2 State 3 State k Figure 1: Block ladder structure of two-stage stochastic linear optimization. 2 Reformulation of the Basic Model The basic model format is: cT x + VAL = minimumx,y s.t. fT y Ax Bx + x≥0 3 = b Dy = d y≥0, whose optimal objective value is denoted by VAL. We re-write this as: VAL = minimumx cT x + z(x) s.t. Ax = b x≥0, where: P2 : z(x) = minimumy fT y s.t. Dy = d − Bx y≥0. We call this problem “P2” because it represents the stage-2 decision problem, once the stage-1 variables x have been chosen. Applying duality, we can also construct z(x) through the dual of P2, which we denote by D2: D2 : z(x) = maximump pT (d − Bx) s.t. DT p ≤ f . The feasible region of D2 is the set: � � D2 := p | DT p ≤ f , whose extreme points and extreme rays can be enumerated: p1 , . . . , pI are the extreme points of P2 , and 4 r1 , . . . , rJ are the extreme rays of P2 . If we use an algorithm to solve D2, exactly one of two cases will arise: D2 is unbounded from above, or D2 has an optimal solution. In the ﬁrst case, the algorithm will ﬁnd that D2 is unbounded from above, and it will return one of the extreme rays r¯ = rj for some j, with the property that (rj )T (d − Bx) > 0 in which case z(x) = +∞ . In the second case, the algorithm will solve D2, and it will return one of the extreme points p̄ = pi for some i as well as the optimal objective function value z(x) which will satisfy: z(x) = (pi )T (d − Bx) = max (pk )T (d − Bx) . k=1,...,I Therefore we can re-write D2 again as the following problem: D2 : z(x) = minimumz z s.t. (pi )T (d − Bx) ≤ z (rj )T (d − Bx) ≤ 0 i = 1, . . . , I j = 1, . . . , J . This problem simply states that the value of D2 is just the dual objective function evaluated at the best extreme point of the dual feasible region, provided that the objective function makes an obtuse angle with every extreme ray of the dual feasible region. 5 We now take this reformulation of z(x) and place it in the orginal problem: FMP : cT x + z VAL = minimumx,z s.t. Ax = b x≥0 (pi )T (d − Bx) ≤ z i = 1, . . . , I (rj )T (d − Bx) ≤ 0 j = 1, . . . , J . We call this problem FMP for the full master problem. Comparing FMP to the original version of the problem, we see that we have eliminated the variables y from the problem, we have added a single scalar variable z, and we have also added a generically extremely large number of constraints. 3 Delayed Constraint Generation The idea behind delayed constraint generation is to try to solve FMP using only a small subset of the constraints, and to check afterward whether any of the non-included constraints are violated. Consider the restricted master problem RMP composed of only k of the extreme point / extreme ray constraints from FMP: RMPk : VALk cT x + z = minimumx,z s.t. Ax = b x≥0 (pi )T (d − Bx) ≤ z i = 1, . . . , k − l (rj )T (d − Bx) ≤ 0 j = 1, . . . , l . 6 We solve this problem, obtaining the optimal objective value VALk and ¯ z̄ that is optimal for this problem. We ﬁrst observe that the solution x, k VAL is a lower bound on VAL, i.e.: VALk ≤ VAL . In order to check whether the solution x̄, z̄ is optimal for the full master problem, we must check whether x, ¯ z̄ violates any of the non-included constraints. We check for this by solving the linear optimization problem: Q(¯ x) : = maximump pT (d − B x) ¯ = minimumy s.t. DT p ≤ f fT y s.t. Dy = d − B x ¯ y≥0. If Q(¯ x) is unbounded, the algorithm will return an extreme ray r¯ = rj for some j, where x̄ will satisfy: (rj )T (d − B x̄) > 0 . We therefore have found that x̄ has violated the constraint: (rj )T (d − Bx) ≤ 0 , and so we add this constraint to RMP and re-solve RMP. If Q(x̄) has an optimal solution, the algorithm will return an optimal extreme point p̄ = pi for some i, as well as the optimal solution ȳ to the minimization problem of Q(x̄). 7 Now notice ﬁrst from Q(¯ x) that (¯ x, y) ¯ is feasible for the original problem. Therefore, if UB is a previously computed upper bound on VAL, we can update the upper bound and our “best” feasible solution as follows: ¯ + f T ȳ ≤ UB, then BESTSOL ← (¯ x, y) ¯ . If cT x ¯ + f T ȳ} . UB ← min{UB, cT x Also, notice that if: ¯ > z¯ , (pi )T (d − B x) then we have found that x̄, z̄ has violated the constraint: (pi )T (d − Bx) ≤ z , and so we add this constraint to RMP and re-solve RMP. If (pi )T (d − B x̄) ≤ z̄, then we have: max (pk )T (d − B x) ¯ = (pi )T (d − B x) ¯ ≤ z¯ , k=1,...,I and so x, ¯ z̄ is feasible and hence optimal for FMP, and so (¯ x, y) ¯ is optimal for the original problem, and so we terminate the algorithm. By our construction of upper and lower bounds on VAL, we can also terminate the algorithm whenever: UB − VALk ≤ ε for a pre-speciﬁed tolerance ε. 8 3.1 Delayed Constraint Generation Algorithm Here is a formal description of the algorithm just described. 0. Set LB= −∞ and UB= +∞. 1. A typical iteration starts with the relaxed master problem RMPk , in which only k constraints of the full master problem FMP are included. An optimal solution x̄, z̄ to the relaxed master problem is computed. We update the lower bound on VAL: LB ← VALk . 2. Solve the subproblem Q(x̄): Q(¯ x) : maximump pT (d − B x) ¯ = minimumy s.t. DT p ≤ f fT y s.t. Dy = d − B x̄ y≥0. 3. If Q(x̄) has an optimal solution, let p̄ and ȳ be the primal and dual ¯ z̄ satisﬁes all constraints ¯ ≤ z̄, then x, solutions of Q(¯ x). If p¯T (d − B x) of FMP, and so x, ¯ z̄ is optimal for FMP. Furthermore, this implies that (¯ x, y) ¯ is optimal for the original problem. If p¯T (d − B x) ¯ > z̄, then we add the constraint: p̄T (d − Bx) ≤ z to the restricted master problem RMP. Update the upper bound and the best solution: 9 If cT x ¯ + f T ȳ ≤ UB, then BESTSOL ← (¯ x, ȳ) UB ← min{UB, cT x ¯ + f T ȳ} If UB − LB ≤ ε, then terminate. Otherwise, return to step 1. 4. If Q(¯ x) is unbounded from above, let r¯ be the extreme ray generated by the algorithm that solves Q(¯ x). Add the constraint: r̄T (d − Bx) ≤ 0 to the restricted master problem RMP, and return to step 1. This algorithm is known formally as Benders’ decomposition. The Benders’ decomposition method was developed in 1962 [2], and is described in many sources on large-scale optimization and stochastic programming. A general treatment of this method can be found in [3, 4]. The algorithm can be initialized by ﬁrst computing any x̄ that is feasible for the ﬁrst stage, that is, A¯ ¯ ≥ 0, and by choosing the value of z̄ x = b, x to be −∞ (or any suitably small number). This will force (¯ x, z) ¯ to violate some constraints of FMP. Notice that the linear optimization problems solved in Benders’ decomposition are “small” relative to the size of the original problem. We need to solve RMPk once every outer iteration and we need to solve Q(x̄) once every outer iteration. Previous optimal solutions of RMPk and Q(x̄) can serve as a starting basis for the current versions of RMPk and Q(x̄) if we use the simplex method. Because interior-point methods do not have the capability of starting from a previous optimal solution, we should think of Benders’ decomposition as a method to use in conjunction with the simplex method primarily. Incidentally, one can show that Benders’ decomposition method is the same as Dantzig-Wolfe decomposition applied to the dual problem. 10 4 Benders’ Decomposition for Problems with BlockLadder Structure Suppose that our problem has the complex block-ladder structure of the following optimization model that we see for two-stage optimization under uncertainty: VAL = minimize x, y1 , . . . , yK s.t. cT x + α1 f1T y1 + α2 f2T y2 + · · · Ty + αK fK K Ax = b B1 x + = d1 D 1 y1 B2 x + D 2 y2 .. . = d2 .. .. . . BK x + D K yK x, y1 , y2 , . . . , yK ≥ 0 . whose optimal objective value is denoted by VAL. We re-write this as: VAL = minimumx cT x + s.t. K � ω=1 αω zω (x) Ax x≥0, where: 11 = b = dK P2ω : zω (x) = minimumyω fωT yω s.t. Dω y ω = dω − Bω x yω ≥ 0 . We call this problem “P2ω ” because it represents the stage-2 decision problem under scenario ω, once the stage-1 variables x have been chosen. Applying duality, we can also construct zω (x) through the dual of P2ω , which we denote by D2ω : D2ω : zω (x) = maximumpω pTω (dω − Bω x) s.t. DωT pω ≤ fω . The feasible region of D2ω is the set: � � D2ω := pω | DωT pω ≤ fω , whose extreme points and extreme rays can be enumerated: pω 1 , . . . , pIωω are the extreme points of D2ω , and rω 1 , . . . , rωJω are the extreme rays of D2ω . If we use an algorithm to solve D2ω , exactly one of two cases will arise: D2ω is unbounded from above, or D2ω has an optimal solution. In the ﬁrst case, the algorithm will ﬁnd that D2ω is unbounded from above, and it will return one of the extreme rays r¯ω = rωj for some j, with the property that 12 (rωj )T (dω − Bω x) > 0 in which case zω (x) = +∞ . In the second case, the algorithm will solve D2ω , and it will return one of the extreme points p̄ω = piω for some i as well as the optimal objective function value zω (x) which will satisfy: zω (x) = (piω )T (dω − Bω x) = max (pkω )T (dω − Bω x) . k=1,...,Iω Therefore we can re-write D2ω again as the following problem: D2ω : zω (x) = minimumzω zω s.t. (piω )T (dω − Bω x) ≤ zω i = 1, . . . , Iω (rωj )T (dω − Bω x) ≤ j = 1, . . . , Jω . 0 This problem simply states that the value of D2ω is just the dual objective function evaluated at the best extreme point of the dual feasible region, provided that the objective function makes an obtuse angle with every extreme ray of the dual feasible region. We now take this reformulation of zω (x) and place it in the orginal problem: 13 cT x + FMP : VAL = minimum x, z1 , . . . , zK s.t. K � ω=1 αω zω Ax = b x≥0 (piω )T (dω − Bω x) ≤ zω (rωj )T (dω − Bω x) ≤ i = 1, . . . , Iω , ω = 1, . . . , K 0 j = 1, . . . , Jω , ω = 1, . . . , K . We call this problem FMP for the full master problem. Comparing FMP to the original version of the problem, we see that we have eliminated the vector variables y1 , . . . , yK from the problem, we have added K scalar variables zω for ω = 1, . . . , K, and we have also added a generically extremely large number of constraints. As in the basic model, we consider the restricted master problem RMP composed of only k of the extreme point / extreme ray constraints from FMP: RMPk : VALk = minimum x, z1 , . . . , zK s.t. cT x + K � ω=1 αω zω Ax = b x≥0 (piω )T (dω − Bω x) ≤ zω (rωj )T (dω − Bω x) ≤ for some i and ω 0 for some i and ω , where there are a total of k of the inequality constraints. We solve this problem, obtaining the optimal objective value VALk and the solution x, ¯ z¯1 , . . . , z̄K 14 that is optimal for this problem. We ﬁrst observe that VALk is a lower bound on VAL, i.e.: VALk ≤ VAL . In order to check whether the solution x, ¯ z¯1 , . . . , z̄K is optimal for the full master problem, we must check whether x, ¯ z¯1 , . . . , z̄K violates any of the non-included constraints. We check for this by solving the K linear optimization problems: x) : = maximumpω Qω (¯ pTω (dω − Bω x) ¯ = minimumyω s.t. DωT pω ≤ fω fωT yω s.t. Dω yω = dω − Bω x ¯ yω ≥ 0 . If Qω (¯ x) is unbounded, the algorithm will return an extreme ray r¯ω = rωj for some j, where x̄ will satisfy: ¯ >0. (rωj )T (dω − Bω x) We therefore have found that x̄ has violated the constraint: (rωj )T (dω − Bω x) ≤ 0 , and so we add this constraint to RMP. If Qω (x̄) has an optimal solution, the algorithm will return an optimal extreme point p¯ω = piω for some i, as well as the optimal solution y¯ω of the minimization problem of Qω (x̄). If (piω )T (dω − Bω x) ¯ > z̄ω , 15 then we have found that x, ¯ z¯ω has violated the constraint: (piω )T (dω − Bω x) ≤ zω , and so we add this constraint to RMP. If Qω (x̄) has ﬁnite optimal objective function values for all ω = 1, . . . , K, then x, ¯ y¯1 , . . . , ȳK satisﬁes all of the constraints of the original problem. Therefore, if UB is a previously computed upper bound on VAL, we can update the upper bound and our “best” feasible solution as follows: If cT x ¯+ K � αω fω T ȳω ≤ UB, then BESTSOL ← (¯ x, y¯1 , . . . , ȳK ) . ω=1 T ¯+ UB ← min{UB, c x K � αω fω T ȳω } . ω=1 ¯ ≤ z̄ω for all ω = 1, . . . , K, then we have: If (piω )T (dω − Bω x) ¯ = (piω )T (dω − Bω x) ¯ ≤ z̄ω max (plω )T (dω − Bω x) l=1,...,Iω for all ω = 1, . . . , K , ¯ z¯1 , . . . , z̄K and so x, ¯ z¯1 , . . . , z̄K satisﬁes all of the constraints in FMP. Therefore, x, is feasible and hence optimal for FMP, and so (¯ x, y¯1 , . . . , ȳK ) is optimal for the original problem, and we can terminate the algorithm. Also, by our construction of upper and lower bounds on VAL, we can also terminate the algorithm whenever: UB − VALk ≤ ε for a pre-speciﬁed tolerance ε. 16 4.1 The Delayed Constraint Generation Algorithm for Problems with Block-Ladder Structure Here is a formal description of the algorithm just described for problems with block-ladder structure. 0. Set LB= −∞ and UB= +∞. 1. A typical iteration starts with the relaxed master problem RMPk , in which only k constraints of the full master problem FMP are included. An optimal solution x, ¯ z¯1 , . . . , z̄K to the relaxed master problem is computed. We update the lower bound on VAL: LB ← VALk . 2. For ω = 1, . . . , K, solve the subproblem Qω (x̄): x) : = maximumpω Qω (¯ pTω (dω − Bω x) ¯ = minimumyω s.t. DωT pω ≤ fω fωT yω s.t. Dω yω = dω − Bω x ¯ yω ≥ 0 . – If Qω (x̄) is unbounded from above, let r̄ω be the extreme ray generated by the algorithm that solves Q(x̄). Add the constraint: r̄ωT (dω − Bω x) ≤ 0 to the restricted master problem RMP. x) has an optimal solution, let p¯ω and y¯ω be the primal and – If Qω (¯ x). If (piω )T (dω − Bω x) ¯ > z̄ω , then we add dual solutions of Qω (¯ the constraint: 17 p̄Tω (dω − Bω x) ≤ zω to the restricted master problem RMP. 3. If (piω )T (dω − Bω x) ¯ ≤ z̄ω for all ω = 1, . . . , K, then (¯ x, y¯1 , . . . , ȳK ) is optimal for the original problem, and the algorithm terminates. 4. If Qω (x̄) has an optimal solution for all ω = 1, . . . , K, then update the upper bound on VAL as follows: ¯+ If cT x K � αω fω T ȳω ≤ UB, then BESTSOL ← (¯ x, y¯1 , . . . , ȳK ) . ω=1 T UB ← min{UB, c x ¯+ K � αω fω T ȳω } . ω=1 If UB − LB ≤ ε, then terminate. Otherwise, add all of the new constraints to RMP and return to step 1. In this version of Benders’ decomposition, we might add as many as K new constraints per major iteration. 5 The Power Plant Investment Model Recall the power plant investment model described earlier. The full model is: 18 min x,y s.t. 4 � ci xi + αω 5 � 15 5 � � ω=1 i=1 4 � 125 � 0.01fi (ω)hj yijkω i=1 j=1 k=1 ci xi ≤ 10, 000 (Budget constraint) i=1 x4 ≤ 5.0 (Hydroelectric constraint) yijkω ≤ xi 5 � for i = 1, . . . , 4, all j, k, ω yijkω ≥ Djkω for all j, k, ω (Capacity constraints) (Demand constraints) i=1 x ≥ 0, y≥0 (1) where αω is the probability of scenario ω, and Djkω is the power demand in block j and year k under scenario ω. This model has 4 stage-1 variables, and 46,875 stage-2 variables. Furthermore, other than nonnegativity constraints, there are 2 stage-1 constraints, and 375 constraints for each of the 125 scenarios in stage-2. This yields a total of 46, 877 constraints. It is typical for problems of this type to be solved by Benders’ decomposition, which exploits the problem structure by decomposing it into smaller problems. Consider a ﬁxed x which is feasible (i.e., x ≥ 0 and satisﬁes the budget constraint and hydroelectric constraint). Then the optimal second-stage variables yijkω can be determined by solving for each ω the problem: zω (x) � min y subject to 5 � 15 5 � � 0.01fi (ω)hj yijkω i=1 j=1 k=1 yijkω ≤ xi 5 � for i = 1, . . . , 4, all j, k, yijkω ≥ Djkω for all j, k (Capacity constraints) (Demand constraints) i=1 y ≥ 0. (2) 19 Although this means solving 125 problems (one for each scenario), each subproblem has only 375 variables and 375 constraints. (Furthermore, the structure of this problem lends itself to a simple greedy optimal solution.) The dual of the subproblem (2) is: zω (x) � max p,q 5 � 15 4 � � i=1 j=1 k=1 xi pijkω + 15 5 � � Djkω qjkω j=1 k=1 subject to pijkω + qjkω ≤ 0.01fi (ω)hj i = 1, . . . , 4, for all j, k (3) pijkω ≤ 0 qjkω ≥ 0. 6 Computational Results All models were solved in OPLStudio on a Sony Viao Laptop with a Pentium III 750 MHz Processor Running Windows 2000. Table 1 shows the number of iterations taken by the simplex method to solve the original model. Simplex Original Model Iterations 41,334 CPT Time (minutes:seconds) 5:03 Table 1: Iterations to solve the original model. We next solved the problem using Benders’ decomposition, using a duality gap tolerance of ε = 10−2 . We ﬁrst solved the model treating all stage-2 variables and constraints as one large block, using the basic Benders’ decomposition algorithm. We next solved the model treating each of the 125 stage-2 scenarios as a separate block, using the more advanced version of Benders’ decomposition that takes full advantage of the block-ladder structure. The number of master iterations and the number of constraints 20 generated by the two diﬀerent solution methods are shown in Table 2. Notice that while the block-ladder version generates more constraints, it requires fewer major iterations. Table 3 shows the computation times of all of the solution methods. From this table, it is evident that the method that takes full advantage of the block-ladder structure was the most eﬃcient solution method. Master Iterations Generated Constraints Benders’ Decomposition Basic (1 Block) Block-Ladder (125 Blocks) 39 15 39 1,875 Table 2: Master iterations and generated constraints for the two implementations of Benders’ decomposition. No Decomposition Original Model CPU Time (min:sec) Benders’ Decomposition Basic (1 Block) Block-Ladder (125 Blocks) 5:03 3:22 1:27 Table 3: Computation times (min:sec) for the original model and the two implementations of Benders’ decomposition. Figure 2 shows the the upper and lower bounds generated by Benders’ decomposition at each iteration. Notice that there is initial rapid convergence of the bounds to the optimal objective function value, followed by very slow convergence thereafter. This is typical pattern that is observed in decomposition methods in general. Figure 3 shows the duality gap between the upper and lower bounds on a logarithmic scale. Here we see that the convergence rate of the duality gap is roughly linear. 21 50 Block Ladder 125 Blocks Basic − 1 Block Bounds (in billion $) 40 30 20 10 0 0 5 10 15 Iteration Figure 2: Upper and lower bounds for Benders’ decomposition. 7 Exercise The current version of the powerplant planning model uses the strategy of adding one constraint for every scenario in the model at each iteration of the Benders’ decomposition method. This means that k = 125 new constraints are added to the model at each outer iteration of the method. As discussed in class, this strategy outperforms the alternative strategy of adding only k = 1 new constraint at each outer iteration of the model. In this question, you are asked to explore intermediate strategies that might improve computation time for solving the powerplant planning model using Benders’ decomposition. By adding your own control logic in the region indicated in the script ﬁle, bender1.osc, experiment with diﬀerent strategies for limiting and/or controlling the number of new constraints added to the model at each outer iteration. You may also modify the model ﬁles or the data ﬁle if necessary. (One strategy that you might try would be to add a ﬁxed number of 22 5 Block Ladder 125 Blocks Basic − 1 Block 10 4 10 Upper bound − lower bound 3 10 2 10 1 10 0 10 −1 10 −2 10 0 5 10 15 20 Iteration 25 30 35 40 Figure 3: Gap between upper and lower bounds for Benders’ decomposition. constraints, say k = 10, 20, or 30, etc., constraints per outer iteration.) 23