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