FICO
FICO Xpress Optimization Examples Repository
FICO Optimization Community FICO Xpress Optimization Home
Back to examples browserPrevious exampleNext example

MIP start solutions and subtour elimination algorithm for the traveling salesman problem (TSP)

Description
This example shows how to construct and load solutions for the MIP branch-and-bound search. The model version f5touroptcbrandom.mos shows how to add subtour elimination constraints via solver callbacks during the MIP Branch-and-bound search.
  • Model f5touroptcbrandom.mos: several heuristic start solutions are loaded into a MIP model for solving symmetric TSP via subtour elimination constraints that are added during the MIP Branch-and-bound search using the cut management functionality of Xpress Optimizer in the OPTNODE callback. The initial MIP problem statement is incomplete, all symmetry or dominance related features therefore need to be disabled and the PREINTSOL callback is used to reject any heuristic solutions that contain subtours.
  • Model f5tour3.mos: a CP model generates a start solution that is loaded into the Optimizer before the MIP Branch-and-bound search. With the model parameter ALG set to 2 the CP search uses the (rounded) solution values of the LP relaxation as initial target values for its search.
  • Model f5tour4.mos: a CP model is run at the nodes of the Branch-and-bound tree using the current LP relaxation solution as input. If a solution is found, it is loaded into the Optimizer for exploitation by the MIP heuristics.

f5tour_h.zip[download all files]

Source Files
By clicking on a file name, a preview is opened at the bottom of this page.
f5touroptcbrandom.mos[download]
f5tour3.mos[download]
f5tour4.mos[download]
f5tour4m.mos[download]

Data Files
gr96.dat[download]
st70.dat[download]
gr120.dat[download]




f5touroptcbrandom.mos

(!******************************************************
   Mosel Example Problems
   ======================

   file f5touroptcbrandom.mos
   ``````````````````````````
   Problem 9.6: Planning a flight tour
   (using randomly generated data with integer indices)
   
   - Creatig and loading heuristic start solutions -
   - Using optnode callback to add subtour elimination constraints -  
   
   (c) 2015 Fair Isaac Corporation
       author: S. Heipcke, Nov. 2015, rev. Jan. 2023
*******************************************************!)

model "F-5 Tour planning (nodecb)"
 uses "mmxprs", "mmsystem"

 parameters
  NUM=20
  TOL=10E-10
  EXTRA_OUTPUT = false		! Whether to print information about solutions checked and cuts created
 end-parameters 

 forward procedure print_sol
 forward procedure check_sol
 forward procedure two_opt

 forward procedure create_initial_tour(Tour : array(set of mpvar) of real)
 forward procedure create_initial_tour2(Tour : array(set of mpvar) of real)
 forward procedure cb_preintsol(isheur : boolean, cutoff : real)
 forward function cb_optnode : boolean


 declarations
  starttime: real
  CITIES = 1..NUM                  ! Set of cities
  X,Y: array(CITIES) of real
 end-declarations

 starttime:=gettime
 
 setrandseed(3)
 forall(i in CITIES) X(i):=100*random
 forall(i in CITIES) Y(i):=100*random

 declarations
  NCITIES=getsize(CITIES)
  
  DIST: array(CITIES,CITIES) of integer  ! Distance between cities
  NEXTC: array(CITIES) of integer         ! Next city after i in the solution

  fly: array(CITIES,CITIES) of mpvar     ! 1 if flight from i to j

  InitTour, InitTour2 : array(set of mpvar) of real
  MipSolutions : integer
 end-declarations
 
 forall(i,j in CITIES | i<>j) 
  DIST(i,j):=round(sqrt((X(i)-X(j))^2 + (Y(i)-Y(j))^2))

! Objective: total distance
 TotalDist:= sum(i,j in CITIES | i<>j) DIST(i,j)*fly(i,j)

! Visit every city once
 forall(i in CITIES) sum(j in CITIES | i<>j) fly(i,j) = 1
 forall(j in CITIES) sum(i in CITIES | i<>j) fly(i,j) = 1

 forall(i,j in CITIES | i<>j) fly(i,j) is_binary


! Set a callback to check heuristic solutions
 setcallback(XPRS_CB_PREINTSOL, ->cb_preintsol)

! Set a callback to create subtour elimination constraints for integer solutions
 setcallback(XPRS_CB_OPTNODE, ->cb_optnode)

! Prevent dominance or symmetry reductions by the Optimizer
 setparam("XPRS_MIPDUALREDUCTIONS", 0)

! Set some additional tuning/output controls for the solve.
 setparam("XPRS_VERBOSE", true)        ! Enable Optimizer output

