Wagner model for total tardiness minimization

Dear friends,

I hope this message finds you all well.

I am working on a Wagner family model with positional decision variables for the job shop scheduling problem, aiming to minimize total tardiness. For the toy model presented below, the global optimal solution is 6 time units. However, in my modeling, the tardiness is zero. I am unsure how to correctly compute the completion times and tardiness. Could anybody provide some assistance?

Best regards, Victor

Source code:

P = [2 5; 6 4; 2 3; 5 4; 6 7; 4 3; 6 5] # processing times
O = [1 2; 2 1; 2 1; 1 2; 2 1; 1 2; 1 2] # operations
d = [8 15 18 22 25 30 32] # due dates
V = 10000 # big-M
M = size(P,2)
N = size(P,1)

    # Build matrix R[j,l,k]
    R = zeros(Int64, N, N, M)

    for j in 1:N
        for i in 1:M
            maq = O[j,i]
            R[j,i,maq] = 1
        end
    end

    # Decision variables
    @variable(model, x[j in 1:N, h in 1:N, i in 1:M], Bin)
    @variable(model, C[h in 1:N, i in 1:M] >=0)
    @variable(model, T[j in 1:N] >=0)

    # Objective function
    @objective(model, Min, sum(T[j] for j in 1:N))

    # Constraints
    @constraint(model, [j=1:N, k=1:M], sum(x[j,h,k] for h=1:N) == 1)
    @constraint(model, [h=1:N, k=1:M], sum(x[j,h,k] for j=1:N) == 1)
    @constraint(model, [k=1:M], C[1,k] == sum(P[j,k]*x[j,1,k] for j=1:N))
    @constraint(model, [h=2:N,k=1:M],C[h,k] >= C[h-1,k] + sum(P[j,k]*x[j,h,k] for j=1:N))
    @constraint(model, [j=1:N,h=1:N,hh=1:N,l=2:M],
                        sum(R[j,l,k]*C[h,k] for k=1:M) >= sum(R[j,l-1,k]*C[hh,k] for k=1:M) +
                        sum(R[j,l,k]*P[j,k] for k=1:M) - V*(2 - sum(R[j,l,k]*x[j,h,k] for k=1:M) -
                        sum(R[j,l-1,k]*x[j,hh,k] for k=1:M)))
   
    @constraint(model, [j in 1:N, k in 1:M], T[j] >= C[j,k] - sum(d[j]*x[j,h,k] for h=1:N))

Hi @Victor_Neto, I’ve moved your question to the “Optimization (Mathematical)” section.

It’s not immediately obvious to me where the issue is. I assume that you have a typo somewhere in your constraints. If you have the mathematical formulation, I suggest you double check each constraint in details. If you’re still stuck, please post the mathematical formulation of what you are trying to model and we can help translate it into JuMP.

I rewrote your code slightly to make a reproducible example:

using JuMP, HiGHS
P = [2 5; 6 4; 2 3; 5 4; 6 7; 4 3; 6 5]
O = [1 2; 2 1; 2 1; 1 2; 2 1; 1 2; 1 2]
d = [8, 15, 18, 22, 25, 30, 32]
V = 10_000
N, M = size(P)
R = zeros(Int, N, N, M)
for j in 1:N, i in 1:M
    R[j,i,O[j,i]] = 1
end
model = Model(HiGHS.Optimizer)
@variables(model, begin
    x[j in 1:N, h in 1:N, i in 1:M], Bin
    C[h in 1:N, i in 1:M] >= 0
    T[j in 1:N] >= 0
end)
@objective(model, Min, sum(T))
@constraints(model, begin
    [j in 1:N, k in 1:M], sum(x[j,:,k]) == 1
    [h in 1:N, k in 1:M], sum(x[:,h,k]) == 1
    [k in 1:M], C[1,k] == P[:,k]' * x[:,1,k]
    [h in 2:N, k in 1:M], C[h,k] >= C[h-1,k] + P[:,k]' * x[:,h,k]
    [j in 1:N, h in 1:N, hh in 1:N, l in 2:M],
        R[j,l,:]' * C[h,:] >=
            R[j,l-1,:]' * C[hh,:] +
            R[j,l,:]'*P[j,:] -
            V * (2 - R[j,l,:]' * x[j,h,:] - R[j,l-1,:]' * x[j,hh,:])
    [j in 1:N, k in 1:M], T[j] >= C[j,k] - d[j] * sum(x[j,:,k])
end)
optimize!(model)

A useful tip might be to try and visualize the solution that gives 0 tardiness. If you can see what constraints are being violated (e.g., the jobs are overlapping) that should give you a clue that the corresponding constraint has a bug.

I took a closer look. Your tardiness constraint does’t look correct.

[j in 1:N, k in 1:M], T[j] >= C[j,k] - d[j] * sum(x[j,:,k])

but C[j, k] is the completion time of machine k in position j; it isn’t the completion time of job j on machine k. And by your first constraint, sum(x[j,:,k]) is == 1.

You probably want something like:

[j in 1:N, h in 1:N], T[j] >= C[h,O[j,end]] - d[j] - sum(P) * (1 - x[j,h,O[j,end]])

Dear Odow, thank you very much for your help.
The rewritten constraint for the problem worked correctly, thank you very much.
Victor