# Run robust optimization on the DA network
option cplex_options 'timelimit=1200 timing=1 mipdisplay=0 mipgap=0.010 display 1';   
   
#option show_stats 1;    
option presolve_eps 1.e-9;   
option auxfiles rc;  
#option presolve 0;  
   
model robust_model;   
data  DA.txt;   
   
param Totalr;
let Totalr := 0;
let B := 0;
for {(o,d) in D} {
  for {s in S} {
    if r[o,d,s] > B then let B := r[o,d,s];
    let Totalr := Totalr + r[o,d,s];
  }
}
display B;
display Totalr;

#display reach;   
#display nospans;   
#display dist;   
#display RonL;   
#display AonL;   
#----------------Problem SCENARIO---------------------------------------   

problem SCENARIO {s in S}:   
  #  objective
               costs,   
  #  variables
  {p in P}     x[p,s],   
  {n in N}     LTE[n,s],   
  {e in E}     eLTE[e,s],  
  {e in E}     A[e,s],   
  {e in E}     R[e,s],   
  {e in E}     Fiber[e,s],   
  {e in E}     Channel[e,s],   
  {e in E}     Z[e,s],   
  #  constraints
  {(o,d) in D} Demands[o,d,s],   
  {e in E}     xtoZ,   
  {e in E}     LTEsOnEdge1[e,s],   
  {n in N}     LTEsAtNode[n,s],  
  {e in E}     FsOnLink1[e,s],   
  {e in E}     CsOnLink1[e,s],   
  {e in E}     FtoA[e,s],   
  {e in E}     CtoR[e,s],   
  {e in E}     FiberCapacity[e,s];  
   
param LTETotal  {S};   
param eLTETotal {S};  
param   RTotal  {S};   
param   ATotal  {S};   
param EquipCost {S};   
param TotalScenario;
param BudgetSum;
let Budget[3] := 1.0e+30;
let Budget[1] := 0.0;
let TotalScenario := 0;
let BudgetSum     := 0;

for {s in S} {   
  printf "\nSCENARIO MODEL: Scenario %d\n",s;   
  for {j in S} let delta[j] := 0; 
  let delta[s] := 1; 
  problem SCENARIO[s]; 
  solve   SCENARIO[s];   
  let LTETotal[s]  := sum {n in N} LTE[n,s];   
  let   RTotal[s]  := sum {e in E}   R[e,s];   
  let   ATotal[s]  := sum {e in E}   A[e,s];   
  let EquipCost[s] :=    
     costLTE*LTETotal[s] + costR*RTotal[s] + costA*ATotal[s];   
  printf "Equipment Cost for Scenario %d = %lf\n",s, EquipCost[s];   
  let TotalScenario := TotalScenario + 1;
  let BudgetSum     := BudgetSum + EquipCost[s];
  if EquipCost[s] < Budget[3] then  let Budget[3] := EquipCost[s]; 
  if EquipCost[s] > Budget[1] then  let Budget[1] := EquipCost[s]; 
}   
   
let Budget[2] := BudgetSum/TotalScenario;

display LTETotal;  display RTotal; display ATotal;   
param  ELTETotal;   param ERTotal; param  EATotal; param EEqCost; 
let    ELTETotal := sum {s in S} Prob[s]*LTETotal[s]; 
let    ERTotal   := sum {s in S} Prob[s]*RTotal[s]; 
let    EATotal   := sum {s in S} Prob[s]*ATotal[s]; 
let    EEqCost   := sum {s in S} Prob[s]*EquipCost[s]; 
display ELTETotal;  display ERTotal; display EATotal; display EEqCost; 
#display eLTE; 
#display x;   
#display LTE;   
#display A;   
#display R;   
#display Fiber;   
#display Channel;   
#display Z;   
#display eLTE;   
 display Budget;

#---------------Problem MEANVALUE---------------------------------------   

param EVLTETotal {1..3};
param EVRTotal   {1..3};
param EVATotal   {1..3};

for {(o,d) in D} let EVr[o,d] := ceil(sum {s in S} Prob[s]*r[o,d,s]);

