 FICO Xpress Optimization Examples Repository
 FICO Optimization Community FICO Xpress Optimization Home   Outer approximation for quadratic facility ocation: solving different problem types in sequence

Description
A quadratic facility location problem is formulated and solved as an MIQP problem (quadfacloc.mos) and via an outer approximation algorithms that iterates over MIP and QP subproblems (quadfaclocoa.mos).

Source Files

```(!*********************************************************************
Mosel NL examples
=================
`````````````````````

Model of a quadratic facility location problem solved
by outer approximation (master modeL: MIP, submodel: 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
*********************************************************************!)

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 master
LastMast: real               ! Latest obj solution value for master problem
LastX: array(I,J) of real    ! Latest 'x' solution values from subproblem
LastSub: real                ! Latest obj solution value for subproblem
end-declarations

! **** Master model (MIP)
declarations
Master: mpproblem
end-declarations

with Master 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
Fixings: array(I) of linctr
end-declarations

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 Master do
minimize(LinCost)
! writeprob('Master_'+Iter,'l')
if getprobstat<>XPRS_OPT then break; end-if
forall(i in I) LastPoint(i) := getsol(z(i))
LastMast := getobjval
end-do

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

with Master do
! After the first iteration add the linearization term to the objective
if Iter=1 then LinCost+= eta; end-if
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 := LastMast
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

```   