diff --git a/docs/make.jl b/docs/make.jl index f43b4fbe41c..3a1fe17a9a3 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -140,6 +140,7 @@ const _PAGES = [ "tutorials/linear/transp.md", "tutorials/linear/callbacks.md", "tutorials/linear/constraint_programming.md", + "tutorials/linear/mip_duality.md", ], "Nonlinear programs" => [ "tutorials/nonlinear/introduction.md", diff --git a/docs/src/reference/variables.md b/docs/src/reference/variables.md index 313bb84cf2d..19b6b2c11bf 100644 --- a/docs/src/reference/variables.md +++ b/docs/src/reference/variables.md @@ -95,6 +95,7 @@ BinaryRef ```@docs relax_integrality +fix_discrete_variables ``` ## Extensions diff --git a/docs/src/tutorials/linear/mip_duality.jl b/docs/src/tutorials/linear/mip_duality.jl new file mode 100644 index 00000000000..3feb8af099d --- /dev/null +++ b/docs/src/tutorials/linear/mip_duality.jl @@ -0,0 +1,124 @@ +# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors #src +# This Source Code Form is subject to the terms of the Mozilla Public License #src +# v.2.0. If a copy of the MPL was not distributed with this file, You can #src +# obtain one at https://mozilla.org/MPL/2.0/. #src + +# # Computing the duals of a mixed-integer program + +# This tutorial explains how to compute the duals of a mixed-integer linear +# program by fixing the discrete variables to their optimal solution and +# resolving as a linear program. + +# This tutorial uses the following packages: + +using JuMP +import HiGHS + +# ## The model + +# Our example model is the unit commitment example from [Unit commitment](@ref). +# The details are unimportant, other than to note that there are two types of +# continuous variables, `g` and `w`, representing the quantity of generation +# from thermal and wind plants, and a discrete variable `dispatch`, which is +# `1` if plant `i` is operating, and `0` if not. + +# We are interested in the "dual" of the `power_balance` constraint, because it +# represents the marginal price of electricity that consumers should pay for +# their consumption. + +generators = [ + (min = 0.0, max = 1000.0, fixed_cost = 1000.0, variable_cost = 50.0), + (min = 300.0, max = 1000.0, fixed_cost = 0.0, variable_cost = 100.0), +] +N = length(generators) +model = Model(HiGHS.Optimizer) +set_silent(model) +@variables(model, begin + generators[i].min <= g[i = 1:N] <= generators[i].max + 0 <= w <= 200 + dispatch[i = 1:N], Bin +end) +@constraints(model, begin + power_balance, sum(g[i] for i in 1:N) + w == 1500 + [i = 1:N], g[i] <= generators[i].max * dispatch[i] + [i = 1:N], g[i] >= generators[i].min * dispatch[i] +end) +@objective( + model, + Min, + sum( + generators[i].fixed_cost * dispatch[i] + + generators[i].variable_cost * g[i] for i in 1:N + ) +) +print(model) + +# ## Manually fix the variables + +# If we optimize this model, we obtain a [`dual_status`](@ref) of `NO_SOLUTION`: + +optimize!(model) +dual_status(model) + +# This is because HiGHS cannot compute the duals of a mixed-integer program. We +# can work around this problem by fixing the integer variables to their optimal +# solution, relaxing integrality, and re-solving as a linear program. + +discrete_values = value.(dispatch) +fix.(dispatch, discrete_values; force = true) +unset_binary.(dispatch) +print(model) + +# Now if we re-solve the problem, we obtain a `FEASIBLE_POINT` for the dual: + +optimize!(model) +dual_status(model) + +# and a marginal price of electricity of \$100/MWh: + +dual(power_balance) + +# To reset the problem back to a mixed-integer linear program, we need to +# [`unfix`](@ref) and call [`set_binary`](@ref): + +unfix.(dispatch) +set_binary.(dispatch) +print(model) + +# ## Use `fix_discrete_variables` + +# Manually choosing the variables to relax and fix works for our small example, +# but it becomes more difficult in problems with a larger number of binary and +# integer variables. To automate the process we just did manually, JuMP provides +# the [`fix_discrete_variables`](@ref) function: + +optimize!(model) +dual_status(model) + +#- + +undo = fix_discrete_variables(model); + +# Here `undo` is a function that, when called with no arguments, returns the +# model to the original mixed-integer formulation. + +# !!! tip +# Afer calling [`fix_discrete_variables`](@ref), you can set a new solver +# with [`set_optimizer`](@ref) if your mixed-integer solver does not support +# computing a dual solution. + +print(model) + +#- + +optimize!(model) +dual_status(model) + +#- + +dual(power_balance) + +# Finally, call `undo` to revert the reformulation + +undo() +print(model) diff --git a/src/variables.jl b/src/variables.jl index e75119a956e..cda61191e2f 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -1591,7 +1591,67 @@ Subject to x binary ``` """ -function relax_integrality(model::Model) +relax_integrality(model::Model) = _relax_or_fix_integrality(nothing, model) + +""" + fix_discrete_variables([var_value::Function = value,] model::Model) + +Modifies `model` to convert all binary and integer variables to continuous +variables with fixed bounds of `var_value(x)`. + +## Return + +Returns a function that can be called without any arguments to restore the +original model. The behavior of this function is undefined if additional +changes are made to the affected variables in the meantime. + +## Notes + +- An error is thrown if semi-continuous or semi-integer constraints are + present (support may be added for these in the future). +- All other constraints are ignored (left in place). This includes discrete + constraints like SOS and indicator constraints. + +## Example + +```jldoctest; setup=:(using JuMP) +julia> model = Model(); + +julia> @variable(model, x, Bin, start = 1); + +julia> @variable(model, 1 <= y <= 10, Int, start = 2); + +julia> @objective(model, Min, x + y); + +julia> undo_relax = fix_discrete_variables(start_value, model); + +julia> print(model) +Min x + y +Subject to + x = 1.0 + y = 2.0 + +julia> undo_relax() + +julia> print(model) +Min x + y +Subject to + y ≥ 1.0 + y ≤ 10.0 + y integer + x binary +``` +""" +function fix_discrete_variables(var_value::Function, model::Model) + return _relax_or_fix_integrality(var_value, model) +end + +fix_discrete_variables(model::Model) = fix_discrete_variables(value, model) + +function _relax_or_fix_integrality( + var_value::Union{Nothing,Function}, + model::Model, +) if num_constraints(model, VariableRef, MOI.Semicontinuous{Float64}) > 0 error( "Support for relaxing semicontinuous constraints is not " * @@ -1604,18 +1664,17 @@ function relax_integrality(model::Model) "yet implemented.", ) end - - bin_int_constraints = vcat( + discrete_variable_constraints = vcat( all_constraints(model, VariableRef, MOI.ZeroOne), all_constraints(model, VariableRef, MOI.Integer), ) - bin_int_variables = VariableRef.(bin_int_constraints) - info_pre_relaxation = - map(v -> (v, _info_from_variable(v)), bin_int_variables) - # We gather the info first because some solvers perform poorly when you - # interleave queries and changes. See, e.g., - # https://github.com/jump-dev/Gurobi.jl/pull/301. - for (v, info) in info_pre_relaxation + # We gather the info first because we cannot modify-then-query. + info_pre_relaxation = map(VariableRef.(discrete_variable_constraints)) do v + solution = var_value === nothing ? nothing : var_value(v) + return (v, solution, _info_from_variable(v)) + end + # Now we can modify. + for (v, solution, info) in info_pre_relaxation if info.integer unset_integer(v) elseif info.binary @@ -1630,24 +1689,37 @@ function relax_integrality(model::Model) ) end end + if solution !== nothing + fix(v, solution; force = true) + end end function unrelax() - for (v, info) in info_pre_relaxation + for (v, solution, info) in info_pre_relaxation + if solution !== nothing + unfix(v) + end + if info.has_lb + set_lower_bound(v, info.lower_bound) + end + if info.has_ub + set_upper_bound(v, info.upper_bound) + end if info.integer set_integer(v) - elseif info.binary + end + if info.binary set_binary(v) - if !info.has_fix - if info.has_lb - set_lower_bound(v, info.lower_bound) - else - delete_lower_bound(v) - end - if info.has_ub - set_upper_bound(v, info.upper_bound) - else - delete_upper_bound(v) - end + end + # Now a special case: when binary variables are relaxed, we add + # [0, 1] bounds, but only if the variable was not previously fixed + # and we did not provide a fixed value, and a bound did not already + # exist. In this case, delete the new bounds that we added. + if solution === nothing && info.binary && !info.has_fix + if !info.has_lb + delete_lower_bound(v) + end + if !info.has_ub + delete_upper_bound(v) end end end diff --git a/test/test_variable.jl b/test/test_variable.jl index 9dd80eff82f..f7d0d732e1e 100644 --- a/test/test_variable.jl +++ b/test/test_variable.jl @@ -1331,4 +1331,32 @@ function test_Hermitian_PSD_anon() return end +function test_fix_discrete_variables() + le = JuMP._math_symbol(MIME("text/plain"), :leq) + ge = JuMP._math_symbol(MIME("text/plain"), :geq) + eq = JuMP._math_symbol(MIME("text/plain"), :eq) + model = Model() + @variable(model, x, Bin, start = 1) + @variable(model, 1 <= y <= 10, Int, start = 2) + @objective(model, Min, x + y) + @test sprint(print, model) == + "Min x + y\nSubject to\n y $ge 1.0\n y $le 10.0\n y integer\n x binary\n" + undo_relax = fix_discrete_variables(start_value, model) + @test sprint(print, model) == + "Min x + y\nSubject to\n x $eq 1.0\n y $eq 2.0\n" + undo_relax() + @test sprint(print, model) == + "Min x + y\nSubject to\n y $ge 1.0\n y $le 10.0\n y integer\n x binary\n" + return +end + +function test_fix_discrete_variables_value() + model = Model() + @variable(model, x, Bin, start = 1) + @variable(model, 1 <= y <= 10, Int, start = 2) + @objective(model, Min, x + y) + @test_throws OptimizeNotCalled fix_discrete_variables(model) + return +end + end # module TestVariable