############################### Problem Input ######################## 
param SEED; 
param MTSO; 
param TOWER; 
param SUBSCRIBER; # Number of subscriber locations 
param dLo;        # Lower bound of demands at each subscriber location 
param dHi;        # Upper bound of demands at each subscriber location 

param r;          # annual revenue genreated by a channel 
param aLo;        # Lower bound of a[l] 
param aHi;        # Upper bound of a[l] 
param bLo;        # Lower bound of b[k] 
param bHi;        # Upper bound of b[k] 
param cperfoot;   # The annualized cost of providing 
                  # a foot of link from base station to MTSO 
param hperfoot;   # The annualized cost of providing 
                  # a foot of link between two MTSOs 
param SIRmin;     # minimum signal to interference ratio 
param x_mile; 
param y_mile; 
param alpha;
 
data p30;             
display SEED,MTSO,TOWER,SUBSCRIBER,dLo,dHi;
 
param min_dist;   # minimum distance between all data 
let   min_dist := 0;
###############################  Model ###############################  
  
###############################  Sets  ###############################  
 
set L default 1 .. TOWER;   # candidate tower locations  
set M default 1 .. SUBSCRIBER;  # candidate subscriber locations  
set K default 1 .. MTSO;   # candidate MTSO locations,   
printf "\n"; 
#display L,M,K,SEED; 
 
set K0 default {0} union K; # MTSO 0 is the dummy node  
  
set C{M} within L default L;     
                              # C[m] denotes the set of towers that can   
                              # service subscribers at m  
  
set P{L} within M default {}; # P[l] denotes the subscriber locations  
                              # that can be assigned to tower l  
                        
set LV within L default {};   # Towers for which SIRConst is violated  
 
set A within {M, L} default {}; # Active x[m,l]'s (x[m,l] < 0)  
 
  
###############################  Constants  ###############################  
# In this section only list the parameters that appear in the constraints  
# move the other ones to the code section  
 
param d{M};    #default 5     
               # d[m] denotes the number of channels   
               # required to service all subscribers at   
               # location m  
  
param dd {M};  # dd[m]= min(d[m],1+1/SIRmin) 
               # used in constraint BuildTowers 
 
#param r;      # default 42820;   
               # annual revenue generated by a channel  
 
param rho default 0.25; # mandated minimum service requirement # 0 < rho < 1  
  
param a{L};    # default 100000; #145945;   
               # a[l] denotes the cost for constructing  
               # a tower at location l    
               # (this needs to be an annual cost)  
  
param b{K};    # default 145945;   
               # b[k] deontes the annualized cost for  
               # constructing an MTSO at location k  
  
param c{L,K};  # c[l,k] is the annualized cost of providing  
               # a link from base station l to MTSO k  
  
param h{K,K0}; # h[k1,k2] is the annualized cost of providing a link   
               # from MTSO k1 to MTSO k2  
  
#param SIRmin default 0.03125; #default 0.009789;   
                               # minimum signal to interference ratio   
  
param g{M,L};         # g[m,l] denotes the attenuation factor   
                      # between subscriber m and tower l  
  
param G {m in M, l in L, j in C[m]}; 
      # G[m,l,j]=g[m,l]/g[m,j] 
      # used in constraint SIRConst 
 
param TotalDemand;    # total demand;  
  
param beta {L} default 0;    # used in constraint SIRConst    
  
###### Variables ######  
  
var y{L} >=0, <=1;# binary;      # y[l] = 1 if a tower is built at location l;  
                      # and  = 0, otherwise  
var yy {L} binary;  
 
var x{m in M,C[m]} >= 0 ;#integer;     
                      # x[m,l] denotes the number of subscribers at   
                      # location m assigned to tower l  
 
var xx {m in M,C[m]} >=0 integer; 
 
 
var z {K0} binary; # z[k] = 1 if and only if MTSO k is constructed  
  
var s {L,K} >=0,<=1 ;#binary; # s[l,k] = 1 if and only if base station l is  
# connected to MTSO k  
  
var u {k1 in K, k2 in K0} >= 0 ;   
# u[k1,k2] = the flow from MTSO k1 to MTSO k2 
 
