fix isinteger + continuous columns if sp has only continuous vars (#…

…551)

* fix #550

* correct

src/Algorithm/colgen.jl

@@ -235,19 +235,15 @@ set_bestcol_id!(spinfo::SubprobInfo, varid::VarId) = spinfo.bestcol_id = varid
 function insert_cols_in_master!(
  masterform::Formulation, spinfo::SubprobInfo, phase::Int64, spform::Formulation,
 )
- sp_uid = getuid(spform)
  nb_of_gen_col = 0
 
  for sol_id in spinfo.recorded_sol_ids
  nb_of_gen_col += 1
  name = string("MC_", getsortuid(sol_id))
  lb = 0.0
  ub = Inf
- kind = Continuous
  duty = MasterCol
- mc = setcol_from_sp_primalsol!(
- masterform, spform, sol_id, name, duty; lb=lb, ub=ub, kind=kind
- )
+ mc = setcol_from_sp_primalsol!(masterform, spform, sol_id, name, duty; lb=lb, ub=ub)
  if phase == 1
  setcurcost!(masterform, mc, 0.0)
  end

src/MathProg/duties.jl

@@ -19,6 +19,7 @@ mutable struct DwSp <: AbstractSpDuty
  setup_var::Union{VarId, Nothing}
  lower_multiplicity::Int
  upper_multiplicity::Int
+ column_var_kind::VarKind
 end
 
 "A Benders subproblem of a formulation decomposed using Benders."
@@ -47,14 +48,14 @@ struct BendersSp <: AbstractSpDuty end
  AbstractDwSpVar <= Duty{Variable}
  DwSpPricingVar <= AbstractDwSpVar
  DwSpSetupVar <= AbstractDwSpVar
- DwSpPureVar <= AbstractDwSpVar
+ #DwSpPureVar <= AbstractDwSpVar
  DwSpPrimalSol <= AbstractDwSpVar
  AbstractBendSpVar <= Duty{Variable}
  AbstractBendSpSlackMastVar <= AbstractBendSpVar
  BendSpSlackFirstStageVar <= AbstractBendSpSlackMastVar
  BendSpSlackSecondStageCostVar <= AbstractBendSpSlackMastVar
  BendSpSepVar <= AbstractBendSpVar
- BendSpPureVar <= AbstractBendSpVar
+ #BendSpPureVar <= AbstractBendSpVar
  BendSpPrimalSol <= AbstractBendSpVar
 end
 
@@ -108,11 +109,11 @@ function isaStaticDuty(duty::NestedEnum)
  duty <= MasterRepPricingSetupVar ||
  duty <= DwSpPricingVar ||
  duty <= DwSpSetupVar ||
- duty <= DwSpPureVar ||
+ #duty <= DwSpPureVar ||
  duty <= DwSpPrimalSol ||
  duty <= DwSpDualSol ||
  duty <= BendSpSepVar ||
- duty <= BendSpPureVar ||
+ #duty <= BendSpPureVar ||
  duty <= BendSpSlackFirstStageVar ||
  duty <= BendSpSlackSecondStageCostVar ||
  duty <= OriginalConstr ||

src/MathProg/formulation.jl

@@ -301,7 +301,7 @@ end
 
 function setcol_from_sp_primalsol!(
  masterform::Formulation, spform::Formulation, sol_id::VarId, name::String,
- duty::Duty{Variable}; lb::Float64 = 0.0, ub::Float64 = Inf, kind::VarKind = Continuous,
+ duty::Duty{Variable}; lb::Float64 = 0.0, ub::Float64 = Inf,
  inc_val::Float64 = 0.0, is_active::Bool = true, is_explicit::Bool = true,
  moi_index::MoiVarIndex = MoiVarIndex(), custom_data::Union{Nothing, AbstractCustomData} = nothing
 )
