## Bender's Decomposition implementation problem from unbounded constraint in master problem

Problems with modeling
Supakarn Sompasong
User
Posts: 2
Joined: 3 months ago

### Bender's Decomposition implementation problem from unbounded constraint in master problem

Hi everyone

I am currently trying to implement MILP or MIP bender's decomposition GAMS code from Erwin Kalvelagen with 5 continuous variables and 2 binary variables. The problem is I got infeasible equation from unbounded equation in Master problem (* Exec Error at line 482: Equation infeasible due to rhs value) (0 =l= -1 or 1 =l=0) The equation look like this

Code: Select all

``````unboundedcut(unbcutset)..
cutconst(unbcutset)
+sum((I,J),cutcoeff1(unbcutset,I,J)*yij(I,J))
+sum((K,L),cutcoeff2(unbcutset,K,L)*ykl(K,L))+constant =l= 0;``````
I think the constant term is not the foundation of the problem since I was implemented previous model with no binary variables (continuous variables only) and it can run successfully ( 1 iteration UB = LB )

I am really new in this software (1 year experience). what should I do to solve this problem?
Supakarn Sompasong
User
Posts: 2
Joined: 3 months ago

### Re: Bender's Decomposition implementation problem from unbounded constraint in master problem

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'
;
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:

*---------------------------------------------------------------------
* 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:

*---------------------------------------------------------------------
* 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:
* 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);