Hey guys,
I have been struggeling to speed up my optimization for quite a while now. Any help is highly appreciated!
In contrast to most of the research, I deal with entities (=worker) and not with machines. But since they are all the same to me, there is no problem. In my case a job can be assigned to multiple entities and the jobs duration depends on the amount of assigned entities. Furthermore an entity can be of type with exoskeleton or without (for the entire project).
The literature I found for symmetrie breaking does not work for me, since a job can be assigned to more then just one entity (e.g. lexicographically decreasing columns).
Any idea how I can further constrain the decision matrix “x” to get rid of the symmetry?
It is imperative as I am running into major performance issues by increasing the size of the process to be optimized.
# Decision matrix resource allocation [Exo: y/n ;jobs ;Entities]
@variable( bilap, x[1:2,1:jobs,1:Entities], Bin)
# Any entity is of type with exoskeleton or of 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
For those who are interested, here is the entire model/problem:
using vOptGeneric, JuMP, Gurobi, DataFrames, GLPK, VegaLite, PyPlot, Plots
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.1 )
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
global Front=filter(row -> row["Rang"] == 1, Total)
end
happy new year…