Variables are not copied correctly using `copy_model` with `filter_constraints`

[phguo]

Hi,

I am implementing a Benders decomposition. I hope to build the master problem use copy_model with filter_constraints, but the lower_bound is not copied to the new model.

using JuMP
import GLPK

function original_model()
    model = Model(GLPK.Optimizer, add_bridges = false)

    f = [i for i=1:3]
    c = [i for i=1:3]
    y_dim = length(f)
    x_dim = length(c)
    A = randn(3, 3)
    b = [30 for i=1:3]
    B = randn(3, 3)
    D = randn(3, 3)
    d = [30 for i=1:3]

    @variable(model, x[1:x_dim], lower_bound = 0, base_name = "__S__")
    @variable(model, y[1:y_dim], Int, lower_bound = 0, base_name = "__M__")
    @constraint(model, A * y .<= b, base_name = "__M__")
    @constraint(model, B * y + D * x .<= d)
    @expression(model, obj, sum(f[i] * x[i] for i=1:x_dim) + sum(c[i] * y[i] for i=1:y_dim))
    @objective(model, Min, obj)
    return model
end

# The original model.
o_model = original_model()

# The lower_bound of the original model.
for var in all_variables(o_model)
    if occursin("__M__", name(var))
        print(name(var), lower_bound(var), '\n')
    end
end

# Copy constraint with "__M__" in name.
cons_filter(cons::ConstraintRef) = occursin("__M__", name(cons))
m_model, reference_map = copy_model(o_model, filter_constraints=cons_filter)

# The lower_bound of the copied model.
for var in all_variables(m_model)
    if occursin("__M__", name(var))
        print(name(var), lower_bound(var), '\n')
    end
end

The output:

__M__[1]0.0
__M__[2]0.0
__M__[3]0.0
ERROR: Variable __M__[1] does not have a lower bound.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:33
 [2] lower_bound(v::VariableRef)
   @ JuMP ~/.julia/packages/JuMP/0C6kd/src/variables.jl:571
 [3] top-level scope
   @ ~/GitHub/JuBD.jl/src/problem.jl:44

However, when filter_constraints=cons_filter is removed, the lower bounds are copied correctly.

[odow]

If you print out each ref, you'll see that variable bounds have an empty name "", so they fail your comparison:

julia> function cons_filter(cons::ConstraintRef)
           @show cons, name(cons)
           return occursin("__M__", name(cons))
       end
cons_filter (generic function with 1 method)

julia> m_model, reference_map = copy_model(o_model, filter_constraints=cons_filter);
(cons, name(cons)) = (__S__[1] ≥ 0.0, "")
(cons, name(cons)) = (__S__[2] ≥ 0.0, "")
(cons, name(cons)) = (__S__[3] ≥ 0.0, "")
(cons, name(cons)) = (__M__[1] ≥ 0.0, "")
(cons, name(cons)) = (__M__[2] ≥ 0.0, "")
(cons, name(cons)) = (__M__[3] ≥ 0.0, "")
(cons, name(cons)) = (__M__[1] integer, "")
(cons, name(cons)) = (__M__[2] integer, "")
(cons, name(cons)) = (__M__[3] integer, "")
(cons, name(cons)) = (__M__ : 1.432853835876625 __M__[1] + 0.667998045246818 __M__[2] - 1.3386057944346192 __M__[3] ≤ 30.0, "__M__")
(cons, name(cons)) = (__M__ : -0.8646067833456793 __M__[1] + 0.18869417558811097 __M__[2] - 0.35356068773861055 __M__[3] ≤ 30.0, "__M__")
(cons, name(cons)) = (__M__ : 0.1434215900934186 __M__[1] + 0.17988978977901005 __M__[2] - 1.4528201290040352 __M__[3] ≤ 30.0, "__M__")
(cons, name(cons)) = (-1.34717267442682 __S__[1] - 0.4483362319424767 __S__[2] + 0.8137697717052725 __S__[3] - 0.14113468042813018 __M__[1] - 0.8578201592682566 __M__[2] - 0.37916046919022073 __M__[3] ≤ 30.0, "")
(cons, name(cons)) = (0.5331974637213535 __S__[1] - 1.4819971199973525 __S__[2] - 0.8583056394197597 __S__[3] - 0.28090757509492004 __M__[1] + 0.933894320164304 __M__[2] + 0.8947800620933404 __M__[3] ≤ 30.0, "")
(cons, name(cons)) = (1.0415461567128537 __S__[1] + 1.2536740165273854 __S__[2] - 1.838837482913547 __S__[3] + 2.5073324473937584 __M__[1] + 0.26035700398151207 __M__[2] - 0.5824176257694471 __M__[3] ≤ 30.0, "")
(cons, name(cons)) = (__M__ : 1.432853835876625 __M__[1] + 0.667998045246818 __M__[2] - 1.3386057944346192 __M__[3] ≤ 30.0, "__M__")
(cons, name(cons)) = (__M__ : -0.8646067833456793 __M__[1] + 0.18869417558811097 __M__[2] - 0.35356068773861055 __M__[3] ≤ 30.0, "__M__")
(cons, name(cons)) = (__M__ : 0.1434215900934186 __M__[1] + 0.17988978977901005 __M__[2] - 1.4528201290040352 __M__[3] ≤ 30.0, "__M__")
(cons, name(cons)) = (-1.34717267442682 __S__[1] - 0.4483362319424767 __S__[2] + 0.8137697717052725 __S__[3] - 0.14113468042813018 __M__[1] - 0.8578201592682566 __M__[2] - 0.37916046919022073 __M__[3] ≤ 30.0, "")
(cons, name(cons)) = (0.5331974637213535 __S__[1] - 1.4819971199973525 __S__[2] - 0.8583056394197597 __S__[3] - 0.28090757509492004 __M__[1] + 0.933894320164304 __M__[2] + 0.8947800620933404 __M__[3] ≤ 30.0, "")
(cons, name(cons)) = (1.0415461567128537 __S__[1] + 1.2536740165273854 __S__[2] - 1.838837482913547 __S__[3] + 2.5073324473937584 __M__[1] + 0.26035700398151207 __M__[2] - 0.5824176257694471 __M__[3] ≤ 30.0, "")
(cons, name(cons)) = (__S__[1] ≥ 0.0, "")
(cons, name(cons)) = (__S__[2] ≥ 0.0, "")
(cons, name(cons)) = (__S__[3] ≥ 0.0, "")
(cons, name(cons)) = (__M__[1] ≥ 0.0, "")
(cons, name(cons)) = (__M__[2] ≥ 0.0, "")
(cons, name(cons)) = (__M__[3] ≥ 0.0, "")
(cons, name(cons)) = (__M__[1] integer, "")
(cons, name(cons)) = (__M__[2] integer, "")
(cons, name(cons)) = (__M__[3] integer, "")

