Lagrangian Decomposition - GAMS

Problems with modeling
Post Reply
sivudu
User
User
Posts: 6
Joined: 6 years ago

Lagrangian Decomposition - GAMS

Post by sivudu »

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

Code: Select all

(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

Code: Select all

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.

Code: Select all

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

Code: Select all

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.

Code: Select all

X6+6X7<=9 in (DP2)
sub-problem is violated.

where did i made the mistake in Lagrangian problem coding. Please correct the mistake
sivudu
User
User
Posts: 6
Joined: 6 years ago

Re: Lagrangian Decomposition - Constraint Violation - GAMS

Post by sivudu »

I want to know why one constraint in the lagrangian solution is violated. The constraint in (DP2) sub-problem is violated.

The solution, i got is

Code: Select all

X1=0
X2=4
X3=8
X4=0
X5=0
X6=0
X7=2.
This solution violates the constraint

Code: Select all

X6+6X7<=9 
in (DP2) sub-problem.

please guide me.
Post Reply