Facility location model Q&A

Dear Erwin, 

I always follow and read your interesting posts. Thank you so much for your useful posts on your blogs. For the learning purposes, I was trying to code you MIQCP sample on your blog but not successful.  (http://yetanothermathprogrammingconsultant.blogspot.com.au/2018/01/solving-facility-location-problem-as.html). 

I have checked the code, Lp, and your post 10,000 times. Everything seems to be exactly like what you have addressed however not sure why Gams returns the error as 100 nonlinear lines. I see you are too busy and not work as GAMS support. But since it might be a coding problem with your post (or perhaps something is missed), It would be highly appreciated if you could do a favor and double check that. Much obliged.  

The attached GAMS model looks like:

set i /i1*i50/
set j /facility1,facility10/

parameter maxdist /40/;

set c /x,y/

Table dloc(i,c)
           x         y
i1         65        84
i2         8         72
i3         57        28
i4         70        79
i5         82        80
i6         87        95
i7         14        98
i8         28        6
i9         62        22
i10        22        80
i11        31        94
i12        87        74
i13        80        32
i14        92        48
i15        22        74
i16        0         42
i17        80        90
i18        100       56
i19        73        1
i20        82        84
i21        7         15
i22        3         19
i23        99        89
i24        56        94
i25        18        35
i26        8         30
i27        71        19
i28        10        66
i29        7         71
i30        13        66
i31        68        14
i32        95        43
i33        23        50
i34        95        43
i35        4         63
i36        51        13
i37        56        26
i38        41        56
i39        75        72
i40        41        68
i41        57        91
i42        16        63
i43        75        82
i44        44        94
i45        79        32
i46        72        5
i47        4         52
i48        16        20
i49        31        60
i50        92        63
;

scalar M /1000000/

variable
   floc(j,c)
   n

binary variables
  isopen(j)
  assign(i,j)

equations
   distance,assigndemand,closed,numfacilities,order;

distance(i,j)..          (sum(c,dloc(i,c)- floc(j,c)))**2 =l= Maxdist**2 + M*(1-assign(i,j));
assigndemand(i)..         Sum(j,assign(i,j)) =e= 1;
closed(i,j)..             assign(i,j) =l= isopen(j) ;
Numfacilities..           n =e= sum(j,isopen(j));
order(j+1)..              isopen(j) =g= isopen(j+1);

model location /all/

option limcol = 180;
option limrow =180;

solve location minimizing n using MIQCP;

display n.l,floc.l;

I think this GAMS model as written has a few problems that are interestingly enough to point out.

What is this about?

In [1] a simple facility location problem was discussed. It comprised of two models:

1. Find the number of facilities \(n\) needed, so all demand locations are within MAXDIST kilometers from the closest facility.
2. Given that we know \(n\), place these \(n\) facilities such that some total distance measure is minimized. 

The first model is presented as:

What is wrong with the GAMS model in the e-mail?

Quite a few things. Let me try to review them:

1. The set \(j\) only has two elements. This is not sufficient. We should have more elements in \(j\) than the optimal value of \(n\).  The optimal value for this data set is \(n^* = 3\), so just having two elements in \(j\) will cause the model to be infeasible.
2. The value for \(M\) needs to get much more attention. The value needs to be large enough that the distance equation becomes non-binding when \(\mathit{assign}_{i,j}=0\). A simple value would be \[M=\sum_c \left(\max_i \{ \mathit{dloc}_{i,c}\} -\min_i \{\mathit{dloc}_{i,c}\}\right)^2 \] This is the square of the length of the diagonal of the box containing all demand locations. For this data this leads to \(M=19409\). 
   We can assume facilities will be placed in the convex hull of the demand locations (sorry, I don't know how to prove that). That means we can use the maximum pairwise squared distance between demand locations \((i,i')\) as a good value for \(M\): \[M= \max_{i,i'} \sum_c \left(\mathit{dloc}_{i,c}-\mathit{dloc}_{i',c}\right)^2\] This sets \(M=14116\).
   With some effort better values can be derived by using individual values \(M_{i,j}\) instead of one \(M\). Big-M values with some computation and just having assigned some large value, are always a point of attention.
   Using just very big values for big-\(M\) constants, is something I see very often. This can lead to all kind of problems when solving the MIP problem. 
3. The expression (sum(c,dloc(i,c)- floc(j,c)))**2 takes the square of the sum. This is not the same as a sum of squares:\[\sum_i x_i^2 \ne \left(\sum_i x_i\right)^2\] This expression is just malformed.
4. In GAMS x**y is a general non-linear function, and x**2 is not recognized as quadratic function (GAMS could be much smarter here). This will make the model an MINLP model and not a MIQCP. This will lead to the somewhat cryptic error:
   
     --- Executing SBB: elapsed 0:00:00.140
     Could not load data
     Detected 100 general nonlinear rows in model
   
   Of course you can solve it with an MINLP solver. (Confusingly SBB is an MINLP solver).  Of course it is better to use the sqr(x) function to help GAMS. The left-hand side of the distance constraint should be written as: sum(c,sqr(dloc(i,c)- floc(j,c))).
   Side note: if we would try to solve the original expression (sum(c,dloc(i,c)- floc(j,c)))**2 as a MINLP model, we would still be in a heap of problems. The function x**y assumes \(x\ge 0\). This is obviously not the case here. We would see a ton of domain errors. 
5. For an optimal solution, you will need to add: option optcr=0; GAMS has a default relative gap tolerance of 10%. With this options we set this tolerance to zero.

Modeling is not that easy! We need to pay attention to lots of details to get models working correctly.