Hey guys,
I am working on an multi-objective optimization model (cost, strain, time). The challenge is simultaneous resource allocation and scheduling, as there are corresponding dependencies. However, I have developed a methodology to meet all the requirements. Below is the entire code (maybe it will be useful to someone).
My problem is the performance, since I do not have to optimize only four jobs, as in the attached test case, but about 30 jobs. With each additional job the computing time increases dramatically.
Benchmarking Model with 4 Jobs:
Benchmarking Model with 5 Jobs:
With more than 7 jobs I could not solve the problem yet, even after several hours. Of course, I have gone through the Julia Performance Tips, but unfortunately I do not get further.
Unfortunately, the profiling macro keeps crashing the kernel when I try to run this model, so I can’t provide you with the output
My question is do you have any ideas on how to tackle the performance problem.
Or what your approach would be?
-
Symmetry breaking?
-
Further Gurobi Parameters?
-
Necessity of heuristics?
Thank you so much!! Any help highly appreciated
Model:
using vOptGeneric, JuMP, Gurobi, DataFrames, GLPK, VegaLite
function Pareto(Epsilon::Int)
jobs=4
Entities=3
planning_horizon=100
#costs per minute
minute_costs_p_worker=1
minute_costs_p_worker_Exo=1.1
minute_costs_p_foreman=2
minute_costs_p_foreman_Exo=2.1
# Predecessor of job 2 e.g. is job 1
P1=[0]
P2=[1]
P3=[1,2]
P4=[2,3]
P=[P1,P2,P3,P4]
# Delays between jobs / quasi-paralellizable jobs
# x axis: Successor ; y axis: Predecessor --> job 3 can start after 20% of job 2 are done
Lag=[0 0 0 0
0 0 0.8 0
0 0 0 0
0 0 0 0]
# Duration
M=1000 # BIG M: Headcount too low for job
N=10000 # BIG N: Headcount too high for / no further advantage
# Column index is number of allocated persons to job, figure is time needed for job with given headcount
D1 = [10 6 4 N]
D2 = [5 N N N]
D3 = [M 8 5 N]
D4 = [15 12 8 N]
Duration= [D1,D2,D3,D4]
# Strain in KJ/min. per job
S1 = 18.5
S2 = 10
S3 = 15
S4 = 13
Strain=[S1,S2,S3,S4]
# Reduction of the strain through the use of an exoskeleton
R1=0.8
R2=1
R3=0.8
R4=0.8
R=[R1,R2,R3,R4]
Strain_Exo= Strain.*R
# Foreman alloc. mandatory for job
F1 =0
F2 =1
F3 =0
F4 =0
FO=[F1,F2,F3,F4]
bilap = vModel( Gurobi.Optimizer; add_bridges = false)
set_optimizer_attribute(bilap, "TimeLimit", 100)
set_optimizer_attribute(bilap, "MIPFocus", 3)
set_optimizer_attribute(bilap, "Presolve", 2)
set_optimizer_attribute(bilap, "MIPGap", 0.02)
@variable( bilap, x[1:2,1:jobs,1:Entities], Bin) # Decision matrix resource allocation [Exo: y/n ;jobs ;Entities]
@variable( bilap, t[1:jobs], Int) # Start time
@variable( bilap, Aux_Exo[1:Entities,1:2], Int) # Auxiliary variable equipment (exoskeleton y/n per entity).
@variable( bilap, Aux_Dur[1:Entities,1:jobs], Bin) # Auxiliary variable duration
@variable( bilap, Dur_p_job[1:jobs], Int) # Duration per job (after resource allocation Decision)
@variable( bilap, S[1:2,1:jobs,1:Entities]) # Strain / job / Entity in KJ (after resource allocation Decision)
@variable( bilap, Max) # Strain of the most heavily loaded employee in KJ
@variable( bilap, H_Wo[1:Entities-1], Bin) # Headcount worker
@variable( bilap, H_WoEx[1:Entities-1], Bin) # Headcount worker with exoskeleton
@variable( bilap, H_Fo, Bin) # Headcount Foreman
@variable( bilap, H_FoEx, Bin) # Headcount Foreman with exoskeleton
@variable( bilap, H_normal,Int) # Headcount Entities without exoskeleton
@variable( bilap, H_Exos,Int) # Headcount Entities with exoskeleton
@variable( bilap, start[1:planning_horizon,1:jobs], Bin) # Binary sequencing variable, decides when to start job (as soon as value =0 --> start)
@variable( bilap, stop[1:planning_horizon,1:jobs], Bin) # Binary sequencing variable, decides when to stop job (as soon as value =0 --> stop)
@variable( bilap, Aktiv[t=1:planning_horizon,i=1:jobs], Bin) # job i beeing activ at time t? (stop-start)
@variable( bilap, totalduration,Int) # Total duration of the project (Until last job completed)
# Any entity over total duration of the type with exoskeleton or of the type without exoskeleton
@constraint( bilap , [a=1:Entities,o=1:2],Aux_Exo[a,o]==sum(x[o,j,a] for j in 1:jobs))
for a = 1:Entities
@constraint(bilap, Aux_Exo[a,:] in SOS1())
end
# Store allocated entities per job
@expression( bilap , H_Entities[i=1:jobs], sum(x[a,i,p] for a in 1:2, p in 1:Entities))
# Foreman is mandatory for the job (1. and only 1. Entity is foreman)
@constraint( bilap , [j=1:jobs], sum(x[i,j,1] for i in 1:2) >= FO[j])
# Allocate as many entities as possible at maximum
@constraint( bilap , [a=1:jobs], H_Entities[a] <= length(Duration[a]) - count(i->(i== N),Duration[a]))
# Allocate as many entities as needed at minimum
@constraint( bilap , [a=1:jobs], H_Entities[a] >= 1 + count(i->(i== M),Duration[a]))
# Realized Duration
# binary aux. variable can only become 1 once --> Index where this happens determines number of allocated entities
@constraint( bilap , [i=1:jobs],Dur_p_job[i]==sum(Aux_Dur[j,i]*Duration[i][j] for j in 1:Entities))
@constraint( bilap , [i=1:jobs],sum(Aux_Dur[:,i])==1)
@constraint( bilap , [a=1:jobs],H_Entities[a]==sum(Aux_Dur[j,a]*j for j in 1:Entities))
# Scheduling
@constraint( bilap , [i=1:jobs], t[i]>=0)
@constraint( bilap , [i=1:jobs], Dur_p_job[i]>=0)
# Sequence condition + reduction of time in case of quasi-parallelization (Lag)
@constraint( bilap , [j = 2:jobs, h in P[j]], t[h]+Dur_p_job[h]*(1-Lag[h,j]) <= t[j])
# Capacity restrictions
# start and stop are monotonic {0, 1} variables that start at 1 and fall to 0.
# The job starts when start => 0, and it ends when stop => 0.
# If both values are 1, the job has not yet started. If they are both 0, the job is stopped.
# Assign each entity to only one job at any given time
@constraint( bilap , [t=2:planning_horizon, i=1:jobs], start[t-1,i] >= start[t,i])
@constraint( bilap , [t=2:planning_horizon, i=1:jobs], stop[t-1,i] >= stop[t,i])
@constraint( bilap , [t=1:planning_horizon,i=1:jobs], Aktiv[t,i]==(stop[t,i] - start[t,i]))
@constraint( bilap , [i=1:jobs], sum(Aktiv[t,i] for t=1:planning_horizon) == Dur_p_job[i])
@constraint( bilap , [i=1:jobs], sum(start[t,i] for t=1:planning_horizon) == t[i])
@constraint( bilap , [w = 1:Entities, o=1:planning_horizon], sum(Aktiv[o, j] * x[i,j,w] for j = 1:jobs, i=1:2 ) <= 1)
# Collect information about used resources from the decision matrix
@constraint( bilap , [i=1:Entities-1,j=1:jobs], H_Wo[i]>=x[1,j,i+1])
@expression( bilap , Agg_cost_Wo, sum(H_Wo)*minute_costs_p_worker)
@constraint( bilap , [i=1:Entities-1,j=1:jobs], H_WoEx[i]>=x[2,j,i+1])
@expression( bilap , Agg_cost_WoEx, sum(H_WoEx)*minute_costs_p_worker_Exo)
@constraint( bilap , [j=1:jobs],H_Fo>=x[1,j,1])
@expression( bilap , Agg_cost_Fo, H_Fo*minute_costs_p_foreman)
@constraint( bilap , [j=1:jobs],H_FoEx>=x[2,j,1])
@expression( bilap , Agg_cost_FoEx, H_FoEx*minute_costs_p_foreman_Exo)
@constraint( bilap , H_normal== sum(H_Wo) + H_Fo )
@constraint( bilap , H_Exos==sum(H_WoEx) + H_FoEx)
@constraint( bilap , [i=1:jobs], S[1,i,:].==x[1,i,:]*Strain[i]*Dur_p_job[i])
@constraint( bilap , [i=1:jobs], S[2,i,:].==x[2,i,:]*Strain_Exo[i]*Dur_p_job[i])
# Aggregate the strain across the respective Entity and find most loaded employee
@expression( bilap , S_p_E[i=1:Entities], sum(S[:,:,i]))
@constraint( bilap , [i=1:Entities], Max>= S_p_E[i])
# Project duration
# last job terminates project
@constraint( bilap , totalduration==sum(stop[t,jobs] for t=1:planning_horizon))
@constraint( bilap , totalduration==Epsilon)
# Global standstills are excluded during the project duration (no global breaks)
@constraint( bilap , [o=1:Epsilon],sum(Aktiv[o, j] for j in 1:jobs)>=1)
# Relate cost and load to time
@expression( bilap , Cost, Epsilon *(Agg_cost_Wo+Agg_cost_WoEx+Agg_cost_Fo+Agg_cost_FoEx))
@expression( bilap , Strain_max, Max/Epsilon)
@addobjective( bilap ,Min, Cost)
@addobjective( bilap ,Min,Strain_max)
vSolve( bilap, method=:epsilon, step=0.3 )
Y_N = getY_N( bilap )
Cost=[]
Duration=[]
Strain=[]
Allocation=[]
Schedule=[]
Headcount_without_Exo=[]
Headcount_with_Exo=[]
Duration_p_job=[]
starting=[]
for i in 1:length(Y_N)
push!(Cost,Y_N[i][1]);
push!(Strain,Y_N[i][2]);
push!(Duration,Epsilon);
push!(Allocation,value.(x,[i]));
push!(Schedule,value.(Aktiv,[i]));
push!(Headcount_without_Exo,value.(H_normal,[i]));
push!(Headcount_with_Exo,value.(H_Exos,[i]));
push!(Duration_p_job,value.(Dur_p_job,[i]));
push!(starting,value.(t,[i]));
end
return Y_N,Cost, Strain, Duration, Allocation, Schedule, Headcount_without_Exo, Headcount_with_Exo, Duration_p_job, starting
end
function Test()
Cost=[]
Duration=[]
Strain=[]
Allocation=[]
Schedule=[]
Headcount_without_Exo=[]
Headcount_with_Exo=[]
Duration_p_job=[]
starting=[]
u_lbound=10
u_ubound=60
stepsize=1
steps=convert(Int,(u_ubound-u_lbound)/stepsize+1)
epsilon=u_ubound;
while epsilon>=u_lbound
KPI=Pareto(epsilon)
for i in 1:length(KPI[1])
push!(Cost,KPI[2][i]);
push!(Strain,KPI[3][i]);
push!(Duration,KPI[4][i]);
push!(Allocation,KPI[5][i]);
push!(Schedule,KPI[6][i]);
push!(Headcount_without_Exo,KPI[7][i]);
push!(Headcount_with_Exo,KPI[8][i]);
push!(Duration_p_job,KPI[9][i]);
push!(starting,KPI[10][i]);
end
epsilon-=stepsize;
end
Modus=[]
for i=1:length(Headcount_without_Exo)
push!(Modus,string(Int(round(Headcount_without_Exo[i][1])), " w/o Exo; ", Int(round(Headcount_with_Exo[i][1]))," Exo "))
end
Rang=ones(length(Headcount_without_Exo))
for j in 1:length(Headcount_without_Exo)
for i in 1:length(Headcount_without_Exo)
if j!=i
if Duration[j]>=Duration[i] && Strain[j]>=Strain[i] && Cost[j]>=Cost[i]
Rang[j]=2
end
end
end
end
# local pareto optima
Total=DataFrame("Duration"=>Duration,"Cost (€)"=>Cost,"Strain (KJ/min.)"=>Strain,"Rang"=>Rang,"Modus"=>Modus,"Duration per job"=>Duration_p_job,"start"=>starting,"Schedule"=>Schedule,"Allocation"=>Allocation)
# global pareto optima
Front=filter(row -> row["Rang"] == 1, Total)
end
Section of the model output:
