Columns and Constraints generation

Specific Domains Optimization (Mathematical)
, , ,
1

Hello, I’m looking to implement two-stage robust optimization using column and constraint generation. I’ve encountered an issue within while loop related to the optimize!() function when trying to implement the algorithm. Thank you for your help. Here’s the relevant code:

using JuMP
using Gurobi

# Constant creation
f = [400, 414, 326]
a = [18, 25, 20]
C = [22 33 24;
     33 23 30;
     20 25 27]
D = [206 + 40, 274 + 40, 220 + 40]
dl = [206, 274, 220]
du = [40, 40, 40]
k = 0  # Iterative counting

MP = Model(Gurobi.Optimizer)

# Variables
@variable(MP, x[1:3, 1:3] >= 0)
@variable(MP, y[1:length(f)], Bin)
@variable(MP, z[1:length(a)] >= 0)
@variable(MP, d[1:3] >= 0)
@variable(MP, eta >= 0)
@variable(MP, 0<=g[1:3]<=1)

# Objective function
@objective(MP, Min, sum(f[i] * y[i] for i in 1:length(f)) + sum(a[i] * z[i] for i in 1:length(a)) + eta)

# Constraints
@constraint(MP, MP_Cons_1[i in 1:3], z[i] <= 800 * y[i])
@constraint(MP, MP_Cons_2, sum(z[i] for i in 1:3) >= 772)

# Iteration constraints
@constraint(MP, MP_Cons_3[i in 1:3], sum(x[i, j] for j in 1:3) <= z[i])

@constraint(MP, MP_Cons_4[j in 1:3], sum(x[i, j] for i in 1:3) >= d[j])
@constraint(MP, MP_Cons_eta, eta >= sum(x[i, j]*C[i,j] for i in 1:3 for j in 1:3))
# Master-problem uncertainty constraints
@constraint(MP, MP_Cons_uncertainty_1[i in 1:3], d[i] == dl[i] + du[i] * g[i])
@constraint(MP, MP_Cons_uncertainty_2, sum(g[i] for i in 1:3) <= 1.8)
@constraint(MP, MP_Cons_uncertainty_3, sum(g[i] for i in 1:2) <= 1.2)

# Solve the problem
optimize!(MP)

# Get the lower bound
LB = objective_value(MP)

# Display the results if needed
println("Lower Bound: $LB")


SP = Model(optimizer_with_attributes(Gurobi.Optimizer, "OutputFlag" => 1))

# Variables
@variable(SP, x_sub[1:3, 1:3] >= 0)
@variable(SP, d_sub[1:3] >= 0)
@variable(SP, g_sub[1:3], Bin)
@variable(SP, pi[1:6] >= 0)
@variable(SP, v[1:6], Bin)
@variable(SP, w[1:3, 1:3], Bin)
M = 10000

# Objective function
@objective(SP, Max, sum(C[i,j] * x_sub[i, j] for i in 1:3, j in 1:3))
# Constraints
@constraint(SP, SP_Cons_1[i in 1:3], sum(x_sub[i, j] for j in 1:3) <= value.(z[i]))

@constraint(SP, SP_Cons_2[j in 1:3], sum(x_sub[i, j] for i in 1:3) >= d_sub[j])

@constraint(SP, SP_Cons_3[i in 1:3, j in 1:3], - pi[i] + pi[j + 3] <= C[i,j])
# Slack constraints part 1
@constraint(SP, SP_SLACK_CONS_1[i in 1:3], value.(z[i]) - sum(x_sub[i, j] for j in 1:3) <= M * (1 - v[i]))
@constraint(SP, SP_SLACK_CONS_2[j in 1:3], sum(x_sub[i, j] for i in 1:3) - d_sub[j] <= M * (1 - v[j + 3]))
@constraint(SP, SP_SLACK_CONS_3[i in 1:6], pi[i] <= M * v[i])

# Slack constraints part 2
@constraint(SP, SP_SLACK_CONS_4[i in 1:3, j in 1:3], C[i,j] + pi[i] - pi[j + 3] <= M * (1 - w[i, j]))
@constraint(SP, SP_SLACK_CONS_5[i in 1:3, j in 1:3], x_sub[i, j] <= M * w[i, j])

# Uncertainty
@constraint(SP, SP_Cons_uncertainty_1[i in 1:3], d_sub[i] == dl[i] + du[i] * g_sub[i])
@constraint(SP, SP_Cons_uncertainty_2, sum(g_sub[i] for i in 1:3) <= 1.8)
@constraint(SP, SP_Cons_uncertainty_3, sum(g_sub[i] for i in 1:2) <= 1.2)


# Solve the sub-problem
optimize!(SP)

# Get the sub-problem objective value
SP_objval = objective_value(SP)


# Calculate the upper bound
UB = LB - value.(eta) + SP_objval

println("Upper Bound: $UB")

while abs(UB - LB) > 1e-5
    k = k + 1
    # Master-problem
    @variable(MP, x_new[1:3, 1:3] >= 0)
    @constraint(MP, MP_Cons_3_new[i in 1:3], sum(x_new[i, j] for j in 1:3) <= value(z[i]))
    @constraint(MP, MP_Cons_4_new[j in 1:3], sum(x_new[i, j] for i in 1:3) >= value(d_sub[j]))
    @constraint(MP, MP_Cons_eta, eta >= sum(x_new[i, j] * C[i,j] for i in 1:3, j in 1:3))
  optimize!(MP)
    LB = max(LB, objective_value(MP))
    println("Lower Bound: $LB")
    # Sub-problem update
    # Remove old constraints related to z
    for i in 1:3
        delete(MP_Cons_3[i])
        delete(SP_SLACK_CONS_1[i])
    end
    # Add new constraints related to z
    for i in 1:3
        @constraint(MP, MP_Cons_3_new[i], sum(x_sub[i, j] for j in 1:3) <= value(z[i]))
        @constraint(SP, SP_Cons_1[i], sum(x_sub[i, j] for j in 1:3) <= value(z[i]))
        @constraint(SP, SP_SLACK_CONS_1[i], value(z[i]) - sum(x_sub[i, j] for j in 1:3) <= M * (1 - value(v[i])))
    end
   optimize!(SP)
    SP_objval = objective_value(SP)
    UB = LB - value(eta) + SP_objval
    println("Upper Bound: $UB")
