(!******************************************************
Mosel Example Problems
======================
file tspsub.mos
```````````````
Subtour elimination algorithm from f5tour2.mos
(Applications of optimization with Xpress,
Problem 9.6: Planning a flight tour,
second formulation for other data format)
- A subset of nodes may be fixed via the input data -
(c) 2011 Fair Isaac Corporation
author: S. Heipcke, Oct. 2011, rev. Aug 2018
*******************************************************!)
model "TSP"
uses "mmxprs", "mmsystem"
parameters
NUM=1 ! Model instance
PB=0 ! Optimization problem index
LEVEL=1 ! =1: first round, >1: rounds with partial unfix
BIN = "" ! File format (binary)
RMT = "" ! Whether to use remote files
ZIP = "" ! File format (compression)
MAXUF = 8 ! Max. number of nodes to unfix in each subproblem
PRINTDETAIL = true ! Whether to display result detail
end-parameters
writeln("Sub ", NUM, ": ", PB, ", level: ", LEVEL)
forward procedure break_subtour
forward procedure print_sol
forward procedure two_opt
declarations
starttime: real ! Time measure
NODES: range ! Set of nodes
Squares: range ! Subproblems
NODEX,NODEY: array(NODES) of integer ! X/Y coordinates
SNODES: array(NSq:range) of set of integer ! Nodes per subproblem
LEV: array(Squares) of integer ! Iteration round
SBORDER: array(Squares) of integer ! Offset from border
BTYPE: array(Squares) of string ! X:horizontal, Y:vertical neighbors
SQSETS: array(Squares,1..2) of set of integer ! Predecessor subproblems
PREC: array(Squares) of set of integer ! Predecessor subproblems
NodeSet: set of integer ! Subproblem node set
UnfixedSet: set of integer ! Unfixed nodes in this subproblem
SOL: dynamic array(NODES) of integer ! Dynamic to write index
ATemp1, AIndx1: dynamic array(R:set of integer) of integer ! Aux. node data
ATemp2, AIndx2: dynamic array(R2:set of integer) of integer ! Aux. node data
r1,r2: set of integer ! Auxiliary node sets
end-declarations
starttime:=gettime
initializations from BIN+ZIP+RMT+"tspgendata"
Squares NSq
NODEX
NODEY
SNODES
LEV
SBORDER
BTYPE
SQSETS
PREC
end-initializations
initializations from BIN+ZIP+RMT+"tspsolinit"+NUM
SOL
end-initializations
fdelete(RMT+"tspsolinit"+NUM)
if LEV(PB)=1 then
NodeSet:=SNODES(PB)
else
NodeSet:= union(s in SQSETS(PB,1), n in SNODES(s)) {n} +
union(s in SQSETS(PB,2), n in SNODES(s)) {n}
end-if
finalize(NodeSet)
if LEV(PB)>1 then
if BTYPE(PB)="X" then
forall(s in SQSETS(PB,1), n in SNODES(s))
ATemp1(n):=abs(NODEX(n)-SBORDER(PB))
r1:=union(i in R | exists(ATemp1(i))) {i}
qsort(SYS_UP, ATemp1, AIndx1, r1)
ct:=0
forall(ct as counter, i in min(j in r1) j .. max(k in r1) k |
exists(AIndx1(i)) ) do
UnfixedSet+= {AIndx1(i)}
if ct>=MAXUF then break; end-if
end-do
forall(s in SQSETS(PB,2), n in SNODES(s))
ATemp2(n):=abs(NODEX(n)-SBORDER(PB))
r2:=union(i in R2 | exists(ATemp2(i))) {i}
qsort(SYS_UP, ATemp2, AIndx2, r2)
ct:=0
forall(ct as counter, i in min(j in r2) j .. max(k in r2) k |
exists(AIndx2(i)) ) do
UnfixedSet+= {AIndx2(i)}
if ct>=MAXUF then break; end-if
end-do
else
forall(s in SQSETS(PB,1), n in SNODES(s))
ATemp1(n):=abs(NODEY(n)-SBORDER(PB))
r1:=union(i in R | exists(ATemp1(i))) {i}
qsort(SYS_UP, ATemp1, AIndx1, r1)
ct:=0
forall(ct as counter, i in min(j in r1) j .. max(k in r1) k |
exists(AIndx1(i))) do
UnfixedSet+= {AIndx1(i)}
if ct>=MAXUF then break; end-if
end-do
forall(s in SQSETS(PB,2), n in SNODES(s))
ATemp2(n):=abs(NODEY(n)-SBORDER(PB))
r2:=union(i in R2 | exists(ATemp2(i))) {i}
qsort(SYS_UP, ATemp2, AIndx2, r2)
ct:=0
forall(ct as counter, i in min(j in r2) j .. max(k in r2) k |
exists(AIndx2(i)) ) do
UnfixedSet+= {AIndx2(i)}
if ct>=MAXUF then break; end-if
end-do
end-if
else
UnfixedSet:=NodeSet
end-if
MinN:=min(n in UnfixedSet) n
FixedSet:=NodeSet-UnfixedSet
! writeln("Unfixed: ", UnfixedSet, " Fixed: ", FixedSet)
! forall(n in NodeSet) if SOL(n)=0 then writeln(PB, " =0: ", n); end-if
declarations
NNODES=getsize(NodeSet) ! Number of subproblem nodes
DIST: array(NodeSet,NodeSet) of real ! Distance between cities
NEXTC: array(NodeSet) of integer ! Next city after i in the solution
fly: array(NodeSet,NodeSet) of mpvar ! 1 if flight from i to j
end-declarations
forall(i in NodeSet, j in NodeSet)
DIST(i,j):=sqrt((NODEX(i)-NODEX(j))^2 + (NODEY(i)-NODEY(j))^2)
! DIST(i,j):=round(10*(sqrt((NODEX(i)-NODEX(j))^2 + (NODEY(i)-NODEY(j))^2)))/10
! Objective: total distance
TotalDist:= sum(i,j in NodeSet | i<>j) DIST(i,j)*fly(i,j)
! Visit every city once
forall(i in NodeSet) sum(j in NodeSet | i<>j) fly(i,j) = 1
forall(j in NodeSet) sum(i in NodeSet | i<>j) fly(i,j) = 1
forall(i,j in NodeSet | i<>j) fly(i,j) is_binary
! Fix part of the variables
if LEV(PB)>1 then
forall(i in FixedSet | SOL(i) not in UnfixedSet) do
fly(i,SOL(i)) = 1
end-do
end-if
setparam("XPRS_THREADS", 1)
! Solve the problem
minimize(TotalDist)
! Eliminate subtours
break_subtour
! Perform 2-opt after solving partially fixed problems
if LEVEL>1 then two_opt; end-if
! Save solution
forall(n in NodeSet) SOL(n):=NEXTC(n)
if getprobstat<>XPRS_OPT then
writeln("(", gettime-starttime, "s) ", NUM, "-", PB, ": INFEAS ")
else
write("(", gettime-starttime, "s) ", NUM, "-", PB, ": Total distance: ",
getobjval)
if(PRINTDETAIL) then
writeln(" ", NEXTC)
else
writeln
end-if
end-if
initializations to BIN+ZIP+RMT+"tspsol"+NUM
PB
evaluation of array(n in NodeSet) SOL(n) as "SOL"
end-initializations
!-----------------------------------------------------------------
procedure break_subtour
declarations
TOUR,SMALLEST,ALLCITIES: set of integer
TOL=10E-10
end-declarations
forall(i in NodeSet)
forall(j in NodeSet)
if(getsol(fly(i,j))>=1-TOL) then
NEXTC(i):=j
break
end-if
! Get (sub)tour containing city 1
TOUR:={}
first:=MinN
repeat
TOUR+={first}
first:=NEXTC(first)
until first=MinN
size:=getsize(TOUR)
! Find smallest subtour
if size < NNODES then
SMALLEST:=TOUR
if size>2 then
ALLCITIES:=TOUR
forall(i in NodeSet) do
if(i not in ALLCITIES) then
TOUR:={}
first:=i
repeat
TOUR+={first}
first:=NEXTC(first)
until first=i
ALLCITIES+=TOUR
if getsize(TOUR)<size then
SMALLEST:=TOUR
size:=getsize(SMALLEST)
end-if
if size=2 then
break
end-if
end-if
end-do
end-if
!*** Use at least one of the following cuts***
! Add a subtour breaking constraint
sum(i in SMALLEST) fly(i,NEXTC(i)) <= getsize(SMALLEST) - 1
! Optional: Also exclude the inverse subtour
if SMALLEST.size>2 then
sum(i in SMALLEST) fly(NEXTC(i),i) <= getsize(SMALLEST) - 1
end-if
! Optional: Add a stronger subtour elimination cut
sum(i in SMALLEST,j in NodeSet-SMALLEST) fly(i,j) >= 1
! Re-solve the problem
minimize(TotalDist)
break_subtour
end-if
end-procedure
!-----------------------------------------------------------------
! Print the current solution
procedure print_sol
declarations
ALLCITIES: set of integer
end-declarations
writeln("(", gettime-starttime, "s) Total distance: ", getobjval)
ALLCITIES:={}
forall(i in NodeSet) do
if(i not in ALLCITIES) then
write(i)
first:=i
repeat
ALLCITIES+={first}
write(" - ", NEXTC(first))
first:=NEXTC(first)
until first=i
writeln
end-if
end-do
end-procedure
!-----------------------------------------------------------------
! 2-opt algorithm: for all pairs of arcs in a tour, try out the effect
! of exchanging the arcs
procedure two_opt
declarations
NCtemp: array(NodeSet) of integer
end-declarations
if getsize(NodeSet)>3 then
forall(n in NodeSet) do
secondN:=NEXTC(n)
ADist:=DIST(n,NEXTC(n))+DIST(secondN,NEXTC(secondN))
AInd:=secondN
repeat
if DIST(n,NEXTC(n))+DIST(secondN,NEXTC(secondN)) >
DIST(n,secondN)+DIST(NEXTC(n),NEXTC(secondN)) and
ADist > DIST(n,secondN)+DIST(NEXTC(n),NEXTC(secondN)) then
ADist:=DIST(n,secondN)+DIST(NEXTC(n),NEXTC(secondN))
AInd:=secondN
end-if
secondN:=NEXTC(secondN)
until secondN=n
if AInd<>NEXTC(n) then
writeln("2-opt crossover: ", n, "-", NEXTC(n), " and ",
AInd, "-", NEXTC(AInd), " Dist: ", ADist, "<",
DIST(n,NEXTC(n))+DIST(NEXTC(n),NEXTC(NEXTC(n))))
forall(m in NodeSet) NCtemp(m):=-1
i:=NEXTC(n)
nextI:=NEXTC(i)
repeat
NCtemp(nextI):=i
i:=nextI
nextI:=NEXTC(i)
until i=AInd
! Exchange the two arcs
NEXTC(NEXTC(n)):=NEXTC(AInd)
NEXTC(n):=AInd
! Inverse the sense of the tour between NEXTC(n) and AInd
forall(m in NodeSet | NCtemp(m)<>-1) NEXTC(m):=NCtemp(m)
end-if
end-do
end-if
end-procedure
end-model
|
(!*******************************************************
Mosel Example Problems
======================
file tspmain.mos
````````````````
Solving the TSP problem by a series of optimization subproblems.
Master model coordinating all submodels.
We start one Mosel instance per specified node (entries of NODELIST).
We then start K model instances per Mosel instance.
This model calculates random (X,Y) coordinates, organized in
squares containing SN nodes each, arranged in NUMY rows.
The areas per optimization run are calculated upfront.
TSP algorithm:
round 1: solve the TSP associated with each square
round n: for 2 neighbouring areas from the previous run, unfix
the variables closest to the common border and reoptimize
Each submodel instance is sent an optimization problem (set of nodes,
coordinates, and possibly previous results). The results are passed
back via a file (located at the same place as this master model,
no write access to remote instances is required).
Once the result has been displayed, the submodel is restarted
for a new optimization run if any are available.
This master model and also the submodels may run on any platform.
- Submodels running through the Mosel Distributed Framework -
(c) 2011 Fair Isaac Corporation
author: S. Heipcke, Oct. 2011, rev. Aug. 2018
*******************************************************!)
model "TSP (main) no XAD"
uses "mmjobs","mmsystem"
parameters
K = 2 ! Number of submodels per Mosel instance
NUMNODE = 320 ! Instance size
NUM = 0 ! Model ID
CWIDTH =1000 ! Window height (>=NUMPX)
CHEIGHT =1000 ! Window width (>=NUMPX)
SN = 20 ! Size of subproblems (no of nodes)
NUMY = 4 ! Number of rows of squares
XOFF = 100 ! Offset to left canvas border
BOFF = 10 ! Offset around canvas
BIN = "bin:" ! File format (binary)
ZIP = "zlib.deflate:" ! File format (compression)
RMT = "rmt:" ! Whether to use remote files
SUBMODEL = "tspsub" ! Name of the submodel
end-parameters
!!! Configure this list with machines in your local network,
!!! there must be at least 1 entry in this list.
NODELIST:=[""]
declarations
M: integer ! Number of remote Mosel instances
A: range
B: range
INODE: array(B) of string
end-declarations
forall(n in NODELIST, M as counter) INODE(M):=n
A:= 1..M*K
declarations
modPar: array(A) of Model ! Submodels
moselInst: array(B) of Mosel ! Mosel instances
modendct, config, nbrunning: integer ! Counters
ev: Event ! Event message
SQUARE: array(Squares:range,1..4) of integer ! Border coordinates
SOL: array(range) of integer ! Solution values (route)
NODES = 1..NUMNODE ! Set of nodes
NODEX,NODEY: dynamic array(NODES) of integer ! Dynamic to write out index
SNODES: array(NSq:range) of set of integer ! Nodes per subproblem
LEV: array(Squares) of integer ! Iteration round
SBORDER: array(Squares) of integer ! Offset from border
BTYPE: array(Squares) of string ! X:horizontal, Y:vertical neighbors
SQSETS: array(Squares,1..2) of set of integer ! Predecessor subproblems
NumProb: integer ! No. of original subproblems
PREC: array(Squares) of set of integer ! Predecessor subproblems
SQSTATUS: array(Squares) of integer ! Subproblem solve status
ToStart: set of integer ! Data instances to be run
idlemod: list of integer ! List of free submodel instances
starttime: real ! Time measure
end-declarations
starttime:=gettime
!***************** Subroutines ******************
! **** Generate the node coordinates and
! determine node sets for all optimization problems ****
procedure generate_coordinates(seed: integer)
setrandseed(seed)
NumProb:=ceil(NUMNODE/SN)
NumNodes:=floor(NUMNODE/NumProb)
! Positioning of subproblems ("squares")
NUMX:=ceil(NumProb/NUMY)
ct:=0
forall(x in 1..NUMX, y in 1..NUMY, ct as counter) do
SQUARE(ct,1):= round((x-1)*((CWIDTH-2*BOFF)/NUMX)) + BOFF
SQUARE(ct,2):= round(x*((CWIDTH-2*BOFF)/NUMX)) + BOFF
SQUARE(ct,3):= round((y-1)*((CHEIGHT-2*BOFF)/NUMY)) +BOFF
SQUARE(ct,4):= round(y*((CHEIGHT-2*BOFF)/NUMY)) +BOFF
LEV(ct):=1
end-do
! Node coordinates
RestNum:=NUMNODE mod NumProb
ct:=0; lastn:=0
forall(s in 1..NumProb, ct as counter) do
SNODES(s):=lastn+1..lastn+NumNodes+if(ct<=RestNum,1,0)
lastn:=lastn+NumNodes+if(ct<=RestNum,1,0)
forall(n in SNODES(s)) do
NODEX(n) := round(SQUARE(ct,1)+(SQUARE(ct,2)-SQUARE(ct,1))*random)
NODEY(n) := round(SQUARE(ct,3)+(SQUARE(ct,4)-SQUARE(ct,3))*random)
end-do
PREC(s):={}
SQSETS(s,1):={s}
SQSETS(s,2):={}
end-do
finalize(NSq)
! Initialize solution
forall(n in NODES) SOL(n):=n
! Level 2
forall(p in 1..floor(NumProb/2)) do
LEV(NumProb+p):=2
SBORDER(NumProb+p):= SQUARE(p*2-1,4)
BTYPE(NumProb+p):="Y"
SQSETS(NumProb+p,1):={p*2-1}
SQSETS(NumProb+p,2):={p*2}
PREC(NumProb+p):={p*2-1,p*2}
end-do
! Level 3
qs:=Squares.size
forall(p in 1..floor(NumProb/8)) do
LEV(qs+2*p-1):=3
SBORDER(qs+2*p-1):= SQUARE((p-1)*8+1,2)
BTYPE(qs+2*p-1):="X"
SQSETS(qs+2*p-1,1):={(p-1)*8+1,(p-1)*8+2}
SQSETS(qs+2*p-1,2):={(p-1)*8+5,(p-1)*8+6}
PREC(qs+2*p-1):={NumProb+(p-1)*4+1,NumProb+(p-1)*4+3}
LEV(qs+2*p):=3
SBORDER(qs+2*p):= SQUARE((p-1)*8+1,2)
BTYPE(qs+2*p):="X"
SQSETS(qs+2*p,1):={(p-1)*8+3,(p-1)*8+4}
SQSETS(qs+2*p,2):={(p-1)*8+7,p*8}
PREC(qs+(2*p)):={NumProb+(p-1)*4+2,NumProb+p*4}
end-do
! Level 4
ps:=qs
qs:=Squares.size
forall(p in 1..floor(NumProb/8)) do
LEV(qs+p):=4
SBORDER(qs+p):= SQUARE((p-1)*8+2,4)
BTYPE(qs+p):="Y"
SQSETS(qs+p,1):={(p-1)*8+1,(p-1)*8+2,(p-1)*8+5,(p-1)*8+6}
SQSETS(qs+p,2):={(p-1)*8+3,(p-1)*8+4,(p-1)*8+7,p*8}
PREC(qs+p):={ps+2*p-1,ps+2*p}
end-do
! Level 5
l:=3
repeat
l+=1
ps:=qs
qs:=Squares.size
forall(p in 1..floor(NumProb/2^l)) do
LEV(qs+p):=l+1
SBORDER(qs+p):= SQUARE(round((p-1)*2^l+2^(l-1)),2)
BTYPE(qs+p):="X"
SQSETS(qs+p,1):=union(i in round((p-1)*2^l)+1..minlist(NumProb,round((p-1)*2^l+2^(l-1)))) {i}
SQSETS(qs+p,2):=union(i in minlist(NumProb,round((p-1)*2^l+2^(l-1))+1)..minlist(NumProb,round(p*2^l))) {i}
PREC(qs+p):={ps+2*p-1,ps+2*p}
end-do
until ceil(NumProb/2^l)=1
finalize(Squares)
initializations to BIN+ZIP+"tspgendata"
Squares NSq
NODEX
NODEY
SNODES
LEV
SBORDER
BTYPE
SQSETS
PREC
SQUARE
end-initializations
end-procedure
! **** Find the next available area for (re)optimization ****
function get_next_square: integer
returned:=0
forall(s in Squares) do
if SQSTATUS(s)=0 and (and(p in PREC(s)) SQSTATUS(p)>1) then
SQSTATUS(s):=1
returned:=s
break
end-if
end-do
end-function
! **** Retrieve the solution from a run ****
procedure write_box(num: integer, ifall: boolean)
declarations
sol: array(ASet:set of integer) of integer
PB: integer
ALLA: set of integer
end-declarations
initializations from BIN+ZIP+"tspsol"+num
sol as "SOL"
PB
end-initializations
forall(x in ASet) SOL(x):= sol(x)
writeln("Problem ", PB, ":")
! Subtour printout:
ALLA:={}
forall(i in ASet) do
if(i not in ALLA) then
write(i)
first:=i
repeat
ALLA+={first}
write(" - ", SOL(first))
first:=SOL(first)
until first=i
writeln
end-if
end-do
fdelete("tspsol"+num )
SQSTATUS(PB):=2
end-procedure
! **** Solution algorithm ****
procedure solve_and_display
! Start first lot of remote model executions
nbrunning:=0
idlemod:=[]
modendct:=0
forall(n in 1..minlist(M*K,NumProb)) do
nbrunning+=1
SQSTATUS(n):=1
initializations to BIN+ZIP+"tspsolinit"+n
SOL
end-initializations
run(modPar(n), "PB="+n + ",LEVEL=" + LEV(n) + ",NUM=" + n +
",BIN='" + BIN + "',RMT='" + RMT + "',ZIP='" + ZIP + "'" +
",PRINTDETAIL=false")
end-do
! Set of remaining jobs (=data instances)
ToStart:=union(i in minlist(M*K,NumProb)+1..getsize(Squares)) {i}
! Wait for termination and start remaining
while ( modendct<getsize(Squares) ) do
wait
ev:=getnextevent
if getclass(ev)=EVENT_END then
modendct+=1
num:= ev.fromuid
idlemod+=[num]
nbrunning-=1
write_box(num,true)
newnum:=0
ct:=0
while (ToStart<>{} and newnum=0) do
newnum:=get_next_square
if (newnum>0) then
ToStart-={newnum}
num:=getfirst(idlemod)
cuthead(idlemod,1)
nbrunning+=1
initializations to BIN+ZIP+"tspsolinit"+num
SOL
end-initializations
run(modPar(num), "PB="+newnum + ",LEVEL=" + LEV(newnum) +
",NUM=" + num + ",BIN='" + BIN + "',RMT='" + RMT +
"',ZIP='" + ZIP + "'" + ",PRINTDETAIL=false")
else
break
end-if
end-do
end-if
! if nbrunning=0 then break; end-if
end-do
end-procedure
!***************** Main ******************
! Start child nodes
forall(i in B)
if connect(moselInst(i), INODE(i))<>0 then exit(2); end-if
if compile("g",SUBMODEL+".mos")<>0 then exit(3); end-if
forall(j in A) do ! Load models
load(moselInst(j mod M + 1), modPar(j), "rmt:"+SUBMODEL+".bim")
modPar(j).uid:= j
end-do
! Generate data
generate_coordinates(5)
solve_and_display
writeln("Total elapsed time: ", gettime-starttime, "s")
! Cleaning up
forall(i in B) disconnect(moselInst(i))
fdelete(SUBMODEL+".bim")
fdelete("tspgendata")
forall(i in A) fdelete("tspsolinit"+i)
end-model
|
