Performance issues with Multi-Objective Optimization

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:

grafik

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:

I guess a duplicate of Symmetry breaking: scheduling problems with job duration dependent on ressource allocation?

This was originally under the general section; I moved it to the Optimization (Mathematical) section. If you have JuMP or optimization related questions, this is the place to post them. (It’s easy for questions to get missed in the general section.)