Methods
1.0 Decomposition Introduction Consider the XYZ corporation that has 10 departments, each of which have a certain function necessary to the overall productivity of the corporation. The corporation has capability to make 100 different products, but it any particular month, it makes some of them and does not make others. The decision of which products to make is made by the CEO. The CEOs decision is a yes/no decision on each of the 100 products. Of course, the CEOs decision depends, in part, on the productivity of the departments and their capability to make a profit given the decision of which products to make. Each department knows its own particular business very well, and each has developed sophisticated mathematical programs (optimization problems) which provide the best (most profitable) way to use their resources given identification of which products to make. The organization works like this. The CEO makes a tentative decision on which products to make, based on his/her own mathematical program which assumes certain profits from each department based on that decision. He/she then passes that decision to the various departments. Each of the departments uses that information to determine how it is going to operate in order to maximize profitability. Then each department passes that information back to the CEO. Some departments may not be able to be profitable at all with the CEOs selection of products to make: these departments will also communicate this information to the CEO. Once the CEO gets all the information back from the departments, he or she will re-reun their optimization to select the products, likely resulting in a modified choice of products, and then the process will
1
repeat. At some point, the optimization problem solved by the CEO will not change from one iteration to the next. At this point, the CEO will believe the current selection of products is best. This is an example of a multidivisional problem [1, pg. 219]. Such problems involve coordinating the decisions of separate divisions, or departments, of a large organization, when the divisions operate autonomously. Solution of such problems often may be facilitated by separating them into a single master problem and subproblems where the master corresponds to the problem addressed by the CEO and the subproblems correspond to the problems addressed by the various departments. We developed our example where the master problem involved choice of integer variables and the subproblems involved choice of continuous variables. The reason for this is that it conforms to the form of a mixed-integer-programming (MIP) problem, which is the kind of problem we have most recently had interest. However, the master-subproblem relationship may be otherwise. It may also involve decisions on the part of the CEO to directly modify resources for each department. By resources, we mean the right-hand-side of the constraints. Such a scheme is referred to as a resource-directed approach. Alternatively, the master-subproblem relationship may involve decisions on the part of the CEO to indirectly modify resources by charging each department a price for the amount of resources that are used. The CEO would then modify the prices, and the departments would adjust accordingly. Such a scheme is called a price-directed approach. These types of optimization approaches are referred to as decomposition methods.
2
2.0 Connection with optimization: problem structure [2] Recall that our optimization problems were always like this: Minimize f(x) Subject to c1xb1 c2xb2 cmxbm We may place all of the row-vectors ci into a matrix, and all elements bi into a vector b, so that our optimization problem is now: Minimize f(x) Subject to A x b Problems that have special structures in the constraint matrices A are typically more amenable to decomposition methods. Almost all of these structures involve the constraint matrix A being blockangular. A block angular constraint matrix is illustrated in Fig. 1. In this matrix, the yellow-cross-hatched regions represent sub-matrices that contain non-zero elements. The remaining sub-matrices, not yellow-cross-hatched, contain all zeros. We may think of each yellow-cross-hatched as a department. The decision variables x1 are important only to department 1; the decision variables x2 are important only to department 2; and the decision variables x3 are important only to department 3. In this particular structure, we have no need of a CEO at all. All departments are completely independent!
3
x1
b1
x2
b2
x3
b3
Fig. 1: Block-angular structure Fig. 2 represents a structure where the decisions are linked at the CEO level, who must watch out for the entire organizations consumption of resources. In this case, the departments are independent, i.e., they are concerned only with decision on variables for which no other department is concerned, BUT the CEO is concerned with constraints that span across the variables for all departments. And so we refer to this structure as block-angular with linking constraints.
c0 1 c1
2
c2
3
c3
Fig. 2: Block-angular structure with linking constraints
4
This previous situation, characterized by Fig. 2, is not, however, the situation that we originally described. In the situation of Fig. 2, the CEO allocates resources. In the original description, the CEO chose values (1 or 0) for certain variables for which it was assumed would affect each department. The original situation would have different departments linked by variables, then, and not by constraints. The structure of the constraint matrix for this situation is shown in Fig. 3.
x1 b1
x2
b2
x3 x4
b3
Fig. 3: Block-angular structure with linking variables 3.0 Motivation for decomposition methods: solution speed To motivate decomposition methods, we consider introducing security constraints to a familiar problem: the OPF. The OPF may be posed as problem P0. Min f 0 ( x0 , u 0 )
P0
s.t.
hk ( xk , u 0 ) = 0 g k ( xk , u 0 ) g k
max
k=0 k =0
5
where hk(xk,u0)=0 represents the power flow equations and gk(xk,u0)gkmax represents the line-flow constraints. The index k=0 indicates this problem is posed for only the normal condition, i.e., the condition with no contingencies. Denote the number of constraints for this problem as N. Assumption: Lets assume that running time T of the algorithm we use to solve the above problem is proportional to the square of the number of constraints, i.e., N2. For simplicity, we assume the constant of proportionality is 1, so that T=N2. Now lets consider the security-constrained OPF (SCOPF). Its problem statement is given as problem Pc: Min f 0 ( x0 , u 0 )
g k ( xk , u0 ) g k k = 0,1, 2,..., c Notice that there are c contingencies to be addressed in the SCOPF, and that there are a complete new set of constraints for each of these c contingencies. Each set of contingency-related constraints is exactly like the original set of constraints (those for problem P0), except it corresponds to the system with an element removed.
Pc
s.t.
hk ( x k , u 0 ) = 0
max
k = 0,1,2,..., c
So the SCOPF must deal with the original N constraints, and also another set of N constraints for every contingency. Therefore, the total number of constraints for Problem PC is N+cN=(c+1)N. Under our assumption that running time is proportional to the square of the number of constraints, then the running time will be proportional to [(c+1)N]2=(c+1)2N2=(c+1)2T. What does this mean?
6
It means that the running time of the SCOPF is (c+1)2 times the running time of the OPF. So if it takes OPF 1 minute to run, and you want to run SCOPF with 100 contingencies, it will take you 1012 minutes, or 10,201 minutes to run the SCOPF. This is 170 hours, about 1 week!!!! Many systems need to address 1000 contingencies. This would take about 2 years! So this is what you do..
Solve OPF k=0 (normal condition)
Solve OPF k=1 (contingency 1)
Solve OPF k=2 (contingency 2)
Solve OPF k=3 (contingency 3)
Solve OPF k=c (contingency c)
Fig. 4: Decomposition solution strategy The solution strategy first solves the OPF (master problem) and then takes contingency 1 and re-solves the OPF, then contingency 2 and resolves the OPF, and so on (these are the subproblems). For any contingency-OPFs which require a redispatch, relative to the k=0 OPF, an appropriate constraint is generated, and at the end of the cycle these constraints are gathered and applied to the k=0 OPF. Then the k=0 OPF is resolved, and the cycle starts again. Experience has it that such an approach usually requires only 2-3 cycles.
7
Denote the number of cycles as m. Each of the individual problems has only N constraints and therefore requires only T minutes. There are (c+1) individual problems for every cycle. There are m cycles. So the amount of running time is m(c+1)T. If c=100 and m=3, T=1 minute, this approach requires 303 minutes. That would be about 5 hours (instead of 1 week). If c=1000 and m=3, T=1 minute, this approach requires about 50 hours (instead of 2 years). In addition, this approach is easily parallelizable, i.e., each individual OPF problem can be sent to its own CPU. This will save even more time. Fig. 5 compares computing time for a toy system. The comparison is between a full SCOPF, a decomposed SCOPF (DSCOPF), and a decomposed SCOPF where the individual OPF problems have been sent to separate CPUs.
Fig. 5
8
4.0 Benders decomposition J. F. Benders [3] proposed solving a mixed-integer programming problem by partitioning the problem into two parts an integer part and a continuous part. It uses the branch-and-bound method on the integer part and linear programming on the continuous part. The approach is well-characterized by the linking-variable problem illustrated in Fig. 3 where here the linking variables are the integer variables. In the words of A. Geoffrion [4], J.F. Benders devised a clever approach for exploiting the structure of mathematical programming problems with complicating variables (variables which, when temporarily fixed, render the remaining optimization problem considerably more tractable). Note in the below problem statements, all variables except z1* and z2* are considered to be vectors. There are three steps to Benders decomposition for mixed integer programs. The problem can be generally specified as follows:
T T max z1 = c1 x + c 2 w
s.t. Dw e A1 x + A2 w b x, w 0 w integer
Define the master problem and primal subproblem as
Master :
T * max z1 = c 2 w + z 2
Primal subproblem :
T max z 2 = c1 x
s.t. Dw e w0 w integer
s.t. A1 x b A2 w * x0
9
Some comments on duality for linear programs:
1. Number of dual decision variables is number of primal constraints. Number of dual constraints is number of primal decision variables.
2. Coefficients of decision variables in dual objective are righthand-sides of primal constraints. Problem D Problem P min G = 41 + 122 + 183 max F = 3x1 + 5 x2 subject to s.t. x1 4 1 + 33 3 2 x2 12 22 + 23 5 3 x1 + 2 x2 18 1 0, 2 0, 3 0 x1 0, x2 0 14444244443 144 2444 4 3
Primal Problem
Dual Problem
3. Coefficients of decision variables in primal objective are righthand-sides of dual constraints. Problem D Problem P min G = 41 + 122 + 183 max F = 3x1 + 5 x2 subject to s.t. x1 4 1 + 33 3 2 x2 12 22 + 23 5 3 x1 + 2 x2 18 1 0, 2 0, 3 0 x1 0, x2 0 14444244443 144 2444 4 3
Primal Problem
Dual Problem
4. Coefficients of one variable across multiple primal constraints are coefficients of multiple variables in one dual constraint.
Problem P max F = 3 x1 + 5 x2 s.t. x1 4 2 x2 12 3 x1 + 2 x2 18 x1 0, x2 0 144 2444 4 3
Primal Problem
Problem D min G = 41 + 122 + 183 subject to
1 + 33 3 22 + 23 5 1 0, 2 0, 3 0
Dual Problem
14444244443
5. If primal obj. is maximization, dual obj. is minimization. 6. If primal constraints are , dual constraints are .
10
From the above, we can write the dual of our primal subproblem. Dual subproblem : Primal subproblem :
T max z 2 = c1 x
min z 2 = (b A2 w *)
T
s.t. A1 x b A2 w *
s.t. A1T c1
x0 k 0, k Now consider the master problem and dual subproblem together: Master : Dual subproblem : T * max z1 = c 2 w + z 2 T min z 2 = (b A2 w *) s.t.
Dw e w0
s.t.
A1T c1
k 0, k w integer 1. Interdependence: The master depends on the subproblem optimal objective z2*, and the subproblem depends on the master optimal solution w*. Therefore, the solution to each problem depends on the solution obtained in the other problem. 2. Iterative procedure: We will solve the overall problem by iterating between the master and the subproblem. The master will be used to generate a solution w*, given a value (or a guess) for z2*. Then the subproblem will be used to get a new value of z2* and * using the solution w* obtained in the master. This will tell us one very important thing: when we resolve the master, we should constrain z2* to be no bigger than (b-A2w*)T*, i.e., z2*(b-A2w*)T*.
11
3. Upper bound: a. Initial solution: Start the solution procedure by solving the master problem with a guess for z2*. Since the dual problem is going to minimize (lower) z2, lets be safe and guess a large value of z2* at this initial master problem solution. Since this value of z2* is chosen large, we can be sure that the solution to the master, z1*, will be above the actual (overall problem optimal) solution, and so we will consider this solution z1* to be an upper bound on the actual solution. b. Successive solutions: As the iterations proceed, we will add constraints, so that the master problem solution z1*, will continuously decrease to either the actual (overall problem optimal) solution, or it will decrease to value greater than the actual solution. Thus, the value of z1*, obtained from the master problem, serves as an upper bound on the actual (overall problem optimal) solution. 4. Lower bound: The dual problem results in a new value of z2*, and it can then be added to c2Tw* (where w* was obtained from the last master problem solution) to provide another estimate of z1*. Since the dual problem minimizes z2, the solution c2Tw*+z2* will be a lower bound on z1. These constraints are referred to as optimality cuts. 5. Monotonicity: The upper bound decreases monotonically through the iterations, but the lower bound may not. 6. Feasibility: If the dual problem results in an unbounded solution, then it means the primal problem is infeasible. In this case, we must resolve the master problem with more restrictive constraints on w. The associated constraints are called feasibility cuts. 7. Algorithm: In what follows, we specify Q as the set of constraints binding the master program. It will change via the addition of feasibility and optimality cuts as the algorithm proceeds. Initially, Q={wk large number, for all k, z2*M, M large}.
12
Master problem:
T * max z1 = c2 w + z 2
s.t. wk constrained as in Q k
* z 2 M , M is large
wk Sub-problem (dual):
s.t.
integer k
min z 2 = (b A2 w*) T A1T c1
k 1. Solve the master problem using Branch and Bound (or any other integer programming method). Designate the solution as w*. 2. Using the value of w* found in step 1, solve the sub-problem (the dual) which gives z2* and *. There are two possibilities: a. If the solution is unbounded (implying the primal is infeasible), adjoin the most constraining feasibility constraint from (b-A2w)T0 to Q, and go to step 1. b. Otherwise, designate the solution as * and go to step 3. T 3. Compare z1 found in step 1 to c2 w * + z 2 * where z2 * = (b A2 w*)T * found in step 2. There are two possibilities: a. If they are equal (or within of each other), then the solution (w*, *) corresponding to the subproblem dual solution, is optimal and the primal variables x* are found as the dual variables1 within the subproblem. b. If they are not equal, adjoin an optimality constraint to Q T given by z 2 * (b A2 w*) * and go to step 1.
1
k 0,
(Dual variables are the coefficients of the objective function in the final iteration of the simplex method and are provided with the LP solution by a solver like CPLEX. They may be found from the fact that primal and dual objective functions must be equal at the optimum: cTx=(b-A2w*)T*.
13
Step 3 is a check on Benders optimal rule. Benders optimal rule: If (z1*, w*) is the optimal solution to the master problem, and (z2*, *) is the optimal solution to dual subproblem, and if
from master problem T c2 w * + (b A2 w*)T * = 1 3 14 244 2 4 3 z 2 * from subproblem from master problem * z1 {
,
then (z1*, w*, *) is the optimal solution for the complete problem. We will work an example using the formalized nomenclature of the previous summarized steps. Consider the following problem P0 [5]:
max subject to : P0
z1 = 4 x1 + 3x2 + 5w 2 x1 + 3x2 + w 12 2 x1 + x2 + 3w 12 w 20 x1 , x2 , w 0, w integer
Clearly P0 is a mixed integer problem. There are two ways to think about this problem. First way: Lets redefine the objective function as
z1 = 5w + z 2
where
z 2 = 4 x1 + 3x2
so that problem P1P below is equivalent to problem P0 above:
14
P1P
max subject to max z1 = 5w + subject to : w 20 w 0, w integer
z 2 = 4 x1 + 3x2
2 x1 + 3 x2 12 w 2 x1 + x2 12 3w x1 , x2 0
This way is similar to the way J. Bloom described his two-stage generation planning problem in [6]. Second way: Here, we will simply cast the problem into the general form outlined in our three-step procedure. Comparing to our general formulation,
T T max z1 = c1 x + c2 w
s.t. A1 x + A2 w b x, w 0 w
we find that
integer
4 c1 = , 3 2 A1 = 2
c2 = 5 3 1 12 , A2 = , b = 1 3 12
15
Step 1: The master problem is
T * max z1 = c 2 w + z 2
s.t.
* Q : wk constrained k , z 2 M , M is large
wk 0 and integer k
or
* T max z1 = c 2 w + z 2 * max z1 = 5w + z 2
s.t.
* Q : w 20 , z 2 M
s.t.
* Q : w 20, z 2 M
w 0 and integer
w 0 and integer
The solution to this problem is trivial: since the objective function is being maximized, we make w and z2* as large as possible, resulting in w*=20, z2*=M. Step 2: Using the value of w found in the master, get the dual: T 12 1 1 min z 2 = w * 12 3 2 min z 2 = (b A2 w*) T s.t. s.t.
A1T c1
0
2 2 1 4 3 1 3 2 1 , 2 0
Substituting, from step 1, w*=20, the subproblem becomes:
16
12 1 min z 2 = 20 1 = 81 48 2 12 3 2 s.t. 2 2 1 4 3 1 3 2 1 , 2 0 Because the ks are non-negative, all terms in the objective function are negative. Noting the ks are constrained from below, we may make them as large as we like, implying the objective function is unbounded. This occurs because the coefficients in the objective function are negative. The coefficients in the objective function are negative because the master problem, since it was not sufficiently constrained, yielded a poor choice of w. We need to correct this situation, by taking step 2b, which means we will add a feasibility constraint to the master problem. This feasibility constraint is contained in (b-A2w)T0, or
T
12 1 w 12 3
or
T
1 0 2
(12 w)1 + (12 3w)2 0 (12 w) 0, (12 3w) 0
To guarantee the above, without concern for what values of k are chosen, we must make resulting in
12 w,
4w
Clearly, w must be chosen to satisfy w4. This constraint is added to Q, and we repeat step 1.
17
Step 1:
* max z1 = 5w + z 2
s.t.
* Q : w 20, w 4, z 2 M
w 0 and integer The solution is clearly w=4, z2*=M, with z1*=5(4)+M=20+M.
Step 2: Using the value of w=4 found in the master, get the dual. T 12 1 1 min z 2 = w * 12 3 2 min z 2 = (b A2 w*) T s.t. s.t.
A1T c1
0
2 2 1 4 3 1 3 2 1 , 2 0
Substituting, from step 1, w*=4, the subproblem becomes: T 12 1 1 min z 2 = 4 = 81 12 3 2
s.t. 2 2 1 4 3 1 3 2 1 , 2 0
We can use CPLEX LP solver (or any other LP solver) to solve the above, obtaining the solution 1*=0, 2*=3, with objective function value z2*=0. (Intuitively, one observes that minimization of the objective subject to nonnegativity constraint on 1 requires 1=0; then 2 can be anything as long as it satisfies
18
22 4 2 2
2 3
Therefore a finite optimal solution is 1*=0, 2*=3.) Since we have a bounded dual solution, our primal is feasible, and we may proceed to step 3.
T Step 3: Compare z1 found in step 1 to c2 w * + z 2 * where z2 * = (b A2 w*)T * is found in step 2.
In step 1, solution of the master problem resulted in z1*=20+M. In step 2, solution of the subproblem resulted in z2*=0. In both problems, c2=5, and we found (master) or used (sub) w*=4. Benders optimal rule is
from master problem
c w * + (b A2 w*) * = 1 3 14 244 2 4 3
T 2 T z 2 * from subproblem
?
from master problem
* z1 {
Substitution yields:
from master problem
5 4 + {
z 2 * from subproblem
0 {
= 2023 1+ M
from master problem
?
The fact that they are not equal indicates that our solution is not optimal, since it does not satisfy Benders optimal rule. These two problems, the master and the subproblem, are really part of a single problem, and therefore for the single problem to be solved, the solutions to the master and subproblems must be consistent. That is, T when we maximize z1= c 2 w * + z 2 * in the master using a value of z2*, we need to find this value to be the same as the solution that the subproblem gives for z2*. If we do that (since c2w* is the same for both), the objective function from the master problem, z1*, will be the same as the sum of {c2Tw+z2*} where z2* is the objective function from the subproblem.
19
If we find that z2* differs in the master and subproblem, as we have found here, then we impose a constraint in the master based on the answer obtained in the subproblem. The fact that this constraint is imposed in order to obtain optimality makes it an optimality constraint, or in the language of Benders, an optimality cut. We distinguish between a feasibility cut and an optimality cut: Feasibility cut: Takes place as a result of finding an unbounded dual subproblem, which, by duality, implies an infeasible primal subproblem. It means that for the value of w found in the master problem, there is no possible solution in the primal subproblem. We address this by adding a feasibility cut (a constraint on w) to the master problem, where that cut is obtained from (b-A2w)T0. Optimality cut: Takes place as a result of finding that Benders optimal rule is not satisfied, i.e., that
from master problem T c2 w * + (b A2 w*)T * < 1 3 14 244 2 4 3 z 2 * from subproblem from master problem * z1 {
It means that the value of z2* used in the master problem is larger than the value of z2* computed in the subproblem. We address this by adding an optimality cut (a constraint on z2*) to the master problem, where that cut is obtained from z 2 * (b A2 w*)T * Recalling that we are still in Step 3 of our example, and that we found that Benders optimal rule is not satisfied, we need to obtain T the optimality cut from z 2 * (b A2 w*) * . With 1 1 0 12 b = , A2 = , = 3 2 3 12
12 1 0 0 * z 2 w = [12 w 12 3w] = 36 9 w 12 3 3 3
T
Now we return to step 1.
20
Step 1: Adjoin the optimality cut to Q, resulting in the following master problem: max z1 = 5w + (z * ) 2
subject to :
* Q : w 20, w 4, z * M, z 2 36 9 w 2
w 0, w integer
This all-integer program can be solved using a branch and bound algorithm (both CPLEX and Matlab have one), but the solution can be identified using enumeration, since w can only be 0, 1, 2, 3, or 4. For example, letting w=0, we have max z1 = (z * ) 2
* subject to : Q : z * M, z 2 36 2 The solution is recognized immediately, as z2*=36, z1*=36.
Likewise, letting w=1, we have max z1 = 5 + (z * 2
)
* subject to : Q : z * M, z 2 27 2 The solution is recognized immediately, as z2*=27, z1*=32.
Continuing on, we find the complete set of solutions are w=0, z2*=36, z1=36 w=2, z2*=27, z1=32 w=2, z2*=18, z1=28 w=3, z2*=9, z1=24 w=4, z2*=0, z1=20 Since the first one results in maximizing z1, our solution is w*=0, z2*=36, z1*=36.
21
Step 2: Using the value of w found in the master, get the dual: T 12 1 1 min z 2 = w * 12 3 2 min z 2 = (b A2 w*) T s.t. s.t.
A1T c1
0
2 2 1 4 3 1 3 2 1 , 2 0
Substituting, from step 1, w*=0, the subproblem becomes:
12 1 min z 2 = 0 1 = 121 + 12 2 12 3 2 s.t. 2 2 1 4 3 1 3 2 1 , 2 0 We can use CPLEX LP solver (or any other LP solver) to solve the above, obtaining the solution 1*=2, 2*=0, with objective function value z2*=24. Since we have a bounded dual solution, our primal is feasible, and we may proceed to step 3. Brief tutorial for using CPLEX. CPLEX version 10.1.0 resides on an ISU server called pluto. To access it, you will need to logon to pluto. To do that, you will need a telnet and ftp facility. You may find instructions on getting the appropriate telnet and ftp facilities at http://home.eng.iastate.edu/~jdm/ee458/Intro_CPLEX.pdf.
The first thing to do is to construct a file containing the problem. To construct this file, you can use the program called notepad under the accessories selection of the start button in Windows. Once you open notepad, you can immediately save to your local directory under the filename filename.lp. You can choose
22
T
filename to be whatever you want, but you will need the extension lp. To obtain the extension lp when you save, you should do save as and then choose all files. Otherwise, it will assign the suffix .txt to your file. Here is what I typed into the file I called mip.lp maximize 12 x1 + 12 x2 subject to 2 x1 + 2 x2 >= 4 3 x1 + x2 >= 3 x1 >= 0 x2 >= 0 end Once you get the above file onto Pluto, may simply type cplex101
CPLEX> read mip.lp CPLEX> primopt Detailed information on using CPLEX is available at
http://home.eng.iastate.edu/~jdm/ee458/Intro_CPLEX.pdf.
T Step 3: Compare z1 found in step 1 to c2 w * + z 2 * where z2 * = (b A2 w*)T * is found in step 2.
In step 1, solution of the master problem resulted in z1*=36 In step 2, solution of the subproblem resulted in z2*=24. In both problems, c2=5, and we found (master) or used (sub) w*=0.
23
Benders optimal rule is
from master problem
c w * + (b A2 w*) * = 1 3 14 244 2 4 3
T 2 T
z 2 * from subproblem
?
from master problem
* z1 {
Substitution yields:
from master problem
50 + {
z 2 * from subproblem
24 {
=
?
from master problem
36 {
Benders optimal rule is not satisfied, we need to obtain the T optimality cut from z 2 * (b A2 w*) * . With 1 1 2 12 b = , A2 = , = 3 2 0 12
12 1 2 2 * z 2 w = [12 w 12 3w] = 24 2 w 0 12 3 0
T
Now we return to step 1. Step 1: Adjoin the optimality cut to Q, resulting in the following master problem: max z1 = 5w + (z * ) 2
subject to :
* * Q : w 20, w 4, z * M, z 2 36 9w, z 2 24 2w 2
w 0, w integer This all-integer program can be solved using a branch and bound algorithm (both CPLEX and Matlab have one), but the solution can be identified using enumeration, since w can only be 0, 1, 2, 3, or 4. For example, letting w=0, we have max z1 = z * 2
(
)
* * subject to : Q : z * M, z 2 36, z 2 24 2 The solution is recognized immediately, as z2*=24, z1*=24.
Likewise, letting w=1, we have
24
max
z1 = 5 + z * 2
(
)
* * subject to : Q : z * M, z 2 27, z 2 22 2 The solution is recognized immediately, as z2*=22, z1*=27.
Continuing on, we find the complete set of solutions are w=0, z2*=24, z1=24 w=1, z2*=22, z1=27 w=2, z2*=18, z1=28 w=3, z2*=9, z1=24 w=4, z2*=0, z1=20 And so the third one results in maximizing z1, so our solution is w*=2, z2*=18, z1*=28. Step 2: Using the value of w found in the master, get the dual: T 12 1 1 min z 2 = w * 12 3 2 min z 2 = (b A2 w*) T s.t. s.t.
A1T c1
0
2 2 1 4 3 1 3 2 1 , 2 0
Substituting, from step 1, w*=2, the subproblem becomes:
12 1 min z 2 = 2 1 = 101 + 6 2 12 3 2 s.t. 2 2 1 4 3 1 3 2 1 , 2 0 We can use CPLEX LP solver (or any other LP solver) to solve the above, obtaining the solution 1*=0.5, 2*=1.5, with objective
25
T
function value z2*=14. Since we have a bounded dual solution, our primal is feasible, and we may proceed to step 3.
T Step 3: Compare z1 found in step 1 to c2 w * + z 2 * where z2 * = (b A2 w*)T * is found in step 2.
In step 1, solution of the master problem resulted in z1*=28 In step 2, solution of the subproblem resulted in z2*=14. In both problems, c2=5, and we found (master) or used (sub) w*=2.
Benders optimal rule is
from master problem
c w * + (b A2 w*) * = 1 3 14 244 2 4 3
T 2 T z 2 * from subproblem
?
from master problem
* z1 {
Substitution yields:
from master problem
5 2 + {
z 2 * from subproblem
14 {
=
?
from master problem
28 {
Benders optimal rule is not satisfied, we need to obtain the T optimality cut from z 2 * (b A2 w*) * . With 1 1 0.5 12 b = , A2 = , = 3 2 1.5 12
12 1 0.5 0 .5 z w = [12 w 12 3w] = 24 5w 1.5 12 3 1.5
* 2 T
Now we return to step 1.
26
Step 1: Adjoin the optimality cut to Q, resulting in the following master problem:
max subject to : z1 = 5 w + z * 2
(
)
* * * Q : w 20, w 4, z * M, z 2 36 9 w, z 2 24 2 w, z 2 24 5w 2
w 0, w integer
This all-integer program can be solved using a branch and bound algorithm (both CPLEX and Matlab have one), but the solution can be identified using enumeration, since w can only be 0, 1, 2, 3, or 4. Enumerating the solutions to this problem results in w=0: z2*=24, z1*=24 w=1: z2*=19, z1*=24 w=2: z2*=14, z1*=24 w=3: z2*=9, z1*=24 w=4: z2*=0, z1*=20 We see that w=0, 1, 2, and 3 are equally good solutions Steps 2 and 3: for each of these solutions, using the value of w found in the master, get the dual. Then check Benders rule. The general form of the dual is below. T 12 1 1 min z 2 = w * 12 3 2 min z 2 = (b A2 w*) T s.t. s.t.
A1T c1
0
T 2
2 2 1 4 3 1 3 2 1 , 2 0
* c w * + (b A2 w*) * = { z1 2 4 3 Benders optimal rule is 1 3 14 244 from master T from master problem z 2 * from subproblem problem
?
27
w*=0, the subproblem becomes:
12 1 min z 2 = 0 1 = 121 + 2 2 12 3 2 s.t. 2 2 1 4 3 1 3 2 1 , 2 0 Solution from CPLEX is 1=2, 2=0, with objective function value z2*=24.
T
Benders rule:
from master problem
50 + {
z 2 * from subproblem
24 {
=
?
from master problem
24 {
This solution is optimal. Dual variables obtained from CPLEX are x1=6, x2=0. w*=1, the subproblem becomes:
12 1 min z 2 = 1 1 = 111 + 9 2 12 3 2 s.t. 2 2 1 4 3 1 3 2 1 , 2 0 Solution from CPLEX is 1=0.5, 2=1.5, with objective function value z2*=19.
T
Benders rule:
from master problem
5 1 + {
z 2 * from subproblem
19 {
=
?
from master problem
24 {
This solution is optimal. Dual variables obtained from CPLEX are x1=4, x2=1.
28
w*=2, the subproblem becomes:
12 1 min z 2 = 2 1 = 101 + 6 2 12 3 2 s.t. 2 2 1 4 3 1 3 2 1 , 2 0 Solution from CPLEX is 1=0.5, 2=1.5, with objective function value z2*=14.
T
Benders rule:
from master problem
5 2 + {
z 2 * from subproblem
14 {
=
?
from master problem
24 {
This solution is optimal. Dual variables obtained from CPLEX are x1=2, x2=2. w*=3, the subproblem becomes:
12 1 min z 2 = 3 1 = 91 + 3 2 12 3 2 s.t. 2 2 1 4 3 1 3 2 1 , 2 0 Solution from CPLEX is 1=0, 2=3, with objective function value z2*=9.
T
Benders rule:
from master problem
53 + {
z 2 * from subproblem
9 {
=
?
from master {
This problem
24 solution is optimal. Dual variables obtained from CPLEX are x1=0, x2=3.
29
Problem summary: Recall our original problem:
max subject to : P0
z1 = 4 x1 + 3x2 + 5w 2 x1 + 3x2 + w 12 2 x1 + x2 + 3w 12 w 20 x1 , x2 , w 0, w integer
Optimal solutions to this problem result in an objective function value of z1=24 and are: w=0, x1=6, x2=0 w=1, x1=4, x2=1 w=2, x1=2, x2=2 w=3, x1=0, x2=3 Some comments about this problem: 1. It is coincidence that the values of x1 and x2 for the optimal solution also turn out to be integers. 2. The fact that there are multiple solutions is typical of MIP problems. MIP problems are generally non-convex. 5.0 Benders simplifications In the previous section, we studied problems having the following structure:
T T max z1 = c1 x + c 2 w
s.t. Dw e A1 x + A2 w b x, w 0 w integer
30
and we defined the master problem and primal subproblem as
Master :
T max z 2 = c1 x s.t. s.t. Dw e A1 x b A2 w * w0 x0 w integer However, what if our original problem appears as below, which is the same as the original problem except that it does not contain an x in the objective function, although the x still remains in one of the constraints.
T max z1 = c2 w
T * max z1 = c 2 w + z 2
Primal subproblem :
s.t. Dw e A1 x + A2 w b x, w 0 w integer
In this case, the master problem and the primal subproblem become:
Master :
max z2 = ??? s.t. s.t. Dw e A1 x b A2 w * w0 x0 w integer One sees clearly here that the primal subproblem has no z2 to maximize! One way to address this issue is to introduce a vector of non-negative slack variables s having one element for each constraint. We will minimize the sum of these slack variables, so that a non-zero value of this sum indicates the subproblem is
31
T max z1 = c2 w
Primal subproblem :
infeasible. That is, we replace our primal subproblem with a feasibility check subproblem, as follows: Master : Feasibility check subproblem : max z = cT w
1 2
s.t. Dw e w0
min
v = OnesT s
s.t. A1 x s b A2 w *
x 0, s 0 w integer Here, Ones is a column vector of 1s, so that v=OnesTs is the summation of all elements in the column vector s. When v=0, each constraint in A1x-sb-A2w* is satisfied so that A1xb-A2w*, which means the constraints to the original problem are in fact satisfied.
In this case, one observes that if v=0, then the problem is solved since Benders optimality rule will always be satisfied.
from master problem T c2 w * + (b A2 w*)T * = 1 3 14 244 2 4 3 z 2 * from subproblem from master problem * z1 {
Here, z1 is always zero, and the other two terms come from the master problem, therefore if the problem is feasible, it is optimal, and no step 3 is necessary. One question does arise, however, and that is what should be the feasibility cuts returned to the master problem if the feasibility check subproblem results in v>0? The answer to this is stated in [7] and shown in [8] to be v + A2(w* w) < 0 This kind of problem is actually very common. Figure 4, using a SCOPF to motivate decomposition methods for enhancing computational efficiency, is of this type. This is very similar to the so-called simultaneous feasibility test (SFT) of industry.
32
The SFT (Simultaneous Feasibility Test) is widely used in SCED and SCUC [9,10, 11]. SFT is a contingency analysis process. The objective of SFT is to determine violations in all post-contingency states and to produce generic constraints to feed into economic dispatch or unit commitment, where a generic constraint is a transmission constraint formulated using linear sensitivity coefficients/factors. The ED or UC is first solved without considering network constraints and security constraints. The results are sent to perform the security assessment in a typical power flow. If there is an existing violation, the new constraints are generated using the sensitivity coefficients/ factors and are added to the original problem to solve repetitively until no violation exists. The common flowchart is shown in Fig. 6.
Fig. 6 This section has focused on the very common case where the general Benders approach degenerates to a feasibility test problem only, i.e., the optimality test does not need to be done. There are at least three other degenerate forms of Benders: No feasibility problem: In some situations, the optimality problem will be always feasible, and so the feasibility problem is unnecessary.
33
Dual-role feasibility and optimality problem: In some applications, the feasibility and optimality problem can be the same problem.
Reference [7] provides examples of these degenerate forms of Benders decomposition. 6.0 Application of Benders to other Problem Types This section is best communicated by quoting from Geoffrion [4] (highlight added), considered the originator of Generalized Benders. J.F. Benders devised a clever approach for exploiting the structure of mathematical programming problems with complicating variables (variables which, when temporarily fixed, render the remaining optimization problem considerably more tractable). For the class of problems specifically considered by Benders, fixing the values of the complicating variables reduces the given problem to an ordinary linear program, parameterized, of course, by the value of the complicating variables vector. The algorithm he proposed for finding the optimal value of this vector employs a cutting-plane approach for building up adequate representations of (i) the extremal value of the linear program as a function of the parameterizing vector and (ii) the set of values of the parameterizing vector for which the linear program is feasible. Linear programming duality theory was employed to derive the natural families of cuts characterizing these representations, and the parameterized linear program itself is used to generate what are usually deepest cuts for building up the representations. In this paper, Benders' approach is generalized to a broader class of programs in which the parametrized subproblem need no longer be a linear program. Nonlinear convex duality theory is employed to derive the natural families of cuts corresponding to those in Benders' case. The conditions under which such a generalization is possible and appropriate are examined in detail.
34
The spirit of the above quotation is captured by the below modified formulation of our problem. The problem can be generally specified as follows:
T max z1 = f ( x) + c 2 w
s.t. Dw e F ( x) + A2 w b x, w 0 w integer
Define the master problem and primal subproblem as
Master :
T * max z1 = c 2 w + z 2
s.t. Dw e w0 w integer
Primal subproblem : max z 2 = f ( x) s.t. F ( x) b A2 w * x0
The Benders process must be generalized to solve the above problem since the subproblem is a nonlinear program (NLP) rather than a linear program (LP). Geoffrion shows how to do this [4]. In the above problem, w is integer, the master is therefore a linear integer program (LIP); the complete problem is therefore an integer NLP. If Benders can solve this problem, then it will also solve the problem when w is continuous, so that the master is LP and subproblem is NLP. If this is the case, then Benders will also solve the problem where both master and subproblem are LP, which is a very common approach to solving very-large-scale linear programs. Table 1 summarizes the various problems Benders is known to be able to solve. Table 1 Master Subproblem LP NLP
35
ILP
LP
One might ask whether Benders can handle a nonlinear integer program in the master, but it is generally unnecessary to do so since such problems can usually be decomposed to an ILP master with a NLP subproblem. 7.0 Application of Benders to Stochastic Programming For good, but brief overviews of Stochastic Programming, see [12] and [13]. In our example problem, we considered only a single subproblem, as shown below. T T max z1 = c1 x + c 2 w s.t. Dw e A1 x + A2 w b x, w 0 w integer To prepare for our generalization, we rewrite the above in a slightly different form, using slightly different notation: max z1 = cT w + d1T x1 s.t. Dw e B1w + A1 x1 b1 x1 , w 0 w integer Now we are in position to extend our problem statement so that it includes more than a single subproblem, as indicated in the structure provided below.
36
T T max z1 = c T w + d1T x1 + d 2 x 2 + .... + d n x n
s.t. Dw e B1 w + A1 x1 B2 w Bm w xk , w 0 + A2 x 2 O
b1
+ An x n
b2 M bn
In this case, the master problem is
max z1 = c w + zi ( xi )
T i =1 n
s.t. Dw e w0 where zi provide values of the maximization subproblem given by: max z i = d iT xi
s.t. Ai xi bi Bi w xi 0 Note that the constraint matrix for the complete problem appears as: D w e B A x b 1 1 1 1 B2 x2 b2 A2 M O M M Bn An xn bn
37
The constraint matrix shown above, if one only considers D, B1, and A1, has an L-shape, as indicated below. D w e B A x b 1 1 1 1 B2 x2 b2 A2 M O M M Bn An xn bn Consequently, methods to solve these kinds of problems, when they are formulated as stochastic programs, are called L-shaped methods. But what is, exactly, a stochastic program [12]? A stochastic program is an optimization approach to solving decision problems under uncertainty where we make some choices for now (the current period) represented by w, in order to minimize our present costs. After making these choices, event i happens, so that we take recourse2, represented by x, in order to minimize our costs under each event i that could occur in the next period. Our decision must be made a-priori, however, and so we do not know which event will take place, but we do know that each event i will have probability pi. Our goal, then, is to minimize the cost of the decision for now (the current period) plus the expected cost of the later recourse decisions (made in the next period).
2
Recourse is the act of turning or applying to a person or thing for aid.
38
A good application of this problem for power systems is the security-constrained optimal power flow (SCOPF) with corrective action. In this problem, we dispatch generation to minimize costs for the network topology that exists in this 15 minute period. Each unit generation level is a choice, and the complete decision is captured by the vector w. The dispatch costs are represented by cTw. These normal conditions are constrained by the power flow equations and by the branch flow and unit constraints, all of which are captured by Dwe. Any one of i=1,,n contingencies may occur in the next 15 minute period. Given that we are operating at w during this period, each contingency i requires that we take corrective action (modify the dispatch, drop load, or reconfigure the network) specified by xi. The cost of the corrective action for contingency i is diTxi, so that the expected costs over all possible contingencies is pidiTxi. Each contingency scenario is constrained by the post-contingency power flow equations, and by the branch flow and unit constraints, represented by Biw+Aixibi. The dependency on w (the pre-contingency dispatch) occurs as a result of, for example, unit ramp rate limitations. A 2-stage recourse problem is formulated below:
min z1 = c w + pi d iT xi
T i =1 n
s.t. Dw e B1 w + A1 x1 B2 w Bm w x, w 0
39
+ A2 x 2 O
b1
+ An x n
b2 M bn
where pi is the (scalar) probability of event i, and di is the vector of costs associated with taking recourse action xi. Each constraint equation Biw+Aixibi limits the recourse actions that can be taken in response to event i, and depends on the decisions w made for the current period. Formulation of this problem for solution by Benders (the L-shaped method) results in the master problem as
min z1 = c w + z i
T i =1 n
s.t. Dw e w0 where zi is minimized in the subproblem given by: min z i = pi d iT xi s.t. Ai xi bi Bi w xi 0 Note that the first-period decision, w, does not depend on which second-period scenario actually occurs (but does depend on a probabilistic weighting of the various possible futures). This is called the nonanticipativity property. The future is uncertain and so today's decision cannot take advantage of knowledge of the future.
Recourse models can be extended to handle multistage problems, where a decision is made now (in the current period), we wait for some uncertainty to be resolved, and then we make another decision based on what happened. The objective is to minimize the expected costs of all decisions taken. This problem can be appropriately thought of as the coverage of a decision tree, as shown in Fig. 7, where each level of the tree corresponds to another stochastic program.
40
Fig. 7 Multistage stochastic programs have been applied to handle uncertainty in planning problems before. This is a reasonable approach; however, one should be aware that computational requirements increase with number of time periods and number of scenarios (contingencies) per time period. Reference [14] by J. Beasley provides a good, but brief overview of multistage stochastic programming. Reference [15], notes for an entire course, provides a comprehensive treatment of stochastic programming including material on multistage stochastic programming.
41
8.0 Two related problem structures We found the general form of the stochastic program to be T T max z1 = c T w + d1T x1 + d 2 x2 + .... + d n xn s.t. Dw e b1 B1 w + A1 x1 + A2 x2 b2 B2 w O M + An xn bn Bm w
xk , w 0 so that the constraint matrix appears as D w e x b B A 1 1 1 1 B2 x2 b2 A2 M O M M Bn An xn bn If we move the w vector to the bottom of the decision vector, the constraint matrix appears as B1 x1 b1 A1 A2 B2 x 2 b2 O M M M An Bn x n bn D w e Notice that this is the structure that we introduced in Fig. 3 at the beginning of these notes, repeated here for convenience. We referred to this structure as block angular with linking variables. Now we will call it the Benders structure.
42
This means that coupling exists between what would otherwise be independent optimization problems via various parts of the problem occur via variables, in this case w. For example, in the SCOPF with corrective action, the almostindependent optimization problems are the n optimization problems related to the n contingencies. The dependency occurs via the dispatch determined for the existing (no-contingency) condition.
x1
b1
x2
b2
x3 x4
b3
A simpler example harkens back to our illustrations of the CEO and his or her organization. In this case, the CEO must decide on how to allocate a certain number of employees to each department, then each department minimizes costs (or maximizes profits) subject to constraints on financial resources. Recall from our discussion of duality that one way that primal and dual problems are related is that - coefficients of one variable across multiple primal constraints - are coefficients of multiple variables in one dual constraint, as illustrated below.
43
Problem P max F = 3 x1 + 5 x2 s.t. x1 4 2 x2 12 3 x1 + 2 x2 18 x1 0, x2 0 144 2444 4 3
Primal Problem
Problem D min G = 41 + 122 + 183 subject to
1 + 33 3 22 + 23 5 1 0, 2 0, 3 0
Dual Problem
14444244443
This can be more succinctly stated by saying that the dual constraint matrix is the transpose of the primal constraint matrix. You should be able to see, then, that the structure of the dual problem to a problem with the Benders structure looks like Fig. 2, repeated here for convenience.
c0 1 c1
2
c2
3
c3
This structure differs from that of the Benders structure (where the subproblems were linked by variables) in that now, the subproblems are linked by constraints. This problem is actually solved most effectively by another algorithm called Dantzig-Wolfe (DW) decomposition. We will refer to the above structure as the DW structure.
44
9.0 The Dantzig-Wolfe Decomposition Procedure Most of the material from this section is adapted from [5]. We attempt to solve the following problem P.
max subject to z= A00 x 0
T c1 x1 + A01 x1 T + c2 x2 + A02 x 2 T + L + ch xh + L + A0 h x h
= b0 = b1 = b2 M
A11 x1
P
x1 , x 2 , K , x h 0
A22 x 2
O
Ahh x h
= bh
where cj and xj have dimension nj1 Aij has dimension minj bi has dimension mi1 This problem has mi constraints and
i =0 h
n
i =0
h
i
variables.
Note: Constraints 2, 3, ,h can be thought of as constraints on individual departments. Constraint 1 can be thought of as a constraint on total corporate (organizational) resources (sum across all departments). Definition 1: For a linear program (LP), an extreme point is a corner point and represents a possible solution to the LP. The solution to any feasible linear program is always an extreme point. Figure 8 below illustrates 10 extreme points for some linear program.
45
5 4 6
100
3 7
90 80 70 60
8 2 9 1 10
50 40 30 20 10
5
Fig. 8 Observation 1: Each individual subproblem represented by
max z = ciT xi subject to Aii xi = bi
has its own set of extreme points. We refer to these extreme points as the subproblem i extreme points, denoted by xik , k = 1, K , pi where pi is the number of extreme points for subproblem i. Observation 2: Corresponding to each extreme point solution xik , the amount of corporate resources used is A0i xik , an m01 vector. Observation 3: The contribution to the objective function of the extreme point solution is ciT xik , a scalar. Definition 2: A convex combination of extreme points is a point
xik yik , where
k =1 k i
pi
y
k =1
pi
k i
= 1 , so that y ik is the fraction of extreme point
x in the convex combination. Figure 9 below shows a white dot in
the interior of the region illustrating a convex combination of the two extreme points numbered 4 and 9.
46
5 4 6
100
3 7
90 80 70 60
8 2 9 1 10
50 40 30 20 10
5
Fig. 9 Fact 1: Any convex combination of extreme points must be feasible to the problem. This should be self-evident from Fig. 9 and can be understood from the fact that the convex combination is a weighted average of the extreme points and therefore must lie between them (for two extreme points) or interior to them (for multiple extreme points). Fact 2: Any point in the feasible region may be identified by appropriately choosing the y ik . Observation 4a: Since the convex combination of extreme points is a weighted average of those extreme points, then the total resource usage by that convex combination will also be a weighted average of the resource usage of the extreme points, i.e.,
(A
pi k =1
0i
xik y ik . The
)
total resources for the complete problem is the summation over all of the subproblems,
(A
h pi i =1 k =1
0i
xik y ik .
)
47
Observation 4b: Since the convex combination of extreme points is a weighted average of those extreme points, then the contribution to the objective function contribution by that convex combination will also be a weighted average of the objective function contribution of the extreme points, i.e.,
(c
pi k =1
T i
xik y ik . The total objective function can
)
be expressed as the summation over all subproblems,
(c
h pi i =1 k =1
T i
xik y ik .
)
Based on Observations 4a and 4b, we may now transform our optimization problem P as a search over all possible combinations of points within the feasible regions of the subproblems to maximize the total objective function, subject to the constraint that the y ik must sum to 1.0 and must be nonnegative, i.e.,
max z = ciT xik y ik
i =1 k =1 h pi
(
)
P-T
subject to A00 x0 + A0i xik y ik = b0
i =1 k =1
h
pi
(
)
m0 constraints h constraints
y
k =1
pi
k i
= 1,
i = 1, K , h i = 1, K , h,
y ik 0,
k = 1, K , pi
Note that this new problem P-T (P-transformed) has m0+h constraints, in contrast to problem P which has has only
m
i =0
h
i
constraints.
Therefore it has far fewer constraints. However, whereas problem P
n
i =0
h
i
variables, this new problem has as many variables as it
has total number of extreme points across the h subproblems, and so it has a much larger number of variables.
p
i =0
h
i
,
48
The DW decomposition method solves the new problem without explicitly considering all of the variables. Understanding the DW method requires having a background in linear programming so that one is familiar with the revised simplex algorithm. We do not have time to cover this algorithm in this class, but it is standard in any linear programming class to do so. Instead, we provide an economic interpretation to the DW method. In the first constraint of problem P-T, b0 can be thought of as representing shared resources among the various subproblems i=1,,h. Let the first m0 dual variables of problem P-T be contained in the m01 vector . Each of these dual variables provide the change in the objective as the corresponding right-hand-side (a resource) is changed. That is, if b0k is changed by b0k+, then the optimal value of the objective is modified by adding k . Likewise, if the ith subproblem (department) increases its use of resource b0k by (instead of increasing the amount of the resource by ), then we can consider that that subproblem (department) has incurred a charge of k . This charge worsens its contribution to the complete problem P-T objective, and accounting for all shared resources b0k, k=1, m0, the contribution to the objective is ciT xi T A0i xi where A0i xi is the amount of shared resources consumed by the ith subproblem (department). One may think of these dual variables contained in as the prices that each subproblem (department) must pay for use of the corresponding shared resources.
49
Assuming each subproblem i=1,,h represents a different department in the CEOs organization, the DW-method may be interpreted in the following way, paraphrased from [2]: - If each department i worked independently of the others, then each would simply minimize its part of the objective function, i.e., ciT xi . - However, the departments are not independent but are linked by the constraints of using resources shared on a global level. - The right-hand sides b0k, k=1,,m0, are the total amounts of resources to be distributed among the various departments. - The DW method consists of having the CEO make each department pay a unit price, k, for use of each resource k. - Thus, the departments react by including the prices of the resources in its own objective function. In other words, o Each department will look for new activity levels xi which minimize ciT xi T A0i xi . o Each department performs this search by solving the following problem:
max ciT xi T A0i xi xik 0 k subject to Aii xi = bi
- The departments make proposals of activity levels xi back to the CEO, and the CEO then determines the optimal weights y ik for the proposals by solving problem P-T, getting a new set of prices , and the process repeats until all proposals remain the same. Reference [2], p. 346, and [16], pp. 349-350, provide good articulations of the above DW economic interpretation.
50
10.0 Movitating interest in multicommodity flow problems Multicommodity flow problems are optimization problems where one is interested in identifying the least cost method of moving two or more types of quantities (commodities) from certain locations where they are produced to certain locations where they are consumed, and the transport mechanism has the form of a network. Clearly, multicommodity flow problems are of great interest in commodity transportation problems. For example, both the rail and trucking industries extensively utilize multicommodity. Airlines also use it (first-class and economy passengers may be viewed as two different commodities). Although the transport mechanism for electric energy certainly has the form of a network (the transmission grid!), electric power planners have not traditionally been interested in multicommodity flow problems because of two issues: (a) Commodity flow (single or multiple) problems do not respect Kirchoffs Voltage Law (they do respect nodal balance). (b) Transmission networks flow only 1 commodity, electric energy. Issue (a) is one electric planners can live with in some cases when the flow on a branch is directly controllable. This is the case for DC lines and tends to be the case for AC intertie lines between different control areas. So-called network flow methods have therefore been deployed frequently to solve certain kinds of problems in electric planning. For example, Minimal cost network flow has been employed to solve problems combining fuel transportation (rail, natural gas) with electric energy transportation (transmission). Example applications include [17, 18]. A maximum flow algorithm has been employed to solve the multiarea reliability problem for electric power systems. An extensive set of notes on this method is located at www.ee.iastate.edu/~jdm/ee653/U21-inclass.doc, Section U21.7.
51
Issue (b), however, has precluded most electric power system interest in multicommodity flows. That may need to change in the future. Why? The reason is simple. Transportation systems, e.g., rail, truck (or highway), and airline systems have been almost totally decoupled from the electric system. There is fairly clear indication that this decoupled state is coming to an end in that the pluggable hybrid electric vehicle (PHEV) or the pluggable electric vehicle (PEV) is already becoming quite common, and projections are that electric utility grids are likely to see millions of such vehicles plugging into the grid within the next 5-10 years. Although rail in the US is almost entirely powered by on-board diesel generators, both Europe and Japan have utilized electric trains extensively, indicating the technology is available. These electric trains can have speeds significantly higher than existing US rail, which typically travels no faster than 100 mph. Increasing that speed (doubling?) could attract increased passenger travel, resulting in decreased airline use. These observations motivate an interest in consideration of the increased coupling, and thus interdependencies, between the electric system (including the fuel production and delivery systems) and the transportation system. Figure 10 illustrates these various systems. These ideas motivate interest in the development of modeling methods that provide optimized investment plans for electric infrastructure and transportation infrastructure simultaneously. Such modeling is likely to make use of multicommodity transportation flow problems.
52
Fig. 10 11.0 Multicommodity flow problem: formulation Multicommodity flow problems have the quality that the branches (arcs) between nodes may be capacitated, for each commodity flow, and for the sum total of all commodity flows. The below formulation is adapted from [16].
53
Suppose that we have a network with m nodes and n arcs in which there will flow t different commodities. We use the following nomenclature: ui: vector of upper limits on flow for commodity i in the arcs of the network, so that uipq is the upper limit on the flow of commodity i in arc (p,q). u: vector of upper limits on the sum of all commodities flowing in the arcs of the network, so that upq is the upper limit on the sum of all commodity flows in arc (p,q). ci: vector of arc costs in the network for commodity i, so that cipq is the unit cost of commodity i or arc (p,q). bi: vector of supplies (or demands) of commodity i in the network, so that biq is the supply (if biq>0) or demand (if biq<0) of commodity i at node q. A: node-arc incidence matrix of the network. xi: vector of flows of commodity i in the network. The linear programming formulation for the multicommodity minimal cost flow problem is as follows:
min subject to
c
i =1 t i =1
t
T i
xi
i
x
u i = 1, K , t
Axi = bi
0 xi u i i = 1, K , t
Observe that the first constraint is across multiple commodities. The second constraint, actually a set of constraints, is each across a single commodity. 12.0 Multicommodity flow problem: example This example is adapted from [16]. It is a 2-commodity minimal cost flow problem and is illustrated in Fig. 11.
54
(2, 0) (upq,u1pq,u2pq,c1pq,c2pq)
(4, 3, 4, 5, -2) 1 2
(1, -3)
(1, 5, 4, 1, -1) (b1q,b2q) 4 (-2, 0)
(3, 4, 2, 1, 0) (5, 1, 6, 4, 6)
3 (7, 5, 3, 0, 2) (-1, 3)
The constraint matrix for this system is displayed below.
0 0 0 01 0 0 0 0 4 u12 1 0 1 0 0 0 01 0 0 0 x 3 u 23 112 0 01 0 0 0 01 0 0 7 u 34 x123 u 0 01 0 0 0 01 0 1 0 x134 41 0 0 0 01 0 0 0 0 1 5 u 42 x141 b 2 0 0 1 0 0 0 0 0 0 1 x142 11 1 1 0 0 0 1 0 0 0 0 1 = b x 212 12 0 0 0 0 0 0 0 1 b 0 1 1 x 223 13 0 b 0 1 1 1 0 0 0 0 0 2 x 234 14 0 b 0 0 0 0 01 0 0 1 0 x 241 21 b 0 0 0 0 1 1 0 0 1 3 0 x 242 22 b 0 3 0 0 0 0 0 1 1 0 0 23 0 0 0 0 0 0 1 1 1 0 b24 0
The first five equations represent constraints on total commodity flow and are less than or equal to their respective right-hand-sides. The last 8 equations represent nodal balance constraints for the four nodes, commodity 1, and the four nodes, commodity 2. Observe that the structure is block-angular with linking constraints, i.e., this is what we called the GW-structure. We will therefore be able to effectively apply Dantzig-Wolfe decomposition to efficiently solve this problem.
55
[1] F. Hillier and G. Lieberman, Introduction to Operations Research, 4th edition, Holden-Day, Oakland California, 1986. [2] M. Minoux, Mathematical Programming: Theory and Algorithms, Wiley, 1986. [3] J. Benders, Partitioning procedures for solving mixed variables programming problems, Numerische Mathematics, 4, 238-252, 1962. [4] A. M. Geoffrion, Generalized benders decomposition, Journal of Optimization Theory and Applications, vol. 10, no. 4, pp. 237260, Oct. 1972. [5] S. Zionts, Linear and Integer Programming, Prentice-Hall, 1974. [6] J. Bloom, Solving an Electricity Generating Capacity Expansion Planning Problem by Generalized Benders' Decomposition, Operations Research, Vol. 31, No. 1, January-February 1983. [7] Y. Li, Decision making under uncertainty in power system using Benders decomposition, PhD Dissertation, Iowa State University, December 2008. [8] S. Granville, M. V. F. Pereira, G. B. Dantzig, B. Avi-Itzhak, M. Avriel, A. Monticelli, and L. M. V. G. Pinto, Mathematical decomposition techniques for power system, Tech. Rep. 2473-6, EPRI, 1988. [9] D. Streiffert, R. Philbrick, and A. Ott, A mixed integer programming solution for market clearing and reliability analysis, Power Engineering Society General Meeting, 2005. IEEE, pp. 27242731 Vol. 3, June 2005. [10] PJM Interconnection, On line training materials, www.pjm.com/services/training/training.html. [11] ISO New England, On line training materials, www.isone.com/support/training/courses/index.html. [12] www-fp.mcs.anl.gov/otc/Guide/OptWeb/continuous/constrained/stochastic/. [13] http://users.iems.northwestern.edu/~jrbirge/html/dholmes/StoProIntro.html. [14] http://people.brunel.ac.uk/~mastjjb/jeb/or/sp.html [15] http://www.lehigh.edu/~jtl3/teaching/ie495/ [16] M. Bazaraa, J. Jarvis, and H. Sherali, Linear Programming and Network Flows, second edition, Wiley, 1990. [17] A. Quelhas, E. Gil, J. McCalley, and S. Ryan, A Multiperiod Generalized Network Flow Model of the U.S. Integrated Energy System: Part I Model Description, IEEE Transactions on Power Systems, Volume 22, Issue 2, May 2007, pp 829 836. [18] A. Quelhas and J. McCalley, A Multiperiod Generalized Network Flow Model of the U.S. Integrated Energy System: Part II Simulation Results, IEEE Transactions on Power Systems Volume 22, Issue 2, May 2007, pp. 837 844.
56