FICO Xpress Optimization Examples Repository
 FICO Optimization Community FICO Xpress Optimization Home

Outer approximation for quadratic facility location: 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 algorithm that iterates over MIP and QP subproblems (quadfaclocoa.mos).

Source Files
By clicking on a file name, a preview is opened at the bottom of this page.

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


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
*********************************************************************!)

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
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 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
! 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
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

`