Implementation
We are using the following algorithm for modeling and solving this problem:
| Define the MIP machine assignment problem. Define the operations of the CP model. Start the MIP Branch-and-Bound search. At every node of the MIP search: while function generate_cuts returns true re-solve the LP-relaxation Function generate_cuts for all machines m call generate_cut_machine(m) if at least one cut has been generated Return true otherwise Return false Function generate_cut_machine(m) Collect all operations assigned to machine m if more than one operation assigned to m Solve the CP sequencing problem for m if sequencing succeeds Save the solution otherwise Add an infeasibility cut for machine m to the MIP |
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