@@ -322,7 +322,7 @@ function setcol_from_sp_primalsol!(
  cost = cost,
  lb = lb,
  ub = ub,
- kind = kind,
+ kind = spform.duty_data.column_var_kind,
  inc_val = inc_val,
  is_active = is_active,
  is_explicit = is_explicit,
@@ -331,6 +331,9 @@ function setcol_from_sp_primalsol!(
  id = sol_id,
  custom_data = get(spform.manager.primal_sols_custom_data, sol_id, nothing)
  )
+
+ setcurkind!(masterform, mast_col, Continuous)
+
  return mast_col
 end
 
@@ -481,29 +484,36 @@ function _addlocalartvar!(form::Formulation, constr::Constraint, abs_cost::Float
  return
 end
 
+"""
+ enforce_integrality!(formulation)
+
+Set the current kind of each active & explicit variable of the formulation to its perennial kind.
+"""
 function enforce_integrality!(form::Formulation)
- @logmsg LogLevel(-1) string("Enforcing integrality of formulation ", getuid(form))
- for (varid, var) in getvars(form)
- !iscuractive(form, varid) && continue
- !isexplicit(form, varid) && continue
- getcurkind(form, varid) == Integ && continue
- getcurkind(form, varid) == Binary && continue
- if getduty(varid) <= MasterCol || getperenkind(form, varid) != Continuous
- @logmsg LogLevel(-3) string("Setting kind of var ", getname(form, var), " to Integer")
- setcurkind!(form, varid, Integ)
+ for (_, var) in getvars(form)
+ enforce = iscuractive(form, var) && isexplicit(form, var)
+ enforce &= getcurkind(form, var) === Continuous
+ enforce &= getperenkind(form, var) !== Continuous
+ if enforce
+ setcurkind!(form, var, getperenkind(form, var))
  end
  end
  return
 end
 
-function relax_integrality!(form::Formulation) # TODO remove : should be in Algorithm
- @logmsg LogLevel(-1) string("Relaxing integrality of formulation ", getuid(form))
- for (varid, var) in getvars(form)
- !iscuractive(form, varid) && continue
- !isexplicit(form, varid) && continue
- getcurkind(form, var) == Continuous && continue
- @logmsg LogLevel(-3) string("Setting kind of var ", getname(form, var), " to continuous")
- setcurkind!(form, varid, Continuous)
+"""
+ relax_integrality!(formulation)
+
+Set the current kind of each active & explicit integer or binary variable of the formulation
+to continuous.
+"""
+function relax_integrality!(form::Formulation)
+ for (_, var) in getvars(form)
+ relax = iscuractive(form, var) && isexplicit(form, var)
+ relax &= getcurkind(form, var) !== Continuous
+ if relax
+ setcurkind!(form, var, Continuous)
+ end
  end
  return
 end

src/MathProg/solutions.jl

@@ -20,15 +20,13 @@ end
 
 function Base.isinteger(sol::Solution)
  for (vc_id, val) in sol
- #if getperenkind(sol.model, vc_id) != Continuous
- abs(round(val) - val) <= 1e-5 || return false
- #end
+ if getperenkind(sol.model, vc_id) !== Continuous && abs(round(val) - val) > 1e-5
+ return false
+ end
  end
  return true
 end
 
-isfractional(sol::Solution) = !Base.isinteger(sol)
-
 function contains(sol::PrimalSolution, f::Function)
  for (varid, val) in sol
  f(varid) && return true

src/decomposition.jl

@@ -60,7 +60,7 @@ function instantiatesp!(
 )
  spform = create_formulation!(
  env,
- MathProg.DwSp(nothing, 1, 1);
+ MathProg.DwSp(nothing, 1, 1, Integ);
  parent_formulation = masterform,
  obj_sense = getobjsense(masterform)
  )
@@ -138,7 +138,9 @@ function create_side_vars_constrs!(
  @assert length(setupvars) == 1
  setupvar = collect(values(setupvars))[1]
  setuprepvar = clonevar!(origform, masterform, spform, setupvar, MasterRepPricingSetupVar, is_explicit = false)
- # create convexity constraint
+
+ # create convexity constraint & storing information about the convexity constraint
+ # in the duty data of the formulation
  lb_mult = Float64(BD.getlowermultiplicity(ann))
  name = string("sp_lb_", spuid)
  lb_conv_constr = setconstr!(
@@ -156,11 +158,24 @@ function create_side_vars_constrs!(
  kind = Essential, sense = Less, inc_val = 100.0, 
  loc_art_var_abs_cost = env.params.local_art_var_cost
  )
+ masterform.parent_formulation.dw_pricing_sp_ub[spuid] = getid(ub_conv_constr)
+ coefmatrix[getid(ub_conv_constr), getid(setuprepvar)] = 1.0
+
  spform.duty_data.lower_multiplicity = lb_mult
  spform.duty_data.upper_multiplicity = ub_mult
  spform.duty_data.setup_var = getid(setupvar)
- masterform.parent_formulation.dw_pricing_sp_ub[spuid] = getid(ub_conv_constr)
- coefmatrix[getid(ub_conv_constr), getid(setuprepvar)] = 1.0
+
+ # If pricing subproblem variables are continuous, the master columns generated by
+ # the subproblem must have a continuous perenkind.
+ # This piece of information is stored in the duty data of the formulation.
+ continuous_columns = true
+ for (varid, var) in getvars(spform)
+ if getduty(varid) <= DwSpPricingVar && getperenkind(spform, var) !== Continuous
+ continuous_columns = false
+ break
+ end
+ end
+ spform.duty_data.column_var_kind = continuous_columns ? Continuous : Integ
  end
  return
 end

test/issues_tests.jl

@@ -234,6 +234,78 @@ function branching_file_completion()
  @test JuMP.termination_status(model) == MOI.OPTIMAL
 end
 
+# Issue https://github.com/atoptima/Coluna.jl/issues/550
+function continuous_vars_in_sp()
+ # Simple min cost flow problem
+ # 
+ # n1 ------ f[1] -------> n3
+ # \ ^
+ # \ /
+ # -f[2]-> n2 --f[3]--
+ # 
+ # n1: demand = -10.1
+ # n2: demand = 0
+ # n3: demand = 10.1 
+ # f[1]: cost = 0, capacity = 8.5, in mip model only integer flow allowed
+ # f[2]: cost = 50, (capacity = 5 can be activated by removing comment at constraint, line 93)
+ # f[3]: cost = 50
+ #
+ # Correct solution for non-integer f[1]
+ # f[1] = 8.5, f[2] = f[3] = 1.6, cost = 8.5*0 + 1.6*2*50 = 160
+ # Correct solution for integer f[1]
+ # f[1] = 8, f[2] = f[3] = 2.1, cost = 8.5*0 + 2.1*2*50 = 210
+ #
+ function solve_flow_model(f1_integer, coluna)
+ @axis(M, 1:1)
+ model = BlockDecomposition.BlockModel(coluna, direct_model=true)
+ @variable(model, f[1:3, m in M])
+ if f1_integer
+ JuMP.set_integer(f[1, 1])
+ end
+ @constraint(model, n1[m in M], f[1,m] + f[2,m] == 10.1)
+ @constraint(model, n2[m in M], f[2,m] == f[3,m])
+ @constraint(model, n3[m in M], f[1,m] + f[3,m] == 10.1)
+ @constraint(model, cap1, sum(f[1,m] for m in M) <= 8.5)
+ #@JuMP.constraint(model, cap2, sum(f[2,m] for m in M) <= 5)
+ @objective(model, Min, 50 * f[2,1] + 50 * f[3,1])
+
+ @dantzig_wolfe_decomposition(model, decomposition, M)
+
+ subproblems = BlockDecomposition.getsubproblems(decomposition)
+ BlockDecomposition.specify!.(subproblems, lower_multiplicity=1, upper_multiplicity=1)
+
+ optimize!(model)
+
+ if f1_integer
+ @test termination_status(model) == MOI.OPTIMAL
+ @test primal_status(model) == MOI.FEASIBLE_POINT
+ @test objective_value(model) ≈ 210
+ @test value(f[1,1]) ≈ 8
+ @test value(f[2,1]) ≈ 2.1
+ @test value(f[3,1]) ≈ 2.1
+ else
+ @test termination_status(model) == MOI.OPTIMAL
+ @test primal_status(model) == MOI.FEASIBLE_POINT
+ @test objective_value(model) ≈ 160
+ @test value(f[1,1]) ≈ 8.5
+ @test value(f[2,1]) ≈ 1.6
+ @test value(f[3,1]) ≈ 1.6
+ end
+ end
+
+ coluna = JuMP.optimizer_with_attributes(
+ Coluna.Optimizer,
+ "params" => Coluna.Params(
+ solver=Coluna.Algorithm.TreeSearchAlgorithm(),
+ ),
+ "default_optimizer" => GLPK.Optimizer 
+ );
+
+ solve_flow_model(false, coluna)
+ solve_flow_model(true, coluna)
+ return
+end
+
 function test_issues_fixed()
  @testset "no_decomposition" begin
  solve_with_no_decomposition()
@@ -262,6 +334,10 @@ function test_issues_fixed()
  @testset "branching_file_completion" begin
  branching_file_completion()
  end
+
+ @testset "continuous_vars_in_sp" begin
+ continuous_vars_in_sp()
+ end
 end
 
 test_issues_fixed()