# JuMP reduced_cost questions

Hi, can anyone please help me with some questions on `reduced_cost` in JuMP? Thank you!

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

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

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
4. Correct. There is no way to compute a reduced cost for arbitrary `x`.

If youâre looking for sensitivity analysis, see

1 Like

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

``````4-element Vector{Float64}:
1.0
2.0
-3.0
-4.0
``````

Thanks again!

1 Like
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.

1 Like

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.

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/

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?

I was reading the JuMP documentation tutorial on 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?

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

``````min     c'x
L <= Ax <= U
l <=  x <= u
``````

solvers rewrite to

``````min     c'x
Ax - Is == 0
l <=  x      <= u
L <=       s <= U
``````

which simplies to

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

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

1 Like

Thank you @odow ! These make things very clear!

To make sure I understand this, I tried something like:

``````@variable(model, x[1:n])
fix.(x, x_values; force = true)
``````

and print the model, then I saw these âconstraintsâ

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

``````@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:

``````@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:

``````@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?

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!

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:

``````@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)
``````
1 Like

Thank you @odow !

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!

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`.)

1 Like

Great! Thank you @odow !!!

1 Like