(!******************************************************
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. Apr. 2019
*******************************************************!)
model "F-5 Tour planning (nodecb)"
uses "mmxprs", "mmsystem"
parameters
NUM=20
! GLOBAL = false ! Apply global cuts
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 set_optimizer_controls
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 public procedure cb_preintsol(isheur : boolean, cutoff : real)
forward public 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
set_optimizer_controls
! 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
!-----------------------------------------------------------------
! Disable Optimizer features that can reduce the search space. Since
! the Optimizer only sees the initial set of constraints, it could
! make inferences that drop some solutions due to dominance or
! symmetries. These reductions might not be correct when all constraints
! are considered.
!
! So we disable all symmetry or dominance related features.
procedure set_optimizer_controls
declarations
NewPresolveOps : integer
NewMipPresolve : integer
end-declarations
NewPresolveOps := getparam("XPRS_PRESOLVEOPS")
if (bittest(NewPresolveOps, 1) <> 0) then NewPresolveOps -= 1; end-if ! Disable singleton column removal
if (bittest(NewPresolveOps, 8) <> 0) then NewPresolveOps -= 8; end-if ! Disable dual reductions
if (bittest(NewPresolveOps, 64) <> 0) then NewPresolveOps -= 64; end-if ! Disable duplicate column removal
if (bittest(NewPresolveOps, 1024) = 0) then NewPresolveOps += 1024; end-if
! Do not detect and create semi-continuous variables
setparam("XPRS_PRESOLVEOPS", NewPresolveOps)
NewMipPresolve := getparam("XPRS_MIPPRESOLVE")
if (bittest(NewMipPresolve, 16) <> 0) then NewMipPresolve -= 16; end-if ! Disable node-to-node dual reductions
setparam("XPRS_MIPPRESOLVE", NewMipPresolve)
! Disable symmetry detection
setparam("XPRS_SYMMETRY", 0)
end-procedure
!-----------------------------------------------------------------
! 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) ! Mosel5
nexti:=min(i in ToVisit) i
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.
public 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.
public 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
(! if GLOBAL then
storecuts(0,cutid,type,cut,ndx_cuts)
writeln("Stored cuts at node ", node, " -> ", getsize(ndx_cuts))
getcplist(1,1,0, violated_cuts,viol)
writeln("Current violated cuts in the cutpool at node ",
node, " -> ", getsize(violated_cuts),
" (obj. ", getparam("XPRS_LPOBJVAL"), ")")
loadcuts(1,1,violated_cuts)
else !)
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
|
(!******************************************************
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.
(c) 2008 Fair Isaac Corporation
author: S. Heipcke, Dec. 2007, rev. Nov. 2017
*******************************************************!)
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 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
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
|
(!******************************************************
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.
(c) 2008 Fair Isaac Corporation
author: S. Heipcke, Dec. 2007, rev. Nov. 2017
*******************************************************!)
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
end-declarations
! cp_set_solution_callback("print_CP_sol")
setparam("XPRS_VERBOSE", true)
setparam("XPRS_MIPLOG", -25)
setparam("XPRS_GOMCUTS",0)
setparam("XPRS_HEURSEARCHEFFORT", 2) ! Increased local search effort
setparam("XPRS_HEURSEARCHROOTSELECT", 0)
setparam("XPRS_HEURSEARCHTREESELECT", 7)
setparam("XPRS_MAXNODE", 300)
setcallback(XPRS_CB_OPTNODE, "solve_CP")
setcallback(XPRS_CB_INTSOL, "update_CP")
!setparam("XPRS_HEURSTRATEGY", 0)
setparam("XPRS_NODESELECTION", 2) ! Breadth-first search
setparam("XPRS_MIPTHREADS", 1) ! Desactivate 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
if getparam("XPRS_NODES")=1 or getparam("XPRS_NODEDEPTH") MOD 4 = 0 then
writeln("Solving CP")
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)
if cp_find_next_sol then
writeln("(", gettime-starttime, "s) CP solution: ", 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
res:=loadmipsol(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
|
