diff --git a/src/Algorithm/colgen/stabilization.jl b/src/Algorithm/colgen/stabilization.jl index e6ce04d67..6468b9384 100644 --- a/src/Algorithm/colgen/stabilization.jl +++ b/src/Algorithm/colgen/stabilization.jl @@ -141,33 +141,43 @@ function _primal_solution(master::Formulation, generated_columns, is_minimizatio return sparsevec(var_ids, var_vals) end -function _dynamic_alpha_schedule( - stab::ColGenStab, smooth_dual_sol, h, primal_solution, is_minimization -) +function _increase(smooth_dual_sol, cur_stab_center, h, primal_solution, is_minimization) # Calculate the in-sep direction. - in_sep_direction = smooth_dual_sol - stab.cur_stab_center + in_sep_direction = smooth_dual_sol - cur_stab_center in_sep_dir_norm = norm(in_sep_direction) + # if in & sep are the same point, we need to decrease α becase it is the weight of the + # stability center (in) in the formula to compute the sep point. + if iszero(in_sep_dir_norm) + return false + end + # Calculate the subgradient subgradient = h.a - h.A * primal_solution subgradient_norm = norm(subgradient) # we now calculate the angle between the in-sep direction and the subgradient - angle = (transpose(in_sep_direction) * subgradient) / (in_sep_dir_norm * subgradient_norm) + cos_angle = (transpose(in_sep_direction) * subgradient) / (in_sep_dir_norm * subgradient_norm) + if !is_minimization - angle *= -1 + cos_angle *= -1 end + return cos_angle < 1e-12 +end +function _dynamic_alpha_schedule( + α, smooth_dual_sol, cur_stab_center, h, primal_solution, is_minimization +) + increase = _increase(smooth_dual_sol, cur_stab_center, h, primal_solution, is_minimization) # we modify the alpha parameter based on the calculated angle - α = angle > 1e-12 ? f_decr(stab.base_α) : f_incr(stab.base_α) - return α + return increase ? f_incr(α) : f_decr(α) end function ColGen.update_stabilization_after_iter!(stab::ColGenStab, ctx, master, generated_columns, mast_dual_sol) if stab.smooth_factor == 1 is_min = ColGen.is_minimization(ctx) primal_sol = _primal_solution(master, generated_columns, is_min) - α = _dynamic_alpha_schedule(stab, mast_dual_sol, subgradient_helper(ctx), primal_sol, is_min) + α = _dynamic_alpha_schedule(stab.base_α, mast_dual_sol, stab.cur_stab_center, subgradient_helper(ctx), primal_sol, is_min) stab.base_α = α end diff --git a/src/Algorithm/colgen/utils.jl b/src/Algorithm/colgen/utils.jl index 178161adf..726521e42 100644 --- a/src/Algorithm/colgen/utils.jl +++ b/src/Algorithm/colgen/utils.jl @@ -47,7 +47,9 @@ function _submatrix( form::Formulation, keep_constr::Function, keep_var::Function, + m::Function = (form, is_min, constr_id, var_id) -> 1.0 ) + is_min = getobjsense(form) == MinSense matrix = getcoefmatrix(form) constr_ids = ConstrId[] var_ids = VarId[] @@ -55,9 +57,10 @@ function _submatrix( for constr_id in Iterators.filter(keep_constr, Iterators.keys(getconstrs(form))) for (var_id, coeff) in @view matrix[constr_id, :] if keep_var(var_id) + c = m(form, is_min, constr_id, var_id) push!(constr_ids, constr_id) push!(var_ids, var_id) - push!(nz, coeff) + push!(nz, c * coeff) end end end @@ -154,18 +157,34 @@ end function SubgradientCalculationHelper(master) constr_ids = ConstrId[] constr_rhs = Float64[] + + m_rhs = (master, is_min, constr_id) -> begin + constr_sense = getcursense(master, constr_id) + if is_min + return constr_sense == Less ? -1.0 : 1.0 + else + return constr_sense == Greater ? -1.0 : 1.0 + end + end + m_submatrix = (master, is_min, constr_id, var_id) -> begin + m_rhs(master, is_min, constr_id) + end + + is_min = getobjsense(master) == MinSense for (constr_id, constr) in getconstrs(master) if !(getduty(constr_id) <= MasterConvexityConstr) && iscuractive(master, constr) && isexplicit(master, constr) push!(constr_ids, constr_id) - push!(constr_rhs, getcurrhs(master, constr_id)) + push!(constr_rhs, m_rhs(master, is_min, constr_id) * getcurrhs(master, constr_id)) end end + a = sparsevec(constr_ids, constr_rhs, Coluna.MAX_NB_ELEMS) A = _submatrix( master, constr_id -> !(getduty(constr_id) <= MasterConvexityConstr), - var_id -> getduty(var_id) <= MasterPureVar || getduty(var_id) <= MasterRepPricingVar + var_id -> getduty(var_id) <= MasterPureVar || getduty(var_id) <= MasterRepPricingVar, + m_submatrix ) return SubgradientCalculationHelper(a, A) end \ No newline at end of file diff --git a/test/unit/ColGen/colgen_default.jl b/test/unit/ColGen/colgen_default.jl index b9d3bd8b8..032ff1e24 100644 --- a/test/unit/ColGen/colgen_default.jl +++ b/test/unit/ColGen/colgen_default.jl @@ -155,9 +155,9 @@ function test_subgradient_calculation_helper() helper = ClA.SubgradientCalculationHelper(master) @test helper.a[cids["c1"]] == 10 - @test helper.a[cids["c2"]] == 100 + @test helper.a[cids["c2"]] == -100 @test helper.a[cids["c3"]] == 100 - @test helper.a[cids["c4"]] == 5 + @test helper.a[cids["c4"]] == -5 @test helper.A[cids["c1"], vids["x1"]] == 1 @test helper.A[cids["c1"], vids["x2"]] == 1 @@ -167,11 +167,11 @@ function test_subgradient_calculation_helper() @test helper.A[cids["c1"], vids["y3"]] == 1 @test helper.A[cids["c1"], vids["z1"]] == 2 @test helper.A[cids["c1"], vids["z2"]] == 1 - @test helper.A[cids["c2"], vids["x1"]] == 1 - @test helper.A[cids["c2"], vids["x2"]] == 2 - @test helper.A[cids["c2"], vids["y1"]] == 1 - @test helper.A[cids["c2"], vids["y2"]] == 2 - @test helper.A[cids["c2"], vids["z1"]] == 1 + @test helper.A[cids["c2"], vids["x1"]] == -1 + @test helper.A[cids["c2"], vids["x2"]] == -2 + @test helper.A[cids["c2"], vids["y1"]] == -1 + @test helper.A[cids["c2"], vids["y2"]] == -2 + @test helper.A[cids["c2"], vids["z1"]] == -1 @test helper.A[cids["c2"], vids["z2"]] == 0 @test helper.A[cids["c3"], vids["x1"]] == 1 @test helper.A[cids["c3"], vids["x3"]] == 3 @@ -179,10 +179,10 @@ function test_subgradient_calculation_helper() @test helper.A[cids["c3"], vids["y3"]] == 3 @test helper.A[cids["c3"], vids["z1"]] == 0 @test helper.A[cids["c3"], vids["z2"]] == 0 - @test helper.A[cids["c4"], vids["z1"]] == 1 - @test helper.A[cids["c4"], vids["z2"]] == 1 + @test helper.A[cids["c4"], vids["z1"]] == -1 + @test helper.A[cids["c4"], vids["z2"]] == -1 end -register!(unit_tests, "colgen_default", test_subgradient_calculation_helper) +register!(unit_tests, "colgen_default", test_subgradient_calculation_helper; f = true) # All the tests are based on the Generalized Assignment problem. # x_mj = 1 if job j is assigned to machine m diff --git a/test/unit/ColGen/colgen_stabilization.jl b/test/unit/ColGen/colgen_stabilization.jl index 8e1799b06..c676129c3 100644 --- a/test/unit/ColGen/colgen_stabilization.jl +++ b/test/unit/ColGen/colgen_stabilization.jl @@ -64,7 +64,6 @@ form_primal_solution() = """ z2 >= 3 """ - function test_primal_solution() _, master, sps, _, _ = reformfromstring(form_primal_solution()) @@ -113,11 +112,180 @@ function test_primal_solution() end register!(unit_tests, "colgen_stabilization", test_primal_solution; f = true) +form_primal_solution2() = """ + master + min + 3x_11 + 2x_12 + 2x_21 + x_22 + z1 + z2 + 0s1 + 0s2 + s.t. + x_11 + x_12 + x_21 + 2z1 + z2 >= 10 + x_11 + 2x_12 + x_21 + 2x_22 + z1 <= 100 + x_11 + x_22 + == 100 + z1 + z2 >= 5 + s1 >= 1 {MasterConvexityConstr} + s1 <= 2 {MasterConvexityConstr} + s2 >= 0 {MasterConvexityConstr} + s2 <= 3 {MasterConvexityConstr} + + dw_sp + min + x_11 + x_12 + 0s1 + s.t. + x_11 + x_12 >= 10 + + dw_sp + min + x_21 + x_22 + 0s2 + s.t. + x_21 + x_22 >= 10 + + integer + representatives + x_11, x_12, x_21, x_22 + + pure + z1, z2 + + pricing_setup + s1, s2 + + bounds + x_11 >= 0 + x_12 >= 0 + x_21 >= 0 + x_22 >= 1 + z1 >= 0 + z2 >= 3 +""" + +# We consider the master with the following coefficient matrix: +# master_coeff_matrix = [ +# 1 1 1 0 1 1; +# -1 -2 -1 -2 -1 0; +# 1 0 0 1 0 0; # is it correct to handle an "==" constraint like this in subgradient computation? +# 0 0 0 0 1 1; +# ] +# the following rhs: rhs = [10, -100, 100, 5] + +# We consider the primal solution: primal = [1, 2, 12, 12, 0, 0] +# The subgradient is therefore: rhs - master_coeff_matrix * primal = [-5, -59, 87, 5] +# We use the following stability center: stab = [1, 2, 0, 1] + +function _test_angle_primal_sol(master, sp1, sp2) + vids = get_name_to_varids(master) + pool = Coluna.Algorithm.ColumnsSet() + + sol1 = Coluna.MathProg.PrimalSolution( + sp1, + [vids["x_11"], vids["x_12"]], + [1.0, 2.0], + 11.0, + Coluna.MathProg.FEASIBLE_SOL + ) + Coluna.Algorithm.add_primal_sol!(pool.subprob_primal_sols, sol1, false) # improving = true + + sol2 = Coluna.MathProg.PrimalSolution( + sp2, + [vids["x_21"], vids["x_22"]], + [4.0, 4.0], + 13.0, + Coluna.MathProg.FEASIBLE_SOL + ) + Coluna.Algorithm.add_primal_sol!(pool.subprob_primal_sols, sol2, true) + is_minimization = true + primal_sol = Coluna.Algorithm._primal_solution(master, pool, is_minimization) + return primal_sol +end + +function _test_angle_stab_center(master) + cids = get_name_to_constrids(master) + return Coluna.MathProg.DualSolution( + master, + [cids["c1"], cids["c2"], cids["c4"]], + [1.0, 2.0, 1.0], + Coluna.MathProg.VarId[], Float64[], Coluna.MathProg.ActiveBound[], + 0.0, + Coluna.MathProg.FEASIBLE_SOL + ) +end + # Make sure the angle is well computed. -function test_angle() +# Here we test the can where the in and sep points are the same. +# In that case, we should decrease the value of α. +function test_angle_1() + _, master, sps, _, _ = reformfromstring(form_primal_solution2()) + sp1, sp2 = sps[2], sps[1] + cids = get_name_to_constrids(master) + smooth_dual_sol = Coluna.MathProg.DualSolution( + master, + [cids["c1"], cids["c2"], cids["c4"]], + [1.0, 2.0, 1.0], + Coluna.MathProg.VarId[], Float64[], Coluna.MathProg.ActiveBound[], + 0.0, + Coluna.MathProg.FEASIBLE_SOL + ) + cur_stab_center = _test_angle_stab_center(master) + + h = Coluna.Algorithm.SubgradientCalculationHelper(master) + is_minimization = true + primal_sol = _test_angle_primal_sol(master, sp1, sp2) + increase = Coluna.Algorithm._increase(smooth_dual_sol, cur_stab_center, h, primal_sol, is_minimization) + @test increase == false +end +register!(unit_tests, "colgen_stabilization", test_angle_1; f = true) + +# Let's consider the following sep point: sep = [5, 7, 0, 3] +# The direction will be [4, 5, 0, 2] and should lead to a negative cosinus for the angle. +# In that case, we need to increase the value of α. +function test_angle_2() + _, master, sps, _, _ = reformfromstring(form_primal_solution2()) + sp1, sp2 = sps[2], sps[1] + cids = get_name_to_constrids(master) + + smooth_dual_sol = Coluna.MathProg.DualSolution( + master, + [cids["c1"], cids["c2"], cids["c4"]], + [5.0, 7.0, 3.0], + Coluna.MathProg.VarId[], Float64[], Coluna.MathProg.ActiveBound[], + 0.0, + Coluna.MathProg.FEASIBLE_SOL + ) + cur_stab_center = _test_angle_stab_center(master) + + h = Coluna.Algorithm.SubgradientCalculationHelper(master) + is_minimization = true + primal_sol = _test_angle_primal_sol(master, sp1, sp2) + increase = Coluna.Algorithm._increase(smooth_dual_sol, cur_stab_center, h, primal_sol, is_minimization) + @test increase == true end -register!(unit_tests, "colgen_stabilization", test_angle; f = true) +register!(unit_tests, "colgen_stabilization", test_angle_2; f = true) + +# Let's consider the following sep point: sep = [5, 1, 10, 3] +# The direction will be [4, 1, 10, 2] and should lead to a positive cosinus for the angle. +# In that case, we need to decrease the value of α. +function test_angle_3() + _, master, sps, _, _ = reformfromstring(form_primal_solution2()) + sp1, sp2 = sps[2], sps[1] + cids = get_name_to_constrids(master) + + smooth_dual_sol = Coluna.MathProg.DualSolution( + master, + [cids["c1"], cids["c2"], cids["c3"], cids["c4"]], + [5.0, 1.0, 10.0, 3.0], + Coluna.MathProg.VarId[], Float64[], Coluna.MathProg.ActiveBound[], + 0.0, + Coluna.MathProg.FEASIBLE_SOL + ) + cur_stab_center = _test_angle_stab_center(master) + + h = Coluna.Algorithm.SubgradientCalculationHelper(master) + is_minimization = true + primal_sol = _test_angle_primal_sol(master, sp1, sp2) + increase = Coluna.Algorithm._increase(smooth_dual_sol, cur_stab_center, h, primal_sol, is_minimization) + @test increase == false +end +register!(unit_tests, "colgen_stabilization", test_angle_3; f = true) + function test_dynamic_alpha_schedule()