# ----------------------------------------
# LAGRANGIAN RELAXATION FOR
# THE LOCATION-TRANSPORTATION PROBLEM
# ----------------------------------------

printf "\nLP RELAXATION\n\n";

model trnloc2a.mod;
data trnloc2.dat;

option omit_zero_rows 1;
option display_eps .000001;
option solution_round 8;

option solver cplex;
option solver_msg 0;

option relax_integrality 1;
objective Shipping_Cost;

solve;

param LB; param UB;
let LB := Shipping_Cost.val;
let UB := sum {j in CITY} max {i in CITY} ship_cost[i,j];

option relax_integrality 0;

problem LowerBound: Build, Ship, Supply, Limit, Lagrangian;
problem UpperBound: Ship, Supply, Demand, Limit, Shipping_Cost;

let {j in CITY} mult[j] := 0;

param slack {CITY};
param scale  default 1;
param norm;
param step;

param same  default 0;
param same_limit := 3;

param iter_limit := 20;
param LBlog {0..iter_limit};  let LBlog[0] := LB;
param UBlog {0..iter_limit};  let UBlog[0] := UB;
param scalelog {1..iter_limit};
param steplog {1..iter_limit};

for {k in 1..iter_limit} { printf "\nITERATION %d\n\n", k;
   solve LowerBound;

   let {j in CITY} slack[j] := sum {i in CITY} Ship[i,j] - demand[j];

   if Lagrangian > LB + 0.000001 then {
      let LB := Lagrangian;
      let same := 0; }
   else let same := same + 1;

   if same = same_limit then {
      let scale := scale / 2;
      let same := 0;
      };

   let norm := sum {j in CITY} slack[j]^2;
   let step := scale * (UB - Lagrangian) / norm;

   let {j in CITY} mult[j] := mult[j] less step * slack[j];

   if sum {i in CITY} supply[i] * Build[i] 
         >= sum {j in CITY} demand[j] - 1e-8 then { 
      solve UpperBound;
      let UB := min (UB, Shipping_Cost);
      }

   let LBlog[k] := LB;
   let UBlog[k] := UB;
   let scalelog[k] := scale;
   let steplog[k] := step;
   }

printf "\n\n";
display LBlog, UBlog, scalelog, steplog;