(!****************************************************** 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. Mar. 2014 *******************************************************!) 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(Sq: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" 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 break end-do Use at least one of the following Add a subtour breaking constraint in <=getsize(SMALLEST) Optional: Also exclude inverse 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