Initializing help system before first use

MIP start solutions and subtour elimination algorithm for the traveling salesman problem (TSP)
Type: Traveling salesman problem
Rating: 5 (difficult)
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.
File(s): f5tourcbrandom.mos, f5tour3.mos, f5tour4.mos, f5tour4m.mos
Data file(s): gr96.dat, st70.dat, gr120.dat
f5tourcbrandom.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 preintsol callback to add subtour elimination constraints -  
   
   (c) 2015 Fair Isaac Corporation
       author: S. Heipcke, Nov. 2015, rev. Feb. 2024
*******************************************************!)

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

 parameters
  NUM=40
  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)

 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 and create subtour elimination 
! constraints for integer solutions
 setcallback(XPRS_CB_PREINTSOL, ->cb_preintsol)

! Set a callback to display intermediate solutions
 if EXTRA_OUTPUT: setcallback(XPRS_CB_INTSOL, ->print_sol)

! Prevent dominance or symmetry reductions by the Optimizer:
! Since the optimizer does not see the whole model (subtour elimination 
! constraints are only generated on the fly), dual reductions may cut off 
! the optimal solution.
 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. If 'isheur' is true then this is done by calling
! 'rejectintsol'. If 'isheur' is false then the solution came from 
! an integral node. In this case we add at least one violated subtour
! elimination constraint. 
 procedure cb_preintsol(isheur : boolean, cutoff : real)
  declarations
    iCity: integer
    nCitiesVisited,toursize: integer
    SUBTOUR: set of integer
    VISITEDCITIES: set of 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

  ! 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
   ! The tour given by the current solution does not pass through
   ! all the nodes and is thus infeasible.
    if isheur then
    ! Reject infeasible heuristic solution
      if EXTRA_OUTPUT:
        writeln("Rejected heuristic solution (has a ", nCitiesVisited, " city subtour)")
      rejectintsol

    else
    ! Create sub-tour elimination constraints
      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
			
          ! Create the sub-tour elimination constraints as cut(s)
	  toursize:=SUBTOUR.size
          if toursize < NCITIES then
          ! Exclude the subtour (at most toursize-1 of its legs can be used)
            cutid(ncut):= 0
            type(ncut):= CT_LEQ
            cut(ncut):= sum(i in SUBTOUR) fly(i,NEXTC(i)) - (toursize - 1)
            ncut+=1
            if EXTRA_OUTPUT: writeln("Cut ", ncut, ": ", SUBTOUR)
          ! Optional: Also exclude the inverse subtour 
            cutid(ncut):= 0
            type(ncut):= CT_LEQ
            cut(ncut):= sum(i in SUBTOUR) fly(NEXTC(i),i) - (toursize - 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

    end-if
  end-if

 end-procedure

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

! 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
f5tour3.mos
(!******************************************************
   Mosel Example Problems
   ======================

   file f5tour3.mos
   ````````````````
   Tour planning - traveling salesman problem (TSP)
   (See "Applications of optimization with Xpress-MP",
        Section 9.6: Planning a flight tour)

   Combination with CP as start solution does not 
   work with the iterative MIP algorithm in 'f5tour.mos', 
   the MIP model used here therefore contains the complete 
   problem definition.
   
   ALG 1+2: use CP with local-search style enumeration
            to generate a start solution for the MIP search
 
   For small instances (up to approx. 50 cities) pure MIP tends
   to beat the combined approach. Beyond approx. 120 cities
   the problem grows too large for the formulations used here. 

   *** This model cannot be run with a Community Licence 
       for the provided data instances ***

   (c) 2008 Fair Isaac Corporation
       author: S. Heipcke, Dec. 2007, rev. Nov. 2019
*******************************************************!)

model "F-5 Tour planning (MIP + CP)"
 uses "mmxprs", "mmsystem", "kalis"

 parameters
  DATAFILE="Data/gr96.dat"
  ALG=2                ! 0: pure MIP (no start solution)
                       ! 1: CP heuristic for initial solution
                       ! 2: LP as initial solution
 end-parameters 

 forward public procedure print_CP_sol
 forward procedure print_MIP_sol

 declarations
  starttime: real
  CITIES: range                          ! Set of cities
 end-declarations

 starttime:=gettime

 initializations from DATAFILE
  CITIES
 end-initializations

 finalize(CITIES)
 CITIES1:= CITIES-{min(i in CITIES) i}

! **** MIP model ****

 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 
  y: array(CITIES) of mpvar              ! Variables for excluding subtours
 end-declarations

 initializations from DATAFILE
  DIST
 end-initializations

! 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

! Exclude subtours
 forall(i in CITIES, j in CITIES1 | i<>j) y(j) >= y(i) + 1 - NCITIES * (1 - fly(i,j))

! **** CP model ****

 setparam("KALIS_DEFAULT_LB", 0)
 setparam("KALIS_DEFAULT_UB", NCITIES-1)

 declarations
  D: array(CITIES) of integer            ! Auxiliary array
  flyc: array(CITIES) of cpvar           ! Successor of a city 
 end-declarations

 MAXD:= sum(i in CITIES) max(j in CITIES) DIST(i,j)

 setparam("KALIS_DEFAULT_UB", MAXD)

 declarations
  totaldist: cpvar                       ! Total travel distance (objective)
 end-declarations

 MAXD:= max(i,j in CITIES) DIST(i,j)

 setparam("KALIS_DEFAULT_UB", MAXD)

 declarations
  dist: array(CITIES) of cpvar           ! Distance to successor of city i
  fsol: array(CITIES) of integer         ! CP solution values
 end-declarations

! Assign values to 'flyc' variables as to obtain a single cycle,
! 'totaldist' is the distance travelled
 cycle(flyc, totaldist, DIST)

! Element constraints, defining 'dist' variables used in search
! (redundant constraints)
 forall(i in CITIES) do
  forall(j in CITIES) D(j):= DIST(i,j)
  dist(i) = element(D, flyc(i))
 end-do 

! **** Solving: CP start solution ****

 declarations
  CORDER: array(1..NCITIES) of integer   ! Ordered city indices
  dsol: array(CITIES) of real            ! LP solution values
  allsol: array(mpvar) of real           ! MIP start solution
  sdist: integer                         ! CP objective function value
 end-declarations

 cp_set_solution_callback("print_CP_sol")
 setparam("KALIS_MAX_COMPUTATION_TIME", 60)

! Find a first solution
 case ALG of
 1: do    ! Initial solution generated by CP:
   cp_set_branching(assign_var(KALIS_LARGEST_MAX, KALIS_MIN_TO_MAX, dist))
   if not (cp_find_next_sol) then
    writeln("No initial solution found")
    cpsol:=false    
   end-if
 
   forall(i in CITIES) settarget(flyc(i), getsol(flyc(i)))
   forall(i in CITIES) settarget(dist(i), getsol(dist(i)))
   forall(i in CITIES) fsol(i):=flyc(i).sol
   cpsol:=true
   sdist:= getsol(totaldist)
   cp_reset_search
   totaldist <= sdist - maxlist(sdist * 0.001,1)
   cp_set_branching(assign_var(KALIS_LARGEST_MAX, KALIS_NEAREST_VALUE, dist))
   setparam("KALIS_MAX_COMPUTATION_TIME", 10)
  end-do

 2: do    ! Using LP results as target values for CP:
   minimize(XPRS_LIN,TotalDist)       ! Solve the LP-relaxation
   forall(i in CITIES) do
    minj:=i; minv:=0.0
    forall(j in CITIES | i<>j) 
     if getsol(fly(i,j))>minv then
      minj:=j; minv:=getsol(fly(i,j))
     end-if
    dsol(i):= minv
    settarget(flyc(i), minj)
    settarget(dist(i), 
              round(sum(j in CITIES | i<>j) DIST(i,j)*getsol(fly(i,j))) )
   end-do
   qsort(SYS_DOWN, dsol, CORDER)

   forall(i in CITIES)
    Strategy(i):= assign_var(KALIS_LARGEST_MAX, KALIS_NEAREST_VALUE, 
                             {dist(CORDER(i+1))})
   cp_set_branching(Strategy)
  end-do
 end-case

 if ALG<>0 then
! Solve the CP problem: minimize total distance using nearest-neighbor search
! on distance variables
  cpstart:=gettime
  while (gettime-cpstart<5 and cp_find_next_sol) do
   cpsol:=true
   sdist:= getsol(totaldist)
   forall(i in CITIES) settarget(dist(i),getsol(dist(i)))
   forall(i in CITIES) fsol(i):=flyc(i).sol
   cp_reset_search
   totaldist <= sdist - maxlist(sdist * 0.001,1)
   setparam("KALIS_MAX_COMPUTATION_TIME", 10)
  end-do

  if cpsol then                      
   loadprob(TotalDist)             ! Load best CP solution as MIP start sol.
   forall(i in CITIES) allsol(fly(i,fsol(i))):=1
   res:=loadmipsol(allsol)
  else
   writeln("(", gettime-starttime, "s) No CP start solution found.") 
  end-if
 end-if

! **** Solving: MIP ****

 setparam("XPRS_VERBOSE", true)
 setparam("XPRS_MIPLOG", -500)

 setparam("XPRS_GOMCUTS",0)    
 setparam("XPRS_HEURSEARCHEFFORT", 2)      ! Increased local search effort
 setparam("XPRS_HEURSEARCHROOTSELECT", 7)
 setparam("XPRS_HEURSEARCHTREESELECT", 3)
 if ALG<>0 then
  setparam("XPRS_MAXNODE", 1000)
 else
  setparam("XPRS_MAXNODE", 1500)
 end-if

 minimize(TotalDist)                ! Solve the MIP problem
 
 if getparam("XPRS_MIPSOLS") > 0 then
  print_MIP_sol
 else
  writeln("(", gettime-starttime, "s) No MIP solutions.")   
 end-if
 
!-----------------------------------------------------------------
! Print the current CP solution
 public procedure print_CP_sol
  declarations
   ALLCITIES: set of integer
  end-declarations
   
  writeln("(", gettime-starttime, "s) Total distance: ", getsol(totaldist))
  ALLCITIES:={}
  forall(i in CITIES) do
   if(i not in ALLCITIES) then
    write(i)
    first:=i
    repeat
     ALLCITIES+={first}
     write(" - ", flyc(first).sol)
     first:= flyc(first).sol
    until first=i
    writeln 
   end-if
  end-do        
 end-procedure

! Print the current MIP solution
 procedure print_MIP_sol
  declarations
   ALLCITIES: set of integer
  end-declarations
   
  writeln("(", gettime-starttime, "s) Total distance: ", getobjval)
  forall(i in CITIES) 
   forall(j in CITIES)
    if(getsol(fly(i,j))>0.99) then 
     NEXTC(i):=j
     break
    end-if

  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

end-model
f5tour4.mos
(!******************************************************
   Mosel Example Problems
   ======================

   file f5tour4.mos
   ````````````````
   Tour planning - traveling salesman problem (TSP)
   (See "Applications of optimization with Xpress-MP",
        Section 9.6: Planning a flight tour)

   Generate and load heuristic CP solutions during MIP search
   using the 'optimized node' and 'integer solution' callbacks.
   
   For small instances (up to approx. 50 cities) pure MIP tends
   to beat the combined approach. Beyond approx. 120 cities
   the problem grows too large for the formulations used here. 

   *** This model cannot be run with a Community Licence 
       for the provided data instances ***

   (c) 2008 Fair Isaac Corporation
       author: S. Heipcke, Dec. 2007, rev. Mar. 2022
*******************************************************!)

model "F-5 Tour planning (MIP + CP)"
 uses "mmxprs", "mmsystem", "kalis"

 parameters
  DATAFILE="Data/gr96.dat"
 end-parameters 

 forward procedure update_CP
 forward function solve_CP: boolean
 forward procedure print_CP_sol
 forward procedure print_MIP_sol

 declarations
  starttime: real
  CITIES: range                          ! Set of cities
 end-declarations

 starttime:=gettime

 initializations from DATAFILE
  CITIES
 end-initializations

 finalize(CITIES)
 CITIES1:= CITIES-{min(i in CITIES) i}

! **** MIP model ****

 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 
  y: array(CITIES) of mpvar              ! Variables for excluding subtours
 end-declarations

 initializations from DATAFILE
  DIST
 end-initializations

! 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

! Exclude subtours
 forall(i in CITIES, j in CITIES1 | i<>j) y(j) >= y(i) + 1 - NCITIES * (1 - fly(i,j))

! **** CP model ****

 setparam("KALIS_DEFAULT_LB", 0)
 setparam("KALIS_DEFAULT_UB", NCITIES-1)

 declarations
  D: array(CITIES) of integer            ! Auxiliary array
  flyc: array(CITIES) of cpvar           ! Successor of a city 
 end-declarations

 MAXD:= sum(i in CITIES) max(j in CITIES) DIST(i,j)

 setparam("KALIS_DEFAULT_UB", MAXD)

 declarations
  totaldist: cpvar                       ! Total travel distance (objective)
 end-declarations

 MAXD:= max(i,j in CITIES) DIST(i,j)

 setparam("KALIS_DEFAULT_UB", MAXD)

 declarations
  dist: array(CITIES) of cpvar           ! Distance to successor of city i
  fsol: array(CITIES) of integer         ! CP solution values
 end-declarations

! Assign values to 'flyc' variables as to obtain a single cycle,
! 'totaldist' is the distance travelled
 cycle(flyc, totaldist, DIST)

! Element constraints, defining 'dist' variables used in search
! (redundant constraints)
 forall(i in CITIES) do
  forall(j in CITIES) D(j):= DIST(i,j)
  dist(i) = element(D, flyc(i))
 end-do 

! **** Solving ****

 declarations
  CORDER: array(1..NCITIES) of integer   ! Ordered city indices
  dsol: array(CITIES) of real            ! LP solution values
  allsol: array(mpvar) of real           ! MIP start solution
  sdist: integer                         ! CP objective function value
  solct: integer                         ! CP solution counter
 end-declarations

! cp_set_solution_callback("print_CP_sol")

 setparam("XPRS_VERBOSE", true)
 setparam("XPRS_MIPLOG", -25)

 setparam("XPRS_GOMCUTS",0)    
 setparam("XPRS_MAXNODE", 300)

 setcallback(XPRS_CB_OPTNODE, ->solve_CP)
 setcallback(XPRS_CB_INTSOL, ->update_CP)
 
 setparam("XPRS_NODESELECTION", 2)       ! Breadth-first search
 setparam("XPRS_MIPTHREADS", 1)          ! Deactivate parallel MIP
 
 minimize(TotalDist)                     ! Solve the MIP problem
 
 if getparam("XPRS_MIPSOLS") > 0 then
  print_MIP_sol
 else
  writeln("(", gettime-starttime, "s) No MIP solutions.")   
 end-if

!-----------------------------------------------------------------
! Update CP with solutions found by MIP
 procedure update_CP
  totaldist <= getparam("XPRS_LPOBJVAL") - 1
 end-procedure 

! Generate and load CP solutions into MIP
 function solve_CP:boolean
  nd:=getparam("XPRS_NODES")
  if nd=1 or getparam("XPRS_NODEDEPTH") MOD 4 = 0 then
   writeln("Solving CP (node ", nd,")")
   forall(i in CITIES) do
    minj:=i; minv:=0.0
    forall(j in CITIES | i<>j) 
     if getsol(fly(i,j))>minv then
      minj:=j; minv:=getsol(fly(i,j))
     end-if
    dsol(i):= minv
    settarget(flyc(i), minj)
    settarget(dist(i), 
              round(sum(j in CITIES | i<>j) DIST(i,j)*getsol(fly(i,j))) )
   end-do
   qsort(SYS_DOWN, dsol, CORDER)

   forall(i in CITIES)
    Strategy(i):= assign_var(KALIS_LARGEST_MAX, KALIS_NEAREST_VALUE, 
                             {dist(CORDER(i+1))})
   cp_set_branching(Strategy)
   setparam("KALIS_MAX_COMPUTATION_TIME", ceil(NCITIES/10))
   setparam("KALIS_MAX_SOLUTIONS", 1)
!   setparam("KALIS_VERBOSE_LEVEL", 2)    

   if cp_find_next_sol then 
    solct+=1
    writeln("(", gettime-starttime, "s) CP solution ", solct,": ", getsol(totaldist))
    sdist:= getsol(totaldist)
    forall(i in CITIES) fsol(i):=flyc(i).sol
    forall(i in CITIES) allsol(fly(i,fsol(i))):=1

    cp_reset_search

    addmipsol("CP"+solct, allsol)

    forall(i in CITIES) allsol(fly(i,fsol(i))):=0
    totaldist <= sdist - maxlist(sdist * 0.0001,1)
   end-if
  end-if
 end-function
 
!-----------------------------------------------------------------
! Print the current CP solution
 procedure print_CP_sol
  declarations
   ALLCITIES: set of integer
  end-declarations
   
  writeln("(", gettime-starttime, "s) Total distance: ", getsol(totaldist))
  ALLCITIES:={}
  forall(i in CITIES) do
   if(i not in ALLCITIES) then
    write(i)
    first:=i
    repeat
     ALLCITIES+={first}
     write(" - ", flyc(first).sol)
     first:= flyc(first).sol
    until first=i
    writeln 
   end-if
  end-do        
 end-procedure

! Print the current MIP solution
 procedure print_MIP_sol
  declarations
   ALLCITIES: set of integer
  end-declarations
   
  writeln("(", gettime-starttime, "s) Total distance: ", getobjval)
  forall(i in CITIES) 
   forall(j in CITIES)
    if(getsol(fly(i,j))>0.99) then 
     NEXTC(i):=j
     break
    end-if

  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

end-model
f5tour4m.mos
(!******************************************************
   Mosel Example Problems
   ======================

   file f5tour4m.mos
   `````````````````
   Tour planning - traveling salesman problem (TSP)
   (See "Applications of optimization with Xpress-MP",
        Section 9.6: Planning a flight tour)

   Generate and load heuristic CP solutions during MIP search
   using the 'optimized node' and 'integer solution' callbacks.
   - Passing callback functions by name -
   
   For small instances (up to approx. 50 cities) pure MIP tends
   to beat the combined approach. Beyond approx. 120 cities
   the problem grows too large for the formulations used here. 

   *** This model cannot be run with a Community Licence 
       for the provided data instances ***

   (c) 2008 Fair Isaac Corporation
       author: S. Heipcke, Dec. 2007, rev. Jun. 2021
*******************************************************!)

model "F-5 Tour planning (MIP + CP)"
 uses "mmxprs", "mmsystem", "kalis"

 parameters
  DATAFILE="Data/gr96.dat"
 end-parameters 

 forward public procedure update_CP
 forward public function solve_CP: boolean
 forward procedure print_CP_sol
 forward procedure print_MIP_sol

 declarations
  starttime: real
  CITIES: range                          ! Set of cities
 end-declarations

 starttime:=gettime

 initializations from DATAFILE
  CITIES
 end-initializations

 finalize(CITIES)
 CITIES1:= CITIES-{min(i in CITIES) i}

! **** MIP model ****

 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 
  y: array(CITIES) of mpvar              ! Variables for excluding subtours
 end-declarations

 initializations from DATAFILE
  DIST
 end-initializations

! 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

! Exclude subtours
 forall(i in CITIES, j in CITIES1 | i<>j) y(j) >= y(i) + 1 - NCITIES * (1 - fly(i,j))

! **** CP model ****

 setparam("KALIS_DEFAULT_LB", 0)
 setparam("KALIS_DEFAULT_UB", NCITIES-1)

 declarations
  D: array(CITIES) of integer            ! Auxiliary array
  flyc: array(CITIES) of cpvar           ! Successor of a city 
 end-declarations

 MAXD:= sum(i in CITIES) max(j in CITIES) DIST(i,j)

 setparam("KALIS_DEFAULT_UB", MAXD)

 declarations
  totaldist: cpvar                       ! Total travel distance (objective)
 end-declarations

 MAXD:= max(i,j in CITIES) DIST(i,j)

 setparam("KALIS_DEFAULT_UB", MAXD)

 declarations
  dist: array(CITIES) of cpvar           ! Distance to successor of city i
  fsol: array(CITIES) of integer         ! CP solution values
 end-declarations

! Assign values to 'flyc' variables as to obtain a single cycle,
! 'totaldist' is the distance travelled
 cycle(flyc, totaldist, DIST)

! Element constraints, defining 'dist' variables used in search
! (redundant constraints)
 forall(i in CITIES) do
  forall(j in CITIES) D(j):= DIST(i,j)
  dist(i) = element(D, flyc(i))
 end-do 

! **** Solving ****

 declarations
  CORDER: array(1..NCITIES) of integer   ! Ordered city indices
  dsol: array(CITIES) of real            ! LP solution values
  allsol: array(mpvar) of real           ! MIP start solution
  sdist: integer                         ! CP objective function value
  solct: integer                         ! CP solution counter
 end-declarations

! cp_set_solution_callback("print_CP_sol")

 setparam("XPRS_VERBOSE", true)
 setparam("XPRS_MIPLOG", -25)

 setparam("XPRS_GOMCUTS",0)    
 setparam("XPRS_MAXNODE", 300)

 setcallback(XPRS_CB_OPTNODE, "solve_CP")
 setcallback(XPRS_CB_INTSOL, "update_CP")
 
 setparam("XPRS_NODESELECTION", 2)       ! Breadth-first search
 setparam("XPRS_MIPTHREADS", 1)          ! Deactivate parallel MIP
 
 minimize(TotalDist)                     ! Solve the MIP problem
 
 if getparam("XPRS_MIPSOLS") > 0 then
  print_MIP_sol
 else
  writeln("(", gettime-starttime, "s) No MIP solutions.")   
 end-if

!-----------------------------------------------------------------
! Update CP with solutions found by MIP
 public procedure update_CP
  totaldist <= getparam("XPRS_LPOBJVAL") - 1
 end-procedure 

! Generate and load CP solutions into MIP
 public function solve_CP:boolean
  nd:=getparam("XPRS_NODES")
  if nd=1 or getparam("XPRS_NODEDEPTH") MOD 4 = 0 then
   writeln("Solving CP (node ", nd,")")
   forall(i in CITIES) do
    minj:=i; minv:=0.0
    forall(j in CITIES | i<>j) 
     if getsol(fly(i,j))>minv then
      minj:=j; minv:=getsol(fly(i,j))
     end-if
    dsol(i):= minv
    settarget(flyc(i), minj)
    settarget(dist(i), 
              round(sum(j in CITIES | i<>j) DIST(i,j)*getsol(fly(i,j))) )
   end-do
   qsort(SYS_DOWN, dsol, CORDER)

   forall(i in CITIES)
    Strategy(i):= assign_var(KALIS_LARGEST_MAX, KALIS_NEAREST_VALUE, 
                             {dist(CORDER(i+1))})
   cp_set_branching(Strategy)
   setparam("KALIS_MAX_COMPUTATION_TIME", ceil(NCITIES/10))
   setparam("KALIS_MAX_SOLUTIONS", 1)
!   setparam("KALIS_VERBOSE_LEVEL", 2)    

   if cp_find_next_sol then 
    solct+=1
    writeln("(", gettime-starttime, "s) CP solution ", solct,": ", getsol(totaldist))
    sdist:= getsol(totaldist)
    forall(i in CITIES) fsol(i):=flyc(i).sol
    forall(i in CITIES) allsol(fly(i,fsol(i))):=1

    cp_reset_search

    addmipsol("CP"+solct, allsol)

    forall(i in CITIES) allsol(fly(i,fsol(i))):=0
    totaldist <= sdist - maxlist(sdist * 0.0001,1)
   end-if
  end-if
 end-function
 
!-----------------------------------------------------------------
! Print the current CP solution
 procedure print_CP_sol
  declarations
   ALLCITIES: set of integer
  end-declarations
   
  writeln("(", gettime-starttime, "s) Total distance: ", getsol(totaldist))
  ALLCITIES:={}
  forall(i in CITIES) do
   if(i not in ALLCITIES) then
    write(i)
    first:=i
    repeat
     ALLCITIES+={first}
     write(" - ", flyc(first).sol)
     first:= flyc(first).sol
    until first=i
    writeln 
   end-if
  end-do        
 end-procedure

! Print the current MIP solution
 procedure print_MIP_sol
  declarations
   ALLCITIES: set of integer
  end-declarations
   
  writeln("(", gettime-starttime, "s) Total distance: ", getobjval)
  forall(i in CITIES) 
   forall(j in CITIES)
    if(getsol(fly(i,j))>0.99) then 
     NEXTC(i):=j
     break
    end-if

  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

end-model
Parent Topic Mosel language: Examples from the whitepaper 'Hybrid MIP/CP solving with Xpress Optimizer and Xpress Kalis'

© 2001-2025 Fair Isaac Corporation. All rights reserved. This documentation is the property of Fair Isaac Corporation (“FICO”). Receipt or possession of this documentation does not convey rights to disclose, reproduce, make derivative works, use, or allow others to use it except solely for internal evaluation purposes to determine whether to purchase a license to the software described in this documentation, or as otherwise set forth in a written software license agreement between you and FICO (or a FICO affiliate). Use of this documentation and the software described in it must conform strictly to the foregoing permitted uses, and no other use is permitted.