You'd need something like this:

cons_filter(cons::ConstraintRef) = occursin("__M__", name(cons))

function cons_filter(cons::ConstraintRef{Model,<:MOI.ConstraintIndex{MOI.VariableIndex}})
    x = VariableRef(owner_model(cons), MOI.VariableIndex(index(cons).value))
    return occursin("__M__", name(x))
end

m_model, reference_map = copy_model(o_model, filter_constraints=cons_filter)

[odow]

There are two action items here:

1. VariableRef(::ConstraintRef) that returns the variable
2. Look into why the variable constraints were printed twice

[phguo]

Thank you @odow ! It works for me.

I am new to Julia and JuMP. Am I correct that lower bounds are saved as constraints? Would it be better to save lower bounds as variable attributes?

[odow]

I am new to Julia and JuMP

Welcome! What are you trying to do? Copying models with a filter is quite an advanced operation, and not something I would recommend for new users. There is likely a better way to achieve what you are trying to do.

Am I correct that lower bounds are saved as constraints?

Yes.

Would it be better to save lower bounds as variable attributes?

No. JuMP uses this standard form for representing problems: https://jump.dev/JuMP.jl/stable/moi/manual/standard_form

Take a read of https://pubsonline.informs.org/doi/10.1287/ijoc.2021.1067 (Or slightly older pre-print: https://arxiv.org/abs/2002.03447).

[odow]

What are you trying to do?

Oh, I see, you want the user to write the full model, and then you want to decompose based on the name of each constraint.

don't do this. Ask the user to explicitly model the first- and second-stage subproblems, and get them to declare which are the state variables that map between the two subproblems.

Potentially useful links:

• https://github.com/martinbiel/StochasticPrograms.jl
• https://github.com/odow/SDDP.jl
• https://github.com/atoptima/Coluna.jl
• https://jump.dev/JuMP.jl/stable/tutorials/algorithms/benders_decomposition/

[odow]

Look into why the variable constraints were printed twice

This isn't a bug. The constraints get copied once, and then it gets called a second time to see if it should be in the extension dictionary: https://github.com/jump-dev/JuMP.jl/blob/249d96938fafcfee1aca84c2173573fca967064f/src/copy.jl#L165

[phguo]

Thank you @odow! The links you provide are valuable for me.

I will try your suggestion

Ask the user to explicitly model the first- and second-stage subproblems

[odow]

I will try your suggestion

Great. If you have other questions, please post on the community forum: https://discourse.julialang.org/c/domain/opt/13. There are quite a few people who have developed similar things, so we should be able to help if you get stuck.