diff --git a/LICENSE b/LICENSE index 901ac1b5..467a9c8c 100644 --- a/LICENSE +++ b/LICENSE @@ -1,21 +1,22 @@ -MIT License +The BilevelJuMP.jl package is licensed under the MIT "Expat" License: -Copyright (c) 2019 Joaquim Dias Garcia - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. +> Copyright (c) 2019 Joaquim Dias Garcia, and contributors +> +> Permission is hereby granted, free of charge, to any person obtaining a copy +> of this software and associated documentation files (the "Software"), to deal +> in the Software without restriction, including without limitation the rights +> to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +> copies of the Software, and to permit persons to whom the Software is +> furnished to do so, subject to the following conditions: +> +> The above copyright notice and this permission notice shall be included in all +> copies or substantial portions of the Software. +> +> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +> IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +> FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +> AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +> LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +> OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +> SOFTWARE. +> \ No newline at end of file diff --git a/Project.toml b/Project.toml new file mode 100644 index 00000000..a5758265 --- /dev/null +++ b/Project.toml @@ -0,0 +1,24 @@ +name = "BilevelJuMP" +uuid = "485130c0-026e-11ea-0f1a-6992cd14145c" +authors = ["Joaquim Garcia "] +version = "0.1.0" + +[deps] +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +Dualization = "191a621a-6537-11e9-281d-650236a99e60" + +[compat] +julia = "1" +MathOptInterface = "~0.9.0" +JuMP = "~0.20.0" +Dualization = "~0.2.0" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b" +GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6" +MathOptFormat = "f4570300-c277-12e8-125c-4912f86ce65d" + +[targets] +test = ["Test", "Gurobi", "GLPK", "MathOptFormat"] diff --git a/src/BilevelJuMP.jl b/src/BilevelJuMP.jl new file mode 100644 index 00000000..7fabd4df --- /dev/null +++ b/src/BilevelJuMP.jl @@ -0,0 +1,16 @@ +module BilevelJuMP + +using MathOptInterface +const MOI = MathOptInterface +const MOIU = MathOptInterface.Utilities +using JuMP +using Dualization + + +export +BilevelModel, UpperToLower, LowerToUpper, Upper, Lower + +include("moi.jl") +include("jump.jl") + +end \ No newline at end of file diff --git a/src/jump.jl b/src/jump.jl new file mode 100644 index 00000000..4dc0b65a --- /dev/null +++ b/src/jump.jl @@ -0,0 +1,548 @@ + +# The following is largely inspired from JuMP/test/JuMPExtension.jl + +@enum Level BOTH LOWER UPPER + +abstract type AbstractBilevelModel <: JuMP.AbstractModel end + +mutable struct BilevelModel <: AbstractBilevelModel + # Structured data + upper::JuMP.AbstractModel + lower::JuMP.AbstractModel + + # Model data + nextvaridx::Int # Next variable index is nextvaridx+1 + variables::Dict{Int, JuMP.AbstractVariable} # Map varidx -> variable + varnames::Dict{Int, String} # Map varidx -> name + var_level::Dict{Int, Level} + var_upper::Dict{Int, JuMP.AbstractVariableRef} + var_lower::Dict{Int, JuMP.AbstractVariableRef} + + # upper level decisions that are "paramters" of the second level + upper_to_lower_link::Dict{JuMP.AbstractVariableRef, JuMP.AbstractVariableRef} + # lower level decisions that are input to upper + lower_to_upper_link::Dict{JuMP.AbstractVariableRef, JuMP.AbstractVariableRef} + # joint link + link::Dict{JuMP.AbstractVariableRef, JuMP.AbstractVariableRef} + + nextconidx::Int # Next constraint index is nextconidx+1 + constraints::Dict{Int, JuMP.AbstractConstraint} # Map conidx -> variable + connames::Dict{Int, String} # Map varidx -> name + ctr_level::Dict{Int, Level} + ctr_upper::Dict{Int, JuMP.ConstraintRef} + ctr_lower::Dict{Int, JuMP.ConstraintRef} + + upper_objective_sense::MOI.OptimizationSense + upper_objective_function::JuMP.AbstractJuMPScalar + + lower_objective_sense::MOI.OptimizationSense + lower_objective_function::JuMP.AbstractJuMPScalar + + # solution data + solver#::MOI.ModelLike + upper_to_sblm + lower_to_sblm + sblm_to_solver + + objdict::Dict{Symbol, Any} # Same that JuMP.Model's field `objdict` + + function BilevelModel() + + model = new( + JuMP.Model(), + JuMP.Model(), + + # var + 0, Dict{Int, JuMP.AbstractVariable}(), Dict{Int, String}(), # Model Variables + Dict{Int, Level}(), Dict{Int, JuMP.AbstractVariable}(), Dict{Int, JuMP.AbstractVariable}(), + # links + Dict{Int, JuMP.AbstractVariable}(), Dict{Int, JuMP.AbstractVariable}(), + Dict{Int, JuMP.AbstractVariable}(), + #ctr + 0, Dict{Int, JuMP.AbstractConstraint}(), Dict{Int, String}(), # Model Constraints + Dict{Int, Level}(), Dict{Int, JuMP.AbstractConstraint}(), Dict{Int, JuMP.AbstractConstraint}(), + #obj + MOI.FEASIBILITY_SENSE, + zero(JuMP.GenericAffExpr{Float64, BilevelVariableRef}), # Model objective + MOI.FEASIBILITY_SENSE, + zero(JuMP.GenericAffExpr{Float64, BilevelVariableRef}), # Model objective + + nothing, + nothing, + nothing, + nothing, + Dict{Symbol, Any}(), + ) + + return model + end +end + +abstract type InnerBilevelModel <: AbstractBilevelModel end +struct UpperModel <: InnerBilevelModel + m::BilevelModel +end +Upper(m::BilevelModel) = UpperModel(m) +struct LowerModel <: InnerBilevelModel + m::BilevelModel +end +Lower(m::BilevelModel) = LowerModel(m) +bilevel_model(m::InnerBilevelModel) = m.m +mylevel_model(m::UpperModel) = bilevel_model(m).upper +mylevel_model(m::LowerModel) = bilevel_model(m).lower +level(m::LowerModel) = LOWER +level(m::UpperModel) = UPPER +mylevel_ctr_list(m::LowerModel) = bilevel_model(m).ctr_lower +mylevel_ctr_list(m::UpperModel) = bilevel_model(m).ctr_upper +mylevel_var_list(m::LowerModel) = bilevel_model(m).var_lower +mylevel_var_list(m::UpperModel) = bilevel_model(m).var_upper + +# obj + +mylevel_obj_sense(m::LowerModel) = bilevel_model(m).lower_objective_sense +mylevel_obj_function(m::LowerModel) = bilevel_model(m).lower_objective_function +mylevel_obj_sense(m::UpperModel) = bilevel_model(m).upper_objective_sense +mylevel_obj_function(m::UpperModel) = bilevel_model(m).upper_objective_function + +set_mylevel_obj_sense(m::LowerModel, val) = bilevel_model(m).lower_objective_sense = val +set_mylevel_obj_function(m::LowerModel, val) = bilevel_model(m).lower_objective_function = val +set_mylevel_obj_sense(m::UpperModel, val) = bilevel_model(m).upper_objective_sense = val +set_mylevel_obj_function(m::UpperModel, val) = bilevel_model(m).upper_objective_function = val + +# UpperToLower / LowerParameter / ParamterInLower +# LowerToUpper / ArgMin +abstract type BridgeBilevelModel <: AbstractBilevelModel end +struct UpperToLowerModel <: BridgeBilevelModel + m::BilevelModel +end +UpperToLower(m::BilevelModel) = UpperToLowerModel(m) +struct LowerToUpperModel <: BridgeBilevelModel + m::BilevelModel +end +LowerToUpper(m::BilevelModel) = LowerToUpperModel(m) +bilevel_model(m::BridgeBilevelModel) = m.m + +function set_link!(m::UpperToLowerModel, upper::JuMP.AbstractVariableRef, lower::JuMP.AbstractVariableRef) + bilevel_model(m).upper_to_lower_link[upper] = lower + bilevel_model(m).link[upper] = lower + nothing +end +function set_link!(m::LowerToUpperModel, upper::JuMP.AbstractVariableRef, lower::JuMP.AbstractVariableRef) + bilevel_model(m).lower_to_upper_link[lower] = upper + bilevel_model(m).link[upper] = lower + nothing +end + +#### Model #### + +# Variables +struct BilevelVariableRef <: JuMP.AbstractVariableRef + model::BilevelModel # `model` owning the variable + idx::Int # Index in `model.variables` + level::Level +end +mylevel(v::BilevelVariableRef) = v.level +function solver_ref(v::BilevelVariableRef) + m = v.model + if mylevel(v) == LOWER + return m.sblm_to_solver[ + m.lower_to_sblm[JuMP.index(m.var_lower[v.idx])]] + else + return m.sblm_to_solver[ + m.upper_to_sblm[JuMP.index(m.var_upper[v.idx])]] + end +end +Base.broadcastable(v::BilevelVariableRef) = Ref(v) +Base.copy(v::BilevelVariableRef) = v +Base.:(==)(v::BilevelVariableRef, w::BilevelVariableRef) = + v.model === w.model && v.idx == w.idx && v.level == w.level +JuMP.owner_model(v::BilevelVariableRef) = v.model +JuMP.isequal_canonical(v::BilevelVariableRef, w::BilevelVariableRef) = v == w +JuMP.variable_type(::AbstractBilevelModel) = BilevelVariableRef +# add in BOTH levels +function JuMP.add_variable(bb::BridgeBilevelModel, v::JuMP.AbstractVariable, name::String="") + m = bilevel_model(bb) + m.nextvaridx += 1 + vref = BilevelVariableRef(m, m.nextvaridx, BOTH) + v_upper = JuMP.add_variable(m.upper, v, name) + m.var_upper[vref.idx] = v_upper + v_lower = JuMP.add_variable(m.lower, v, name) + m.var_lower[vref.idx] = v_lower + m.var_level[vref.idx] = BOTH + set_link!(bb, v_upper, v_lower) + m.variables[vref.idx] = v + JuMP.set_name(vref, name) + vref +end +function JuMP.add_variable(inner::InnerBilevelModel, v::JuMP.AbstractVariable, name::String="") + m = bilevel_model(inner) + m.nextvaridx += 1 + vref = BilevelVariableRef(m, m.nextvaridx, level(inner)) + v_level = JuMP.add_variable(mylevel_model(inner), v, name) + mylevel_var_list(inner)[vref.idx] = v_level + m.var_level[vref.idx] = level(inner) + m.variables[vref.idx] = v + JuMP.set_name(vref, name) + vref +end +function MOI.delete!(m::AbstractBilevelModel, vref::BilevelVariableRef) + error("No deletion on bilevel models") + delete!(m.variables, vref.idx) + delete!(m.varnames, vref.idx) +end +MOI.is_valid(m::BilevelModel, vref::BilevelVariableRef) = vref.idx in keys(m.variables) +JuMP.num_variables(m::BilevelModel) = length(m.variables) +JuMP.num_variables(m::InnerBilevelModel) = JuMP.num_variables(bilevel_model(m)) + +# ------------------------- +# begin +# Unchanged from StructJuMP +# ------------------------- + +# Internal function +variable_info(vref::BilevelVariableRef) = vref.model.variables[vref.idx].info +function update_variable_info(vref::BilevelVariableRef, info::JuMP.VariableInfo) + vref.model.variables[vref.idx] = JuMP.ScalarVariable(info) +end + +JuMP.has_lower_bound(vref::BilevelVariableRef) = variable_info(vref).has_lb +function JuMP.lower_bound(vref::BilevelVariableRef) + @assert !JuMP.is_fixed(vref) + variable_info(vref).lower_bound +end +function JuMP.set_lower_bound(vref::BilevelVariableRef, lower) + info = variable_info(vref) + update_variable_info(vref, + JuMP.VariableInfo(true, lower, + info.has_ub, info.upper_bound, + info.has_fix, info.fixed_value, + info.has_start, info.start, + info.binary, info.integer)) +end +function JuMP.delete_lower_bound(vref::BilevelVariableRef) + info = variable_info(vref) + update_variable_info(vref, + JuMP.VariableInfo(false, info.lower_bound, + info.has_ub, info.upper_bound, + info.has_fix, info.fixed_value, + info.has_start, info.start, + info.binary, info.integer)) +end +JuMP.has_upper_bound(vref::BilevelVariableRef) = variable_info(vref).has_ub +function JuMP.upper_bound(vref::BilevelVariableRef) + @assert !JuMP.is_fixed(vref) + variable_info(vref).upper_bound +end +function JuMP.set_upper_bound(vref::BilevelVariableRef, upper) + info = variable_info(vref) + update_variable_info(vref, + JuMP.VariableInfo(info.has_lb, info.lower_bound, + true, upper, + info.has_fix, info.fixed_value, + info.has_start, info.start, + info.binary, info.integer)) +end +function JuMP.delete_upper_bound(vref::BilevelVariableRef) + info = variable_info(vref) + update_variable_info(vref, + JuMP.VariableInfo(info.has_lb, info.lower_bound, + false, info.upper_bound, + info.has_fix, info.fixed_value, + info.has_start, info.start, + info.binary, info.integer)) +end +JuMP.is_fixed(vref::BilevelVariableRef) = variable_info(vref).has_fix +JuMP.fix_value(vref::BilevelVariableRef) = variable_info(vref).fixed_value +function JuMP.fix(vref::BilevelVariableRef, value) + info = variable_info(vref) + update_variable_info(vref, + JuMP.VariableInfo(info.has_lb, info.lower_bound, + info.has_ub, info.upper_bound, + true, value, + info.has_start, info.start, + info.binary, info.integer)) +end +function JuMP.unfix(vref::BilevelVariableRef) + info = variable_info(vref) + update_variable_info(vref, + JuMP.VariableInfo(info.has_lb, info.lower_bound, + info.has_ub, info.upper_bound, + false, info.fixed_value, + info.has_start, info.start, + info.binary, info.integer)) +end +JuMP.start_value(vref::BilevelVariableRef) = variable_info(vref).start +function JuMP.set_start_value(vref::BilevelVariableRef, start) + info = variable_info(vref) + update_variable_info(vref, + JuMP.VariableInfo(info.has_lb, info.lower_bound, + info.has_ub, info.upper_bound, + info.has_fix, info.fixed_value, + true, start, + info.binary, info.integer)) +end +JuMP.is_binary(vref::BilevelVariableRef) = variable_info(vref).binary +function JuMP.set_binary(vref::BilevelVariableRef) + @assert !JuMP.is_integer(vref) + info = variable_info(vref) + update_variable_info(vref, + JuMP.VariableInfo(info.has_lb, info.lower_bound, + info.has_ub, info.upper_bound, + info.has_fix, info.fixed_value, + info.has_start, info.start, + true, info.integer)) +end +function JuMP.unset_binary(vref::BilevelVariableRef) + info = variable_info(vref) + update_variable_info(vref, + JuMP.VariableInfo(info.has_lb, info.lower_bound, + info.has_ub, info.upper_bound, + info.has_fix, info.fixed_value, + info.has_start, info.start, + false, info.integer)) +end +JuMP.is_integer(vref::BilevelVariableRef) = variable_info(vref).integer +function JuMP.set_integer(vref::BilevelVariableRef) + @assert !JuMP.is_binary(vref) + info = variable_info(vref) + update_variable_info(vref, + JuMP.VariableInfo(info.has_lb, info.lower_bound, + info.has_ub, info.upper_bound, + info.has_fix, info.fixed_value, + info.has_start, info.start, + info.binary, true)) +end +function JuMP.unset_integer(vref::BilevelVariableRef) + info = variable_info(vref) + update_variable_info(vref, + JuMP.VariableInfo(info.has_lb, info.lower_bound, + info.has_ub, info.upper_bound, + info.has_fix, info.fixed_value, + info.has_start, info.start, + info.binary, false)) +end + +# ------------------------- +# end +# Unchanged from StructJuMP +# ------------------------- + + +# Constraints +struct BilevelConstraintRef + model::BilevelModel # `model` owning the constraint + idx::Int # Index in `model.constraints` +end +JuMP.constraint_type(::AbstractBilevelModel) = BilevelConstraintRef +function JuMP.add_constraint(m::BilevelModel, c::JuMP.AbstractConstraint, name::String="") + error( + "Can't add constraint directly to the bilevel model `m`, "* + "attach the constraint to the upper or lower model "* + "with @constraint(Upper(m), ...) or @constraint(Lower(m), ...)") +end +# function constraint_object(con_ref::ConstraintRef{Model, _MOICON{FuncType, SetType}}) where +# {FuncType <: MOI.AbstractScalarFunction, SetType <: MOI.AbstractScalarSet} +# model = con_ref.model +# f = MOI.get(model, MOI.ConstraintFunction(), con_ref)::FuncType +# s = MOI.get(model, MOI.ConstraintSet(), con_ref)::SetType +# return ScalarConstraint(jump_function(model, f), s) +# end +JuMP.add_constraint(m::UpperModel, c::JuMP.VectorConstraint, name::String="") = +error("no vec ctr") +function JuMP.add_constraint(m::InnerBilevelModel, c::JuMP.ScalarConstraint{F,S}, name::String="") where {F,S} + blm = bilevel_model(m) + blm.nextconidx += 1 + cref = BilevelConstraintRef(blm, blm.nextconidx) + func = JuMP.jump_function(c) + level_func = replace_variables(func, bilevel_model(m), mylevel_model(m), mylevel_var_list(m)) + level_c = JuMP.build_constraint(error, level_func, c.set) + level_cref = JuMP.add_constraint(mylevel_model(m), level_c, name) + blm.ctr_level[cref.idx] = level(m) + mylevel_ctr_list(m)[cref.idx] = level_cref + blm.constraints[cref.idx] = c + JuMP.set_name(cref, name) + cref +end +function MOI.delete!(m::AbstractBilevelModel, cref::BilevelConstraintRef) + error("can't delete") + delete!(m.constraints, cref.idx) + delete!(m.connames, cref.idx) +end +MOI.is_valid(m::BilevelModel, cref::BilevelConstraintRef) = cref.idx in keys(m.constraints) +MOI.is_valid(m::InnerBilevelModel, cref::BilevelConstraintRef) = + MOI.is_valid(bilevel_model(m), cref) && bilevel_model(m).ctr_level[cref.idx] == level(m) +function JuMP.constraint_object(cref::BilevelConstraintRef, F::Type, S::Type) + c = cref.model.constraints[cref.idx] + # `TypeError` should be thrown is `F` and `S` are not correct + # This is needed for the tests in `constraints.jl` + c.func::F + c.set::S + c +end + +# Objective +function JuMP.set_objective(m::InnerBilevelModel, sense::MOI.OptimizationSense, + f::JuMP.AbstractJuMPScalar) + set_mylevel_obj_sense(m, sense) + set_mylevel_obj_function(m, f) + level_f = replace_variables(f, bilevel_model(m), mylevel_model(m), mylevel_var_list(m)) + JuMP.set_objective(mylevel_model(m), sense, level_f) +end +JuMP.objective_sense(m::InnerBilevelModel) = mylevel_obj_sense(m) +JuMP.objective_function_type(m::InnerBilevelModel) = typeof(mylevel_obj_function(m)) +JuMP.objective_function(m::InnerBilevelModel) = mylevel_obj_function(m) +function JuMP.objective_function(m::InnerBilevelModel, FT::Type) + mylevel_obj_function(m) isa FT || error("The objective function is not of type $FT") + mylevel_obj_function(m) +end + +# todo remove +JuMP.objective_sense(m::AbstractBilevelModel) = MOI.FEASIBILITY_SENSE +# end todo remove +JuMP.num_variables(m::AbstractBilevelModel) = JuMP.num_variables(bilevel_model(m)) +JuMP.show_constraints_summary(::Any, ::AbstractBilevelModel) = "no summary" +JuMP.show_backend_summary(::Any, ::AbstractBilevelModel) = "no summary" +JuMP.object_dictionary(m::BilevelModel) = m.objdict +JuMP.object_dictionary(m::AbstractBilevelModel) = JuMP.object_dictionary(bilevel_model(m)) +JuMP.show_objective_function_summary(::IO, ::AbstractBilevelModel) = "no summary" + +bileve_obj_error() = error("There is no objective for BilevelModel use Upper(.) and Lower(.)") + +function JuMP.set_objective(m::BilevelModel, sense::MOI.OptimizationSense, + f::JuMP.AbstractJuMPScalar) + bileve_obj_error() +end +JuMP.objective_sense(m::BilevelModel) = JuMP.objective_sense(m.upper)#bileve_obj_error() +JuMP.objective_function_type(model::BilevelModel) = bileve_obj_error() +JuMP.objective_function(model::BilevelModel) = bileve_obj_error() +function JuMP.objective_function(model::BilevelModel, FT::Type) + bileve_obj_error() +end + +# Names +JuMP.name(vref::BilevelVariableRef) = vref.model.varnames[vref.idx] +function JuMP.set_name(vref::BilevelVariableRef, name::String) + vref.model.varnames[vref.idx] = name +end +JuMP.name(cref::BilevelConstraintRef) = cref.model.connames[cref.idx] +function JuMP.set_name(cref::BilevelConstraintRef, name::String) + cref.model.connames[cref.idx] = name +end + + +# replace variables +function replace_variables(var::BilevelVariableRef, + model::BilevelModel, + inner::JuMP.AbstractModel, + variable_map::Dict{Int, V}) where {V<:JuMP.AbstractVariableRef} + if var.model === model + return variable_map[var.idx] + else + error("A BilevelModel cannot have expression using variables of a BilevelModel different from itself") + end +end +function replace_variables(aff::JuMP.GenericAffExpr{C, BilevelVariableRef}, + model::BilevelModel, + inner::JuMP.AbstractModel, + variable_map::Dict{Int, V}) where {C,V<:JuMP.AbstractVariableRef} + result = JuMP.GenericAffExpr{C, JuMP.VariableRef}(0.0)#zero(aff) + result.constant = aff.constant + for (coef, var) in JuMP.linear_terms(aff) + JuMP.add_to_expression!(result, + coef, + replace_variables(var, model, model, variable_map)) + end + return result +end +# function replace_variables(quad::JuMP.GenericQuadExpr{C, BilevelVariableRef}, +# model::BilevelModel, +# inner::JuMP.AbstractModel, +# variable_map::Dict{Int, V}) where {C,V<:JuMP.AbstractVariableRef} +# error("A BilevelModel cannot have quadratic function") +# end +function replace_variables(quad::C, + model::BilevelModel, + inner::JuMP.AbstractModel, + variable_map::Dict{Int, V}) where {C,V<:JuMP.AbstractVariableRef} + error("A BilevelModel cannot have $(C) function") +end +replace_variables(funcs::Vector, args...) = map(f -> replace_variables(f, args...), funcs) +using MathOptFormat +function print_lp(m, name) + lp_model = MathOptFormat.MOF.Model() + MOI.copy_to(lp_model, m) + MOI.write_to_file(lp_model, name) +end + +JuMP.optimize!(::T) where {T<:AbstractBilevelModel} = + error("cant solve a model of type: $T ") +function JuMP.optimize!(model::BilevelModel, optimizer) + + upper = JuMP.backend(model.upper) + lower = JuMP.backend(model.lower) + + # print_lp(upper, "upper.lp") + # print_lp(lower, "lower.lp") + + moi_upper = JuMP.index.( + collect(values(model.upper_to_lower_link))) + moi_link = JuMP.index(model.link) + + single_blm, upper_to_sblm, lower_to_sblm = build_bilivel(upper, lower, moi_link, moi_upper) + + solver = MOI.Bridges.full_bridge_optimizer(optimizer, Float64) + # MOI.empty!(solver) + # MOI.is_empty(solver) + sblm_to_solver = MOI.copy_to(solver, single_blm) + + MOI.optimize!(solver) + + # map from bridged to single_blm + # map from single_blm to upper & lowel + # map from upper & lower to model + model.solver = solver#::MOI.ModelLike + model.upper_to_sblm = upper_to_sblm + model.lower_to_sblm = lower_to_sblm + model.sblm_to_solver = sblm_to_solver + + # feasible = JuMP.primal_status(single_blm)# == MOI.FEASIBLE_POINT + # termination = JuMP.termination_status(single_blm)# == MOI.OPTIMAL + # objective_value = NaN + # if feasible == MOI.FEASIBLE_POINT + # objective_value = JuMP.objective_value(single_blm) + # end + # variable_value = Dict{JuMP.VariableRef, Float64}() + return nothing +end + +function JuMP.index(d::Dict) + ret = Dict{VI,VI}() + # sizehint!(ret, length(d)) + for (k,v) in d + ret[JuMP.index(k)] = JuMP.index(v) + end + return ret +end + +function JuMP.value(v::BilevelVariableRef)::Float64 + m = owner_model(v) + solver = m.solver + ref = solver_ref(v) + return MOI.get(solver, MOI.VariablePrimal(), ref) +end +function JuMP.value(v::BilevelConstraintRef)::Float64 + error("value of BilevelConstraintRef not enabled") +end +function JuMP.dual(v::BilevelConstraintRef)::Float64 + error("dual of BilevelConstraintRef not enabled") +end +function JuMP.primal_status(model::BilevelModel) + return MOI.get(model.solver, MOI.PrimalStatus()) +end +JuMP.dual_status(model::BilevelModel) = error("dual status cant be queried for BilevelModel") +function JuMP.termination_status(model::BilevelModel) + return MOI.get(model.solver, MOI.TerminationStatus()) +end +function JuMP.objective_value(model::BilevelModel) + return MOI.get(model.solver, MOI.ObjectiveValue()) +end \ No newline at end of file diff --git a/src/moi.jl b/src/moi.jl new file mode 100644 index 00000000..8e1a4f02 --- /dev/null +++ b/src/moi.jl @@ -0,0 +1,188 @@ +const SVF = MOI.SingleVariable +const VVF = MOI.VectorOfVariables +const SAF{T} = MOI.ScalarAffineFunction{T} +const VAF{T} = MOI.VectorAffineFunction{T} + +const VI = MOI.VariableIndex +const CI = MOI.ConstraintIndex + +struct Complement#{M1 <: MOI.ModelLike, M2 <: MOI.ModelLike, F, S} + # primal::M1 + func#::F + set#::S + # dual::M2 + variable#::VI +end + +function get_canonical_complements(primal_model, primal_dual_map) + map = primal_dual_map.primal_con_dual_var + out = Complement[] + T = Float64 + for ci in keys(map) + func = MOI.get(primal_model, MOI.ConstraintFunction(), ci) + set = MOI.get(primal_model, MOI.ConstraintSet(), ci) + i = 1 + func.constant = Dualization.set_dot(i, set, T)*Dualization.get_scalar_term(primal_model, i, ci) + # todo - set dot on function + con = Complement(func, set_with_zero(set), map[ci][1]) + push!(out, con) + end + return out +end + +function set_with_zero(set::Union{MOI.LessThan{T}, + MOI.GreaterThan{T}, MOI.EqualTo{T}}) where T + return typeof(set)(0.0) +end +function set_with_zero(set) + return copy(set) +end + +function build_bilivel(upper::MOI.ModelLike, lower::MOI.ModelLike, + link::Dict{VI,VI}, upper_variables::Vector{VI}) + + # Start with an empty problem + mode = MOIU.AUTOMATIC + m = MOIU.CachingOptimizer(MOIU.Model{Float64}(), mode) + + # add the first level + copy_names = false + # key are from src, value are from dest + upper_idxmap = MOIU.default_copy_to(m, upper, copy_names) + pass_names(m, upper, upper_idxmap) + + # append the second level primal + lower_idxmap = MOIU.IndexMap() # + + for i in keys(upper_idxmap.varmap) + lower_idxmap[link[i]] = upper_idxmap[i] + end + + append_to(m, lower, lower_idxmap, copy_names) + pass_names(m, lower, lower_idxmap) + + + # dualize the second level + dual_problem = dualize(lower, + dual_names = DualNames("dual_","dual_"), + variable_parameters = upper_variables, + ignore_objective = true) + lower_dual = dual_problem.dual_model + lower_primal_dual_map = dual_problem.primal_dual_map + + # appende the second level dual + lower_dual_idxmap = MOIU.IndexMap() + + append_to(m, lower_dual, lower_dual_idxmap, copy_names) + pass_names(m, lower_dual, lower_dual_idxmap) + + # complete KKT + # 1 - primal dual equality (quadratic equality constraint) + # 1a - no slacks + # 1b - use slack + # 2 - complementarity + # 2a - actual complementarity constraints + # 2b - SOS1 (slack * duals) + # 2c - Big-M formulation (Fortuny-Amat and McCarl, 1981) + # 3 - NLP (y = argmin_y problem) + + if true # SOS1 + # futurely add slacks to primal model + comps = get_canonical_complements(lower, lower_primal_dual_map) + for comp in comps + # create a constrained slack - depends on set + add_complement_constraint(m, comp, lower_idxmap, lower_dual_idxmap) + # add directly to the primal equality + # (will require set change - loses dual mapping) + end + + end + + return m, upper_idxmap, lower_idxmap#, lower_primal_dual_map, lower_dual_idxmap +end + +function add_complement_constraint(m, comp::Complement, idxmap_primal, idxmap_dual) + f = comp.func + s = comp.set + v = comp.variable + T = Float64 + slack, slack_in_set = MOI.add_constrained_variable(m, s) + f_dest = MOIU.map_indices.(Ref(idxmap_primal), f) + new_f = MOIU.operate(-, T, f_dest, MOI.SingleVariable(slack)) + equality = MOIU.normalize_and_add_constraint(m, new_f, MOI.EqualTo(zero(T))) + + dual = idxmap_dual[v] + c1 = MOI.add_constraint(m, + MOI.VectorOfVariables([slack, dual]), + MOI.SOS1([1.0, 2.0])) + + nm = MOI.get(m, MOI.VariableName(), dual) + MOI.set(m, MOI.VariableName(), slack, "slk($(nm))") + MOI.set(m, MOI.ConstraintName(), c1, "c_slk($(nm))") + MOI.set(m, MOI.ConstraintName(), slack_in_set, "b_slk($(nm))") + + return slack, slack_in_set, equality, c1 +end + +function pass_names(dest, src, map) + for vi in MOI.get(src, MOI.ListOfVariableIndices()) + name = MOI.get(src, MOI.VariableName(), vi) + MOI.set(dest, MOI.VariableName(), map[vi], name) + end + for (F,S) in MOI.get(src, MOI.ListOfConstraints()) + for con in MOI.get(src, MOI.ListOfConstraintIndices{F,S}()) + name = MOI.get(src, MOI.ConstraintName(), con) + MOI.set(dest, MOI.ConstraintName(), map[con], name) + end + end +end + +function append_to(dest::MOI.ModelLike, src::MOI.ModelLike, idxmap, copy_names::Bool) + # MOI.empty!(dest) + + # idxmap = MOIU.IndexMap() + + vis_src = MOI.get(src, MOI.ListOfVariableIndices()) + constraint_types = MOI.get(src, MOI.ListOfConstraints()) + single_variable_types = [S for (F, S) in constraint_types + if F == MOI.SingleVariable] + vector_of_variables_types = [S for (F, S) in constraint_types + if F == MOI.VectorOfVariables] + + # The `NLPBlock` assumes that the order of variables does not change (#849) + if MOI.NLPBlock() in MOI.get(src, MOI.ListOfModelAttributesSet()) + error("not nlp for now") + vector_of_variables_not_added = [ + MOI.get(src, MOI.ListOfConstraintIndices{MOI.VectorOfVariables, S}()) + for S in vector_of_variables_types + ] + single_variable_not_added = [ + MOI.get(src, MOI.ListOfConstraintIndices{MOI.SingleVariable, S}()) + for S in single_variable_types + ] + else + vector_of_variables_not_added = [ + MOIU.copy_vector_of_variables(dest, src, idxmap, S) + for S in vector_of_variables_types + ] + single_variable_not_added = [ + MOIU.copy_single_variable(dest, src, idxmap, S) + for S in single_variable_types + ] + end + + MOIU.copy_free_variables(dest, idxmap, vis_src, MOI.add_variables) + + # Copy variable attributes + # MOIU.pass_attributes(dest, src, copy_names, idxmap, vis_src) + + # Copy model attributes + # pass_attributes(dest, src, copy_names, idxmap) + + # Copy constraints + MOIU.pass_constraints(dest, src, copy_names, idxmap, + single_variable_types, single_variable_not_added, + vector_of_variables_types, vector_of_variables_not_added) + + return idxmap +end \ No newline at end of file diff --git a/src/other.jl b/src/other.jl new file mode 100644 index 00000000..e132da98 --- /dev/null +++ b/src/other.jl @@ -0,0 +1,56 @@ + +function split_variables(saf::MOI.ScalarAffineFunction{T}, + variable_parameters::Vector{VI}) where T + remove = get_indices_variables(saf, variable_parameters) + new_saf = MOI.ScalarAffineFunction{T}(saf.terms[remove], 0.0) + old_saf = copy(saf) + deleteat!(old_saf, remove) + return old_saf, new_saf +end + +function get_parametric_constant(saf::MOI.ScalarAffineFunction{T}, + variable_parameters::Vector{VI}) where T + terms = get_indices_variables(saf, variable_parameters) + MOI.ScalarAffineFunction{T}(saf.terms[terms], 0.0) +end + +function get_parametric_constants(primal_model::MOI.ModelLike, + # primal_dual_map::PrimalDualMap{T}, dual_names::DualNames, + con_types::Vector{Tuple{DataType, DataType}}, variable_parameters::Vector{VI}) where T + # todo + param_functions = Dict() + for (F, S) in con_types + for ci in MOI.get(primal_model, MOI.ListOfConstraintIndices{F,S}()) # Constraints of type {F, S} + func = MOI.get(primal_model, MOI.ConstraintFunction(), ci) + set = MOI.get(primal_model, MOI.ConstraintSet(), ci) + i = 1 # for scalar affine + val = set_dot(i, set, T)*get_scalar_term(primal_model, i, ci) + + # todo - requires a setdot here + param_functions[ci] = get_parametric_constant(func, variable_parameters) + + + param_functions[ci].constant = val + end + end + return param_functions +end + +function get_canonical_functions(primal_model::MOI.ModelLike) + # todo + T = Float64 + con_types = MOI.get(primal_model, MOI.ListOfConstraints()) + functions = Dict() + for (F, S) in con_types + for ci in MOI.get(primal_model, MOI.ListOfConstraintIndices{F,S}()) # Constraints of type {F, S} + func = MOI.get(primal_model, MOI.ConstraintFunction(), ci) + set = MOI.get(primal_model, MOI.ConstraintSet(), ci) + i = 1 # for scalar affine + func.constant = set_dot(i, set, T)*get_scalar_term(primal_model, i, ci) + + # todo - requires a setdot here + functions[ci] = func + end + end + return functions +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 00000000..38f6f97a --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,343 @@ +using BilevelJuMP +using Test, Gurobi, MathOptInterface, JuMP, Dualization +using MathOptFormat + +const MOI = MathOptInterface +const MOIU = MathOptInterface.Utilities +const MOIB = MathOptInterface.Bridges +const MOIT = MathOptInterface.Test + +# original problem +# min -4x -3y +# s.t. +# 2x + y <= 4 +# x +2y <= 4 +# x >= 0 +# y >= 0 +# +# sol x = y = 4/3 +# obj = 28/3 + +# bilevel modifications +# min -4x -3y +# s.t. +# y = argmin_y y +# 2x + y <= 4 +# x +2y <= 4 +# x >= 0 +# y >= 0 +# +# sol: x = 2, y = 0 +# obj_upper = -8 +# obj_lower = 0 + +@testset "Simple LP" begin + optimizer = Gurobi.Optimizer() + bridged = MOIB.full_bridge_optimizer(optimizer, Float64) + MOI.empty!(bridged) + @test MOI.is_empty(bridged) + + # add 10 variables - only diagonal is relevant + X = MOI.add_variables(bridged, 2) + + c1 = MOIU.normalize_and_add_constraint(bridged, + MOI.ScalarAffineFunction([ + MOI.ScalarAffineTerm(2.0, X[1]), + MOI.ScalarAffineTerm(1.0, X[2]) + ], 0.0), MOI.EqualTo(4.0)) + + c2 = MOIU.normalize_and_add_constraint(bridged, + MOI.ScalarAffineFunction([ + MOI.ScalarAffineTerm(1.0, X[1]), + MOI.ScalarAffineTerm(2.0, X[2]) + ], 0.0), MOI.EqualTo(4.0)) + + b1 = MOIU.normalize_and_add_constraint(bridged, + MOI.ScalarAffineFunction([ + MOI.ScalarAffineTerm(1.0, X[1]) + ], 0.0), MOI.GreaterThan(0.0)) + + b2 = MOIU.normalize_and_add_constraint(bridged, + MOI.ScalarAffineFunction([ + MOI.ScalarAffineTerm(1.0, X[2]) + ], 0.0), MOI.GreaterThan(0.0)) + + MOI.set(bridged, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], [X[1], X[2]]), 0.0) + ) + MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.optimize!(bridged) + + obj = MOI.get(bridged, MOI.ObjectiveValue()) + + @test obj ≈ -9.33333 atol = 1e-2 + + Xr = MOI.get(bridged, MOI.VariablePrimal(), X) + + @test Xr ≈ [1.3333, 1.3333] atol = 1e-2 + +end + +@testset "Simple BLP" begin + optimizer = Gurobi.Optimizer() + bridged = MOIB.full_bridge_optimizer(optimizer, Float64) + MOI.empty!(bridged) + @test MOI.is_empty(bridged) + + X = MOI.add_variables(bridged, 2) + + MOI.set(bridged, MOI.VariableName(), X[1], "x_u") + MOI.set(bridged, MOI.VariableName(), X[2], "y_u") + + MOI.set(bridged, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + MOI.ScalarAffineFunction( + MOI.ScalarAffineTerm.([-4.0, -3.0], [X[1], X[2]]), 0.0) + ) + MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE) + + optimizer2 = Gurobi.Optimizer() + bridged2 = MOIB.full_bridge_optimizer(optimizer2, Float64) + MOI.empty!(bridged2) + @test MOI.is_empty(bridged2) + + X2 = MOI.add_variables(bridged2, 2) + MOI.set(bridged2, MOI.VariableName(), X2[1], "x_l") + MOI.set(bridged2, MOI.VariableName(), X2[2], "y_l") + + + c1 = MOIU.normalize_and_add_constraint(bridged2, + MOI.ScalarAffineFunction([ + MOI.ScalarAffineTerm(2.0, X2[1]), + MOI.ScalarAffineTerm(1.0, X2[2]) + ], 0.0), MOI.LessThan(4.0)) + MOI.set(bridged2, MOI.ConstraintName(), c1, "l_c1") + + c2 = MOIU.normalize_and_add_constraint(bridged2, + MOI.ScalarAffineFunction([ + MOI.ScalarAffineTerm(1.0, X2[1]), + MOI.ScalarAffineTerm(2.0, X2[2]) + ], 0.0), MOI.LessThan(4.0)) + MOI.set(bridged2, MOI.ConstraintName(), c2, "l_c2") + + b1 = MOIU.normalize_and_add_constraint(bridged2, + MOI.ScalarAffineFunction([ + MOI.ScalarAffineTerm(1.0, X2[1]) + ], 0.0), MOI.GreaterThan(0.0)) + MOI.set(bridged2, MOI.ConstraintName(), b1, "l_b1") + + + b2 = MOIU.normalize_and_add_constraint(bridged2, + MOI.ScalarAffineFunction([ + MOI.ScalarAffineTerm(1.0, X2[2]) + ], 0.0), MOI.GreaterThan(0.0)) + MOI.set(bridged2, MOI.ConstraintName(), b2, "l_b12") + + MOI.set(bridged2, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + MOI.ScalarAffineFunction( + MOI.ScalarAffineTerm.([1.0], [X2[2]]), 0.0) + ) + MOI.set(bridged2, MOI.ObjectiveSense(), MOI.MIN_SENSE) + + + links = Dict(X[1] => X2[1], X[2] => X2[2]) + + parametric = [X2[1]] + + blp, _, _ = BilevelJuMP.build_bilivel(bridged, bridged2, links, parametric) + + + # MOI.optimize!(bridged) + + # obj = MOI.get(bridged, MOI.ObjectiveValue()) + + # @test obj ≈ -9.33333 atol = 1e-2 + + # Xr = MOI.get(bridged, MOI.VariablePrimal(), X) + + # @test Xr ≈ [1.3333, 1.3333] atol = 1e-2 + + # lp_model = MathOptFormat.MOF.Model() + # MOI.copy_to(lp_model, blp) + # MOI.write_to_file(lp_model, "my_model.LP") + # @show pwd() + + optimizer = Gurobi.Optimizer() + bridged3 = MOIB.full_bridge_optimizer(optimizer, Float64) + MOI.empty!(bridged3) + @test MOI.is_empty(bridged3) + MOI.copy_to(bridged3, blp) + + MOI.optimize!(bridged3) + +end + +@testset "Simple BLP2" begin + optimizer = Gurobi.Optimizer() + bridged = MOIB.full_bridge_optimizer(optimizer, Float64) + MOI.empty!(bridged) + @test MOI.is_empty(bridged) + + x = MOI.add_variable(bridged) + y = MOI.add_variable(bridged) + + MOI.set(bridged, MOI.VariableName(), x, "x_u") + MOI.set(bridged, MOI.VariableName(), y, "y_u") + + MOI.set(bridged, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + MOI.ScalarAffineFunction( + MOI.ScalarAffineTerm.([1.0, 3.0], [x, y]), 0.0) + ) + MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE) + + optimizer2 = Gurobi.Optimizer() + bridged2 = MOIB.full_bridge_optimizer(optimizer2, Float64) + MOI.empty!(bridged2) + @test MOI.is_empty(bridged2) + + x2 = MOI.add_variable(bridged2) + y2 = MOI.add_variable(bridged2) + + MOI.set(bridged2, MOI.VariableName(), x2, "x_l") + MOI.set(bridged2, MOI.VariableName(), y2, "y_l") + + + c1 = MOIU.normalize_and_add_constraint(bridged2, + MOI.ScalarAffineFunction([ + MOI.ScalarAffineTerm(1.0, x2), + MOI.ScalarAffineTerm(1.0, y2) + ], 0.0), MOI.LessThan(8.0)) + MOI.set(bridged2, MOI.ConstraintName(), c1, "l_c1") + + c2 = MOIU.normalize_and_add_constraint(bridged2, + MOI.ScalarAffineFunction([ + MOI.ScalarAffineTerm(1.0, x2), + MOI.ScalarAffineTerm(4.0, y2) + ], 0.0), MOI.GreaterThan(8.0)) + MOI.set(bridged2, MOI.ConstraintName(), c2, "l_c2") + + c3 = MOIU.normalize_and_add_constraint(bridged2, + MOI.ScalarAffineFunction([ + MOI.ScalarAffineTerm(1.0, x2), + MOI.ScalarAffineTerm(2.0, y2) + ], 0.0), MOI.LessThan(13.0)) + MOI.set(bridged2, MOI.ConstraintName(), c3, "l_c3") + + b1 = MOIU.normalize_and_add_constraint(bridged2, + MOI.ScalarAffineFunction([ + MOI.ScalarAffineTerm(1.0, x2) + ], 0.0), MOI.GreaterThan(1.0)) + MOI.set(bridged2, MOI.ConstraintName(), b1, "l_b1") + + + b2 = MOIU.normalize_and_add_constraint(bridged2, + MOI.ScalarAffineFunction([ + MOI.ScalarAffineTerm(1.0, x2) + ], 0.0), MOI.LessThan(6.0)) + MOI.set(bridged2, MOI.ConstraintName(), b2, "l_b12") + + MOI.set(bridged2, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + MOI.ScalarAffineFunction( + MOI.ScalarAffineTerm.([-1.0], [y2]), 0.0) + ) + MOI.set(bridged2, MOI.ObjectiveSense(), MOI.MIN_SENSE) + + + links = Dict(x => x2, y => y2) + + parametric = [x2] + + blp, _, _ = BilevelJuMP.build_bilivel(bridged, bridged2, links, parametric) + + + # MOI.optimize!(bridged) + + # obj = MOI.get(bridged, MOI.ObjectiveValue()) + + # @test obj ≈ -9.33333 atol = 1e-2 + + # Xr = MOI.get(bridged, MOI.VariablePrimal(), X) + + # @test Xr ≈ [1.3333, 1.3333] atol = 1e-2 + + # lp_model = MathOptFormat.MOF.Model() + # MOI.copy_to(lp_model, blp) + # MOI.write_to_file(lp_model, "my_model.LP") + # @show pwd() + + optimizer = Gurobi.Optimizer() + bridged3 = MOIB.full_bridge_optimizer(optimizer, Float64) + MOI.empty!(bridged3) + @test MOI.is_empty(bridged3) + MOI.copy_to(bridged3, blp) + + MOI.optimize!(bridged3) + +end + +@testset "Simple BLP JuMP" begin + + model = BilevelModel() + + @variable(UpperToLower(model), x) + @variable(LowerToUpper(model), y) + + @objective(Upper(model), Min, -4x -3y) + + @objective(Lower(model), Min, y) + + @constraint(Lower(model), 2x+y <= 4) + @constraint(Lower(model), x+2y <= 4) + @constraint(Lower(model), x >= 0) + @constraint(Lower(model), y >= 0) + + optimize!(model, Gurobi.Optimizer()) + + primal_status(model) + + termination_status(model) + + @test objective_value(model) == -8 + + @test value(x) == 2 + @test value(y) == 0 + +end + +@testset "Simple BLP2 JuMP" begin + + model = BilevelModel() + + @variable(UpperToLower(model), x) + @variable(LowerToUpper(model), y) + + @objective(Upper(model), Min, x + 3y) + + @objective(Lower(model), Min, -y) + + @constraint(Lower(model), x+y <= 8) + @constraint(Lower(model), x+4y >= 8) + @constraint(Lower(model), x+2y <= 13) + @constraint(Lower(model), x >= 1) + @constraint(Lower(model), x <= 6) + @constraint(Lower(model), y >= 0) + + optimize!(model, Gurobi.Optimizer()) + + primal_status(model) + + termination_status(model) + + @test objective_value(model) == 12 + + @test value(x) == 6 + @test value(y) == 2 + +end + +# alternative sintax +# @constraint(model, lower, [y] \in ArgMin([x]))