Optimization using GLPK

I am trying to find an optimal number of employees to satisfy a scheduling problem, this is my code:

using JuMP
using GLPK

# Constants and parameters
Wjk = [19 16 22 22 22 22 22; 19 16 22 22 22 22 22; 14 11 16 16 16 16 16]  # Workload requirement
Ajk = 8  # Some constant value for Ajk
M = 395  # Total number of shifts

# Create a JuMP model
model = Model(GLPK.Optimizer)

# Decision variables for shift allocation
@variable(model, i >= N , Int)  # Number of employees
@variable(model, 1 <= j <= 3, Int)  # Number of shifts per day
@variable(model, 1 <= k <= 7, Int)  # Number of days per week
@variable(model, x[1:N, 1:3, 1:7] >= 0, Int)  # Decision variable for shift allocation


# Objective function
@objective(model, Min, sum(x[i, j, k] for i in 1:N, j in 1:3, k in 1:7))


# Objective function
@objective(model, Min, sum(x[i, j, k] for i in 1:N for j in 1:3 for k in 1:7))


# Constraints
for i in 1:N
    for k in 1:7
        @constraint(model, sum(x[i, j, k] for j in 1:3) <= 2)  # Constraints (8) - (10)
        @constraint(model, sum(x[i, j, k]) + sum(x[i, 2, mod(k + 1, 7) + 1] for j in 2:3) <= 2)  # Constraint (9)
        @constraint(model, sum(x[i, 3, k] for j in 1:3) + sum(x[i, 1, mod(k + 1, 7) + 1]) +
                    sum(x[i, 2, (mod(k + 1, 7)) + 1]) <= 2)  # Constraint (10)
    end
end


for j in 1:3
    for k in 1:7
        @constraint(model, sum(x[i, j, k] for i in 1:N) >= Wjk[k, j])  # Constraint (11)
        @constraint(model, sum(x[i, j, k] for i in 1:N) <= Wjk[k, j] + Ajk)  # Constraint (12)
    end
end

for i in 1:N
    @constraint(model, sum(x[i, j, ((i - 1) % 7) + 1] + x[i, j, i % 7 + 1] for j in 1:3) == 0)  # Constraint (13)
end

for i in 1:N
    for k in 1:5
        @constraint(model, sum(x[i, j, (i + k) % 7 + 1] for j in 1:3) >= 1)  # Constraint (14)
    end
end

@constraint(model, sum(x[i, j, k] for i in 1:N for j in 1:3 for k in 1:7) == M)

# Solve the optimization problem
optimize!(model)

# Display the results
println("Optimal Number of Employees: ", value(N))
println("Optimal Shift Allocation:")
for i in 1:value(N)
    for j in 1:3
        for k in 1:7
            if value(x[i, j, k]) > 0.5
                println("Employee $i, Shift $j on Day $k")
            end
        end
    end
end

When this runs I keep getting an error:

> ERROR: MethodError: Cannot `convert` an object of type VariableRef to an object of type Float64
> Closest candidates are:
>   convert(::Type{T}, ::Base.TwicePrecision) where T<:Number at twiceprecision.jl:250
>   convert(::Type{T}, ::AbstractChar) where T<:Number at char.jl:180
>   convert(::Type{T}, ::CartesianIndex{1}) where T<:Number at multidimensional.jl:136
>   ...
> Stacktrace:
>  [1] MathOptInterface.GreaterThan{Float64}(lower::VariableRef)
>    @ MathOptInterface C:\Users\doriiido\.julia\packages\MathOptInterface\IiXiU\src\sets.jl:171
>  [2] _moi_constrain_variable(moi_backend::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{GLPK.Optimize                                                                                                                                          er}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}, index::MathOptInterface.VariableIndex, info::VariableInfo{VariableRef, Float64, Float64, Float64}, #unused#::Type{Float64})
>    @ JuMP C:\Users\doriiido\.julia\packages\JuMP\ToPd2\src\variables.jl:1754
>  [3] _moi_add_variable(moi_backend::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{GLPK.Optimizer}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}, model::Model, v::ScalarVariable{VariableRef, Float64, Float64, Float64}, name::String)
>    @ JuMP C:\Users\doriiido\.julia\packages\JuMP\ToPd2\src\variables.jl:1737
>  [4] add_variable(model::Model, v::ScalarVariable{VariableRef, Float64, Float64, Float64}, name::String)
>    @ JuMP C:\Users\doriiido\.julia\packages\JuMP\ToPd2\src\variables.jl:1726
>  [5] macro expansion
>    @ C:\Users\doriiido\.julia\packages\JuMP\ToPd2\src\macros.jl:1213 [inlined]
>  [6] top-level scope
>    @ Untitled-2:13

