(!*********************************************************************
   Mosel NL examples
   =================
   file quadfacloc.mos
   ```````````````````
   Quadratic Facility Location model

   Model of quadratic facility location problem. 
   The goal is to assign facilities to customers to minimize cost. 
   
   Basic version (formulation as MIQP)

   (c) 2018 Fair Isaac Corporation
       author: Z.Csizmadia, Sven Leyffer, Aug. 2018
*********************************************************************!)

model QuadFacLoc
 uses "mmxprs", "mmnl", "mmsystem"

 declarations
   Iter: integer
 end-declarations

 declarations
   I = 1..5                      ! Set of facilities
   J = 1..10                     ! Set of customers

   FacX,FacY: array(I) of real   ! X/Y coordinates of facility location
   CustX,CustY: array(J) of real ! X/Y coordinates of customer location
   C: array(I) of real           ! Fixed cost of opening facility at location
   Q: array(I,J) of real         ! Quadr. cost of meeting demand j from facility i
 end-declarations

 ! **** Generate random data:
 ! **** Cost coefficients and coordinates of all facility + customer locations
 setrandseed(42)
 SCALE:=100
 forall(i in I) do
   C(i) := SCALE*100*random
   FacX(i) := SCALE*random
   FacY(i) := SCALE*random
 end-do

 forall(j in J) do
   CustX(j) := SCALE*random
   CustY(j) := SCALE*random
 end-do

 forall(i in I, j in J)        ! Formula used by Lee to create instances
   Q(i,j) := 50 * sqrt( (FacX(i) - CustX(j))*(FacX(i) - CustX(j)) + 
                        (FacY(i) - CustY(j))*(FacY(i) - CustY(j)))

 ! **** Optimization model formulation (MIQP) ****
 declarations
   x: array(I,J) of mpvar       ! Percentage of customer j served by facility i
   z: array(I) of mpvar         ! =1, iff facility i built
   QuadCost: nlctr

   BasicNodes: integer
 end-declarations

  ! Quadratic objective function
  QuadCost:= sum(i in I) C(i)*z(i) + sum(i in I, j in J) Q(i,j)*x(i,j)^2;

  ! Variable upper bound (only serve customers from open facilities)
  forall(i in I, j in J) x(i,j) <= z(i)

  ! Meet all customers' demands
  forall(j in J) sum(i in I) x(i,j) = 1

  forall(i in I) z(i) is_binary

  setparam("xprs_verbose",true) 
  starttime:=gettime
  minimize(QuadCost)
  BasicNodes:= getparam("XPRS_NODES")
  writeln(gettime-starttime, "sec  Solution: ", getobjval, ". Basic formulation - nodes: ", BasicNodes)
  
 ! Display the solution
  write("Open         Customer locations\nfacility")
  forall(j in J) write(formattext("%6d",j))
  forall(i in I | z(i).sol>0) do
    write("\n  ", textfmt(i,-8))
    forall(j in J) write(formattext("%5.1f%%",x(i,j).sol*100))
  end-do
  writeln

end-model