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