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

   file d1wagon.mos
   ````````````````
   Load balancing of train wagons
   
   (c) 2002 Dash Associates
       author: S. Heipcke, Mar. 2002
*******************************************************!)

model "D-1 Wagon load balancing"
 uses "mmxprs"

 forward function solve_heur:real
 forward procedure shell_sort(A:array(range) of integer)

 declarations   
  BOXES = 1..16                      ! Set of boxes
  WAGONS = 1..3                      ! Set of wagons

  WEIGHT: array(BOXES) of integer    ! Box weights
  WMAX: integer                      ! Weight limit per wagon

  load: array(BOXES,WAGONS) of mpvar ! 1 if box loaded on wagon, 0 otherwise
  maxweight: mpvar                   ! Weight of the heaviest wagon load
 end-declarations

 initializations from 'd1wagon.dat'
  WEIGHT WMAX
 end-initializations

! Solve the problem heuristically and terminate the program if the
! heuristic solution is good enough
 if solve_heur<=WMAX then
  writeln("Heuristic solution fits capacity limits")
  exit(0)
 end-if

! Every box into one wagon
 forall(b in BOXES) sum(w in WAGONS) load(b,w) = 1
 
! Limit the weight loaded into every wagon 
 forall(w in WAGONS) sum(b in BOXES) WEIGHT(b)*load(b,w) <= maxweight
 
! Bounds on maximum weight
 maxweight <= WMAX
 maxweight >= ceil((sum(b in BOXES) WEIGHT(b))/3)

 forall(b in BOXES,w in WAGONS) load(b,w) is_binary

! Alternative to lower bound on maxweight: adapt the optimizer cutoff value 
! setparam("XPRS_MIPADDCUTOFF",-0.99999)

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

! Minimize the heaviest load
 minimize(maxweight)
 
! Solution printing
 writeln("Optimal solution:\n Max weight: ", getobjval)
 forall(w in WAGONS) do
  write(" ", w, ":")
  forall(b in BOXES) write(if(getsol(load(b,w))=1, " "+b, ""))
  writeln(" (total weight: ", getsol(sum(b in BOXES) WEIGHT(b)*load(b,w)), ")")
 end-do

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

(! LPT (Longest processing time) heuristic: 
   One at a time place the heaviest unassigned box onto the wagon with 
   the least load
!)
   
 function solve_heur:real
  declarations
   ORDERW: array(BOXES) of integer      ! Box weights in decreasing order
   Load: array(WAGONS,range) of integer ! Boxes loaded onto the wagons 
   CurWeight: array(WAGONS) of integer  ! Current weight of wagon loads 
   CurNum: array(WAGONS) of integer     ! Current number of boxes per wagon
  end-declarations

 ! Copy the box weights into array ORDERW and sort it in decreasing order
  forall(b in BOXES) ORDERW(b):=WEIGHT(b)
  shell_sort(ORDERW)

 ! Distribute the loads to the wagons using the LPT heuristic 
  forall(b in BOXES) do
   v:=1                                 ! Find wagon with the smallest load
   forall(w in WAGONS) v:=if(CurWeight(v)<CurWeight(w), v, w)     
   CurNum(v)+=1                         ! Increase the counter of boxes on v
   Load(v,CurNum(v)):=b                 ! Add box to the wagon
   CurWeight(v)+=ORDERW(b)              ! Update current weight of the wagon
  end-do
 
  returned:= max(w in WAGONS) CurWeight(w)  ! Return the solution value

 ! Solution printing
  writeln("Heuristic solution:\n Max weight: ", returned)
  forall(w in WAGONS) do
   write(" ", w, ":")
   forall(b in 1..CurNum(w)) write(" ", Load(w,b))
   writeln(" (total weight: ", CurWeight(w), ")")
  end-do
  
 end-function

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

(! Sort an array in decreasing order using a Shell sort method:
   * First sort, by straight insertion, small groups of numbers. 
   * Next, combine several small groups and sort them (possibly 
     repeat this step). 
   * Finally, sort the whole list of numbers.

   The spacings between the numbers of groups sorted during each 
   pass through the data are called the increments. A good choice
   is the sequence which can be generated by the recurrence
   i(1)=1, i(k+1)=3i(k)+1, k=1,2,...
   
   The implementation assumes that the entries of the array to sort
   have the indices 1,...,N
!)   

 procedure shell_sort(A:array(range) of integer)
  N:=getsize(A)
  inc:=1                                ! Determine the starting increment
  repeat                         
   inc:=3*inc+1
  until (inc>N)  
 
  repeat                                ! Loop over the partial sorts
   inc:=inc div 3
   forall(i in inc+1..N) do             ! Outer loop of straight insertion
    v:=A(i)
    j:=i
    while (A(j-inc)<v) do               ! Inner loop of straight insertion
     A(j):=A(j-inc)
     j -= inc
     if j<=inc then break; end-if
    end-do
    A(j):= v     
   end-do  
  until (inc<=1)
 end-procedure

end-model