! Create an initial tour as a starting solution
 create_initial_tour(InitTour)
 create_initial_tour2(InitTour2)
 
! Start the solve of the problem and stop after the initial LP solution...
 minimize(XPRS_LPSTOP, TotalDist)
! ... load our initial tours ...
 addmipsol("InitTour", InitTour)
 addmipsol("InitTour2", InitTour2)
! ... and continue.
 minimize(XPRS_CONT, TotalDist)

 MipSolutions := getparam("XPRS_MIPSOLS")
 if (MipSolutions > 0) then
  print_sol
 end-if


!-----------------------------------------------------------------
! Create an initial tour as a starting solution for the MIP search.
! We simply visit every city in sequence.
 procedure create_initial_tour(Tour : array(set of mpvar) of real)	
  forall(i,j in CITIES | i<>j) Tour(fly(i,j)) := 0
  forall(i in CITIES | i>1) Tour(fly(i-1,i)) := 1
  Tour(fly(NCITIES,1)) := 1
 end-procedure

! This start solution first uses a greedy algorithm (nearest neighbour),
! followed by a 2-opt step to find possible local improvements
 procedure create_initial_tour2(Tour : array(set of mpvar) of real)
  declarations
   ToVisit,Visited:set of integer
   first,lasti,nexti: integer
  end-declarations
  forall(i,j in CITIES | i<>j) Tour(fly(i,j)) := 0
  ToVisit:=CITIES
  first:=min(i in CITIES) i
  ToVisit-={first}; Visited+={first}
  lasti:=first
  while (ToVisit<>{}) do
   nexti:=getelt(ToVisit)
   forall(i in ToVisit)
    if DIST(lasti,i)<DIST(lasti,nexti) then
     nexti:=i
    end-if
   NEXTC(lasti):=nexti
   ToVisit-={nexti}; Visited+={nexti}
   lasti:=nexti
  end-do
  NEXTC(lasti):=first

!  writeln("(", gettime-starttime, "s) Total distance 2a: ", sum(i in 1..NCITIES) DIST(i, NEXTC(i)))
  print_sol
  two_opt
  print_sol
!  writeln("(", gettime-starttime, "s) Total distance 2b: ", sum(i in 1..NCITIES) DIST(i, NEXTC(i)))
  forall(i in CITIES) Tour(fly(i,NEXTC(i))) := 1
 end-procedure

!-----------------------------------------------------------------
! Callback function triggered when the Optimizer has a new integer
! solution
!
! We use this function to reject any heuristic solutions that
! contain subtours.
 procedure cb_preintsol(isheur : boolean, cutoff : real)
  declarations
   iCity: integer
   nCitiesVisited: integer
  end-declarations

  ! Get the tour(s)
  forall(i in CITIES) 
   forall(j in CITIES)
    if(getsol(fly(i,j))>=1-TOL) then 
     NEXTC(i):=j
     break
    end-if

  ! Verify that we have a complete tour
  nCitiesVisited := 1;
  iCity := 1
  while (NEXTC(iCity) <> 1) do
   iCity := NEXTC(iCity)
   nCitiesVisited += 1
  end-do
	
  if (nCitiesVisited < NCITIES) then
   if isheur then
    if EXTRA_OUTPUT then
     writeln("Rejected heuristic solution (has a ", nCitiesVisited, " city subtour)")
    end-if
    rejectintsol
   else
    writeln("INVALID INTEGER SOLUTION!!!")
   end-if
  end-if

 end-procedure