var w {k1 in K, k2 in K0} binary; 
# If u[k1,k2]>0, w[k1,k2]=1 
# Otherwise, w[k1,k2]=0 
 
###### Constraints ######  
 
subject to IP1 {l in L}: y[l] = yy[l]; 
 
subject to IP2 {m in M, l in C[m]}: x[m,l] = xx[m,l]; 
 
subject to W0 {j in K,k in K0: j<>k}: 
   2*w[j,k] <= z[j] + z[k]; 
 
subject to W2: sum {j in K} w[j,0] >=1; 
 
subject to SIRConst {l in L}:  
 sum{m in M, j in C[m]} G[m,l,j]*x[m,j] <= 1+1/SIRmin+(1-y[l])*beta[l];  
  
subject to NewSIRConst {l in L}:   
 sum{m in M: l in C[m]} x[m,l] <=  1+1/SIRmin;  
  
subject to BuildTowers {m in M,l in C[m]}:  
  x[m,l] <= dd[m]*y[l];  
  
subject to Demand {m in M}:  
  sum{l in C[m]} x[m,l] <= d[m];  
  
subject to newcuts {m in M, l in C[m], j in C[m]: g[m,l] < g[m,j]}:  
   x[m,l] <= d[m]*(1 - y[j]);  
  
# Connect each base station to an MTSO  
subject to BSMTSO {l in L}: sum {k in K} s[l,k] = y[l];  
  
subject to s_and_z {l in L, k in K}: s[l,k] <= z[k];  
  
# Capacity constraint for MTSOs  
subject to towers_per_MTSO {k in K}: sum{l in L} s[l,k] <= alpha * z[k];  
  
subject to flow_balance {k in K}:  
sum{j in K diff {k}} (u[k,j] - u[j,k]) + u[k,0] = z[k];  
  
subject to sink: sum{k in K} u[k,0] = sum{k in K} z[k];  
 
# Always use the dummy node  
subject to use_dummy_node: z[0] = 1;  
  
subject to u_and_w {j in K,k in K0 diff {j}}: u[j,k]<=card(K)*w[j,k];     
 
###### Objective ######  
  
maximize Profit: r*sum{m in M, l in C[m]} x[m,l] -    
                   sum{l in L} a[l]*y[l] -  
                   sum{k in K} b[k]*z[k] -  
                   sum{l in L, k in K} c[l,k] * s[l,k] -  
                   sum{j in K, k in K0 diff {j}} h[j,k]* w[j,k];  
  
#### Problem definitions ###### 
problem UB: 
y,x,z,s,u,w, 
W0,W2,SIRConst,NewSIRConst,BuildTowers,Demand, 
newcuts,BSMTSO,s_and_z,towers_per_MTSO,flow_balance, 
sink,use_dummy_node,u_and_w,Profit; 
 
problem LB: 
y,x,z,s,u,w, 
yy,xx, 
IP1,IP2, 
W0,W2,SIRConst,NewSIRConst,BuildTowers,Demand, 
newcuts,BSMTSO,s_and_z,towers_per_MTSO,flow_balance, 
sink,use_dummy_node,u_and_w,Profit; 
 
 
###### Calculated Data ######  
  
param lhs {L};       # left-hand side of SIRConst[l]  
param rhs {L};       # right-hand side of SIRConst[l]  
  
param inf_tol default 0.000001;  
param error default 0;  
param max_error default 0;  
param max_lhs default 0;  
param max_rhs default 0;  
param most_violated_SIR default 0;  
  
param bsx{L}; # (bsx[l],bsy[l]) denotes the (x,y) co-ordinate  
param bsy{L};   # for base station l (units are kilometers)  
  
  
param max_g_ratio default 0;  
  
param N default 660;  # total number of channels available  
                      # required for successful operation  
  
param f default 2000; # frequency  
param Hb default 10;  # height of base station  
param Hm default 1;   # height of mobile station  
  
  
param sx{M};       # (sx[m],sy[m]) denotes the (x,y) co-ordinate  
param sy{M};       # for subscriber m  
  
param mtsox{K0};   # (mtsox[k],mtsoy[k]) denotes the (x,y) co-ordinate  
param mtsoy{K0};   # for MTSO k  
  
