Implementation

We are using the following algorithm for modeling and solving this problem:

The implementation of this model is split into two Mosel models: the first, sched_main.mos, contains the MIP master problem and the definition of the cut generation algorithm. The second model, sched_sub.mos, implements the CP single machine sequencing model.

The first part of the master model sets up the data arrays, compiles and loads the CP submodel, calls subroutines for the model definition and problem solving, and finally produces some summary result output. We have defined the filename of the data file as a parameter to be able to change the name of the data file at the execution of the model without having to change the model source. Correspondingly, all data, including the sizes of index sets, are read in from file. At first, we read in only the values of NP and NM. Subsequently, when declaring the sets and arrays that make use of these values, NP and NM are known and the arrays are created as fixed arrays. Otherwise, if their indexing sets are not known, these arrays would automatically be declared as dynamic arrays and for all but arrays of basic types (real, integer, etc.) we have to create their entries explicitly.

model "Schedule (MIP + CP) master problem"
 uses "mmsystem", "mmxprs", "mmjobs"

 parameters
  DATAFILE = "Data/sched_3_12.dat"
  VERBOSE = 1
 end-parameters

 forward procedure define_MIP_model
 forward procedure setup_cutmanager
 forward public function generate_cuts: boolean
 forward public procedure print_solution

 declarations
  NP: integer                             ! Number of operations (products)
  NM: integer                             ! Number of machines
 end-declarations

 initializations from DATAFILE
  NP NM
 end-initializations

 declarations
  PRODS = 1..NP                           ! Set of products
  MACH = 1..NM                            ! Set of machines
	
  REL: array(PRODS) of integer            ! Release dates of orders
  DUE: array(PRODS) of integer            ! Due dates of orders
  MAX_LOAD: integer                       ! max_p DUE(p) - min_p REL(p)
  COST: array(PRODS,MACH) of integer      ! Processing cost of products
  DUR: array(PRODS,MACH) of integer       ! Processing times of products
  starttime: real                         ! Measure program execution time
  ctcut: integer                          ! Counter for cuts
  solstart: array(PRODS) of integer
                                          ! **** MIP model:
  use: array(PRODS,MACH) of mpvar         ! 1 if p uses machine m, otherwise 0
  Cost: linctr                            ! Objective function

  totsolve,totCP: real                    ! Time measurement
  ctrun: integer                          ! Counter of CP runs
  CPmodel: Model                          ! Reference to the CP sequencing model
  ev: Event                               ! Event
  EVENT_SOLVED=2                          ! Event codes sent by submodels
  EVENT_FAILED=3
 end-declarations

 ! Read data from file
 initializations from DATAFILE
  REL DUE COST DUR
 end-initializations

! **** Problem definition ****
 define_MIP_model                         ! Definition of the MIP model
 res:=compile("sched_sub.mos")            ! Compile the CP model
 load(CPmodel, "sched_sub.bim")           ! Load the CP model

! **** Solution algorithm ****
 starttime:= gettime
 setup_cutmanager                         ! Settings for the MIP search

 totsolve:= 0.0
 initializations to "raw:"
  totsolve as "shmem:solvetime"
  REL as "shmem:REL" DUE as "shmem:DUE"
 end-initializations

 minimize(Cost)                           ! Solve the problem

 writeln("Number of cuts generated: ", ctcut)
 writeln("(", gettime-starttime, "sec) Best solution value: ", getobjval)
 initializations from "raw:"
  totsolve as "shmem:solvetime"
 end-initializations
 writeln("Total CP solve time: ", totsolve)
 writeln("Total CP time: ", totCP)
 writeln("CP runs: ", ctrun)

