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

   Model of a quadratic facility location problem solved
   by outer approximation (main problem: MIP, subproblem: QP).
   The goal is to assign facilities to customers to minimize cost. 
   
   Original AMPL model by Sven Leyffer, see
   https://www.researchgate.net/publication/2827082_Solving_Mixed_Integer_Nonlinear_Programs_by_Outer_Approximation
 
   (c) 2018 Fair Isaac Corporation
       author: Z.Csizmadia, Sven Leyffer, Aug. 2018, rev. Apr.2021
*********************************************************************!)

model QuadFacLocOA
 uses "mmxprs", "mmnl", "mmsystem"

 parameters
   ILIM = 100
 end-parameters

 forward procedure printsol

 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
   UBD,LBD: real                 ! Upper/lower bound on solution
 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 ****
 declarations
   eta: mpvar                   ! Linearization of quad. objective terms
   x: array(I,J) of mpvar       ! Percentage of customer j served by facility j
   z: array(I) of mpvar         ! =1, iff facility i built
   LastPoint: array(I) of real  ! Latest discrete solution values from main
   LastMain: real               ! Latest obj solution value for main problem
   LastX: array(I,J) of real    ! Latest 'x' solution values from subproblem
   LastSub: real                ! Latest obj solution value for subproblem
 end-declarations

 ! **** Main problem (MIP)
 declarations
   MainProb: mpproblem
 end-declarations

 with MainProb do
   declarations
     LinCost: linctr
   end-declarations

   eta is_free

  ! Linearized objective function
   LinCost:= sum(i in I) C(i)*z(i)
! After the initial iteration, the objective is changed to this form:
!   LinCost:= sum(i in I) C(i)*z(i) + eta

  ! 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)
 end-do

 ! **** Submodel (continuous QP)
 declarations
   SubNLP: mpproblem
 end-declarations

 with SubNLP do
   declarations
     QuadCost: nlctr
     Fixings: array(I) of linctr
   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

!   setparam("xprs_verbose",true)
 end-do

 ! **** Outer approximation algorithm
 UBD := 1E20; LBD := -1E20
 Iter:= 0; starttime:=gettime
 while (LBD < UBD-(SCALE*1E-2) and Iter <= ILIM) do
   Iter += 1
   with MainProb do
     minimize(LinCost)
     ! writeprob('MainProb_'+Iter,'l')
     if getprobstat<>XPRS_OPT then break; end-if
     forall(i in I) LastPoint(i) := getsol(z(i))
     LastMain := getobjval
   end-do

   with SubNLP do
     forall(i in I) Fixings(i):= z(i) = LastPoint(i)  ! Fix discrete variables
     minimize(QuadCost)
     ! writeprob('SubNLP_'+Iter,'l')
     forall(i in I, j in J) LastX(i,j) := getsol(x(i,j))
     LastSub := getobjval
   end-do

   with MainProb do
    ! After the first iteration add the linearization term to the objective
     if Iter=1 then LinCost+= eta; end-if
    ! Linearization of quadratic part
     eta >= 
       sum(i in I, j in J) Q(i,j)*LastX(i,j)^2 +    ! Objective value
       sum(i in I, j in J) 2*Q(i,j)*LastX(i,j)*(x(i,j) - LastX(i,j))  ! Gradient
   end-do

   UBD := minlist( UBD, LastSub )
   LBD := LastMain
   writeln(gettime-starttime, "sec Iter = ", Iter, ": LBD = ", LBD, "  UBD = ", UBD)
! Uncomment this line to display intermediate solutions in detail:
!   printsol
 end-do
 if Iter > ILIM then writeln("Iteration limit reached"); end-if

 printsol

 ! **** Display the solution ****
 procedure printsol
   write("Open         Customer locations\nfacility")
   forall(j in J) write(formattext("%6d",j))
   forall(i in I | LastPoint(i)>0) do
     write("\n  ", textfmt(i,-8))
     forall(j in J) write(formattext("%5.1f%%",LastX(i,j).sol*100))
   end-do
   writeln
 end-procedure

end-model