Combining CP and MIP

Application problems often combine different subproblems that are solved better with one or another solver, making the complete problem difficult or unmanageable for a single solver. A typical example is production planning and scheduling applications. The long-term (planning) aspects are usually more easily handled by LP solvers whereas the short-term (scheduling) subproblems are better suited for CP solvers. Solving these two parts completely independent of each other may lead to infeasible scheduling subproblems or plans that do not correspond to the reality of production. A possible solution to this dilemma is to iteratively solve LP planning problems and CP scheduling problems, until a feasible schedule for the planned quantities is obtained.

Another method of combined MIP-CP problem solving that provides a tighter integration of the two techniques consists of solving CP subproblems for generating cuts at the nodes of a MIP Branch-and-Bound search. This technique has already been applied successfully to several large-scale planning and scheduling applications by PSA and BASF (see for instance the description of hybrid MIP-CP algorithms implemented with Mosel in [BP03] and [Sad04]). This type of combination is more technical than sequential CP and LP/MIP solving since it requires the developer of the algorithm to interact with the MIP search at every node.

Cut generation algorithms can be implemented with the help of the Xpress Optimizer callbacks (see the `Mosel Language Reference Manual' for the definition of callbacks with Mosel and the `Xpress Optimizer Reference Manual' for an explanation of the Optimizer callback functions).

The original description of the example in this section was published in [JG99]. A prototype implementation was developed by N.  Pisaruk in the context of the EU-project LISCOS.

Example: Machine assignment and sequencing

We need to produce 12 products on a set of three machines. Each machine may produce all of the products but processing times and costs vary (Table Machine-dependent production costs and durations). Furthermore, for every product we are given its release and due dates (Table Release dates and due dates of products). We wish to determine a production plan for all products that minimizes the total production cost.

Model formulation

We are going to represent this problem by two subproblems: the machine assignment problem and the sequencing of operations on machines. The former is implemented by a MIP model; the latter is formulated as a CP (single machine) problem.

MIP model

Let COST_pm denote the production cost and DUR_pm the processing time of product p (p ∈ PRODS, the set of products) on machine m (m ∈ MACH, the set of machines).

To formulate the machine assignment problem we introduce binary variables use_pm that take the value 1 if product p is produced on machine m and zero otherwise. The objective function is then given as

The assignment constraints expressing that each order needs exactly one machine for processing it are defined as follows:

In addition to these constraints that already fully state the problem we may define some additional constraints to strengthen the LP relaxation. All production takes place between the earliest release date and the latest due date. If we denote the length of this interval by MAX_LOAD, we may formulate the following valid inequalities expressing that the total processing time of products assigned to a machine cannot exceed MAX_LOAD:

CP model

Once the set of operations assigned to machine m, ProdMach_m (ProdMach_m ⊆ PRODS), is known, we obtain the following sequencing problem for this machine:

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

Results

The best solution produced for the data set sched_3_12 is the following :

Cost: 92
1 -> 3:  2 - 15  [2, 32]
2 -> 3: 15 - 23  [4, 33]
3 -> 2: 15 - 32  [5, 36]
4 -> 1: 24 - 30  [7, 37]
5 -> 2: 32 - 38  [9, 39]
6 -> 2:  0 -  3  [0, 34]
7 -> 1:  3 - 13  [3, 30]
8 -> 1: 16 - 24  [6, 26]
9 -> 3: 23 - 36  [11, 36]
10 -> 1: 30 - 38  [2, 38]
11 -> 2:  3 - 15  [3, 31]
12 -> 1: 13 - 16  [4, 22]

A total of 1604 cuts are added to the MIP problem by 2691 CP model runs and the Branch-and-Bound search explores 12295 nodes. Optimality is proven within a few seconds on a Pentium IV PC.

It is possible to implement this problem entirely either with Xpress Optimizer or with Xpress Kalis. However, already for this three machines – 12 jobs instance the problem is extremely hard for either technique on its own. With CP it is difficult to prove optimality and with MIP the formulation of the disjunctions makes the definition of a large number of binary variables necessary (roughly in the order of number_of_machines · number_of_products²) which makes the problem impracticable to deal with.

Parallel solving of CP subproblems

Instead of solving the CP single-machine subproblems at every MIP node sequentially, we can modify our Mosel models to solve the subproblems in parallel—especially when working on a multiprocessor machine this may speed up the cut generation process and hence shorten the total run time. We modify the algorithm of Section Implementation as follows:

The modified version of the function generate_cuts looks as follows. For the full example code the reader is referred to the set of User Guide examples provided with the Xpress Kalis distribution (files sched_mainp.mos and sched_subp.mos).

! Collect the operations assigned to machine m
 procedure products_on_machine(m: integer)

  NumOp(m):=0
  forall(p in PRODS) do
   val:=getsol(use(p,m))
   if (! not isintegral(use(p,m)) !) (val > 0 and val < 1) then
    NumOp(m):=0
    break
   elif val>0.5 then
    NumOp(m)+=1
    OpMach(m,NumOp(m)):= p
   end-if
  end-do

 end-procedure

!-----------------------------------------------------------------
! Add an infeasibility cut for machine m to the MIP problem
 procedure add_cut_machine(m: integer)

  Cut:= sum(p in 1..NumOp(m)) use(OpMach(m,p),m) - (NumOp(m)-1)
  if VERBOSE > 1 then
   write(m,": ")
   forall(p in 1..NumOp(m)) write(OpMach(m,p), " ")
   writeln(" <= ", NumOp(m)-1)
  end-if
  addcut(1, CT_LEQ, Cut)

 end-procedure

The implementation of the CP submodels remains largely unchanged, with the exception of the labels employed for passing data via shared memory: we append the machine index to every data item to be able to distinguish between the data used by the different subproblems running in parallel.

For the data set sched_3_12.dat we have observed only a few percent decrease of the total running time on a dual processor machine using the parallel implementation: in many nodes only a single CP subproblem is solved and if there are several subproblems to be solved their execution may be of quite different length. For instances with a larger number of machines the parallelization is likely to show more effect.