(!****************************************************** Mosel Example Problems ====================== file tsp.mos ```````````` TYPE: Symmetric traveling salesman problem DIFFICULTY: 5 FEATURES: MIP problem, loop over problem solving, TSP subtour elimination algorithm; procedure for generating additional constraints, recursive subroutine calls, working with sets, `forall-do', `repeat-until', `getsize', `not', graphical representation of solutions DESCRIPTION: A flight tour starts from airport 1, visits all the other airports and then comes back to the starting point. The distances between the airports are symmetric. In which order should the airports be visited to minimize the total distance covered? FURTHER INFO: `Applications of optimization with Xpress-MP', Section 11,5 `Planning a flight tour', Section 7,5 `Paint production' (c) 2008 Fair Isaac Corporation author: S. Heipcke, 2002, rev. Sep. 2017 *******************************************************!) model "Tour planning" uses "mmxprs", "mmsvg" forward procedure break_subtour forward procedure print_sol forward procedure draw_sol declarations NCITIES = 7 CITIES = 1..NCITIES ! Cities DIST: array(CITIES,CITIES) of integer ! Distance between cities NEXTC: array(CITIES) of integer ! Next city after i in the solution X,Y: array(CITIES) of integer ! x-y-coordinates for sol. drawing fly: array(CITIES,CITIES) of mpvar ! 1 if flight from i to j end-declarations initializations from 'tsp.dat' DIST [X,Y] as 'POS' end-initializations forall(i,j in CITIES | i<j) DIST(j,i):=DIST(i,j) Objective: total distance TotalDist:=sum(i,j in CITIES>j) DIST(i,j)*fly(i,j) ! Visit every city once forall(i in CITIES) OneVisitI(i):= sum(j in CITIES | i<>j) fly(i,j) = 1 forall(j in CITIES) OneVisitJ(j):= sum(i in CITIES | i<>j) fly(i,j) = 1 forall(i,j in CITIES | i<>j) fly(i,j) is_binary ! Solve the problem minimize(TotalDist) ! Eliminate subtours break_subtour svgwaitclose("Close browser window to terminate model execution.", 1) !----------------------------------------------------------------- procedure break_subtour declarations TOUR,SMALLEST,ALLCITIES: set of integer end-declarations forall(i in CITIES) NEXTC(i):= integer(round(getsol(sum(j in CITIES) j*fly(i,j) ))) ! Print the current solution print_sol draw_sol ! Get (sub)tour containing city 1 TOUR:={} first:=1 repeat TOUR+={first} first:=NEXTC(first) until first=1 size:=getsize(TOUR) ! Find smallest subtour if size < NCITIES then SMALLEST:=TOUR if size>2 then ALLCITIES:=TOUR forall(i in CITIES) 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 break end-do Add a subtour breaking constraint in <=getsize(SMALLEST) also use these additional constraints Also exclude the inverse SMALLEST.size>2 then sum(i in SMALLEST) fly(NEXTC(i),i) <= getsize(SMALLEST) - 1 end-if ! A stronger subtour elimination constraint sum(i in SMALLEST,j in CITIES-SMALLEST) fly(i,j) >= 1 !) ! Re-solve the problem minimize(TotalDist) ! Closing the graphical display will interrupt the algorithm if not svgclosing then break_subtour; end-if end-if end-procedure !----------------------------------------------------------------- ! Print the current solution procedure print_sol declarations ALLCITIES: set of integer end-declarations writeln("Total distance: ", getobjval) 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 !----------------------------------------------------------------- ! Draw the current solution procedure draw_sol declarations ALLCITIES: set of integer end-declarations svgerase minx:=min(i in CITIES) X(i); miny:= min(i in CITIES) Y(i) ! svgsetgraphviewbox(minx-10, miny-10, ! max(i in CITIES) X(i)+20-minx, max(i in CITIES) Y(i)+20-miny) ! Draw the tour(s) ALLCITIES:={}; ct:=0 forall(i in CITIES) do if(i not in ALLCITIES) then svgaddgroup("T"+ct, "Tour "+(ct+1)) first:=i repeat ALLCITIES+={first} svgaddline(X(first), Y(first), X(NEXTC(first)), Y(NEXTC(first))) first:=NEXTC(first) until first=i ct+=1 end-if end-do ! Draw the cities svgaddgroup("C", "Cities", SVG_BLACK) forall(i in CITIES) svgaddtext(X(i), Y(i), text(i)) svgaddtext(minx-5,miny-5, "Tour length: "+getobjval) svgsetgraphscale(5) svgrefresh ! Uncomment to pause at every iteration: if ct>1 then svgpause; end-if end-procedure end-model