problem MEANVALUE {b in {1..3}}:
  #  objective
                       EVcosts, 
  #  variables
                       EVETerm,
  {p in P}             EVx[p],   
  {n in N}             EVLTE[n],   
  {e in E}             EVeLTE[e],  
  {e in E}             EVA[e],   
  {e in E}             EVR[e],   
  {e in E}             EVFiber[e],   
  {e in E}             EVChannel[e],   
  {e in E}             EVZ[e],   
  #  constraints
                       EVETermConst,
                       EVBudget[b],
  {(o,d) in D}         EVDemandConst[o,d],   
  {e in E}             EVxtoZ,   
  {e in E}             EVLTEsOnEdge1[e],   
  {n in N}             EVLTEsAtNode[n],  
  {e in E}             EVFsOnLink1[e],   
  {e in E}             EVCsOnLink1[e],   
  {e in E}             EVFtoA[e],   
  {e in E}             EVCtoR[e],   
  {e in E}             EVFiberCapacity[e];  

#---------------Problem STOCHASTIC---------------------------------------   
param SPLTETotal {1..3};
param SPRTotal   {1..3};
param SPATotal   {1..3};

param ScenarioDemand {S};
for {s in S}
  let   ScenarioDemand[s] := sum{(o,d) in D} r[o,d,s];
param TotalDemand;
let   TotalDemand := sum{s in S} ScenarioDemand[s];
param TotalEquipCost;
let   TotalEquipCost := sum {s in S} EquipCost[s];
let   SPPenalty := 5.0*TotalEquipCost/TotalDemand;
display SPPenalty;

problem STOCHASTIC {b in {1..3}}:
  #  objective
                       SPcosts, 
  #  variables
                       SPETerm,
  {p in P}             SPx[p],   
  {n in N}             SPLTE[n],   
  {e in E}             SPeLTE[e],  
  {e in E}             SPA[e],   
  {e in E}             SPR[e],   
  {e in E}             SPFiber[e],   
  {e in E}             SPChannel[e],   
  {e in E}             SPZ[e],   
                       SPZTerm,
  {(o,d) in D,s in S}  SPzp[o,d,s],
  {(o,d) in D,s in S}  SPzm[o,d,s],
  #  constraints
                       SPETermConst,
                       SPBudget[b],
  {(o,d) in D,s in S}  SPDemands[o,d,s],   
  {e in E}             SPxtoZ,   
  {e in E}             SPLTEsOnEdge1[e],   
  {n in N}             SPLTEsAtNode[n],  
  {e in E}             SPFsOnLink1[e],   
  {e in E}             SPCsOnLink1[e],   
  {e in E}             SPFtoA[e],   
  {e in E}             SPCtoR[e],   
  {e in E}             SPFiberCapacity[e],  
                       SPZTermConst;

#---------------Problem WORSTCASE----------------------------------------   
param WCLTETotal {1..3};
param WCRTotal   {1..3};
param WCATotal   {1..3};

problem WORSTCASE {b in {1..3}}:
  #  objective
                       WCcosts, 
  #  variables
                       WCETerm,
  {p in P}             WCx[p],   
  {n in N}             WCLTE[n],   
  {e in E}             WCeLTE[e],  
  {e in E}             WCA[e],   
  {e in E}             WCR[e],   
  {e in E}             WCFiber[e],   
  {e in E}             WCChannel[e],   
  {e in E}             WCZ[e],   
                       WCMAXS,
  {(o,d) in D,s in S}  WCzp[o,d,s],
  {(o,d) in D,s in S}  WCzm[o,d,s],
  #  constraints
                       WCETermConst,
                       WCBudget[b],
  {(o,d) in D,s in S}  WCDemands[o,d,s],   
  {e in E}             WCxtoZ,   
  {e in E}             WCLTEsOnEdge1[e],   
  {n in N}             WCLTEsAtNode[n],  
  {e in E}             WCFsOnLink1[e],   
  {e in E}             WCCsOnLink1[e],   
  {e in E}             WCFtoA[e],   
  {e in E}             WCCtoR[e],   
  {e in E}             WCFiberCapacity[e],  
                       WCMAXSConst;

#---------------Problem ROBUSTOPTIMIZATION-------------------------------   
param ROLTETotal {1..3};
param RORTotal   {1..3};
param ROATotal   {1..3};

let PW := {1,2,3,4};
for {pw in PW} {
  let slopep[pw] := (2*pw-1)*K;
  let slopem[pw] := (2*pw-1)*G;
}

problem ROBUSTOPTIMIZATION {b in {1..3}}:
  #  objective
                                 Regret, 
  #  variables
                                 ROETerm,
  {p in P}                       ROx[p],   
  {n in N}                       ROLTE[n],   
  {e in E}                       ROeLTE[e],  
  {e in E}                       ROA[e],   
  {e in E}                       ROR[e],   
  {e in E}                       ROFiber[e],   
  {e in E}                       ROChannel[e],   
  {e in E}                       ROZ[e],   
  {(o,d) in D, s in S}           ROzp[o,d,s],
  {(o,d) in D, s in S}           ROzm[o,d,s],
  {(o,d) in D, s in S, pw in PW} ROzppw[o,d,s,pw],
  {(o,d) in D, s in S, pw in PW} ROzmpw[o,d,s,pw],
  #  constraints
                                 ROETermConst,
                                 ROBudget[b],
  {(o,d) in D, s in S}           RODemands[o,d,s],   
  {e in E}                       ROxtoZ,   
  {e in E}                       ROLTEsOnEdge1[e],   
  {n in N}                       ROLTEsAtNode[n],  
  {e in E}                       ROFsOnLink1[e],   
  {e in E}                       ROCsOnLink1[e],   
  {e in E}                       ROFtoA[e],   
  {e in E}                       ROCtoR[e],   
  {e in E}                       ROFiberCapacity[e],  
  {(o,d) in D, s in S}           Piecewise1[o,d,s],
  {(o,d) in D, s in S}           Piecewise2[o,d,s];

#---------------Problem UNDERPROVISION----------------------------------   

#problem UNDERPROVISION {s in S}:
  #  objective 
#              under,
  #  variables
# {e in E}     eLTE[e,s],
# {e in E}     ROeLTE[e],
# {e in E}     uLTE[e,s],
# {e in E}     oLTE[e,s],
#              uLTETotal[s], 
#              oLTETotal[s],
# {e in E}     R[e,s],
# {e in E}     ROR[e],
# {e in E}     uR[e,s],
# {e in E}     oR[e,s],
#              uRTotal[s], 
#              oRTotal[s],
# {e in E}     A[e,s],
# {e in E}     ROA[e],
# {e in E}     uA[e,s],
# {e in E}     oA[e,s],
#              uATotal[s], 
#              oATotal[s],
#  #  constraints
# {e in E}     uoLTE[e,s],
# {e in E}     uoR[e,s],
# {e in E}     uoA[e,s],
#              uLTETotalConst[s], 
#              oLTETotalConst[s],
#              uRTotalConst[s], 
#              oRTotalConst[s],
#              uATotalConst[s], 
#              oATotalConst[s];

param EVUnrouted {1..3};
param SPUnrouted {1..3};
param WCUnrouted {1..3};
param ROUnrouted {1..3};
param EVRegret   {1..3};
param SPRegret   {1..3};
param WCRegret   {1..3};
param RORegret   {1..3};
param SumEVx;
param SumSPx;
param SumWCx;
param SumROx;
param EVz;
param SPz;
param WCz;
param ROz;
for {k in {1..3}} {
  let EVUnrouted[k] := 0;
  let SPUnrouted[k] := 0;
  let WCUnrouted[k] := 0;
  let ROUnrouted[k] := 0;
  let EVRegret[k]   := 0;
  let SPRegret[k]   := 0;
  let WCRegret[k]   := 0;
  let RORegret[k]   := 0;
}

for {k in {1..3} } {
  problem MEANVALUE[k];
  unfix EVETerm;
  unfix EVx;  
  unfix EVLTE; 
  unfix EVeLTE; 
  unfix EVA;   
  unfix EVR;
  unfix EVFiber;   
  unfix EVChannel;   
  unfix EVZ;

  printf"\nMEANVALUE MODEL:  Budget #%d\n",k;
  solve MEANVALUE[k];
  display EVETerm;

  let EVLTETotal[k]  := sum {n in N} EVLTE[n];
  let EVRTotal[k]    := sum {e in E} EVR[e];
  let EVATotal[k]    := sum {e in E} EVA[e];
  problem STOCHASTIC[k];
  unfix SPETerm;
  unfix SPx;  
  unfix SPLTE; 
  unfix SPeLTE; 
  unfix SPA;   
  unfix SPR;
  unfix SPFiber;   
  unfix SPChannel;   
  unfix SPZ;
  unfix SPZTerm;
  unfix SPzp;
  unfix SPzm;

  printf"\nSTOCHASTIC MODEL:  Budget #%d\n",k;
  solve STOCHASTIC[k];
  display SPETerm;

  let SPLTETotal[k]  := sum {n in N} SPLTE[n];
  let SPRTotal[k]    := sum {e in E} SPR[e];
  let SPATotal[k]    := sum {e in E} SPA[e];

  problem WORSTCASE[k];
  unfix WCETerm;
  unfix WCx;  
  unfix WCLTE; 
  unfix WCeLTE; 
  unfix WCA;   
  unfix WCR;
  unfix WCFiber;   
  unfix WCChannel;   
  unfix WCZ;
  unfix WCMAXS;
  unfix WCzp;
  unfix WCzm;

  printf"\nWORSTCASE MODEL:  Budget #%d\n",k;
  solve WORSTCASE[k];
  display WCETerm;

  let WCLTETotal[k]  := sum {n in N} WCLTE[n];
  let WCRTotal[k]    := sum {e in E} WCR[e];
  let WCATotal[k]    := sum {e in E} WCA[e];

  problem ROBUSTOPTIMIZATION[k];
  unfix ROETerm;
  unfix ROx;  
  unfix ROLTE; 
  unfix ROeLTE; 
  unfix ROA;   
  unfix ROR;
  unfix ROFiber;   
  unfix ROChannel;   
  unfix ROZ;
  unfix ROzp;
  unfix ROzm;
  unfix ROzppw;
  unfix ROzmpw;

  printf"\nROBUSTOPTIMIZATION MODEL:  Budget #%d\n",k;
  solve ROBUSTOPTIMIZATION[k];
  display ROETerm;

  let ROLTETotal[k]  := sum {n in N} ROLTE[n];
  let RORTotal[k]    := sum {e in E} ROR[e];
  let ROATotal[k]    := sum {e in E} ROA[e];

  for {(o,d) in D} {
    let SumEVx := sum{p in J[o,d]} EVx[p];
    let SumSPx := sum{p in J[o,d]} SPx[p];
    let SumWCx := sum{p in J[o,d]} WCx[p];
    let SumROx := sum{p in J[o,d]} ROx[p];
    for {s in S} {
      let EVz := r[o,d,s] - SumEVx;
      let SPz := r[o,d,s] - SumSPx;
      let WCz := r[o,d,s] - SumWCx;
      let ROz := r[o,d,s] - SumROx;

      if EVz > 0 then {
        let EVUnrouted[k] := EVUnrouted[k] + EVz;
        let EVRegret[k]   := EVRegret[k] + Prob[s] *
          <<H, 2*H, 3*H; K, 3*K, 5*K, 7*K>> EVz;
      }
      else {
        let EVz := -EVz;
        let EVRegret[k] := EVRegret[k] + Prob[s] *
          <<H, 2*H, 3*H; G, 3*G, 5*G, 7*G>> EVz;
      }

      if SPz > 0 then {
        let SPUnrouted[k] := SPUnrouted[k] + SPz;
        let SPRegret[k]   := SPRegret[k] + Prob[s] *
          <<H, 2*H, 3*H; K, 3*K, 5*K, 7*K>> SPz;
      }
      else {
        let SPz := -SPz;
        let SPRegret[k] := SPRegret[k] + Prob[s] *
          <<H, 2*H, 3*H; G, 3*G, 5*G, 7*G>> SPz;
      }

      if WCz > 0 then {
        let WCUnrouted[k] := WCUnrouted[k] + WCz;
        let WCRegret[k] := WCRegret[k] + Prob[s] *
          <<H, 2*H, 3*H; K, 3*K, 5*K, 7*K>> WCz;
      }
      else {
        let WCz := -WCz;
        let WCRegret[k] := WCRegret[k] + Prob[s] *
          <<H, 2*H, 3*H; G, 3*G, 5*G, 7*G>> WCz;
      }

      if ROz > 0 then {
        let ROUnrouted[k] := ROUnrouted[k] + ROz;
        let RORegret[k] := RORegret[k] + Prob[s] *
          <<H, 2*H, 3*H; K, 3*K, 5*K, 7*K>> ROz;
      }
      else {
        let ROz := -ROz;
        let RORegret[k] := RORegret[k] + Prob[s] *
          <<H, 2*H, 3*H; G, 3*G, 5*G, 7*G>> ROz;
      }
    }
  }

# for {s in S} {
#   printf"\nUNDERPROVISION MODEL: Scenario %d Budget #%d\n",s, k;
#   for {j in S} let delta[j] := 0;
#   let delta[s] := 1;
#   problem UNDERPROVISION[s];
#   fix   eLTE;
#   fix   ROeLTE;
#   fix   R;
#   fix   ROR;
#   fix   A;
#   fix   ROA;
#   unfix uLTE;
#   unfix oLTE;
#   unfix uR;
#   unfix oR;
#   unfix uA;
#   unfix oA;
#   unfix uLTETotal;
#   unfix oLTETotal;
#   unfix uRTotal;
#   unfix oRTotal;
#   unfix uATotal;
#   unfix oATotal;
#   solve   UNDERPROVISION[s];
#}

# display uLTETotal;
# param   uLTEGrandTotal;
# let     uLTEGrandTotal := sum {s in S} uLTETotal[s];
# display uLTEGrandTotal;
# display oLTETotal;
# param   oLTEGrandTotal;
# let     oLTEGrandTotal := sum {s in S} oLTETotal[s];
# display oLTEGrandTotal;
# display uRTotal;
# param   uRGrandTotal;
# let     uRGrandTotal   := sum {s in S} uRTotal[s];
# display uRGrandTotal;
# display oRTotal;
# param   oRGrandTotal;
# let     oRGrandTotal   := sum {s in S} oRTotal[s];
# display oRGrandTotal;
# display uATotal;
# param   uAGrandTotal;
# let     uAGrandTotal   := sum {s in S} uATotal[s];
# display uAGrandTotal;
# display oATotal;
# param   oAGrandTotal;
# let     oAGrandTotal   := sum {s in S} oATotal[s];
# display oAGrandTotal;

}

# Scale the regret
for {k in {1..3} } {
  let EVUnrouted[k] := 100*EVUnrouted[k]/Totalr;
  let SPUnrouted[k] := 100*SPUnrouted[k]/Totalr;
  let WCUnrouted[k] := 100*WCUnrouted[k]/Totalr;
  let ROUnrouted[k] := 100*ROUnrouted[k]/Totalr;
  let EVRegret[k]   := EVRegret[k]/RORegret[k];
  let SPRegret[k]   := SPRegret[k]/RORegret[k];
  let WCRegret[k]   := WCRegret[k]/RORegret[k];
}

display EVLTETotal;
display EVRTotal;
display EVATotal;

display SPLTETotal;
display SPRTotal;
display SPATotal;

display ROLTETotal;
display RORTotal;
display ROATotal;

display WCLTETotal;
display WCRTotal;
display WCATotal;

display EVUnrouted;
display SPUnrouted;
display WCUnrouted;
display ROUnrouted;

display EVRegret;
display SPRegret;
display WCRegret;
display RORegret;

quit;