param dist{M,L};      # dist[m,l] denotes the distance from m to l  
param X;  
param Y;  
param DemandSatisfied;  
param Totaly;  
#display f;  display Hb;  display Hm;  
param logf  := log10(f);  
param logHb := log10(Hb);  
param Const1 := 69.55 + 26.16*logf - 13.82*logHb;  
param CF     := (1.1*logf - 0.7)*Hm - 1.56*logf + 0.8;  
param Const2 := 44.9 - 6.55*logHb;  
param Au;  
param ColSum;  
  
param corrected default 0; # True if the solution was corrected by the  
                           # post-processing routine.  
param mip1_profit default 0;  
param mip2_profit default 0;  
  
################################ Code #########################################  
set GG; # the set of all generated data 
set GL; # the set of generated tower locations 
set GM; # the set of generated subscriber locations 
set GK; # the set of generated MTSO locations 
set GK0; # the union of GK and gateway MTSO 0 
 
set known_w within {K,K0}; # Number of fixed w before solving 
 
let GK:={1..MTSO}; 
let GK0:={0} union GK; 
let GL:= {MTSO+1..MTSO+TOWER}; 
let GM:= {MTSO+TOWER+1..MTSO+TOWER+SUBSCRIBER}; 
let GG:= GK0 union GL union GM; 
#display GK0,GL,GM; 
 
param gx {GG} default 0; 
param gy {GG} default 0; 
 
param besty {L}; 
param bestx {m in M,C[m]}; 
param bestz {K0}; 
param bests {L,K}; 
param bestu {K,K0}; 
param bestw {K,K0}; 
 
# Initialize random number sequence 
param U; 
for {i in {1..SEED}}{ 
 let U:=Uniform(0.0,100.0); 
} 
  
param use_branching_rule default 1;  
  
# The following parameter determines whether or not we use constraint set   
# (23) from "Base Station Location and Service Assingments in W-CDMA Networks" by Kalvenes, Kennington, and Olinick  
 
param use_new_cuts default 0;  
  
# The following parameter determines whether or not we use constraint set   
# (24) from "Base Station Location and Service Assingments in W-CDMA Networks" by Kalvenes, Kennington, and Olinick  
param use_NewSIR default 1;  
  
#printf "|L| = %d\n", card(L);  
#printf "|M| = %d\n", card(M);  
#printf "|K| = %d\n", card(K);  
let alpha := min(alpha, card(L));
#printf "alpha = %d\n\n",alpha;  
  
# Reset C[m] to be empty  
let {m in M} C[m] := {1..TOWER};  
  
 
# Generate random MTS0 locations  
param cum; # cumulative counter 
let cum:=0; 
param radius; 
param angle; 
let angle:= 2*3.1415926/MTSO; 
param angle1; 
param angle2; 
param angle3; 
param counter; 
let counter:=-1; 
 
for {gk in GK0}{ 
   #display gk; 
   let counter:=counter+1; 
 
   if counter==0 then { 
     let gx[gk]:= Uniform(6.0,7.5); 
     let gy[gk]:= Uniform(3.75,4.75); 
     let mtsox[gk]:= gx[gk]; 
     let mtsoy[gk]:= gy[gk]; 
   }#if 
   else { 
     let angle1:= 0 + (counter-1)*angle; 
     let angle2:= angle1 + angle; 
     #display angle1,angle2; 
     let angle3:= Uniform(angle1+0.3,angle2-0.3); 
     let radius:= Uniform(1.80,4.00); 
     #display angle3,radius; 
     #display cos(angle3),sin(angle3); 
     let gx[gk] := mtsox[0]+radius*cos(angle3); 
     let gy[gk] := mtsoy[0]+radius*sin(angle3); 
     let mtsox[gk]:= gx[gk]; 
     let mtsoy[gk]:= gy[gk]; 
   } # else 
}# for 
 
#for {k in K0} {  
#   printf "MTSO %d: x = %f\ty = %f\n",k,mtsox[k],mtsoy[k];  
#} # for {k in K0}  
#printf "\n\n";  
 
# Generate random base station locations  
#printf "\n\n"; 
let cum:=0; 
 
