Hey guys,
I am working on an multiobjective 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 / quasiparalellizable 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:Entities1], Bin) # Headcount worker
@variable( bilap, H_WoEx[1:Entities1], 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? (stopstart)
@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 quasiparallelization (Lag)
@constraint( bilap , [j = 2:jobs, h in P[j]], t[h]+Dur_p_job[h]*(1Lag[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[t1,i] >= start[t,i])
@constraint( bilap , [t=2:planning_horizon, i=1:jobs], stop[t1,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:Entities1,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:Entities1,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_uboundu_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: