# Columns and Constraints generation

**URL:** <https://discourse.julialang.org/t/columns-and-constraints-generation/103813>\
**Category:** Optimization (Mathematical)\
**Tags:** question, jump, optimization, gurobi\
**Created:** [September 13, 2023, 1:45pm UTC](https://discourse.julialang.org/t/columns-and-constraints-generation/103813 "2023-09-13T13:45:14Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![Salma](https://avatars.discourse-cdn.com/v4/letter/s/fbc32d/32.png) [@Salma](https://discourse.julialang.org/u/Salma)\
**Post date:** [September 13, 2023, 1:45pm UTC](https://discourse.julialang.org/t/columns-and-constraints-generation/103813/1 "2023-09-13T13:45:14Z")

</div>

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:

```julia

```

```julia
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

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [September 14, 2023, 1:57am UTC](https://discourse.julialang.org/t/columns-and-constraints-generation/103813/2 "2023-09-14T01:57:09Z")

</div>

Hi @Salma, welcome to the forum.

The key part of the output is this warning:

```julia
┌ 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](https://jump.dev/JuMP.jl/stable/manual/solutions/#OptimizeNotCalled-errors)

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.

```julia
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

```

---

<div class="post-metadata">

**Author:** ![Salma](https://avatars.discourse-cdn.com/v4/letter/s/fbc32d/32.png) [@Salma](https://discourse.julialang.org/u/Salma)\
**Post date:** [September 17, 2023, 12:15pm UTC](https://discourse.julialang.org/t/columns-and-constraints-generation/103813/3 "2023-09-17T12:15:25Z")

</div>

Dear @odow thank you for your help.

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [February 3, 2024, 2:51am UTC](https://discourse.julialang.org/t/columns-and-constraints-generation/103813/4 "2024-02-03T02:51:19Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [February 3, 2024, 3:00am UTC](https://discourse.julialang.org/t/columns-and-constraints-generation/103813/5 "2024-02-03T03:00:20Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [February 3, 2024, 3:00am UTC](https://discourse.julialang.org/t/columns-and-constraints-generation/103813/6 "2024-02-03T03:00:53Z")

</div>

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](https://github.com/atoptima/Coluna.jl)

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [February 3, 2024, 3:01am UTC](https://discourse.julialang.org/t/columns-and-constraints-generation/103813/7 "2024-02-03T03:01:39Z")

</div>

> [@WalterMadelim](#):
>
> I get an optimal cost at 30536, which is different from that in the paper (33680).  
> I’m now inspecting, is there anything wrong.

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

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [February 3, 2024, 3:04am UTC](https://discourse.julialang.org/t/columns-and-constraints-generation/103813/8 "2024-02-03T03:04:07Z")

</div>

[small\_case\_in\_BoZeng2012](https://github.com/WalterMadelim/p8f.jl/blob/main/src/CCG2012.jl)  
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.

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [February 3, 2024, 3:06am UTC](https://discourse.julialang.org/t/columns-and-constraints-generation/103813/9 "2024-02-03T03:06:46Z")

</div>

updated

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [February 3, 2024, 3:13am UTC](https://discourse.julialang.org/t/columns-and-constraints-generation/103813/10 "2024-02-03T03:13:38Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [February 3, 2024, 9:15pm UTC](https://discourse.julialang.org/t/columns-and-constraints-generation/103813/11 "2024-02-03T21:15:24Z")

</div>

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 😄

---

<div class="post-metadata">

**Author:** ![Salma](https://avatars.discourse-cdn.com/v4/letter/s/fbc32d/32.png) [@Salma](https://discourse.julialang.org/u/Salma)\
**Post date:** [February 3, 2024, 9:24pm UTC](https://discourse.julialang.org/t/columns-and-constraints-generation/103813/12 "2024-02-03T21:24:25Z")

</div>

Hey. Here’s the code that works : 🙂

#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)

---

<div class="post-metadata">

**Author:** ![Salma](https://avatars.discourse-cdn.com/v4/letter/s/fbc32d/32.png) [@Salma](https://discourse.julialang.org/u/Salma)\
**Post date:** [February 3, 2024, 9:36pm UTC](https://discourse.julialang.org/t/columns-and-constraints-generation/103813/13 "2024-02-03T21:36:49Z")

</div>

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](https://doi.org/10.1016/j.orl.2013.05.003)

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [February 9, 2024, 1:16pm UTC](https://discourse.julialang.org/t/columns-and-constraints-generation/103813/14 "2024-02-09T13:16:40Z")

</div>

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.