I am not sure how to correct this

The error is coming from these lines shown above. You can safely comment them out since the set indices are not decision variables of the problem.

Thank you, but after doing so this is shown:

ERROR: MethodError: no method matching (::Colon)(::Int64, ::VariableRef)
Closest candidates are:
  (::Colon)(::T, ::Any, ::T) where T<:Real at range.jl:41
  (::Colon)(::A, ::Any, ::C) where {A<:Real, C<:Real} at range.jl:10
  (::Colon)(::T, ::Any, ::T) where T at range.jl:40
  ...
Stacktrace:
 [1] macro expansion
   @ C:\Users\doriiido\.julia\packages\JuMP\ToPd2\src\macros.jl:1213 [inlined]
 [2] top-level scope
   @ c:\Users\doriiido\Downloads\COMP6925\Q1-Scheduling_Optimization.jl:16

N cannot be a decision variable.

Here’s how I’d write your model. It’s currently infeasible because I don’t know if I got the constraints correct (what is sum(x) == M doing?), but it should point you in the right direction:

using JuMP
using HiGHS
Wjk = [
    19 16 22 22 22 22 22
    19 16 22 22 22 22 22
    14 11 16 16 16 16 16
]
Ajk = 8
M = 395
N = 100  # Cannot be a decision variable
model = Model(HiGHS.Optimizer)
@variable(model, x[1:N, 1:3, 1:7] >= 0, Int)
@objective(model, Min, sum(x))
@constraints(model, begin
    # Constraints (8) - (10)
    [i in 1:N, k in 1:7], sum(x[i,:,k]) <= 2
    # Constraint (9)
    [i in 1:N, k in 1:7], x[i,1,k] + sum(x[i,2,(k+1)%7+1] for j in 2:3) <= 2
    # Constraint (10)
    [i in 1:N, k in 1:7], x[i,3,k] + x[i,1,(k+1)%7+1] + x[i,2,(k+1)%7+1] <= 2
    # Constraint (11)
    [j in 1:3, k in 1:7], sum(x[:,j,k]) >= Wjk[j,k]
    # Constraint (12)
    [j in 1:3, k in 1:7], sum(x[:,j,k]) <= Wjk[j,k] + Ajk
    # Constraint (13)
    [i in 1:N], sum(x[i,j,(i-1)%7+1] + x[i,j,i%7+1] for j in 1:3) == 0
    # Constraint (14)
    [i in 1:N, k in 1:5], sum(x[i,j,(i+k)%7+1] for j in 1:3) >= 1
    sum(x) == M
end)
optimize!(model)
println("Number of Employees: ", N)
println("Optimal Shift Allocation:")
for i in 1:N, j in 1:3, k in 1:7
    if value(x[i, j, k]) > 0.5
        println("Employee $i, Shift $j on Day $k")
    end
end
1 Like

These are the mathematical equations taken from Hosein et al, 2019

sum(x) == M I added as a constraint to limit the max value of x to the sum of all shifts for all days, but I guess that the Wjk matrix would already achieve that

I don’t think you have the constraints quite right.