for {gl in GL}{ 
   #display gl; 
   let gx[gl] := Uniform(0.0,13.5); 
   let gy[gl] := Uniform(0.0,8.5); 
  
   repeat until (cum == gl){ 
    let cum:= 0; 
    for {gg in GG: gg < gl }{ 
      #display gg; 
      let X:= gx[gg] - gx[gl]; 
      let Y:= gy[gg] - gy[gl]; 
      if (sqrt(X*X+Y*Y) >= (min_dist*1.609)) then { 
       let cum:= cum + 1; 
      } # if then 
      else { 
       #printf "\nregenerating ..\n"; 
       #display gl,gx[gl],gy[gl]; 
       #display gg,gx[gg],gx[gg]; 
       let cum:=0; 
       let gx[gl] := Uniform(0.0,13.5); 
       let gy[gl] := Uniform(0.0,8.5); 
       #display gx[gl],gy[gl]; 
       break; 
      }#else 
    }# for {gg 
   }# repeat  
   let bsx[gl-MTSO] := gx[gl];  
   let bsy[gl-MTSO] := gy[gl];  
   #printf "bsx[%d] = %f\n",gl-MTSO,gx[gl]; 
   #printf "bsy[%d] = %f\n",gl-MTSO,gy[gl]; 
}# for {gl 
 
#for {l in L}{ 
#   printf "Base Station %d: x = %f\ty = %f\n",l,bsx[l],bsy[l];  
#} # for {l in L}  
#printf "\n\n";  
  
# Gather distribution density of towers 
param towers_q1 default 0; # number of towers in 1st quadrant 
param towers_q2 default 0; # number of towers in 2nd quadrant 
param towers_q3 default 0; # number of towers in 3rd quadrant 
param towers_q4 default 0; # number of towers in 4th quadrant 
param Prob; 
 
let towers_q1:= card({l in L: bsx[l] > mtsox[0] and bsy[l] > mtsoy[0]}); 
let towers_q2:= card({l in L: bsx[l] < mtsox[0] and bsy[l] > mtsoy[0]}); 
let towers_q3:= card({l in L: bsx[l] < mtsox[0] and bsy[l] < mtsoy[0]}); 
let towers_q4:= card({l in L: bsx[l] > mtsox[0] and bsy[l] < mtsoy[0]}); 
 
# Generate random subscriber data  
 
#printf "\n\n"; 
 
for {m in M}{ 
   let sx[m] := Uniform(0,13.5); 
   let sy[m] := Uniform(0,8.5); 
}# for {m 
 
for {m in M} {  
   let d[m] := floor(Uniform(dLo,dHi));  
#   printf "Subscriber %d: x = %f\ty = %f\td = %d\n",m,sx[m],sy[m],d[m];  
} # for {m in M}  
let TotalDemand := sum{m in M} d[m];  
#printf "\n\n";  
printf"TotalDemand = %lf\n",TotalDemand;
 
for {m in M} { 
   let dd[m]:= min(d[m],1+1/SIRmin); 
}  
#display dd; 
 
for {m in M, l in L} {  
   let X         := bsx[l] - sx[m];  
   let Y         := bsy[l] - sy[m];  
   let dist[m,l] := sqrt(X*X + Y*Y);  
} # for {m in M, l in L}  
 
# Reference for Hata path-loss models  
# L. Hanzo, P. Cherriman, J. Streit. 2001. Wireless Video   
# Communications. IEEE Press. The Institute of Electrical   
# and Electronic Engineers, Inc. NY, NY.  
  
for {m in M, l in L} {  
   let Au := Const1 - CF + Const2*log10(dist[m,l]);  
   let Au := Au/10;  
   let g[m,l] := 1.0/(10.0**Au);  
   if g[m,l] < 1.0e-15 then  
     let C[m] := C[m] diff {l};  
} # for {m in M, l in L}  
   
#display C; 
param AveCardOfC;
let   AveCardOfC := (sum{m in M} card(C[m]))/card(M);
printf"AveCardOfC = %lf\n",AveCardOfC;

let {m in M, l in L, j in C[m]} G[m,l,j]:= g[m,l]/g[m,j]; 
 
