I wanted to solve a L.P.P using lagrangian decomposition. The problem consists of 3 sub-problems(DP1,DP2,DP3) and 3 coupling constraints.
The 3 sub-problems are
(DP3) Max X1+2X2+X3
s.t X1+X2<=6
X1+2X2<=8
X3<=10
X3<=2X2
X1,X2,X3>=0
(DP1) Max 5X4+2X5
s.t 2X4+X5<=9
X4,X5>=0
(DP2) Max 3X6+X7
s.t X6+6X7<=9
X6,X7>=0.
coupling constraints are
X1 = X4+X6
X2 =2X7
2X2=X3+X5.
This type of problems are observed in S & O.P (sales and operations planning) in supply chain management.
I have written the lagrangian problem in GAMS as follows.
Positive variables x1,x2,x3,x4,x5,x6,x7;
variables Zmaster,Zsub1,Zsub2,Z1sub1,Z2sub2,Zrelmaster;
parameter w1,w2,w3,w4,w5,w6,w7,Zfeas1,Zfeas2;
equations eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,obj,objsub1,objsub2,obj1sub1,obj2sub2,relmasterobj;
eq1.. x1 + x2 =l=6;
eq2.. x1 + 2*x2=l=8;
eq3.. x3 =l=10;
eq4.. x3 - 2*x2=l=0;
eq5.. x1 -x4 - x6=e=0;
eq6.. x2 -2*x7=e=0;
eq7.. 2*x2 -x3 -x5=e=0;
obj.. Zmaster=e= x1 +2*x2 +x3;
relmasterobj.. Zrelmaster =e= x1+2*x2+x3+w5*(x1-x4-x6)+w6*(x2-2*x7)+w7*(2*x2-x3-x5);
model master/obj,eq1,eq2,eq3,eq4,eq5,eq6,eq7/ ;
model relmaster/relmasterobj,eq1,eq2,eq3,eq4,eq5,eq6,eq7/;
model feasmaster/obj,eq1,eq2,eq3,eq4/;
solve master maximizing Zmaster using rmip;
parameters w,s51,s52,s5,s6,s7,norm5,norm6,norm7,LB,Zlr,improv1,improv2,improv3
,Zlr1,Zlr2,
Zlrmaster ;
parameter tol/1e-4/;
parameter step5/na/
step6/na/
step7/na/
step51/na/
step52/na/
step53/na/
step61/na/
step62/na/
step71/na/
step72/na/;
parameter alpha/1/;
parameter Zfeas/16/;
parameter LB/0/;
parameter Zlbest1 /0/;
parameter Zlbest2/0/;
parameter Zlbestmaster/0/;
set iter/iter1*iter20/;
parameter report(iter,*);
scalars status /0/;
scalar reset/5/
target1
target2
target3
count;
parameters srep(iter,*);
parameters wrep(iter,*) ;
parameters steprep(iter,*);
w1=eq1.m ;
w2 = eq2.m;
w3=eq3.m;
w4=eq4.m;
w5=eq5.m;
w6=eq6.m;
w7=eq7.m;
parameter improv,Zfeasmaster;
display w1,w2,w3,w4,w5,w6,w7;
eq8.. 2*x4 + x5=l=9;
objsub1.. Zsub1=e=5*x4 + 2*x5 + w5*(x1 -x4-x6) + w7*(2*x2- x3 -x5);
model sub1/eq8,objsub1,eq5,eq7,eq6/;
obj1sub1.. Z1sub1=e=5*x4 + 2*x3;
model feas1/eq8,obj1sub1,eq5,eq6,eq7/;
eq9.. x6 + 6*x7=l=9;
objsub2.. Zsub2=e=3*x6 + x7 + w5*(x1 -x4 -x6) + w6*(x2 - 2*x7);
obj2sub2.. Z2sub2=e=3*x6 + x7 ;
model feas2/eq9,obj2sub2,eq5,eq6,eq7/;
model sub2/eq9,objsub2,eq5,eq7,eq6/;
solve feas1 maximizing Z1sub1 using mip;
Zfeas1 = Z1sub1.l;
solve feas2 maximizing Z2sub2 using mip;
Zfeas2 = Z2sub2.l;
solve feasmaster maximizing Zmaster using mip;
Zfeasmaster = Zmaster.l;
option mip=cplex
rmip=cplex;
file results writes iteration report /solution/;
put results 'solvers used:RMIP = 'system.rmip/
' MIP = 'system.mip/;
solve master maximizing Zmaster using rmip;
put /'RMIP objective value = ',Zmaster.l:12:6/;
if(master.modelstat = %modelstat.optimal%,
status = %modelstat.optimal%
else
abort '***relaxed MIP not optimal',
' no subgradient iterations');
count=1;
loop(iter$(status=1),
Zlr1=0;
Zlr2=0;
Zlrmaster=0;
solve sub1 maximizing Zsub1 using mip;
Zlr1=Zlr1+Zsub1.l;
solve sub2 maximizing Zsub2 using mip;
Zlr2 = Zlr2 + Zsub2.l;
solve relmaster maximizing Zrelmaster using mip;
Zlrmaster = Zlrmaster + Zrelmaster.l;
improv1 =0;
improv2=0;
improv3=0;
improv1$(Zlr1 > Zlbest1 )=1;
Zlbest1 = max(Zlbest1,Zlr1);
improv2$(Zlr2 > Zlbest2) =1;
Zlbest2 = max(Zlbest2,Zlr2);
improv3$ (Zlrmaster > Zlbestmaster) =1;
Zlbestmaster = max(Zlbestmaster,Zlrmaster);
Zlr =x1.l +2*x2.l +x3.l ;
s5 = x1.l - x4.l - x6.l;
s6 = x2.l - 2*x7.l;
s7 = 2*x2.l - x3.l - x5.l;
norm5 = sqr(s5);
norm6 = sqr(s6);
norm7 = sqr(s7);
status$((norm5<tol)and (norm6<tol) and (norm7<tol)) =2;
status$(abs(Zlbest1 - Zfeas1)< 1e-4 and (Zlbest2 - Zfeas2) and (Zlbestmaster - Zfeasmaster) < 1e-4)=3;
status$((sub1.modelstat<>%modelstat.optimal%) or (sub2.modelstat<>%modelstat.optimal%))=4;
put results /iter.tl,Zlr1:16:6,Zlr2:16:6,Zlbest1:16:6,Zlbest2:16:6,norm5:16:6,norm6:16:6,norm7:16:6,
abs(Zlbest1-Zfeas1):16:6,abs(Zlbest2-Zfeas2):16:6;
if((status=2),
put // "subgr.method has converged, status =",status:5:0//;
put //"last solution found is optimal for IP problem"//;
);
if((status=3),
put //"subgr.method has converged, status = " , status:5:0//;
);
if((status=4),
put //"something wrong with the last lagrangian subproblem"//;
put //"status = ",status:5:0//;
);
report(iter,'Zlr1') = Zlr1;
report(iter,'Zlr2') = Zlr2;
report(iter,'Zlbest1') = Zlbest1;
report(iter,'Zlbest2') = Zlbest2;
report(iter,'Zlbestmaster') = Zlbestmaster;
report(iter,'norm5')= norm5;
report(iter,'norm6')=norm6;
report(iter,'norm7')=norm7;
srep(iter,'s5')=s5;
srep(iter,'s6')=s6;
srep(iter,'s7')=s7;
if(status=1,
target1 = (Zlbest1+Zfeas1)/2;
target2 = (Zlbest2+Zfeas2)/2;
target3 = (Zlbestmaster+Zfeasmaster)/2;
step51 = (alpha*(target1 -Zlr1)/norm5) $(norm5>tol);
step52 = (alpha*(target2-Zlr2)/norm5) $(norm5>tol);
step53 = (alpha*(target3-Zlrmaster)/norm5) $(norm5>tol);
step5 = max(step51,step52,step53);
step61 = (alpha*(target2-Zlr2)/norm6) $(norm6>tol);
step62 = (alpha*(target3-Zlrmaster)/norm6) $(norm6>tol);
step6 = max(step61,step62);
step71 = (alpha*(target1 -Zlr1)/norm7) $(norm7>tol);
step72 = (alpha*(target3-Zlrmaster)/norm7) $(norm7>tol);
step7 = max(step71,step72);
steprep(iter,'step5')= step5;
steprep(iter,'step6')=step6;
steprep(iter,'step7')=step7;
wrep(iter,'w5')=w5;
wrep(iter,'w6')= w6;
wrep(iter,'w7')=w7;
w5 = w5+step5*s5;
w6 = w6+step6*s6;
w7 = w7+step7*s7;
if(count>reset,
alpha = alpha/2;
count =1;
else if(improv1 or improv2,
count =1;
else
count = count+1;
)
)
)
);
display s5,s6,s7,LB,Zlr,x1.l,x2.l,x3.l,x4.l,x5.l,x6.l,x7.l,w5,w6,w7,Zfeas1,Zfeas2,Zfeasmaster;
display report,wrep,srep,steprep,Zlr1,Zlr2,norm5,norm6,norm7;
put //"best lagrangian bound = " , Zlbest1:10:5,Zlbest2:10:5;
;
I got the solution as
X1=0
X2=4
X3=8
X4=0
X5=0
X6=0
X7=2.
All the constraints are not satisfied by this solution. one constraint is violated.
X6+6X7<=9 in (DP2)
sub-problem is violated.
where did i made the mistake in Lagrangian problem coding. Please correct the mistake