Here is the full code I try to implemented with 2 binary variables (more than 2 nodes i to j to k)
Code: Select all
set i 'sources' /i1*i4/;
set j 'demands' /j1*j3/;
set k customer /k1*k3/
parameter supply(i) /
i1 10
i2 30
i3 40
i4 20
/;
parameter demand(j) /
j1 20
j2 50
j3 30
/;
parameter customer(k) /
k1 25
k2 30
k3 45
/;
table c1(i,j) 'variable cost i to j'
j1 j2 j3
i1 2.0 3.0 4.0
i2 3.0 2.0 1.0
i3 1.0 4.0 3.0
i4 4.0 5.0 2.0
;
table c2(j,k) 'variable cost j to k'
k1 k2 k3
j1 2.0 3.0 4.0
j2 3.0 2.0 1.0
j3 1.0 4.0 3.0
;
table f1(i,j) 'fixed cost i to j'
j1 j2 j3
i1 10.0 30.0 20.0
i2 10.0 30.0 20.0
i3 10.0 30.0 20.0
i4 10.0 30.0 20.0
;
table f2(j,k) 'fixed cost j to k'
k1 k2 k3
j1 10.0 30.0 20.0
j2 10.0 30.0 20.0
j3 10.0 30.0 20.0
;
parameter icap_ij(i,j) 'tight upperbounds for x(i,j)';
icap_ij(i,j) = min(supply(i),demand(j));
parameter icap_jk(j,k) 'tight upperbounds for x(i,j)';
icap_jk(j,k) = min(demand(j),customer(k));
*--------------------------------------------------------------------
* standard MIP problem formulation
*--------------------------------------------------------------------
variables
cost 'objective variable'
xij(i,j) 'shipments i to j'
xjk(j,k) 'shipments j to k'
yij(i,j) 'on-off indicator for link'
yjk(j,k) 'on-off indicator for link'
;
positive variable xij,xjk;
binary variable yij,yjk;
equations
obj 'objective'
cap(i) 'capacity constraint'
equal1(j)
equal2(j)
dem(j) 'demand equation'
cus(k) 'demand equation'
xyij(i,j) 'y=0 => x=0'
xyjk(j,k) 'y=0 => x=0'
;
obj.. cost =e= sum((i,j), f1(i,j)*yij(i,j) + c1(i,j)*xij(i,j))+sum((j,k), f2(j,k)*yjk(j,k) + c2(j,k)*xjk(j,k));
cap(i).. sum(j, xij(i,j)) =l= supply(i);
equal1(j).. sum(i,xij(i,j)) =g= sum(k,xjk(j,k));
equal2(j).. sum(i,xij(i,j)) =l= sum(k,xjk(j,k));
dem(j).. sum(k, xjk(j,k)) =l= demand(j);
cus(k).. sum(j, xjk(j,k)) =g= customer(k);
xyij(i,j).. xij(i,j) =l= icap_ij(i,j)*yij(i,j);
xyjk(j,k).. xjk(j,k) =l= icap_jk(j,k)*yjk(j,k);
display "--------------------- standard MIP formulation----------------------";
option optcr=0;
option limrow=0;
option limcol=0;
model fscp /obj,cap,equal1,equal2,dem,cus,xyij,xyjk/;
solve fscp minimizing cost using mip;
*---------------------------------------------------------------------
* Benders Decomposition Initialization
*---------------------------------------------------------------------
display "--------------------- BENDERS ALGORITHM ----------------------------";
scalar UB 'upperbound' /INF/;
scalar LB 'lowerbound' /-INF/;
yij.l(i,j) = 1;
yjk.l(j,k) = 1;
*---------------------------------------------------------------------
* Benders Subproblem
*---------------------------------------------------------------------
variable z 'objective variable';
positive variables
Du1(i) 'duals for capacity constraint'
Du2(j) 'duals for equal1 constraint'
Du3(j) 'duals for equal2 constraint'
Du4(j) 'duals for demand constraint'
Du5(k) 'duals for customer constraint'
Du6(i,j) 'duals for xyij constraint'
Du7(j,k) 'duals for xyij constraint'
;
equations
subobj 'objective'
subconstr1(i,j) 'dual constraint'
subconstr2(j,k) 'dual constraint'
;
* to detect unbounded subproblem
scalar unbounded /1.0e6/;
z.up = unbounded;
subobj.. z =e= sum(i, -supply(i)*Du1(i))+ sum(j, 0*(Du2(j)+Du3(j))) + sum(j,-demand(j)*Du4(j)) + sum(k,customer(k)*Du5(k))
+ sum((i,j), -icap_ij(i,j)*yij.l(i,j)*Du6(i,j))+ sum((j,k), -icap_jk(j,k)*yjk.l(j,k)*Du7(j,k))
;
subconstr1(i,j).. -Du1(i)+Du2(j)-Du3(j)-Du6(i,j) =l= c1(i,j);
subconstr2(j,k).. -Du2(j)+Du3(j)-Du4(j)+Du5(k)-Du7(j,k) =l= c2(j,k);
model subproblem /subobj, subconstr1, subconstr2/;
* reduce output to listing file:
subproblem.solprint=2;
* speed up by keeping GAMS in memory:
subproblem.solvelink=2;
*---------------------------------------------------------------------
* Benders Modified Subproblem to find unbounded ray
*---------------------------------------------------------------------
variable dummy 'dummy objective variable';
equations
modifiedsubobj 'objective'
modifiedsubconstr1(i,j) 'dual constraint'
modifiedsubconstr2(j,k) 'dual constraint'
edummy;
;
modifiedsubobj..
sum(i, -supply(i)*Du1(i))+ sum(j, 0*(Du2(j)+Du3(j))) + sum(j,-demand(j)*Du4(j)) + sum(k,customer(k)*Du5(k))
+ sum((i,j), -icap_ij(i,j)*yij.l(i,j)*Du6(i,j))+ sum((j,k), -icap_jk(j,k)*yjk.l(j,k)*Du7(j,k)) =e= 1;
modifiedsubconstr1(i,j)..
-Du1(i)+Du2(j)-Du3(j)-Du6(i,j) =l= 0;
modifiedsubconstr2(j,k)..
-Du2(j)+Du3(j)-Du4(j)+Du5(k)-Du7(j,k) =l= 0;
edummy.. dummy =e= 0;
model modifiedsubproblem /modifiedsubobj, modifiedsubconstr1, modifiedsubconstr2, edummy/;
* reduce output to listing file:
modifiedsubproblem.solprint=2;
* speed up by keeping GAMS in memory:
modifiedsubproblem.solvelink=2;
*---------------------------------------------------------------------
* Benders Relaxed Master Problem
*---------------------------------------------------------------------
set iter /iter1*iter50/;
set cutset(iter) 'dynamic set';
cutset(iter)=no;
set unbcutset(iter) 'dynamic set';
unbcutset(iter)=no;
variable z0 'relaxed master objective variable';
equations
cut(iter) 'Benders cut for optimal subproblem'
unboundedcut(iter) 'Benders cut for unbounded subproblem'
;
parameters
cutconst(iter) 'constant term in cuts'
cutcoeff1(iter,i,j)
cutcoeff2(iter,j,k)
;
cut(cutset).. z0 =g= sum((i,j), f1(i,j)*yij(i,j))+sum((j,k), f2(j,k)*yjk(j,k))
+ cutconst(cutset)
+ sum((i,j), cutcoeff1(cutset,i,j)*yij(i,j))
+ sum((j,k), cutcoeff2(cutset,j,k)*yjk(j,k));
unboundedcut(unbcutset)..
cutconst(unbcutset)
+ sum((i,j), cutcoeff1(unbcutset,i,j)*yij(i,j))
+ sum((j,k), cutcoeff2(unbcutset,j,k)*yjk(j,k)) =l= 0;
model master /cut,unboundedcut/;
* reduce output to listing file:
master.solprint=2;
* speed up by keeping GAMS in memory:
master.solvelink=2;
* solve to optimality
master.optcr=0;
*---------------------------------------------------------------------
* Benders Algorithm
*---------------------------------------------------------------------
scalar converged /0/;
scalar iteration;
scalar bound;
parameter yijbest(i,j);
parameter yjkbest(j,k);
parameter log(iter,*) 'logging info';
loop(iter$(not converged),
*
* solve Benders subproblem
*
solve subproblem maximizing z using lp;
*
* check results.
*
abort$(subproblem.modelstat>=2) "Subproblem not solved to optimality";
*
* was subproblem unbounded?
*
if (z.l+1 < unbounded,
*
* no, so update upperbound
*
bound = sum((i,j), f1(i,j)*yij.l(i,j)) + sum((j,k), f2(j,k)*yjk.l(j,k)) + z.l;
if (bound < UB,
UB = bound;
yijbest(i,j) = yij.l(i,j);
yjkbest(j,k) = yjk.l(j,k);
display yijbest,yjkbest;
);
*
* and add Benders' cut to Relaxed Master
*
cutset(iter) = yes;
else
*
* solve modified subproblem
*
solve modifiedsubproblem maximizing dummy using lp;
*
* check results.
*
abort$(modifiedsubproblem.modelstat>=2)
"Modified subproblem not solved to optimality";
*
* and add Benders' cut to Relaxed Master
*
unbcutset(iter) = yes;
);
*
* cut data
*
cutconst(iter) = sum(i, -supply(i)*Du1.l(i))+ sum(j, 0*(Du2.l(j)+Du3.l(j))) + sum(j,-demand(j)*Du4.l(j)) + sum(k,customer(k)*Du5.l(k));
cutcoeff1(iter,i,j) = -icap_ij(i,j)*yij.l(i,j)*Du6.l(i,j) ;
cutcoeff2(iter,j,k) = -icap_jk(j,k)*yjk.l(j,k)*Du7.l(j,k);
*
* solve Relaxed Master Problem
*
option optcr=0;
solve master minimizing z0 using mip;
*
* check results.
*
abort$(master.modelstat=4) "Relaxed Master is infeasible";
abort$(master.modelstat>=2) "Masterproblem not solved to optimality";
*
* update lowerbound
*
LB = z0.l;
log(iter,'LB') = LB;
log(iter,'UB') = UB;
iteration = ord(iter);
display iteration,LB,UB;
converged$( (UB-LB) < 0.1 ) = 1;
display$converged "Converged";
);
display log;
abort$(not converged) "No convergence";
*
* recover solution
*
yij.fx(i,j) = yijbest(i,j);
yjk.fx(j,k) = yjkbest(j,k);
fscp.solvelink=2;
fscp.solprint=2;
solve fscp minimizing cost using rmip;
abort$(fscp.modelstat<>1) "final lp not solved to optimality";
display "Benders solution",yij.l,xij.l,yjk.l,xjk.l,cost.l;
and its got an infeasible equation from unboudedcut equation in master problem (Equation infeasible due to rhs value 1 <= 0) what should I do? or this implementation valid only 1 binary variables yij(i,j) only like in original document?