# Calculate the base station to MTSO connection costs.  
# Let c[l,k] be proportional to the distance from l to k  
for {l in L, k in K} {  
   let X         := bsx[l] - mtsox[k];  
   let Y         := bsy[l] - mtsoy[k];  
   let c[l,k]    := cperfoot*(3281*sqrt(X*X + Y*Y)); # 1 km = 3281 ft  
} # for {l in L, k in K}  
#display c;  
  
# Calculate the MTSO to MTSO connection costs.  
# Let h[k1,k2] be proportional to the distance from k1 to k2  
for {k1 in K, k2 in K0} {  
   let X         := mtsox[k1] - mtsox[k2];  
   let Y         := mtsoy[k1] - mtsoy[k2];  
   let h[k1,k2]  := hperfoot*(3281*sqrt(X*X + Y*Y)); # 1 km = 3281 ft  
} # {k1 in K, k2 in K0}  
#display h;  
  
# Calculate the Tower construction costs. 
for {l in L} { 
   let a[l]:= Uniform(aLo,aHi); 
} 
#display a;   
 
# Calculate the MTSO construction costs. 
for {k in K} { 
   let b[k]:= Uniform(bLo,bHi); 
} 
#display b; 
 
 
#param cplex_time default 0; 
 
problem UB; 
option cplex_options 'timing =1 timelimit=7200'; 
printf"\n\nUpper Bound\n\n";
solve; 
#let cplex_time:= cplex_time + _solve_user_time; 
let Totaly := sum {l in L} y[l]; 
let DemandSatisfied := 100.*(sum {m in M, l in C[m]} x[m,l])/TotalDemand; 
#display y; 
#display z; 
#display x; 
#for {m in M} {
#  for {l in C[m]} {
#    if x[m,l] > 0.0 then printf"x[%d,%d] = %lf\n",m,l,x[m,l];
#  }
#}
#display s; 
#for {l in L} {
#  for {k in K} {
#    if s[l,k] > 0.0 then printf"s[%d,%d] = %lf\n",l,k,s[l,k];
#  }
#}
#display u; 
#display w; 
printf"\nNumber of MTSOs Selected         = %lf\n",sum {k in K} z[k]; 
printf"\nNumber Of Base Stations Selected = %lf\n",Totaly; 
#printf"\nTotal Demand                     = %lf\n",TotalDemand; 
#printf"\nDemand Serviced                  = %lf\n",sum {m in M, l in C[m]} x[m,l]; 
printf"\nDemand Satisfied                 = %7.2lf\%\n\n",DemandSatisfied; 
display Profit; 
#printf "\nTotal CPLEX time = %f\n\n",cplex_time; 
param UBProfit;
let   UBProfit := Profit;
 
 
#let cplex_time:= 0; 
problem LB; 
option cplex_options 'timing =1 timelimit=7200 mipgap=0.05'; 
printf"\n\nLower Bound\n\n";
solve; 
#let cplex_time:= cplex_time + _solve_user_time; 
let Totaly := sum {l in L} y[l]; 
let DemandSatisfied := 100.*(sum {m in M, l in C[m]} x[m,l])/TotalDemand; 
#display y; 
#display z; 
##display x; 
#for {m in M} {
#  for {l in C[m]} {
#    if x[m,l] > 0.0 then printf"x[%d,%d] = %lf\n",m,l,x[m,l];
#  }
#}
#display s; 
#for {l in L} {
#  for {k in K} {
#    if s[l,k] > 0.0 then printf"s[%d,%d] = %lf\n",l,k,s[l,k];
#  }
#}
#display u; 
#display w; 
printf"\nNumber of MTSOs Selected         = %d\n",sum {k in K} z[k]; 
printf"\nNumber Of Base Stations Selected = %7d\n",Totaly; 
#printf"\nTotal Demand                     = %f\n",TotalDemand; 
#printf"\nDemand Serviced                  = %f\n",sum {m in M, l in C[m]} x[m,l]; 
printf"\nDemand Satisfied                 = %7.2lf\%\n\n",DemandSatisfied; 
display Profit; 
param   OptGap;
let     OptGap := 100*(UBProfit-Profit)/Profit;
printf"Optimality Gap = %7.2lf\n\n",OptGap;
#printf "\nTotal CPLEX time = %f\n",cplex_time; 
#include draw.txt; 
quit;