Add fix_discrete_variables for computing the dual of a MIP

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",

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 `relax_integrality`
+
+# 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
+# a `fix` keyword argument to [`relax_integrality`](@ref) function:
+
+optimize!(model)
+dual_status(model)
+
+#-
+
+undo = relax_integrality(model; fix = value);
+
+# Here `undo` is a function that, when called with no arguments, returns the
+# model to the original mixed-integer formulation.
+
+# !!! tip
+#     Afer calling [`relax_integrality`](@ref), you can set a new solver with
+#     [`set_optimizer`](@ref) if your mixed-integer solver does not support
+#     computing a dual solutio.
+
+print(model)
+
+#-
+
+optimize!(model)
+dual_status(model)
+
+#-
+
+dual(power_balance)
+
+# Finally, call `undo` to revert the reformulation
+
+undo()
+print(model)

src/variables.jl

@@ -1543,10 +1543,22 @@ function _info_from_variable(v::VariableRef)
 end
 
 """
-    relax_integrality(model::Model)
+    relax_integrality(model::Model; fix::Union{Nothing,Function} = nothing)
 
 Modifies `model` to "relax" all binary and integrality constraints on
-variables. Specifically,
+variables. In addition, if `fix` is passed a function, then the discrete
+variables are fixed to the value of `fix(x)`.
+
+To fix the model at the last optimal solution, use
+`relax_integrality(model; fix = value)`.
+
+## 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
 
 - Binary constraints are deleted, and variable bounds are tightened if
   necessary to ensure the variable is constrained to the interval ``[0, 1]``.
@@ -1556,17 +1568,14 @@ variables. Specifically,
 - All other constraints are ignored (left in place). This includes discrete
   constraints like SOS and indicator constraints.
 
-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.
+## Example
 
-# Example
 ```jldoctest; setup=:(using JuMP)
 julia> model = Model();
 
-julia> @variable(model, x, Bin);
+julia> @variable(model, x, Bin, start = 1);
 
-julia> @variable(model, 1 <= y <= 10, Int);
+julia> @variable(model, 1 <= y <= 10, Int, start = 2);
 
 julia> @objective(model, Min, x + y);
 
@@ -1582,6 +1591,24 @@ Subject to
 
 julia> undo_relax()
 
+julia> print(model)
+Min x + y
+Subject to
+ y ≥ 1.0
+ y ≤ 10.0
+ y integer
+ x binary
+
+julia> undo_relax = relax_integrality(model; fix = start_value);
+
+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
@@ -1591,7 +1618,7 @@ Subject to
  x binary
 ```
 """
-function relax_integrality(model::Model)
+function relax_integrality(model::Model; fix::Union{Nothing,Function} = nothing)
     if num_constraints(model, VariableRef, MOI.Semicontinuous{Float64}) > 0
         error(
             "Support for relaxing semicontinuous constraints is not " *
@@ -1604,18 +1631,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 = fix === nothing ? nothing : fix(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 +1656,37 @@ function relax_integrality(model::Model)
                 )
             end
         end
+        if solution !== nothing
+            JuMP.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

test/test_variable.jl

@@ -1331,4 +1331,32 @@ function test_Hermitian_PSD_anon()
     return
 end
 
+function test_relax_integrality_fix()
+    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 = relax_integrality(model; fix = start_value)
+    @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_relax_integrality_fix_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 relax_integrality(model; fix = value)
+    return
+end
+
 end  # module TestVariable