end
2

Hi @Salma, welcome to the forum.

The key part of the output is this warning:

┌ Warning: The model has been modified since the last call to `optimize!` (or `optimize!` has not been called yet). If you are iteratively querying solution information and modifying a model, query all the results first, then modify the model.
└ @ JuMP ~/.julia/dev/JuMP/src/optimizer_interface.jl:694

See this section of the docs: Solutions · JuMP

I don’t think your algorithm is quite right. But here’s how I would write it for now. I left a few TODO which need fixing before it’ll work.

using JuMP
using Gurobi

f = [400, 414, 326]
a = [18, 25, 20]
C = [22 33 24; 33 23 30; 20 25 27]
D = [206 + 40, 274 + 40, 220 + 40]
dl = [206, 274, 220]
du = [40, 40, 40]
k = 0

MP = Model(Gurobi.Optimizer)
@variables(MP, begin
    x[1:3, 1:3] >= 0
    y[1:length(f)], Bin
    z[1:length(a)] >= 0
    d[1:3] >= 0
    eta >= 0
    0 <= g[1:3] <= 1
end)
@objective(MP, Min, f' * y + a' * z + eta)
@constraints(MP, begin
    [i in 1:3], z[i] <= 800 * y[i]
    sum(z) >= 772
    MP_Cons_3[i in 1:3], sum(x[i, :]) <= z[i]
    [j in 1:3], sum(x[:, j]) >= d[j]
    eta >= sum(C .* x)
    [i in 1:3], d[i] == dl[i] + du[i] * g[i]
    sum(g[i] for i in 1:3) <= 1.8
    sum(g[i] for i in 1:2) <= 1.2
end)
optimize!(MP)
LB = objective_value(MP)
println("Lower Bound: $LB")

z_star = value.(z)

M = 10_000
SP = Model(Gurobi.Optimizer)
set_silent(SP)
@variables(SP, begin
    x_sub[1:3, 1:3] >= 0
    d_sub[1:3] >= 0
    g_sub[1:3], Bin
    pi[1:6] >= 0
    v[1:6], Bin
    w[1:3, 1:3], Bin
end)
@objective(SP, Max, sum(C .* x_sub))
@constraints(SP, begin
    [i in 1:3], sum(x_sub[i, :]) <= z_star[i]
    [j in 1:3], sum(x_sub[:, j]) >= d_sub[j]
    [i in 1:3, j in 1:3], -pi[i] + pi[j + 3] <= C[i,j]
    SP_SLACK_CONS_1[i in 1:3], z_star[i] - sum(x_sub[i, :]) <= M * (1 - v[i])
    [j in 1:3], sum(x_sub[:, j]) - d_sub[j] <= M * (1 - v[j + 3])
    [i in 1:6], pi[i] <= M * v[i]
    [i in 1:3, j in 1:3], C[i,j] + pi[i] - pi[j + 3] <= M * (1 - w[i, j])
    [i in 1:3, j in 1:3], x_sub[i, j] <= M * w[i, j]
    [i in 1:3], d_sub[i] == dl[i] + du[i] * g_sub[i]
    sum(g_sub[i] for i in 1:3) <= 1.8
    sum(g_sub[i] for i in 1:2) <= 1.2
end)
UB = Inf
while abs(UB - LB) > 1e-5
    optimize!(SP)
    SP_objval = objective_value(SP)
    UB = LB - value(eta) + SP_objval
    println("Upper Bound: $UB")
    # New: query values before modification
    d_sub_star = value.(SP)
    k = k + 1
    x_new = @variable(MP, [1:3, 1:3] >= 0)
    # TODO: do you mean z[i] here, or the previous value of `z-star`?
    @constraint(MP, [i in 1:3], sum(x_new[i, :]) <= value(z[i]))
    @constraint(MP, [j in 1:3], sum(x_new[:, j]) >= d_sub_star[j])
    @constraint(MP, eta >= sum(C .* x_new))
    optimize!(MP)
    LB = max(LB, objective_value(MP))
    println("Lower Bound: $LB")
    # New: query values before modification
    z_star = value.(z)
    v_star = value.(v)
    # What's going on here?
    delete.(MP_Cons_3)
    unregister(MP, :MP_Cons_3)
    delete.(SP_SLACK_CONS_1)
    unregister(SP, :SP_SLACK_CONS_1)
    # TODO: which variables do you mean? `x_sub` is from SP
    @constraint(MP, MP_Cons_3[i in 1:3], sum(x_sub[i, j] for j in 1:3) <= value(z[i]))
    @constraints(SP, begin
        [i in 1:3], sum(x_sub[i, j] for j in 1:3) <= z_star[i]
        SP_SLACK_CONS_1[i in 1:3], z_star[i] - sum(x_sub[i, j] for j in 1:3) <= M * (1 - v_star[i])
    end)
end
2 Likes
3

Dear @odow thank you for your help.

1 Like