Vehicle routing

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

   file vrp.mos
   ````````````
   TYPE:         Vehicle routing problem (VRP)
   DIFFICULTY:   4
   FEATURES:     MIP problem, formulation of constraints to eliminate
                 inadmissible subtours, definition of model cuts,
                 selection  with `|', algorithm for printing the tours, 
                 graphical solution representation
   DESCRIPTION:  A transporter has to deliver heating oil from a refinery 
                 to a certain number of clients. The distances between the
                 clients and the refinery and the demands in liters for 
                 the different sites are given. The transport company uses 
                 tankers with a capacity of 39000 liters for the deliveries. 
                 Determine the tours for delivering to all clients that 
                 minimize the total number of kilometers driven.
   FURTHER INFO: `Applications of optimization with Xpress-MP', 
                 Section 10.4 `Heating oil delivery'
   
   (c) 2008 Fair Isaac Corporation
       author: S. Heipcke, 2002, rev. Sep. 2017
*******************************************************!)

model "Heating oil delivery"
 uses "mmxprs", "mmsvg"

 declarations
  NS = 7
  SITES = 1..NS                       ! Set of locations, 1=refinery
  CLIENTS = 2..NS

  DEM: array(SITES) of integer        ! Demands
  DIST: array(SITES,SITES) of integer ! Distances between locations
  CAP: integer                        ! Lorry capacity

  prec: array(SITES,SITES) of mpvar   ! 1 if i immediately precedes j, 
                                      ! 0 otherwise 
  quant: array(CLIENTS) of mpvar      ! Quantity delivered up to i
 end-declarations

 initializations from 'vrp.dat'
  DEM DIST CAP
 end-initializations

! Objective: total distance driven
 Length:= sum(i,j in SITES | i<>j) DIST(i,j)*prec(i,j)

! Enter and leave every city only once (except the depot)
 forall(j in CLIENTS) EnterOnce(j):= sum(i in SITES| i<>j) prec(i,j) = 1
 forall(i in CLIENTS) LeaveOnce(i):= sum(j in SITES| i<>j) prec(i,j) = 1

! If i is the first client of a tour, then quant(i)=DEM(i)
 forall(i in CLIENTS) QuantFirst(i):= quant(i) <= CAP + (DEM(i)-CAP)*prec(1,i)

! If j comes just after i in a tour, then quant(j) is greater than the 
! quantity delivered during the tour up to i plus the quantity to be 
! delivered at j (to avoid loops and keep capacity limit of the tanker)
 forall(i,j in CLIENTS| i<>j) 
  QuantSucc(i,j):= quant(j) >= quant(i) + DEM(j) - CAP +
                             CAP*prec(i,j) + (CAP-DEM(j)-DEM(i))*prec(j,i)

! Additional constraints:
! If i is not the first client of a tour, quant(i) is larger than the sum 
! of the quantities to deliver to i and to his predecessor on the tour
(!
 declarations
  modcut: array(CLIENTS) of linctr
 end-declarations

 forall(i in CLIENTS) do
  modcut(i):= quant(i) >= DEM(i) + sum(j in SITES| i<>j) DEM(j)*prec(j,i)
  setmodcut(modcut(i))
 end-do
!)

 forall(i in CLIENTS) do
  quant(i) <= CAP
  quant(i) >= DEM(i)
 end-do
 
 forall(i,j in SITES | i<>j) prec(i,j) is_binary

! Uncomment the following line to see the Optimizer log
! setparam("XPRS_VERBOSE",true)

! Solve the problem
 minimize(Length)

! Solution printing
 writeln("Total distance: ", getobjval)
 forall(i in CLIENTS)
  if(getsol(prec(1,i))>0) then
   ct:=DEM(i)
   writeln(1, " -> ", i)
   p:=i
   while(p<>1) do
    n:= integer(round(sum(j in SITES) j*getsol(prec(p,j))))
    writeln(p, " -> ", n)
    ct+=DEM(n)
    p:=n
   end-do
   writeln("Quantity delivered: ", ct)
  end-if

! Solution drawing
 declarations
  X,Y: array(SITES) of integer        ! x-y-coordinates of sites
 end-declarations
 
 initializations from 'vrp.dat'
  [X,Y] as 'POS'
 end-initializations

 minx:=min(i in SITES)X(i)-15; miny:=min(i in SITES)Y(i)-15
 svgsetgraphviewbox(0, 0,
         max(i in SITES)X(i)+15, max(i in SITES)Y(i)+25)
 svgsetgraphscale(2)
 svgsetgraphpointsize(2)

 svgaddgroup("Sites", "Clients", SVG_BLACK)
 forall(i in CLIENTS) do
  svgaddpoint(X(i), Y(i))
  svgaddtext(X(i)+0.5, Y(i)+1, text(i-1))
 end-do
 svgaddgroup("Ref", "Refinery", SVG_BROWN)
 svgaddpoint(X(1), Y(1))
 svgaddtext(X(1)+0.5, Y(1)-5, "Refinery")

 svgaddgroup("Routes", "Delivery routes")
 forall(i in CLIENTS)
  if(getsol(prec(1,i))>0) then
   svgaddarrow(X(1), Y(1), X(i), Y(i))
   p:=i
   while(p<>1) do
    n:= integer(round(sum(j in SITES) j*getsol(prec(p,j))))
    svgaddarrow(X(p), Y(p), X(n), Y(n))
    p:=n
   end-do
  end-if

 svgsave("vrp.svg")
 svgrefresh
 svgwaitclose("Close browser window to terminate model execution.", 1)
end-model