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.

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

If you’re looking for sensitivity analysis, see

1 Like

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

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.

Here’s your original model:


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