# Gurobi fails to invoke the primal simplex method when doing column generation?

**URL:** <https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675>\
**Category:** Optimization (Mathematical)\
**Tags:** jump, gurobi\
**Created:** [November 5, 2025, 11:01am UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675 "2025-11-05T11:01:15Z")\
**Posts on this page:** 19\
**Page:** 1

<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:** [November 5, 2025, 11:01am UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/1 "2025-11-05T11:01:15Z")

</div>

I have a refreshing discovery today, based on my acquaintance [Column generation · JuMP](https://jump.dev/JuMP.jl/dev/tutorials/algorithms/cutting_stock_column_generation/#Column-generation)

wherein, the only thing I changed was the solver (I used Gurobi instead of HiGHS).

This is the master LP’s re-solving logging I acquired when I set `Method` to `0`

```julia-auto
LP warm-start: use basis

Iteration Objective Primal Inf. Dual Inf. Time
       0 3.6848477e+02 0.000000e+00 5.526316e-01 0s
       2 3.5932687e+02 0.000000e+00 0.000000e+00 0s

```

By comparison, this is the counterpart when `Method === 1` (the default setting, just the same as the doc).

```julia-auto
LP warm-start: use basis

Iteration Objective Primal Inf. Dual Inf. Time
       0 -3.3333333e+29 1.333333e+30 3.333333e-01 0s
       7 2.9900000e+02 0.000000e+00 0.000000e+00 0s

```

Okay I believe that everyone would find the former a lot more natural at least in terms of numerics—the primal simplex method which you need to switch to _manually_—unfortunately this is not the default foolproof setting taught by JuMP.

So the conclusion is:

- JuMP’s doc is adopting a default setting that doesn’t attain its expected performance in practice.

**PS** The motivation that yielded my idea today, was that I was wondering why the dual simplex method is the default. My guess is: modern solvers are row generation based, so the spirit of column generation was opposed to that. And lucky enough—I discovered this outcome—affirming my surmise.

> **Code**
>
> ```julia-auto
> using JuMP
> import DataFrames
> import Gurobi
> import SparseArrays
> const GRB_ENV = Gurobi.Env();
> struct Piece
> w::Float64
> d::Int
> end
> struct Data
> pieces::Vector{Piece}
> W::Float64
> end
> function get_data()
> data = [
> 75.0 38
> 70.0 41
> 68.4 34
> 65.5 23
> 59.6 18
> 53.8 33
> 53.0 36
> 51.0 41
> 50.2 35
> 32.2 37
> 30.8 44
> 29.8 49
> 20.1 37
> 16.2 36
> 14.5 42
> 11.0 33
> 8.6 47
> 8.2 35
> 6.6 49
> 5.1 42
> ]
> return Data([Piece(data[i, 1], data[i, 2]) for i in axes(data, 1)], 100.0)
> end
> function solve_pricing(data::Data, π::Vector{Float64})
> I = length(π)
> model = Model(() -> Gurobi.Optimizer(GRB_ENV))
> set_silent(model)
> @variable(model, y[1:I] >= 0, Int)
> @constraint(model, sum(data.pieces[i].w * y[i] for i in 1:I) <= data.W)
> @objective(model, Max, sum(π[i] * y[i] for i in 1:I))
> optimize!(model)
> assert_is_solved_and_feasible(model)
> number_of_rolls_saved = objective_value(model)
> if number_of_rolls_saved > 1 + 1e-8
> # Benefit of pattern is more than the cost of a new roll plus some
> # tolerance
> return SparseArrays.sparse(round.(Int, value.(y)))
> end
> return nothing
> end
> data = get_data();
> I = length(data.pieces)
> J = 1_000 # Some large number
> 
> patterns = map(1:I) do i
> n_pieces = floor(Int, data.W / data.pieces[i].w)
> return SparseArrays.sparsevec([i], [n_pieces], I)
> end
> 
> model = Model(() -> Gurobi.Optimizer(GRB_ENV))
> # JuMP.set_attribute(model, "Method", 1)
> @variable(model, x[1:length(patterns)] >= 0)
> @objective(model, Min, sum(x))
> @constraint(model, demand[i in 1:I], patterns[i]' * x >= data.pieces[i].d)
> optimize!(model)
> assert_is_solved_and_feasible(model)
> solution_summary(model)
> 
> while true
> # Solve the linear relaxation
> optimize!(model)
> assert_is_solved_and_feasible(model; dual = true)
> # Obtain a new dual vector
> π = dual.(demand)
> # Solve the pricing problem
> new_pattern = solve_pricing(data, π)
> # Stop iterating if there is no new pattern
> if new_pattern === nothing
> @info "No new patterns, terminating the algorithm."
> break
> end
> push!(patterns, new_pattern)
> # Create a new column
> push!(x, @variable(model, lower_bound = 0))
> # Update the objective coefficient of the new column
> set_objective_coefficient(model, x[end], 1.0)
> # Update the non-zeros in the coefficient matrix
> for (i, count) in zip(SparseArrays.findnz(new_pattern)...)
> set_normalized_coefficient(demand[i], x[end], count)
> end
> println("Found new pattern. Total patterns = $(length(patterns))")
> end
> 
> ```

**PS2** I further did a test, the results indicate that the default setting (dual simplex method) is indeed the best go-to algorithm for a CTPLN master LP, where I used the (Barrier+noCrossover) as a baseline.

```julia-auto
# Barrier+noCrossover
Row │J cg_time cg_rgap decen_ub decen_rgap decen_time Kverm Kvermu KverM msttimem msttimemu msttimeM 
     │Int64 Int64 Float64 Float64 Float64 Float64 Int64 Float64 Int64 Float64 Float64 Float64  
─────┼────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │16192 300 0.0528173 2.05448e5 0.0499749 317.78 1 2.193 3 0.1394 1.41158 3.92962 
   2 │16192 900 0.0173547 2.01733e5 0.016572 961.554 1 3.21189 5 0.1394 2.48585 7.85246 
   3 │16192 2700 0.00391167 2.00416e5 0.00391474 2817.17 1 4.58745 7 0.1394 4.71001 21.289   
# Dual simplex (default)
 Row │J cg_time cg_rgap decen_ub decen_rgap decen_time Kverm Kvermu KverM msttimem msttimemu msttimeM 
     │Int64 Int64 Float64 Float64 Float64 Float64 Int64 Float64 Int64 Float64 Float64 Float64  
─────┼───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │16192 300 0.0350881 204397.0 0.0338806 332.166 1 2.48141 3 0.00010705 1.14491 7.80573 
   2 │16192 900 0.0130845 2.01481e5 0.0126484 988.064 1 3.48413 5 0.00010705 2.26494 16.9965  
   3 │16192 2700 0.00283829 2.00259e5 0.00277348 2802.26 1 4.87148 8 0.00010705 4.66875 81.4494  

```

---

<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:** [November 5, 2025, 7:02pm UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/2 "2025-11-05T19:02:26Z")

</div>

> [@WalterMadelim](#):
>
> JuMP’s doc is adopting a default setting that doesn’t attain its expected performance in practice.

JuMP uses the default settings chosen by Gurobi.

---

<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:** [November 6, 2025, 12:15am UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/3 "2025-11-06T00:15:14Z")

</div>

This makes sense for a one-shot solve. But not very in this case, where the master LP is built-and-solved incrementally. I’m curious if this behavior (cannot invoke primal simplex) is Gurobi-exclusive. Do you have some interests to take a look at other common solvers? e.g. Highs, Cplex etc. (Or anyone has a different solver has an interest to take a look?)

---

<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:** [November 6, 2025, 2:14am UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/4 "2025-11-06T02:14:50Z")

</div>

What’s your question with the logging?

When you set `Method = 0`, you get primal simplex.

When you set `Method = 1`, you get dual simplex.

When you set `Method = -1` or you don’t set, Gurobi chooses.

Here’s the doc: [Parameter Reference - Gurobi Optimizer Reference Manual](https://docs.gurobi.com/projects/optimizer/en/current/reference/parameters.html#method)

---

<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:** [November 6, 2025, 2:35am UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/5 "2025-11-06T02:35:57Z")

</div>

There are 3 known decomposition techniques in our field:

1. dual decomposition
2. DW decomposition
3. Benders decomposition

In the DW decomposition, _column generation_ is used.  
In the Benders (or dual) decomposition, _row generation_ (cut generation) is used.

They actually has some sort of dual relation.

The current conclusion is: the dual simplex method is suitable for the row generation case, whereas the primal simplex method is suitable for the column generation case. (You can, e.g., ask some chatbot to confim this.)

In a previous topic I’ve mentioned that I favor cut generation (cutting plane method) over column generation. Because cutting plane method is more intuitive and more widely studied. And modern solvers like Gurobi supports cut generation better. This is why the dual simplex method is typically the default. (And thus the results in my #1 post).

The standard method to do column generation is just specified in the JuMP’s column generation doc, which by itself is correct. But I found that Gurobi is not invoking the fitting `Method`. I think this is a mismatch. If you are aware of any stuff from Gurobi (I remember there was one bro here), I think the current situation is worth reporting.

---

<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:** [November 6, 2025, 2:40am UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/6 "2025-11-06T02:40:50Z")

</div>

> The standard method to do column generation is just specified in the JuMP’s column generation doc, which by itself is correct

What do you mean by this? We don’t specify anything.

> But I found that Gurobi is not invoking the fitting `Method`

What evidence do you have for this? Looks okay to me. Note that `Method = 1` is not the default. The default is `Method = -1`

---

<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:** [November 6, 2025, 2:52am UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/7 "2025-11-06T02:52:17Z")

</div>

> [@WalterMadelim](#):
>
> ```julia-auto
> Iteration Objective Primal Inf. Dual Inf. Time
> 0 3.6848477e+02 0.000000e+00 5.526316e-01 0s
> 2 3.5932687e+02 0.000000e+00 0.000000e+00 0s
> 
> ```

The primal simplex method features maintaining a primal feasible point, so the “Primal Inf.” column is typically a zero column, whereas the “Dual Inf.” decrease to zero gradually. (The above logging belongs to the primal simplex’s logging).

Whereas the dual simplex method features maintaining a dual feasible point, so the “Dual Inf.” column is typically a zero column, whereas the “Primal Inf.” decrases to zero gradually.

* * *

> [@odow](#):
>
> What do you mean by this? We don’t specify anything.

You call the standard column generation API, i.e. `set_normalized_coefficint` (although I’ve never written my code with that). That API pertains to the Gurobi’s standard API of modifying a coefficient of a constraint. So I assume it is strongly connected to the usage of column generation.  
Because when we do column generation, we need to modify coefficient of constraints and objective (IIRC), and add variables (i.e. new columns). (By comparison, the cut generation is a lot more simpler in conception, where we only need to add a constraint).

> [@odow](#):
>
> Note that `Method = 1` is not the default. The default is `Method = -1`

You can observe the simplex’s logging, in most cases the simplex logging is identified as the dual simplex’s logging (you can, e.g. ask a chatbot, let it tell you which method corresponds to a particular simplex logging). The dual simplex method is the de-facto default algorithm unless the scale of the LP becomes large.

---

<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:** [November 6, 2025, 3:00am UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/8 "2025-11-06T03:00:56Z")

</div>

What’s the question? I’m a bit confused. So far my answers are “yes, and?”

---

<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:** [November 6, 2025, 3:05am UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/9 "2025-11-06T03:05:41Z")

</div>

Following the JuMP’s doc about column generation, the user’s intention is to perform a column generation process on the master LP.

The expected behavior for Gurobi should be to invoke the primal simplex method (Method = 0).  
But Gurobi cannot identify the user’s intention and using Method = 1 instead. This is a mismatch and affects due performance.

The Gurobi optimizer works well if user do cut generation only,. (e.g. I’m with this style). But the current behavior is not friendly to people who wants to do column generation, unless we explicitly let them know there is an “Method = 0” manual option.

---

<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:** [November 6, 2025, 3:49am UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/10 "2025-11-06T03:49:57Z")

</div>

> [@WalterMadelim](#):
>
> The expected behavior for Gurobi should be to invoke the primal simplex method

This is just not the case. Gurobi is a black-box. It may decide to do anything, so long as it solves the problem.

Regardless, this is not a problem with JuMP, Julia, or even Gurobi.jl. It’s an internal algorithmic choice made by Gurobi.

---

<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:** [November 6, 2025, 5:13am UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/11 "2025-11-06T05:13:39Z")

</div>

😊I’m not averse to it. After all, CG is not my go-to style.

---

<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:** [November 8, 2025, 3:16pm UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/12 "2025-11-08T15:16:19Z")

</div>

😊I’m convinced that @odow 's idea is credible. You always make me rethink so I can make more sound conclusions.

The crux here is I find that it appears the dual simplex method in Gurobi is (probably) way more efficient that its primal simplex method counterpart (may not be true).  
**Edit** : seems to be non-deterministic.

---

<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:** [November 9, 2025, 12:09am UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/13 "2025-11-09T00:09:26Z")

</div>

I do manage to generate a random test case wherein:

- By manually switching to the _primal simplex_ method, we get a consistent **2.5** x performance boost.
- By manually switching to the _barrier+NoCrossover_ method, we get a consistent **21** x performance boost.

**Test:** create a julia text file `s.jl`

#### s.jl

```julia-auto
const (M,) = map(a -> parse(Int, a), ARGS)
import JuMP, Gurobi
import Statistics
import Random; Random.seed!(hash(1))
const N = 499;
function get_warm_model()
    m = JuMP.Model(Gurobi.Optimizer);
    JuMP.set_attribute(m, "Threads", 1)
    M === 2 && JuMP.set_attribute(m, "Crossover", 0)
    M !== 1 && JuMP.set_attribute(m, "Method", M) 
    JuMP.@variable(m, rand(-7:-3) <= y[1:N] <= rand(3:7));
    JuMP.@variable(m, rand(-9:-5) <= x[1:N] <= rand(5:9));
    JuMP.@constraint(m, [1:7*N], rand(-9:0.017:9, N)'x + rand(-9:0.017:9, N)'y <= rand(0:9))
    JuMP.@objective(m, Min, rand(-9:0.017:9, N)'x + rand(-9:0.017:9, N)'y);
    JuMP.optimize!(m) # do an initial solve to warm up
    m
end
const m = get_warm_model();
const R = 20;
const t = Vector{Float64}(undef, R);
const v = Vector{Float64}(undef, R);
function test()
    for r = 1:R
        JuMP.set_objective_coefficient(m, m[:x], rand(-9:0.017:9, N))
        JuMP.set_objective_coefficient(m, m[:y], rand(-9:0.017:9, N))
        JuMP.optimize!(m)
        JuMP.termination_status(m) === JuMP.OPTIMAL || error()
        t[r] = JuMP.solve_time(m)
        v[r] = JuMP.objective_value(m)
    end
    @show t
    @show v
end
test()

```

Then create 3 different new shells, type `julia --threads=1,1 s.jl 0`, `julia --threads=1,1 s.jl 1`, `julia --threads=1,1 s.jl 2`.

#### Result

(each row corresponds to a re-solve with updated objective coefficients)

```julia-auto
julia> df = DataFrames.DataFrame((;tDefault, tPrimSimp, tBarr0Xov, objDefault, objPrimSimp, objBarr0Xov))
20×6 DataFrame
 Row │ tDefault tPrimSimp tBarr0Xov objDefault objPrimSimp objBarr0Xov 
     │ Float64 Float64 Float64 Float64 Float64 Float64     
─────┼──────────────────────────────────────────────────────────────────────
   1 │ 61.4885 24.4547 2.8746 -119.818 -119.818 -119.818
   2 │ 60.4172 24.7372 3.00952 -125.017 -125.017 -125.017
   3 │ 62.5353 25.6472 3.09385 -126.465 -126.465 -126.465
   4 │ 64.9553 25.3034 2.97567 -123.466 -123.466 -123.466
   5 │ 66.5784 25.3754 3.01711 -121.95 -121.95 -121.95
   6 │ 66.3094 25.1885 2.89824 -127.11 -127.11 -127.11
   7 │ 64.453 25.6088 2.87135 -125.203 -125.203 -125.203
   8 │ 60.3566 24.4784 3.02875 -114.185 -114.185 -114.185
   9 │ 60.6781 25.9133 2.99297 -118.85 -118.85 -118.85
  10 │ 61.2646 24.1942 2.97787 -134.872 -134.872 -134.872
  11 │ 62.5482 26.3902 2.99578 -125.065 -125.065 -125.065
  12 │ 64.1979 22.1599 2.92946 -126.156 -126.156 -126.156
  13 │ 63.3141 25.0833 3.17004 -124.922 -124.922 -124.922
  14 │ 62.7033 25.6 2.89424 -125.382 -125.382 -125.382
  15 │ 64.822 25.0511 3.01604 -120.071 -120.071 -120.071
  16 │ 61.9234 25.1266 3.03656 -128.452 -128.452 -128.452
  17 │ 63.1282 24.291 2.87391 -104.031 -104.031 -104.031
  18 │ 60.2861 24.7982 3.12961 -125.406 -125.406 -125.406
  19 │ 58.2621 23.2958 2.89144 -130.608 -130.608 -130.608
  20 │ 62.6119 24.6109 2.85591 -120.871 -120.871 -120.871

julia> DataFrames.describe(df)
6×7 DataFrame
 Row │ variable mean min median max nmissing eltype   
     │ Symbol Float64 Float64 Float64 Float64 Int64 DataType 
─────┼─────────────────────────────────────────────────────────────────────────────────
   1 │ tDefault 62.6417 58.2621 62.5801 66.5784 0 Float64
   2 │ tPrimSimp 24.8654 22.1599 25.0672 26.3902 0 Float64
   3 │ tBarr0Xov 2.97665 2.85591 2.98542 3.17004 0 Float64
   4 │ objDefault -123.395 -134.872 -125.041 -104.031 0 Float64
   5 │ objPrimSimp -123.395 -134.872 -125.041 -104.031 0 Float64
   6 │ objBarr0Xov -123.395 -134.872 -125.041 -104.031 0 Float64

```

One more interesting detail is that:  
The auto’s Logging contains a Warning “Warning: Markowitz tolerance tightened to 0.03125”,

---

<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:** [November 9, 2025, 2:51am UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/14 "2025-11-09T02:51:16Z")

</div>

> [@odow](#):
>
> Regardless, this is not a problem with JuMP, Julia, or even [Gurobi.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/Gurobi). It’s an internal algorithmic choice made by Gurobi.

You’re right. But we could add a page in JuMP’s doc to display the findings in my #13 post.  
Isn’t it thrilling that opting to another method makes our code run **21x faster**?

---

<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 23, 2026, 7:06am UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/15 "2026-02-23T07:06:59Z")

</div>

> [@WalterMadelim](#):
>
> You call the standard column generation API, i.e. `set_normalized_coefficint` (although I’ve never written my code with that). That API pertains to the Gurobi’s standard API of modifying a coefficient of a constraint. So I assume it is strongly connected to the usage of column generation.

It appears that I was wrong—the real native API for column generation is presumably [Model Creation and Modification - Gurobi Optimizer Reference Manual](https://docs.gurobi.com/projects/optimizer/en/current/reference/c/model.html#c.GRBaddvar)

```c
int GRBaddvar(GRBmodel *model, int numnz, int *vind, double *vval, double obj, double lb, double ub, char vtype, const char *varname)

```

, which is a C-API that julia users almost surely shouldn’t call.

Therefore it would be intrinsically awkward to write a column generation algorithm in JuMP (with the Gurobi solver), I assume.

---

<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 24, 2026, 1:47am UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/16 "2026-02-24T01:47:59Z")

</div>

The JuMP developers have discussed this in the past.

We have explicitly decided not to support the native APIs for column generation. Our reasoning is:

1. that some solvers provide an explicit API for adding coefficients when the variable is added
2. that doing so is probably faster than adding a variable and then modifying constraint coefficients one-by-one
3. that supporting this API would add significant complexity to JuMP and MathOptInterface
4. that few people are solving column generation problems
5. that the add-then-modify approach works okay enough in practice

In practice, any one implementing column generation has a lot of other issues to worry about, like coding the column generator, or doing some sort of branch-and-price, that the cost of modifying some coefficients in the master problem is irrelevant. In this case a simple to use API is more important than performance.

---

<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:** [April 30, 2026, 2:58am UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/17 "2026-04-30T02:58:39Z")

</div>

Cutting plane generation and column generation are dual to each other (primal side column generation amounts to dual side row generation). Both are wonderful methods. In some situations writing code in cut generation is more elegant, while in some other situations it’s reversed.

Today I find it elegant to formulate the “π-ascending” problem as a column generation problem. (This problem is used in Benders decomposition to generate stronger lagrangian cuts via tweaking its slope.)

It really looks way more elegant than the original concave maximization where people spent much effort to add “proximal”, “bundle”, “level” tricks.

My code is

```julia-auto
include("src/Settings.jl");
import Gurobi, JuMP, Random; Random.seed!(2);
function build_master(m)
    JuMP.@variable(m, t) # theta
    JuMP.@variable(m, 0 <= y <= 1) # their first-stage variable is symbolized by "y"
    # For generating trial points, we relax the Int constr initially
    JuMP.@objective(m, Min, 0) # their model has no 1st-stage cost; We don't price theta at iter 0.
end;
function build_subproblem(m)
    JuMP.@variable(m, z, Bin) # copy variable
    JuMP.@variable(m, x >= 0) # 2nd-stage inner variable
    JuMP.@constraints(m, begin
        x + 15z >= 8
        3x + 10z >= 13
        x + 10z >= 7
        2x - 10z >= -1
        2x - 70z >= -49
    end)
    JuMP.@objective(m, Min, x) # we always generate optimality cuts (no time for feasibility cuts)
end;
function build_pas(m) # [pi ascending]: we opt to solve it from its dual side, which is a column generation with MIN sense
    JuMP.@variable(m, y) # to be fixed at y_check
    JuMP.@constraint(m, 0 == y)
    JuMP.@constraint(m, 0 == 1)
    # JuMP.@objective(m, Min, 0)
    JuMP.set_attribute(m, "DualReductions", 0) # to get definite INFEASIBLE code (i.e. 3)
end;
function column_generation(m, y_check, y_col, Qy_col)
    Settings.setxdblattrelement(m, 0, "LB", y_check)
    Settings.setxdblattrelement(m, 0, "UB", y_check)
    ter = Settings.opt_and_ter(m)
    Cd = Cdouble[y_col,1] # an aux container
    if ter == 3 # INFEASIBLE
        Gurobi.GRBaddvar(m.o, 1+1, Cint[0,1], Cd, Qy_col, 0., 1e100, Cchar('C'), C_NULL) == 0 || error()
        return NaN
    end
    ter == 2 || error() # OPTIMAL
    Gurobi.GRBgetdblattrarray(m.o, "Pi", 0, 1, Cd) == 0 || error()
    return Cd[1] # revised slope of the Lagrangian cut
end;
function get_check_and_hat_initial(mst, sub)
    y_check = Settings.getxdblattrelement(mst, 1, "X")
    Settings.setxdblattrelement(sub, 0, "Obj", 0.)
    Settings.setxdblattrelement(sub, 0, "LB", y_check)
    Settings.setxdblattrelement(sub, 0, "UB", y_check)
    Settings.setxcharattrelement(sub, 0, "VType", Cchar('C'))
    Settings.opt_ass_opt(sub)
    p_hat = Settings.getxdblattrelement(sub, 0, "RC")
    y_check, p_hat
end;
function get_subLB_and_col(sub, p_hat)
    Settings.setxdblattrelement(sub, 0, "Obj", -1.0 * p_hat);
    Settings.setxdblattrelement(sub, 0, "LB", 0.0);
    Settings.setxdblattrelement(sub, 0, "UB", 1.0);
    Settings.setxcharattrelement(sub, 0, "VType", Cchar('B'));
    Settings.opt_ass_opt(sub);
    subLB = Settings.getmodeldblattr(sub, "ObjBound")
    y_col = Settings.getxdblattrelement(sub, 0, "X")
    Qy_col = Settings.getxdblattrelement(sub, 1, "X")
    subLB, y_col, Qy_col
end;
function get_mst_sub_pas()
    genv = Settings.Env();mst = Settings.Model(genv);o = mst.moi_backend
    build_master(mst)
    mst = (m = mst, o = o, refi = Ref{Cint}(-99), refd = Ref{Cdouble}(NaN))
    genv = Settings.Env();sub = Settings.Model(genv);o = sub.moi_backend
    build_subproblem(sub)
    sub = (m = sub, o = o, refi = Ref{Cint}(-99), refd = Ref{Cdouble}(NaN))
    genv = Settings.Env();pas = Settings.Model(genv);o = pas.moi_backend
    build_pas(pas)
    pas = (m = pas, o = o, refi = Ref{Cint}(-99), refd = Ref{Cdouble}(NaN))
    mst, sub, pas
end;
mst, sub, pas = get_mst_sub_pas();
Settings.opt_ass_opt(mst);

# iter 1: we have an integer trial point, thus add a strengthed Benders cut
y_check, p_hat = get_check_and_hat_initial(mst, sub) # (0.0, -15.0)
subLB, y_col, Qy_col = get_subLB_and_col(sub, p_hat) # (8.0, 0.0, 8.0)
column_generation(pas, y_check, y_col, Qy_col) # NaN
JuMP.@constraint(mst.m, mst.m[:t] >= p_hat * mst.m[:y] + subLB); # add optimality cut to master
Settings.setxdblattrelement(mst, 0, "Obj", 1.0); # activate theta in master
Settings.opt_ass_opt(mst);
# iter 2: we still have an integer trial point, thus add a strengthed Benders cut
y_check, p_hat = get_check_and_hat_initial(mst, sub) # (1.0, 35.0)
subLB, y_col, Qy_col = get_subLB_and_col(sub, p_hat) # (-24.5, 1.0, 10.5)
column_generation(pas, y_check, y_col, Qy_col) # NaN
JuMP.@constraint(mst.m, mst.m[:t] >= p_hat * mst.m[:y] + subLB); # add optimality cut to master
Settings.opt_ass_opt(mst);
# iter 3: we have a fractional trial point, thus we try to do π-ascending so we might be able to generate a stronger cut
y_check, p_hat = get_check_and_hat_initial(mst, sub) # (0.65, 5.0)

julia> column_generation(pas, y_check, NaN, NaN) == 2.5 && @info "We've made a successful update to the dual vector π"
[ Info: We've made a successful update to the dual vector π

```

Notice that finally I fairly easily find the desired slope `2.5`, comapred with the original subgradient method.

 ![VX%%X0E60Y@I%9GAM94O](https://global.discourse-cdn.com/julialang/original/3X/3/9/39d3816e4899b80a49cea1ada44cc8ca6c95221b.png)

 ![{(R1H3P87BO(D3KVOITW](https://global.discourse-cdn.com/julialang/original/3X/f/1/f1a600d98be5197bcf8292374858d1e500c69ced.png)

---

<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:** [April 30, 2026, 8:27am UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/18 "2026-04-30T08:27:14Z")

</div>

I’m a bit more curious here:

Suppose I have a LP (at least one side is feasible, so strong duality holds).  
The primal formulation (min\_sense) has a fixed number of constraints (let’s say `100`), while I may do column generation, there might be exponentially many columns (i.e. new decision variables). Let’s say we are going to add `10^7` new decision variables on the fly.

If I were, instead, encoding the same problem with the dual formulation, there would be a fixed number of decision variables (`100` dual variables), while I need to generate `10^7` new constraints on the fly.

**Which way would be more efficient?**

in terms of, like, encoding size, and solution speed (here we are iteratively solve new LPs when new row/column is added).

I have no idea on this question. But it is indeed valuable to figure out a definite answer. It may depends on the algorithm that the solver uses (as is known, Gurobi’s default in this case is dual simplex algorithm).

@odow any ideas?

AI tools suggest me to do column generation (with the primal-side min formulation). But I intuitively think the cut generation would not be slow—Gurobi just need to cut off the previous trial point iteratively, with dual simplex method.

---

<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:** [April 30, 2026, 9:17pm UTC](https://discourse.julialang.org/t/gurobi-fails-to-invoke-the-primal-simplex-method-when-doing-column-generation/133675/19 "2026-04-30T21:17:20Z")

</div>

Let’s stick to Julia/JuMP related questions. The trade-offs between Benders and CG are complex and problem dependent.