Solving the TSP problem by a series of optimization subproblems

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

   *** Not intended to be run standalone - run from tspmain.mos ***
   
   (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