| |||||||||

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 quadfaclocoa.mos (!********************************************************************* 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 | |||||||||

© Copyright 2021 Fair Isaac Corporation. |