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
4

Hi, Dowson. I’m also studying this C&CG algorithm right now (and the same case study as the one Salma provided). Do you have some ideas where can I find a related package or code snippet about this algorithm, in julia? I wrote one myself, and I’d like to make some comparisons, because I’m not 100% sure that my program is correct.

5

Did you accomplish this small case? What’s your final result?
I get an optimal cost at 30536, which is different from that in the paper (33680).
I’m now inspecting, is there anything wrong.

6

Hi @WalterMadelim, welcome to the forum.

I don’t know if there are any off-the-shelf implementations of row and column generation algorithms.

The closest thing might be GitHub - atoptima/Coluna.jl: Branch-and-Price-and-Cut in Julia

7

Perhaps start a new post with a reproducible example of your code.

8

small_case_in_BoZeng2012
C&CG algorithm

Because new users is limited to comment (<= 3 replies). I update as follows:

I know how the data in Table 1 comes now.
In the initialization phase, I’ve just used a gurobi-default initialization scenario (which happens to be the most advantageous case).
If just to change the initialization scenario, we can recover precisely those 4 numbers (w.r.t. 2 iterations) in Table 1.

I might know where the problem is. This algorithm is a bit involved.

1 Like
9

updated

10

The possible explanation: the subproblem is a Max-Min problem originally. For a given 1st-stage decision y, it may happens that there is a scenario u such that the resulting value is ∞. In this case, we need to incorporate this significant u into the master problem at next iteration.

But if we cast the original Max-Min problem to a single-level Max problem with the KKT condition, we might miss these significant u’s which might have rendered the 2nd-decision phase infeasible.

11

Is there an actionable question about JuMP here?

If you’re trying to reproduce results from a paper, then you likely need to code the algorithm exactly the same. If you get different results, then either you have a mistake in your code, your reformulation is not applicable, or there is a typo in the original paper.

and the same case study as the one Salma provided
I get an optimal cost at 30536, which is different from that in the paper (33680).
I know how the data in Table 1 comes now.

I don’t know what case study, paper, or table you are referring to :smile:

12

Hey. Here’s the code that works : :slight_smile:

#MP:

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)
set_silent(MP)
@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
[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”)

#Recourse:
z_star = value.(z)
M = 10000
SP = Model(Gurobi.Optimizer)
set_silent(SP)
@variables(SP, begin
x_sub[1:3, 1:3] >= 0;
d_sub[1:3] >= 0;
0<=g_sub[1:3]<=1;
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
SP_Cons_1[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)
optimize!(SP)
println(“objective:”, objective_value(SP))

#CCG

UB = Inf
while abs(UB - LB) > 1e-5
optimize!(SP)
SP_objval = objective_value(SP)
println(“SP_objval: $SP_objval”)
UB = LB - value(eta) + SP_objval
println(“Upper Bound: $UB”)
# New: query values before modification
d_sub_star = value.(d_sub)
println(“d_sub=”,value.(d_sub_star))
unregister(MP,:x_new)
@variable(MP, x_new[1:3, 1:3] >= 0)
k = k + 1
@constraint(MP, [i in 1:3], sum(x_new[i, :]) <= 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)
unregister(SP,:x_new)
delete.(SP, SP_Cons_1)
unregister(SP,:SP_Cons_1)
delete.(SP,SP_SLACK_CONS_1)
unregister(SP,:SP_SLACK_CONS_1)
@constraints(SP, begin
SP_Cons_1[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[i])
end)
end
println(“Iteration finished! We found the optimal solution!”)
println(“Final Objective:{0}”,UB)

1 Like
13

I have implemeted the CCG using strong duality. I can share it too.

The paper is “Solving two-stage robust optimization problems using a
column-and-constraint generation method” Redirecting

14

Thank you both. Don’t worry, I’ve figured it out, and have revised the code in my repository as the link I’ve mentioned.

1 Like