!-----------------------------------------------------------------
! Callback function triggered when a node relaxation has been
! solved to optimality and the Optimizer is about to branch or
! declare the node integer feasible.
!
! If the node is about to be declared integer feasible, we need
! to check for subtours and add at least one violated subtour
! elimination constraint.
 function cb_optnode: boolean
  declarations
   SUBTOUR: set of integer
   VISITEDCITIES: set of integer		
   iCity: integer		
   MipInfeas: integer
   ncut:integer
   CutRange: range                        ! Counter for cuts
   cut: array(CutRange) of linctr         ! Cuts
   cutid: array(CutRange) of integer      ! Cut type identification
   type: array(CutRange) of integer       ! Cut constraint type
  end-declarations

  returned:=false                     ! Call this function once per node

  MipInfeas := getparam("XPRS_MIPINFEAS");
  if (MipInfeas = 0) then
  ! There are no fractionals left in the current node solution.
  ! Create a sub-tour elimination constraints if we do not
  ! have a complete tour.
   ncut:= 0

   forall(i in CITIES) 
    forall(j in CITIES)
     if(getsol(fly(i,j))>=1-TOL) then 
      NEXTC(i):=j
      break
     end-if

   VISITEDCITIES := {}
   forall (iFirstCity in CITIES)
    if not iFirstCity in VISITEDCITIES then	
     SUBTOUR := {}
     iCity := iFirstCity
     while (iCity not in SUBTOUR) do
      SUBTOUR += {iCity}
      iCity := NEXTC(iCity)
     end-do
			
     if (getsize(SUBTOUR) < NCITIES) then
     ! Create the sub-tour elimination cut(s)
      cutid(ncut):= 0
      type(ncut):= CT_LEQ
      cut(ncut):= sum(i in SUBTOUR) fly(i,NEXTC(i)) - (getsize(SUBTOUR) - 1)
      ncut+=1
      if EXTRA_OUTPUT then writeln("Cut ", ncut, ": ", SUBTOUR); end-if
     ! Optional: Also exclude the inverse subtour
      cutid(ncut):= 0
      type(ncut):= CT_LEQ
      cut(ncut):= sum(i in SUBTOUR) fly(NEXTC(i),i) - (getsize(SUBTOUR) - 1)
      ncut+=1!)
     ! Optional: Add a stronger subtour elimination cut
      cutid(ncut):= 0
      type(ncut):= CT_GEQ
      cut(ncut) := sum(i in SUBTOUR, j in (CITIES)-SUBTOUR) fly(i,j) - 1
      ncut+=1
     end-if	
     VISITEDCITIES += SUBTOUR
    end-if

   ! Add cuts to the problem
   if(ncut>0) then    
      addcuts(cutid, type, cut);
      num_cuts_added+=ncut
      if(getparam("XPRS_VERBOSE")=true and EXTRA_OUTPUT) then
        node:=getparam("XPRS_NODES")
        depth:=getparam("XPRS_NODEDEPTH")
        writeln("Cuts added : ", ncut, " (depth ", depth, ", node ", 
                node, ", obj. ", getparam("XPRS_LPOBJVAL"), ") count=",
                num_cuts_added)
      end-if   
    end-if

!   returned:=true                   ! Repeat until no new cuts generated
 
  end-if      ! MIPINFEAS

 end-function

!-----------------------------------------------------------------

! Print the current solution
 procedure print_sol 
  declarations
   ALLCITIES: set of integer
   CalcDistance: real
  end-declarations

  ! Get the final tour
   forall(i in CITIES) 
    forall(j in CITIES)
     if(getsol(fly(i,j))>=1-TOL) then 
      NEXTC(i):=j
      break
     end-if

  ! Recalculate the length of the tour
  CalcDistance := sum(i in 1..NCITIES) DIST(i, NEXTC(i))

  writeln("(", gettime-starttime, "s) Total distance: ", CalcDistance)
  ALLCITIES:={}
  forall(i in CITIES) 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

!-----------------------------------------------------------------

! Check whether the current solution is a single tour and if yes, print it
 procedure check_sol
  declarations
   ALLCITIES: set of integer
   ntour: integer
   TOURLIST: array(TR:range) of list of integer
   TOL=10E-5
  end-declarations
  
  ! Get the successor solution value of every city
  forall(i in CITIES) 
   forall(j in CITIES | getsol(fly(i,j))>=1-TOL) do 
    NEXTC(i):=j
    break
   end-do

  ! Calculate tour(s)
  ALLCITIES:={}
  forall(i in CITIES) do
   if(i not in ALLCITIES) then
    ntour+=1
    TOURLIST(ntour)+=[i]
    first:=i
    repeat
     ALLCITIES+={first}
     TOURLIST(ntour)+=[NEXTC(first)]
     first:=NEXTC(first)
    until first=i
   end-if
  end-do

  ! Print tour output
  if TR.size>1 then
   writeln("(", gettime-starttime, "s) ", TR.size, " subtours: ", getobjval)
  else 
   writeln("(", gettime-starttime, "s) Total distance: ", getobjval)
   forall(t in TR) do
     write(gethead(TOURLIST(t),1))
     cuthead(TOURLIST(t),1)
     forall(i in TOURLIST(t)) write(" - ", i)
     writeln      
   end-do
  end-if
 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(CITIES) of integer
     secondN,AInd: integer
     ADist: real
   end-declarations

   if getsize(CITIES)>3 then
     forall(n in CITIES) 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 CITIES) 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 CITIES | NCtemp(m)<>-1) NEXTC(m):=NCtemp(m)
       end-if
         
     end-do
   end-if
 
 end-procedure

end-model
Back to examples browserPrevious exampleNext example
© Copyright 2023 Fair Isaac Corporation.