The MIP model corresponds closely to the mathematical model that we have seen in the previous section.

 procedure define_MIP_model

 ! Objective: total processing cost
  Cost:= sum(p in PRODS, m in MACH) COST(p,m) * use(p,m)

 ! Each order needs exactly one machine for processing
  forall(p in PRODS) sum(m in MACH) use(p,m) = 1

 ! Valid inequalities for strengthening the LP relaxation
  MAX_LOAD:= max(p in PRODS) DUE(p) - min(p in PRODS) REL(p)
  forall(m in MACH) sum(p in PRODS) DUR(p,m) * use(p,m) <= MAX_LOAD

  forall(p in PRODS, m in MACH) use(p,m) is_binary

 end-procedure

The cut generation callback function generate_cuts is called at least once per MIP node. For every machine, it checks whether the assigned operations can be scheduled or whether an infeasibility cut needs to be added. If any cuts have been added, the LP relaxation needs to be re-solved and the cut generation function will be called again, until no more cuts are added. It is important to set and re-set the values of XPRS_solutionfile as shown in our example at the beginning and end of this function if it accesses Xpress Optimizer solution values.
The function generate_cut_machine first collects all tasks that have been assigned to the given machine m into the set ProdMach by calling the procedure products_on_machine. If there are still unassigned tasks the returned set is empty, otherwise, if the set has more than one element it tries to solve the sequencing subproblem (function solve_CP_problem). If this problem cannot be solved, then the function adds a cut to the MIP problem that makes the current assignment of operations to this machine infeasible.

 procedure products_on_machine(m: integer, ProdMach: set of integer)

  forall(p in PRODS) do
   val:=getsol(use(p,m))
   if (val > 0 and val < 1) then
   ProdMach:={}
    break
   elif val>0.5 then
    ProdMach+={p}
   end-if
  end-do

 end-procedure

!-----------------------------------------------------------------
! Generate a cut for machine m if the sequencing subproblem is infeasible
 function generate_cut_machine(m: integer): boolean
  declarations
   ProdMach: set of integer
  end-declarations

 ! Collect the operations assigned to machine m
  products_on_machine(m, ProdMach)

 ! Solve the sequencing problem (CP model): if solved, save the solution,
 ! otherwise add an infeasibility cut to the MIP problem
  size:= getsize(ProdMach)
  returned:= false
  if (size>1) then
   if not solve_CP_problem(m, ProdMach, 1) then
    Cut:= sum(p in ProdMach) use(p,m) - (size-1)
    if VERBOSE > 2 then
     writeln(m,": ", ProdMach, " <= ", size-1)
    end-if
    addcut(1, CT_LEQ, Cut)
    returned:= true
   end-if
  end-if

 end-function

!-----------------------------------------------------------------
! Cut generation callback function
 public function generate_cuts: boolean
  returned:=false; ctcutold:=ctcut

  setparam("XPRS_solutionfile", 0)
  forall(m in MACH) do
   if generate_cut_machine(m) then
    returned:=true                    ! Call function again for this node
    ctcut+=1
   end-if
  end-do
  setparam("XPRS_solutionfile", 1)
  if returned and VERBOSE>1 then
   writeln("Node ", getparam("XPRS_NODES"), ": ", ctcut-ctcutold,
           " cut(s) added")
  end-if

 end-function

The solving of the CP model is started from the function solve_CP_problem that writes out the necessary data to shared memory and starts the execution of the submodel contained in the file sched_sub.mos.

 function solve_CP_problem(m: integer, ProdMach: set of integer,
                           mode: integer): boolean
  declarations
   DURm: dynamic array(range) of integer
   sol: dynamic array(range) of integer
   solvetime: real
  end-declarations

 ! Data for CP model
  forall(p in ProdMach) DURm(p):= DUR(p,m)
  initializations to "raw:"
   ProdMach as "shmem:ProdMach"
   DURm as "shmem:DURm"
  end-initializations

 ! Solve the problem and retrieve the solution if it is feasible
  startsolve:= gettime
  returned:= false
  if(getstatus(CPmodel)=RT_RUNNING) then
    fflush
    writeln("CPmodel is running")
    fflush
    exit(1)
  end-if

  ctrun+=1
  run(CPmodel, "NP=" + NP + ",VERBOSE=" + VERBOSE + ",MODE=" + mode)
  wait                                  ! Wait for a message from the submodel
  ev:= getnextevent                     ! Retrieve the event
  if getclass(ev)=EVENT_SOLVED then
   returned:= true
   if mode = 2 then
    initializations from "raw:"
     sol as "shmem:solstart"
    end-initializations
    forall(p in ProdMach) solstart(p):=sol(p)
   end-if
  elif getclass(ev)<>EVENT_FAILED then
   writeln("Problem with Kalis")
   exit(2)
  end-if
  wait
  dropnextevent                         ! Ignore "submodel finished" event
  totCP+= (gettime-startsolve)
 end-function