using JuMP
using HiGHS
Wjk = [
    19 16 22 22 22 22 22
    19 16 22 22 22 22 22
    14 11 16 16 16 16 16
]
J, K = size(Wjk)
Ajk = 8
N = 395  # TODO
model = Model(HiGHS.Optimizer)
@variable(model, x[1:N, 1:J, 1:K], Bin)
@objective(model, Min, sum(x))
@constraints(model, begin
    # Constraint (8)
    [i in 1:N, k in 1:K], x[i,1,k] + x[i,2,k] + x[i,3,k] <= 2
    # Constraint (9)
    [i in 1:N, k in 1:K], x[i,2,k] + x[i,3,k] + x[i,1,k%7+1] <= 2
    # Constraint (10)
    [i in 1:N, k in 1:K], x[i,3,k] + x[i,1,k%7+1] + x[i,2,k%7+1] <= 2
    # Constraint (11)
    [j in 1:J, k in 1:K], sum(x[:,j,k]) >= Wjk[j,k]
    # Constraint (12)
    [j in 1:J, k in 1:K], sum(x[:,j,k]) <= Wjk[j,k] + Ajk
    # Constraint (13)
    [i in 1:N], sum(x[i,j,(i-1)%7+1] + x[i,j,i%7+1] for j in 1:J) == 0
    # Constraint (14)
    [i in 1:N, k in 1:5], sum(x[i,j,(i+k)%7+1] for j in 1:J) >= 1
end)
optimize!(model)
println("Number of Employees: ", N)
println("Optimal Shift Allocation:")
for i in 1:N, j in 1:J, k in 1:K
    if value(x[i, j, k]) > 0.5
        println("Employee $i, Shift $j on Day $k")
    end
end

But this problem is infeasible because of (14). It doesn’t really make sense. Why is employee i related to day k?

My understanding of it is that this particular constraint 14 is that each employee must work at least one shift per day, excluding the 2 consecutive days off (constraint 13), with the modulus of i+k meaning that for employee i=1, the sum of shifts (j) on day k and k+i must be at least 1

Yes, but it doesn’t make sense to enforce that for all i if N is a constant.

I think you need a formulation like this, which let’s you choose the number of employees:

using JuMP
using HiGHS
Wjk = [
    19 16 22 22 22 22 22
    19 16 22 22 22 22 22
    14 11 16 16 16 16 16
]
J, K = size(Wjk)
Ajk = 8
N = 100  # Upper bound on number of expected employees
model = Model(HiGHS.Optimizer)
@variable(model, x[1:N, 1:J, 1:K], Bin)
@variable(model, y[1:N], Bin)
@objective(model, Min, sum(x))
@constraints(model, begin
    # y[i] is 1 if employee i is used
    [i in 1:N], sum(x[i, :, :]) <= J * K * y[i]
    # Ordering on y[i]. If we don't use y[i-1], then we can't use y[i]
    [i in 2:N], y[i-1] >= y[i]
    # Constraint (8)
    [i in 1:N, k in 1:K], x[i,1,k] + x[i,2,k] + x[i,3,k] <= 2
    # Constraint (9)
    [i in 1:N, k in 1:K], x[i,2,k] + x[i,3,k] + x[i,1,k%7+1] <= 2
    # Constraint (10)
    [i in 1:N, k in 1:K], x[i,3,k] + x[i,1,k%7+1] + x[i,2,k%7+1] <= 2
    # Constraint (11)
    [j in 1:J, k in 1:K], sum(x[:,j,k]) >= Wjk[j,k]
    # Constraint (12)
    [j in 1:J, k in 1:K], sum(x[:,j,k]) <= Wjk[j,k] + Ajk
    # Constraint (13)
    [i in 1:N], sum(x[i,j,(i-1)%7+1] + x[i,j,i%7+1] for j in 1:J) == 0
    # Constraint (14)
    [i in 1:N, k in 1:5], sum(x[i,j,(i+k)%7+1] for j in 1:J) >= y[i]
end)
optimize!(model)
println("Number of Employees: ", N)
println("Optimal Shift Allocation:")
for i in 1:N, j in 1:J, k in 1:K
    if value(x[i, j, k]) > 0.5
        println("Employee $i, Shift $j on Day $k")
    end
end
1 Like

I’m just noticing this so I’m not sure if it was addressed already, but the outcome I’m trying to achieve it to minimize the amount of employees required to fill the shift requirements stated in Wjk, but the model returns whatever the upper bound is set at.

Change your objective to @objective(model, Min, sum(y))