# Symmetry breaking: scheduling problems with job duration dependent on ressource allocation

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=
P2=
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)

vSolve( bilap, method=:epsilon, step=0.1 )
Y_N = getY_N( bilap )

Cost=[]
Duration=[]
Strain=[]
Allocation=[]
Schedule=[]
Duration_p_job=[]
starting=[]
for i in 1:length(Y_N)
push!(Cost,Y_N[i]);
push!(Strain,Y_N[i]);
push!(Duration,Epsilon);
push!(Allocation,value.(x,[i]));
push!(Schedule,value.(Aktiv,[i]));
push!(Duration_p_job,value.(Dur_p_job,[i]));
push!(starting,value.(t,[i]));
end
end

function Test()
Cost=[]
Duration=[]
Strain=[]
Allocation=[]
Schedule=[]
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)
push!(Cost,KPI[i]);
push!(Strain,KPI[i]);
push!(Duration,KPI[i]);
push!(Allocation,KPI[i]);
push!(Schedule,KPI[i]);
push!(Duration_p_job,KPI[i]);
push!(starting,KPI[i]);
end
epsilon-=stepsize;
end

Modus=[]
end

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



I know nothing about the problem, but these are all arrays of type Any , and this is probably very bad for performance. Try using, for example, Float64[] instead.

2 Likes

and this is probably very bad for performance. Try using, for example, Float64[] instead.

Although a reasonable suggestion in general, this is not the problem in this particular case. The time is spent I the optimization problem, not pushing to Julia arrays.

Any idea how I can further constrain the decision matrix “x” to get rid of the symmetry?

I’d start by simplifying the problem. Get rid of the multi objective, and focus on a single objective problem. Does that solve quickly? Is the problem the multi objective solve or the single objective solve?

There are a few things you could improve.

Start by adding bounds to every integer variable. Instead of @constraint( bilap , [i=1:jobs], t[i]>=0), do @variable( bilap, t[1:jobs] >= 0, Int). The former adds a new row to the constraint matrix. The latter modifies the variable bound. You should also add a good upper bound to every variable.

There are also constraints like this:

@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


If I understand, you’re trying to say that for each a, only one o can have a non-zero sum of x[o, :, a], and since x is binary, this is really saying that the other x must be zero. This is quite weak because you have an auxiliary variable integer Aux_Exo with no bounds, and you’re using an SOS1 constraint on the linear sum. But because the x are binary, we know that the sum can be zero iff every element is zero. Gurobi will probably presolve that out, but it might cause some issues.

So I think an equivalent, and slightly stronger formulation is:

    # For each a, only one o can be non-zero
@variable(bilap, z[a=1:Entities, o=1:2], Bin)
@constraint(bilap, [a=1:Entities], sum(z[a, o] for o in 1:2) == 1)
@constraint(bilap, [o=1:2, j=1:jobs, a=1:Entities], x[o, j, a] <= z[a, o])


There are probably a bunch of other examples, but I didn’t look through in detail. In general though, the monotonic start-stop formulation is a good one.

1 Like

Thanks for these tips@odow! They have already helped me a lot.
I have now broken the problem down into sub-problems and found that the scheduling logic in particular is problematic when scaling.
When running the attached model with the small case (4 jobs) it takes only 3-4 sec. But using the real case parameters (as in code below) results in hours of processing time (have not solved yet).
There are 147.456 possible duration per job combinations. Could this be the problem and if so, do you have any idea how to solve it?

any suggestions are greatly appreciated happy new year…

Sub-problem (The focus is solely on minimizing overall duration):

using vOptGeneric, JuMP, Gurobi, DataFrames, GLPK, VegaLite, PyPlot

jobs=22
Entities=3
planning_horizon=1000

#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	=	
P2	=	
P3	=	
P4	=	
P5	=	[2,4]
P6	=	
P7	=	[4,6]
P8	=	[6,5,7]
P9	=	
P10	=	
P11	=	[4,2,10]
P12	=	[8,11]
P13	=	
P14	=	
P15	=	[14,4]
P16	=	[9,12,15]
P17	=	
P18	=	[17,13]
P19	=	
P20	=	[9,13,12]
P21	=	
P22	=	[18,19,20,21]
P=[P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15,P16,P17,P18,P19,P20,P21,P22]

# 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	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0.8	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0.8	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0.8	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0.8	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	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	=[	0	N	N	N	]
D2	=[	35	N	N	N	]
D3	=[	M	35	25	N	]
D4	=[	M	30	25	N	]
D5	=[	18	N	N	N	]
D6	=[	15	9	N	N	]
D7	=[	10	6	N	N	]
D8	=[	80	60	50	N	]
D9	=[	30	20	N	N	]
D10	=[	60	50	N	N	]
D11	=[	M	40	30	N	]
D12	=[	M	60	50	N	]
D13	=[	40	30	N	N	]
D14	=[	20	10	8	N	]
D15	=[	M	50	40	N	]
D16	=[	M	60	50	N	]
D17	=[	20	14	N	N	]
D18	=[	M	30	25	N	]
D19	=[	6	3	N	N	]
D20	=[	10	N	N	N	]
D21	=[	10	N	N	N	]
D22	=[	1	N	N	N	]
Duration=[D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,D13,D14,D15,D16,D17,D18,D19,D20,D21,D22]