We complete the MIP model with settings for the cut manager and the definition of the integer solution callback. The Mosel comparison tolerance is set to a slightly larger value than the tolerance applied by Xpress Optimizer. It is important to switch the LP presolve off since we interfere with the matrix during the execution of the algorithm (alternatively, it is possible to fine-tune presolve to use only non-destructive algorithms). Sufficiently large space for cuts and cut coefficients should be reserved in the matrix. We also enable output printing by the Optimizer and choose among different MIP log frequencies (depending on model parameter VERBOSE.

 procedure setup_cutmanager
  setparam("XPRS_CUTSTRATEGY", 0)           ! Disable automatic cuts
  feastol:= getparam("XPRS_FEASTOL")        ! Get Optimizer zero tolerance
  setparam("zerotol", feastol * 10)         ! Set comparison tolerance of Mosel
  setparam("XPRS_PRESOLVE", 0)              ! Disable presolve
  setparam("XPRS_MIPPRESOLVE", 0)           ! Disable MIP presolve
  command("KEEPARTIFICIALS=0")              ! No global red. cost fixing
  setparam("XPRS_SBBEST", 0)                ! Turn strong branching off
  setparam("XPRS_HEURSTRATEGY", 0)          ! Disable MIP heuristics
  setparam("XPRS_EXTRAROWS", 10000)         ! Reserve space for cuts
  setparam("XPRS_EXTRAELEMS", NP*30000)     ! ... and cut coefficients
  setcallback(XPRS_CB_CUTMGR, "generate_cuts")  ! Define the cut mgr. callback
  setcallback(XPRS_CB_INTSOL, "print_solution") ! Define the integer sol. cb.
  setparam("XPRS_COLORDER", 2)
  case VERBOSE of
  1: do
      setparam("XPRS_VERBOSE", true)
      setparam("XPRS_MIPLOG", -200)
     end-do
  2: do
      setparam("XPRS_VERBOSE", true)
      setparam("XPRS_MIPLOG", -100)
     end-do
  3: do                                     ! Detailed MIP output
      setparam("XPRS_VERBOSE", true)
      setparam("XPRS_MIPLOG", 3)
     end-do
  end-case

 end-procedure

The definition of the integer solution callback is, in parts, similar to the function generate_cut_machine. To obtain a detailed solution output we need to re-solve all CP subproblems, this time with run MODE two, meaning that the CP model writes its solution information to shared memory.

 public procedure print_solution
  declarations
   ProdMach: set of integer
  end-declarations

  writeln("(",gettime-starttime, "sec) Solution ",
          getparam("XPRS_MIPSOLS"), ": Cost: ", getsol(Cost))

  if VERBOSE > 1 then
   forall(p in PRODS) do
    forall(m in MACH) write(getsol(use(p,m))," ")
    writeln
   end-do
  end-if

  if VERBOSE > 0 then
   forall(m in MACH) do
    ProdMach:= {}

  ! Collect the operations assigned to machine m
    products_on_machine(m, ProdMach)

    Size:= getsize(ProdMach)
    if Size > 1 then
   ! (Re)solve the CP sequencing problem
     if not solve_CP_problem(m, ProdMach, 2) then
      writeln("Something wrong here: ", m, ProdMach)
     end-if
    elif Size=1 then
     elem:=min(p in ProdMach) p
     solstart(elem):=REL(elem)
    end-if
   end-do

 ! Print out the result
   forall(p in PRODS) do
    msol:=sum(m in MACH) m*getsol(use(p,m))
    writeln(p, " -> ", msol,": ", strfmt(solstart(p),2), " - ",
            strfmt(DUR(p,round(msol))+solstart(p),2), "  [",
            REL(p), ", ", DUE(p), "]")
   end-do
   writeln("Time: ", gettime - starttime, "sec")
   writeln
   fflush
  end-if
 end-procedure

The following code listing shows the complete CP submodel. At every execution, the set of tasks assigned to one machine and the corresponding durations are read from shared memory. The disjunctions between pairs of tasks are posted explicitly to be able to stop the addition of constraints if an infeasibility is detected during the definition of the problem. The search stops at the first feasible solution. If a solution was found, it is passed back to the master model if the model parameter MODE has the value two. In every case, after termination of the CP search the submodel sends a solution status event back to the master model.

model "Schedule (MIP + CP) CP subproblem"
 uses "kalis", "mmjobs" , "mmsystem"

 parameters
  VERBOSE = 1
  NP = 12                                 ! Number of products
  MODE = 1                                ! 1 - decide feasibility
                                          ! 2 - return complete solution
 end-parameters

 startsolve:= gettime

 declarations
  PRODS = 1..NP                           ! Set of products
  ProdMach: set of integer
 end-declarations

 initializations from "raw:"
  ProdMach as "shmem:ProdMach"
 end-initializations

 finalize(ProdMach)  	

 declarations
  REL: array(PRODS) of integer            ! Release dates of orders
  DUE: array(PRODS) of integer            ! Due dates of orders
  DURm: array(ProdMach) of integer        ! Processing times on machine m
  solstart: array(ProdMach) of integer    ! Solution values for start times

  start: array(ProdMach) of cpvar         ! Start times of tasks
  Disj: array(range) of cpctr             ! Disjunctive constraints
  Strategy: array(range) of cpbranching   ! Enumeration strategy
  EVENT_SOLVED=2                          ! Event codes sent by submodels
  EVENT_FAILED=3
  solvetime: real
 end-declarations

 initializations from "raw:"
  DURm as "shmem:DURm" REL as "shmem:REL" DUE as "shmem:DUE"
 end-initializations

! Bounds on start times
 forall(p in ProdMach) setdomain(start(p), REL(p), DUE(p)-DURm(p))

! Disjunctive constraint
 ct:= 1
 forall(p,q in ProdMach| p<q) do
  Disj(ct):= start(p) + DURm(p) <= start(q) or start(q) + DURm(q) <= start(p)
  ct+= 1
 end-do

! Post disjunctions to the solver
 nDisj:= ct; j:=1; res:= true
 while (res and j<nDisj) do
  res:= cp_post(Disj(j))
  j+=1
 end-do

! Solve the problem
 if res then
  Strategy(1):= settle_disjunction(Disj)
  Strategy(2):= assign_and_forbid(KALIS_SMALLEST_DOMAIN, KALIS_MIN_TO_MAX,
                                 start)
  cp_set_branching(Strategy)
  res:= cp_find_next_sol
 end-if

! Pass solution to master problem
 if res then
  forall(p in ProdMach) solstart(p):= getsol(start(p))
  if MODE=2 then
   initializations to "raw:"
    solstart as "shmem:solstart"
   end-initializations
  end-if
  send(EVENT_SOLVED,0)
 else
  send(EVENT_FAILED,0)
 end-if

! Update total running time measurement
 initializations from "raw:"
  solvetime as "shmem:solvetime"
 end-initializations
 solvetime+= gettime-startsolve
 initializations to "raw:"
  solvetime as "shmem:solvetime"
 end-initializations

end-model