# Euclidean single facility location problem

# Objective Function:  convex
# Constraint Functions: none

var eps >= 1.0e-8;

param m := 200;		# number of existing facilities
param n1 := 5;
param n2 := 5;
param n := n1*n2;		# number of new facilities

param a {1..m, 1..2};   # coordinates of existing facility
param w {1..m, 1..n};	# weights associated with old-new connections
param v {1..n, 1..n};	# weights associated with new-new connections

var x {1..n, 1..2};

minimize sumEucl: 
  eps + 
  sum {i in 1..m, j in 1..n} 
      w[i,j]*sqrt( eps^2 + sum {k in 1..2} (x[j,k] - a[i,k])^2 )
  +
  sum {j in 1..n, jj in 1..n: j<jj}  
      v[j,jj]*sqrt( eps^2 + sum {k in 1..2} (x[j,k] - x[jj,k])^2 ) ;

let eps := 1;

let {i in 1..m, k in 1..2}          a[i,k]  := Uniform01();
let {j in 1..n, jj in 1..n: j < jj} v[j,jj] := 0.2;

let {j1 in 1..n1, j2 in 1..n2} x[j1+n1*(j2-1),1] := (j1-0.5)/n1;
let {j1 in 1..n1, j2 in 1..n2} x[j1+n1*(j2-1),2] := (j2-0.5)/n2;

let {i in 1..m, j in 1..n}          
    w[i,j] := (if abs(a[i,1]-x[j,1]) <= 1/(2*n1) 
	       && abs(a[i,2]-x[j,2]) <= 1/(2*n2) 
	       then 0.95 else 
		   (if abs(a[i,1]-x[j,1]) <= 2/(2*n1)
		    && abs(a[i,2]-x[j,2]) <= 2/(2*n2)
		    then 0.05 else 0
		   )
		);

solve;

display x;

let eps := 0;
display sumEucl;

printf {i in 1..m, j in 1..n: w[i,j] > 0.5}: "%10f %10f \n%10f %10f \n\n",
    a[i,1], a[i,2], x[j,1], x[j,2] > emfl_vareps.out1;

printf {j in 1..n, jj in 1..n: j < jj}: "%10f %10f \n%10f %10f \n\n",
    x[j,1], x[j,2], x[jj,1], x[jj,2] > emfl_vareps.out2;