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.