# JuMP reduced\_cost questions
**URL:** https://discourse.julialang.org/t/jump-reduced-cost-questions/111484
**Category:** Optimization (Mathematical)
**Tags:** gurobi
**Created:** [March 11, 2024, 8:02pm UTC](https://discourse.julialang.org/t/jump-reduced-cost-questions/111484 "2024-03-11T20:02:58Z")
**Posts on this page:** 11
**Page:** 1
### Author:  [@Ken\_Lin](https://discourse.julialang.org/u/Ken_Lin)
#### Post date: [March 11, 2024, 8:02pm UTC](https://discourse.julialang.org/t/jump-reduced-cost-questions/111484/1 "2024-03-11T20:02:58Z")
Hi, can anyone please help me with some questions on `reduced_cost` in JuMP? Thank you!
```julia
toym = Model(() -> Gurobi.Optimizer(GRB_ENV))
@variable(toym, 0 <= toyx[1:4] <= 1)
@constraint(toym, toyx[1]+toyx[2] <= 1)
@constraint(toym, toyx[3]+toyx[4] == 1)
@constraint(toym, toyx[2]+toyx[4] >= 1)
@objective(toym, Min, toyx[1]+2*toyx[2]-3*toyx[3]-4*toyx[4])
optimize!(toym)
```
The optimal objective value is `-4`, and the optimal solution is `[0, 0, 0, 1]`.
Calling `reduced_cost.(toyx)` gives me
```julia
4-element Vector{Float64}:
1.0
2.0
0.0
-1.0
```
My questions on this:
1. Does reduced cost make sense for linear program in arbitrary forms?
2. I don’t understand why `toyx[3]` has reduced cost `0`? If `toyx[3]` increases from 0 to 1, the corresponding solution is `[0,1,1,0]` objective changes to -1, so the reduced cost should be `3`, correct?
3. Why is the reduced cost of `toyx[4]` equal to `-1`? I notice this problem is not in standard form, but `toyx[4]` would be a basic variable in the corresponding standard form for this problem, and should have reduced cost `0`, right?
4. `reduced_cost(x)` only works for an optimal solution `x`, correct? Is there a way to get the reduced cost for any given `x`?
I understand the reduced cost concept by the following definition from wikipedia
`In linear programming, reduced cost, or opportunity cost, is the amount by which an objective function coefficient would have to improve (so increase for maximization problem, decrease for minimization problem) before it would be possible for a corresponding variable to assume a positive value in the optimal solution.`
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [March 12, 2024, 3:57am UTC](https://discourse.julialang.org/t/jump-reduced-cost-questions/111484/2 "2024-03-12T03:57:56Z")
Hi @Ken_Lin, hopefully some answers:
1. Yes
2. The reduced cost is not the change in objective if the variable was increased by one unit, but the change in the objective with respect to an infinitesimal change in the variable that is then scaled to one unit. Think of it like `df/dx` instead of `f(x+1) - f(x)`. Presumably, a small increase in `toyx[3]` leads to a change in the other variables that cancels out the change in objective. (If `toyx[3]` increases, then by the second constraint, `toyx[4]` must decrease, etc)
3. `toys[4]` is a non-basic variable at its upper bound of `1`. You can query basis information using [Basis matrices · JuMP](https://jump.dev/JuMP.jl/dev/tutorials/linear/basis/)
4. Correct. There is no way to compute a reduced cost for arbitrary `x`.
If you’re looking for sensitivity analysis, see
> **[Sensitivity analysis of a linear program · JuMP](https://jump.dev/JuMP.jl/stable/tutorials/linear/lp_sensitivity/)**
>
> Documentation for JuMP.
---
### Author:  [@Ken\_Lin](https://discourse.julialang.org/u/Ken_Lin)
#### Post date: [March 12, 2024, 5:13pm UTC](https://discourse.julialang.org/t/jump-reduced-cost-questions/111484/3 "2024-03-12T17:13:56Z")
Thank you @odow for your answers, they help a lot!
May I ask a few more things please?
1. Sorry I might have got this wrong: suppose `toyx[3]` now equals `a > 0`, then because of the last two constraints, `toyx[4] = 1-a`, `toyx[2] >= a`, and therefore the objective value is at least `2a-3a-4(1-a) = 3a-4`. Compared with the original optimal value `-4`, isn’t the change always `3a`, which is nonzero?
2. I see! I add slack variables and make the example into standard form (`Ax=b, x>=0`), then the reduced cost of `toyx[4]` becomes `0`. Given these two models are equivalent, with the same optimal solution, why does the same variable can have different reduced costs? intuitively I think the same change in `toyx[4]` should lead to the same change in the objective values of the two models, what is wrong with this intuition?
3. Suppose I would like to compute the reduced cost for a given vector `toyx_values`, I tried adding `fix.(toyx, toyx_values; force = true)`, then `optimize!(toym)`, then `reduced_cost.(toyx)` seems to work, this looks like a strange approach, could you please let me know if this makes sense?
For example, I added `fix.(toyx, [1,0,0,1]; force = true)`, and `reduced_cost.(toyx)` returns
```julia
4-element Vector{Float64}:
1.0
2.0
-3.0
-4.0
```
Thanks again!
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [March 12, 2024, 9:25pm UTC](https://discourse.julialang.org/t/jump-reduced-cost-questions/111484/4 "2024-03-12T21:25:38Z")
1. Perhaps this is a poor choice for a toy model, because it has a degenerate solution. `toyx[3]` is basic, so it has a reduced cost of `0`, but it is also at it’s lower bound.
2. Because `toys[4]` is non-basic at its upper bound. When you switch standard forms, it is now a basic variable.
3. Now it is just returning you their objective coefficients. Fixing is equivalent to setting the lower and upper bounds to the value, so they become the non-basic constraints. Why do you need the reduced cost of an arbitrary vector? What result do you expect to get?
> intuitively I think the same change in `toyx[4]` should lead to the same change in the objective values
This is one _interpretation_ of a reduced cost, but it is not how they are calculated. It depends on the formulation and what the basic variables are, etc.
As a comment, I think a lot of your questions stem from the fact that Simplex/reduced costs etc are often taught for `Ax = b, x >= 0`, but Gurobi etc do not solve with this formulation, but with lower and upper variable bounds.
---
### Author:  [@Ken\_Lin](https://discourse.julialang.org/u/Ken_Lin)
#### Post date: [March 13, 2024, 5:18pm UTC](https://discourse.julialang.org/t/jump-reduced-cost-questions/111484/5 "2024-03-13T17:18:25Z")
Thanks @odow for your help, that’s a good comment! You are right, I got confused because I learned the concepts with `Ax=b, x>=0`, and I’m not sure how the solvers do the calculation.
> [@odow](#):
>
> Now it is just returning you their objective coefficients
JuMP documentation has a tutorial on `Computing the duals of a mixed-integer program `,
[https://jump.dev/JuMP.jl/v1.13/tutorials/linear/mip\_duality/](https://jump.dev/JuMP.jl/v1.13/tutorials/linear/mip_duality/)
where they fixed the binary variables to some value, solve the model, and compute the dual solution, so I thought it might make sense in JuMP to fix everything and compute things like dual solutions or reduced costs… Perhaps I understand the tutorial wrong?
> [@odow](#):
>
> Why do you need the reduced cost of an arbitrary vector? What result do you expect to get?
I was reading the JuMP documentation tutorial on column generation,
[https://jump.dev/JuMP.jl/v1.13/tutorials/algorithms/cutting\_stock\_column\_generation/](https://jump.dev/JuMP.jl/v1.13/tutorials/algorithms/cutting_stock_column_generation/)
they solve a relaxation, and at the optimal solution they generate new variables (columns) based on the value of the dual variables, so I was thinking if I could generate new columns not only at the optimal solution, but also other vectors? Do you think this is possible?
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [March 13, 2024, 8:42pm UTC](https://discourse.julialang.org/t/jump-reduced-cost-questions/111484/6 "2024-03-13T20:42:34Z")
> I’m not sure how the solvers do the calculation.
This isn’t a very rigorous explanation, but I hope if gives you the gist.
Given
```julia
min c'x
L <= Ax <= U
l <= x <= u
```
solvers rewrite to
```julia
min c'x
Ax - Is == 0
l <= x <= u
L <= s <= U
```
which simplies to
```julia
min d'y
By == 0
e <= y <= f
```
Then, instead of thinking of varaibles being basic if non-zero, and non-basic if
zero, variables can be non-basic at either their lower or upper bounds.
Here’s your original model:
```Julia
julia> using JuMP, Gurobi
julia> function main_1()
model = Model(Gurobi.Optimizer)
set_silent(model)
@variable(model, 0 <= x[1:4] <= 1)
@constraint(model, x[1]+x[2] <= 1)
@constraint(model, x[3]+x[4] == 1)
@constraint(model, x[2]+x[4] >= 1)
@objective(model, Min, x[1]+2*x[2]-3*x[3]-4*x[4])
optimize!(model)
return value.(x), reduced_cost.(x)
end
main_1 (generic function with 1 method)
julia> main_1()
([0.0, 0.0, 0.0, 1.0], [1.0, 2.0, 0.0, -1.0])
```
Rewriting to matrix form
```Julia
julia> function main_2()
l, u = zeros(4), ones(4)
L, U = [-Inf, 1, 1], [1, 1, Inf]
c = [1, 2, -3, -4]
A = [1 1 0 0; 0 0 1 1; 0 1 0 1]
model = Model(Gurobi.Optimizer)
set_silent(model)
@variable(model, l[i] <= x[i = 1:4] <= u[i])
@objective(model, Min, c' * x)
@constraint(model, L .<= A * x .<= U)
optimize!(model)
return value.(x), reduced_cost.(x)
end
main_2 (generic function with 1 method)
julia> main_2()
([0.0, 0.0, 0.0, 1.0], [1.0, 2.0, 0.0, -1.0])
```
Adding `s`
```Julia
julia> function main_3()
l, u = zeros(4), ones(4)
L, U = [-Inf, 1, 1], [1, 1, Inf]
c = [1, 2, -3, -4]
A = [1 1 0 0; 0 0 1 1; 0 1 0 1]
I = [1 0 0; 0 1 0; 0 0 1]
model = Model(Gurobi.Optimizer)
set_silent(model)
@variable(model, l[i] <= x[i = 1:4] <= u[i])
@variable(model, L[i] <= s[i = 1:3] <= U[i])
@objective(model, Min, c' * x)
@constraint(model, A * x - I * s .== 0)
optimize!(model)
return value.(x), reduced_cost.(x)
end
main_3 (generic function with 1 method)
julia> main_3()
([0.0, 0.0, 0.0, 1.0], [1.0, 2.0, 0.0, -1.0])
```
To `y`
```Julia
julia> function main_4()
l, u = zeros(4), ones(4)
L, U = [-Inf, 1, 1], [1, 1, Inf]
c = [1, 2, -3, -4]
A = [1 1 0 0; 0 0 1 1; 0 1 0 1]
I = [1 0 0; 0 1 0; 0 0 1]
d = [c; zeros(3)]
e, f = [l; L], [u; U]
B = [A -I]
model = Model(Gurobi.Optimizer)
set_silent(model)
@variable(model, e[i] <= y[i = 1:7] <= f[i])
@objective(model, Min, d' * y)
@constraint(model, B * y .== 0)
optimize!(model)
x = y[1:4]
return value.(x), reduced_cost.(x)
end
main_4 (generic function with 1 method)
julia> main_4()
([0.0, 0.0, 0.0, 1.0], [1.0, 2.0, 0.0, -1.0])
```
Now look at optimal and dual solution:
```Julia
julia> function main_5()
l, u = zeros(4), ones(4)
L, U = [-Inf, 1, 1], [1, 1, Inf]
c = [1, 2, -3, -4]
A = [1 1 0 0; 0 0 1 1; 0 1 0 1]
I = [1 0 0; 0 1 0; 0 0 1]
d = [c; zeros(3)]
e, f = [l; L], [u; U]
B = [A -I]
model = Model(Gurobi.Optimizer)
set_silent(model)
@variable(model, e[i] <= y[i = 1:7] <= f[i])
@objective(model, Min, d' * y)
@constraint(model, cons, B * y .== 0)
optimize!(model)
# Reproduce the primal solution by looking at the basis
basis = get_attribute.(y, MOI.VariableBasisStatus())
l_indices = basis .== MOI.NONBASIC_AT_LOWER
u_indices = basis .== MOI.NONBASIC_AT_UPPER
b_indices = basis .== MOI.BASIC
# Storage for optimal solution
y_star = fill(0.0, length(y))
# Some columns are fixed
y_star[l_indices] .= e[l_indices]
y_star[u_indices] .= f[u_indices]
# Move fixed colums in By == 0 to the right-hand side to create
# B[:, b_indices] y[b_indices] == b
b = -B[:, l_indices] * e[l_indices] + -B[:, u_indices] * f[l_indices]
# Basic column solutions
y_star[b_indices] .= B[:, b_indices] \ b
# Check we get optimal solution
@assert value.(y) ≈ y_star
# Reproduce the reduced costs with "c - A'x", which is now `d - B'y` for us
# This follows pretty similarly.
π = dual.(cons)
rc = d .- B' * π
@assert reduced_cost.(y) ≈ rc
end
main_5 (generic function with 1 method)
julia> main_5()
```
> compute the dual solution
Mixed-integer programs don’t have a dual solution. This is a fake “approximation” of a dual that people use sometimes. It mostly works when there are many continuous variables and few integer variables. We fix _only_ the integer variables, so it’s really like solving an LP with some choices fixed (like which generator to turn on). It wouldn’t work if we also fixed continuous variables.
The reduced costs come from the dual solution. If you fix your variable values, a dual solution still exists, but it’s probably not very related to the dual solution of your original model.
It also hints as why we can’t get a reduced cost for an arbitrary point: we’d need a primal-dual feasible solution for an arbitrary point, which solvers don’t easily return.
> I was thinking if I could generate new columns not only at the optimal solution, but also other vectors?
In general, people don’t do this because the top-level problem is simple to solve, and most time is spent finding new columns in the pricing subproblem.
Perhaps what you are getting at is this: you don’t need to solve the LP or get reduced costs etc to create columns. You might have a different heuristic that can look at a solution and create new columns. It doesn’t really matter, so long as you have a way to find new columns that lead to a better solution.
---
### Author:  [@Ken\_Lin](https://discourse.julialang.org/u/Ken_Lin)
#### Post date: [March 15, 2024, 2:46am UTC](https://discourse.julialang.org/t/jump-reduced-cost-questions/111484/7 "2024-03-15T02:46:17Z")
Thank you @odow ! These make things very clear!
> [@odow](#):
>
> Fixing is equivalent to setting the lower and upper bounds to the value
To make sure I understand this, I tried something like:
```julia
@variable(model, x[1:n])
fix.(x, x_values; force = true)
```
and print the model, then I saw these “constraints”
```julia
x[1] = x_values[1]
x[2] = x_values[2]
...
x[n] = x_values[n]
```
What you are saying is that `fix` is the same as
```julia
@variable(model, x_values[i] <= x[i=1:n] <= x_values[i])
```
correct? But how is it different from `@constraint(model, x .== x_values)`?
Again in the JuMP column generation tutorial, they add the new column as a new variable in the model, and update the objective and constraints. However, in terms of indexing, the new column will always be `x[end]`. In my model it is much more convenient to do this with all the index fixed at the beginning, would it be possible? That is, say I have variables `x[1:n]`, the initial set of selected columns `i_selected`, and the set of remaining columns `i_not_selected`, such that `i_selected` and `i_not_selected` are disjoint and their union is `1:n`. Are the following two codes equivalent? For the same constraint in `cons`, do I get the same dual solution in both codes?
Code 1:
```julia
@variable(model, l <= x[1:n] <= u)
fix.(x[i_not_selected], zeros(length(i_not_selected)); force = true)
@constraint(model, cons, L .<= A * x .<= U)
@objective(model, Min, c' * x)
optimize!(model)
```
Code 2:
```julia
@variable(model, l <= x[1:length(i_selected)] <= u)
@constraint(model, cons, L .<= A[:, i_selected] * x .<= U)
@objective(model, Min, c[i_selected]' * x)
optimize!(model)
```
If they are equivalent and give the same dual solutions, it is more convenient to use Code 1 in my case, where the indexing of `x[1:n]` is fixed, and when I generate a new column `i`, I add `i` to `i_selected` and remove `i` from `i_not_selected`.
By the way, do the solvers (like Gurobi) always convert things into the below form you mention (min, equality constraints, variable upper and lower bounds), and use this form to solve the problem? Similarly there are many equivalent formulations of the dual, then what formulation of the dual do they use?
> [@odow](#):
>
> which simplies to
>
> ```julia
> min d'y
> By == 0
> e <= y <= f
>
> ```
As a related question, for some `j`, `@constraint(model, con, L[j] <= A[j, :]' * x) <= U[j]` is actually two constraints `L[j] <= A[j, :]' * x` and `A[j, :]' * x <= U[j]`, correct? But `dual(con)` only returns a single value? Does it convert it into some equality constraint as in your example, and return the dual for that equality?
Thanks again!
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [March 15, 2024, 2:57am UTC](https://discourse.julialang.org/t/jump-reduced-cost-questions/111484/8 "2024-03-15T02:57:27Z")
> [@Ken\_Lin](#):
>
> What you are saying is that `fix` is the same as
>
> ```julia
> @variable(model, x_values[i] <= x[i=1:n] <= x_values[i])
>
> ```
>
> correct? But how is it different from `@constraint(model, x .== x_values)`?
Correct.
`@constraint(model, x .== x_values)` This adds new rows to the constraint matrix. `fix` modifies the variable bounds.
> In my model it is much more convenient to do this with all the index
You don’t have to add to the end. That’s just the data structure we use.
> Are the following two codes equivalent?
Yes
> `fix.(x[i_not_selected], zeros(length(i_not_selected)); force = true)`
You can do `fix.(x[i_not_selected], 0; force = true)`
> always convert things into the below form you mention (min, equality constraints, variable upper and lower bounds), and use this form to solve the problem
No. It again depends on the solver, and what things it supports, and what algorithm it uses to solve the problem.
> But `dual(con)` only returns a single value?
It returns the dual of the active bound (if any). The dual of the side that is non-binding is `0`.
Stepping back though, you shouldn’t need to worry about these details. Use your original formulation (`main_1`), and use the reduced costs that come out of it. JuMP handles the reformulations, you shouldn’t have to. And you don’t need to use the `L .<= A * x .<= U` formulation.
You can get the dual vector from the original formulation as:
```julia
@constraint(model, c1, x[1]+x[2] <= 1)
@constraint(model, c2, x[3]+x[4] == 1)
@constraint(model, c3, x[2]+x[4] >= 1)
cons = [c1, c2, c3]
optimize!(model)
dual.(cons)
```
---
### Author:  [@Ken\_Lin](https://discourse.julialang.org/u/Ken_Lin)
#### Post date: [March 15, 2024, 10:53pm UTC](https://discourse.julialang.org/t/jump-reduced-cost-questions/111484/9 "2024-03-15T22:53:22Z")
Thank you @odow !
> [@odow](#):
>
> Stepping back though, you shouldn’t need to worry about these details. Use your original formulation (`main_1`), and use the reduced costs that come out of it. JuMP handles the reformulations, you shouldn’t have to. And you don’t need to use the `L .<= A * x .<= U` formulation.
Are you suggesting this: In general, for any formulation, I can always understand JuMP’s outputs (like all the results and also dual solutions, reduced costs,…) in terms of the original formulations I give to JuMP? If this is the case it would be nice! It is indeed difficult if some of the results JuMP returns only make sense for JuMP’s reformulation, not my original ones, then I would have to know what reformulations they solve. Thanks!
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [March 17, 2024, 8:11pm UTC](https://discourse.julialang.org/t/jump-reduced-cost-questions/111484/10 "2024-03-17T20:11:43Z")
> in terms of the original formulations I give to JuMP?
Yes.
I bought up the reformulations only to show how the solvers arrive at the values. (Your question was why some specific reduced costs were `0` or `2`.)
---
### Author:  [@Ken\_Lin](https://discourse.julialang.org/u/Ken_Lin)
#### Post date: [March 20, 2024, 2:31pm UTC](https://discourse.julialang.org/t/jump-reduced-cost-questions/111484/11 "2024-03-20T14:31:20Z")
Great! Thank you @odow !!!