# Foreman alloc. mandatory for job
F1	=	0
F2	=	1
F3	=	0
F4	=	0
F5	=	1
F6	=	0
F7	=	0
F8	=	1
F9	=	0
F10	=	0
F11	=	0
F12	=	1
F13	=	0
F14	=	0
F15	=	0
F16	=	1
F17	=	0
F18	=	0
F19	=	0
F20	=	1
F21	=	1
F22	=	0
FO=[F1,F2,F3,F4,F5,F6,F7,F8,F9,F10,F11,F12,F13,F14,F15,F16,F17,F18,F19,F20,F21,F22]

bilap = Model( 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, 1000>=t[1:jobs]>=0, 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, 100>=Dur_p_job[1:jobs]>=0, Int)                              # Duration per job (after resource allocation Decision)
@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)

# Any entity over total duration of the type with exoskeleton or of the type without exoskeleton
@constraint( bilap , [a=1:Entities],sum(Aux_Exo[a,o] for o in 1:2)==1)
@constraint(bilap, [o=1:2, j=1:jobs, a=1:Entities], x[o, j, a] <= Aux_Exo[a, o])

# 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
# 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])
# last index of list is last job
@expression( bilap , totalduration,sum(stop[t,jobs] for t=1:planning_horizon))

# 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)

@objective( bilap ,Min, totalduration)

optimize!(bilap)
value.(totalduration)


Your other option is to consider a column generation/set covering model. I simplified things by removing the exoskeleton stuff (I don’t understand what it’s doing?). But this should point you in the right direction. It’s a nice model, so I should write this up as a tutorial for the documentation.

using JuMP, Gurobi

N_FOREMAN = 1
N_WORKERS = 2

P = [
Int[],
Int[],
,
[1, 3],
,
[3, 5],
[4, 5, 6],
,
,
[1, 3, 9],
[7, 10],
,
,
[13, 3],
[8, 11, 14],
,
[12, 16],
,
[8, 11, 12],
,
[17, 18, 19, 20],
]

Duration = [
35 missing missing
missing 35 25
missing 30 25
18 missing missing
15 9 missing
10 6 missing
80 60 50
30 20 missing
60 50 missing
missing 40 30
missing 60 50
40 30 missing
20 10 8
missing 50 40
missing 60 50
20 14 missing
missing 30 25
6 3 missing
10 missing missing
10 missing missing
1 missing missing
]

# 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	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0.8	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0.8	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0.8	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0.8	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
]

N_JOBS = length(P)
T = sum(minimum(coalesce.(Duration[j, :], 10_000)) for j in 1:N_JOBS)

# Foreman alloc. mandatory for job
F = [0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0]

model = Model(Gurobi.Optimizer)
x = [VariableRef[] for j in 1:N_JOBS]
@variable(model, 1 <= start_times[1:N_JOBS] <= T)
@variable(model, 1 <= finish_times[1:N_JOBS] <= T)
@variable(model, make_span)
@constraints(model, begin
c_n_workers[1:T], 0 <= N_WORKERS
c_n_foreman[1:T], 0 <= N_FOREMAN
c_jobs[j = 1:N_JOBS], 0 == 1
c_foreman[j = 1:N_JOBS], 0 >= F[j]
c_start_time[j = 1:N_JOBS], 0 == start_times[j]
c_finish_time[j = 1:N_JOBS], 0 == finish_times[j]
# Precedence
[j = 1:N_JOBS, k = P[j]], start_times[j] >= finish_times[k] + 1 - Lag[j, k] * (finish_times[k] - start_times[k])
make_span .>= finish_times
end)
@objective(model, Min, make_span)

model::Model,
job_id::Int,
start_time::Int,
n_workers::Int,
n_foreman::Int,
)
c = @variable(
model,
binary = true,
base_name = "x[$job_id,$start_time,$n_workers,$n_foreman]",
)
push!(x[job_id], c)
set_normalized_coefficient(c_jobs[job_id], c, 1.0)
duration = Duration[job_id, n_workers + n_foreman]
if n_workers > 0
for t in start_time:(start_time + duration - 1)
set_normalized_coefficient(c_n_workers[t], c, n_workers)
end
end
if n_foreman > 0
for t in start_time:(start_time + duration - 1)
set_normalized_coefficient(c_n_foreman[t], c, n_foreman)
end
set_normalized_coefficient(c_foreman[job_id], c, 1.0)
end
set_normalized_coefficient(c_start_time[job_id], c, start_time)
set_normalized_coefficient(c_finish_time[job_id], c, start_time + duration - 1)
return
end

# For small problems, we can enumerate very possible column
n = 0
for job_id in 1:N_JOBS, n_foreman in 0:N_FOREMAN, n_workers in 0:N_WORKERS
n_total = n_foreman + n_workers
if n_total == 0 || ismissing(Duration[job_id, n_total])
continue
end
duration = Duration[job_id, n_total]
# Could improve this by computing an earliest start date instead of assuming
# start_time = 1 for every job.
for start_time in 1:(T - duration)