From 8ad3e58c6cf9185f64ca55414616d8701b629cbc Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Mon, 13 Jul 2020 17:10:13 +0200 Subject: [PATCH] [CP-SAT]: Fix overflow in feasibility pump; improve rounding heuristics in feasibility pump; implement better algorithm to provide explanation of infeasibility; work on max hitting set implementation for core based search --- ortools/sat/cp_model.proto | 2 + ortools/sat/cp_model_lns.cc | 13 ++- ortools/sat/cp_model_presolve.cc | 6 ++ ortools/sat/cp_model_solver.cc | 59 ++++++++--- ortools/sat/feasibility_pump.cc | 41 ++++--- ortools/sat/feasibility_pump.h | 7 +- ortools/sat/optimization.cc | 177 ++++++++++++++++++------------- ortools/sat/optimization.h | 9 +- ortools/sat/presolve_context.h | 9 ++ ortools/sat/synchronization.cc | 8 ++ ortools/sat/synchronization.h | 4 + 11 files changed, 215 insertions(+), 120 deletions(-) diff --git a/ortools/sat/cp_model.proto b/ortools/sat/cp_model.proto index 3d6245ece5..06abe9eb94 100644 --- a/ortools/sat/cp_model.proto +++ b/ortools/sat/cp_model.proto @@ -621,6 +621,8 @@ message CpSolverResponse { // If you really want a minimal subset, a possible way to get one is by // changing your model to minimize the number of assumptions at false, but // this is likely an harder problem to solve. + // + // TODO(user): Allows for returning multiple core at once. repeated int32 sufficient_assumptions_for_infeasibility = 23; // This will be true iff the solver was asked to find all solutions to a diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc index 6ca2f6ccc4..0f2270a5ac 100644 --- a/ortools/sat/cp_model_lns.cc +++ b/ortools/sat/cp_model_lns.cc @@ -21,9 +21,11 @@ #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/cp_model_loader.h" #include "ortools/sat/cp_model_utils.h" +#include "ortools/sat/integer.h" #include "ortools/sat/linear_programming_constraint.h" #include "ortools/sat/rins.h" #include "ortools/sat/synchronization.h" +#include "ortools/util/saturated_arithmetic.h" namespace operations_research { namespace sat { @@ -285,8 +287,10 @@ void NeighborhoodGenerator::Synchronize() { // This might not be a final solution, but it does work ok for now. const IntegerValue best_objective_improvement = IsRelaxationGenerator() - ? (data.new_objective_bound - data.initial_best_objective_bound) - : (data.initial_best_objective - data.new_objective); + ? IntegerValue(CapSub(data.new_objective_bound.value(), + data.initial_best_objective_bound.value())) + : IntegerValue(CapSub(data.initial_best_objective.value(), + data.new_objective.value())); if (best_objective_improvement > 0) { num_consecutive_non_improving_calls_ = 0; } else { @@ -630,7 +634,8 @@ Neighborhood RelaxationInducedNeighborhoodGenerator::Generate( (incomplete_solutions_ != nullptr && incomplete_solutions_->HasNewSolution()); - if (!lp_solution_available && !relaxation_solution_available) { + if (!lp_solution_available && !relaxation_solution_available && + !incomplete_solution_available) { return neighborhood; } @@ -651,7 +656,7 @@ Neighborhood RelaxationInducedNeighborhoodGenerator::Generate( neighborhood.source_info = incomplete_solution_available ? "incomplete" : "lp"; } else { - CHECK(relaxation_solution_available); + CHECK(relaxation_solution_available || incomplete_solution_available); rins_neighborhood = GetRINSNeighborhood( response_manager_, relaxation_solutions_, /*lp_solutions=*/nullptr, incomplete_solutions_, random); diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index c82c5071f7..8e503f6a81 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -4649,6 +4649,7 @@ bool CpModelPresolver::Presolve() { if (context_->ModelIsUnsat()) break; } context_->UpdateNewConstraintsVariableUsage(); + context_->RegisterVariablesUsedInAssumptions(); DCHECK(context_->ConstraintVariableUsageIsConsistent()); // Main propagation loop. @@ -4917,6 +4918,11 @@ void ApplyVariableMapping(const std::vector& mapping, } } + // Remap the assumptions. + for (int& mutable_ref : *proto->mutable_assumptions()) { + mapping_function(&mutable_ref); + } + // Remap the search decision heuristic. // Note that we delete any heuristic related to a removed variable. for (DecisionStrategyProto& strategy : *proto->mutable_search_strategy()) { diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 70c89f1128..b3d43b38d4 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -1476,11 +1476,13 @@ void SolveLoadedCpModel(const CpModelProto& model_proto, // Reconfigure search heuristic if it was changed. ConfigureSearchHeuristics(model); + const auto& mapping = *model->GetOrCreate(); SatSolver::Status status; const SatParameters& parameters = *model->GetOrCreate(); if (!model_proto.has_objective()) { while (true) { - status = ResetAndSolveIntegerProblem(/*assumptions=*/{}, model); + status = ResetAndSolveIntegerProblem( + mapping.Literals(model_proto.assumptions()), model); if (status != SatSolver::Status::FEASIBLE) break; solution_observer(); if (!parameters.enumerate_all_solutions()) break; @@ -1490,6 +1492,24 @@ void SolveLoadedCpModel(const CpModelProto& model_proto, shared_response_manager->NotifyThatImprovingProblemIsInfeasible( solution_info); } + if (status == SatSolver::ASSUMPTIONS_UNSAT) { + shared_response_manager->NotifyThatImprovingProblemIsInfeasible( + solution_info); + + // Extract a good subset of assumptions and add it to the response. + auto* sat_solver = model->GetOrCreate(); + std::vector core = sat_solver->GetLastIncompatibleDecisions(); + MinimizeCoreWithPropagation(sat_solver, &core); + std::vector core_in_proto_format; + for (const Literal l : core) { + core_in_proto_format.push_back( + mapping.GetProtoVariableFromBooleanVariable(l.Variable())); + if (!l.IsPositive()) { + core_in_proto_format.back() = NegatedRef(core_in_proto_format.back()); + } + } + shared_response_manager->AddUnsatCore(core_in_proto_format); + } } else { // Optimization problem. const auto& objective = *model->GetOrCreate(); @@ -1501,8 +1521,7 @@ void SolveLoadedCpModel(const CpModelProto& model_proto, // shouldn't be too hard to fix. if (parameters.optimize_with_max_hs()) { status = MinimizeWithHittingSetAndLazyEncoding( - objective_var, objective.vars, objective.coeffs, solution_observer, - model); + objective, solution_observer, model); } else { status = model->Mutable()->Optimize(); } @@ -1551,8 +1570,9 @@ void QuickSolveWithHint(const CpModelProto& model_proto, // Solve decision problem. ConfigureSearchHeuristics(model); - const SatSolver::Status status = - ResetAndSolveIntegerProblem(/*assumptions=*/{}, model); + const auto& mapping = *model->GetOrCreate(); + const SatSolver::Status status = ResetAndSolveIntegerProblem( + mapping.Literals(model_proto.assumptions()), model); const std::string& solution_info = model->Name(); if (status == SatSolver::Status::FEASIBLE) { @@ -2798,15 +2818,22 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) { auto context = absl::make_unique(&new_cp_model_proto, &mapping_proto); - // Load the assumptions if any. For now we just fix the variables. - // - // TODO(user): Handle this properly. - context->InitializeNewDomains(); - for (const int ref : model_proto.assumptions()) { - if (!context->SetLiteralToTrue(ref)) { - final_response.set_status(CpSolverStatus::INFEASIBLE); - final_response.add_sufficient_assumptions_for_infeasibility(ref); - return final_response; + bool degraded_assumptions_support = false; + if (params.num_search_workers() > 1 || model_proto.has_objective()) { + // For the case where the assumptions are currently not supported, we just + // assume they are fixed, and will always report all of them in the UNSAT + // core if the problem turn out to be UNSAT. + // + // If the mode is not degraded, we will hopefully report a small subset + // in case there is no feasible solution under these assumptions. + degraded_assumptions_support = true; + context->InitializeNewDomains(); + for (const int ref : model_proto.assumptions()) { + if (!context->SetLiteralToTrue(ref)) { + final_response.set_status(CpSolverStatus::INFEASIBLE); + final_response.add_sufficient_assumptions_for_infeasibility(ref); + return final_response; + } } } @@ -2982,8 +3009,8 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) { model_proto, std::vector(final_response.solution().begin(), final_response.solution().end()))); } - // TODO(user): Handle this properly. - if (final_response.status() == CpSolverStatus::INFEASIBLE) { + if (degraded_assumptions_support && + final_response.status() == CpSolverStatus::INFEASIBLE) { // For now, just pass in all assumptions. *final_response.mutable_sufficient_assumptions_for_infeasibility() = model_proto.assumptions(); diff --git a/ortools/sat/feasibility_pump.cc b/ortools/sat/feasibility_pump.cc index 9a3aef4605..0a2ace5cc0 100644 --- a/ortools/sat/feasibility_pump.cc +++ b/ortools/sat/feasibility_pump.cc @@ -46,7 +46,9 @@ FeasibilityPump::FeasibilityPump(Model* model) mapping_(model->Get()) { // Tweak the default parameters to make the solve incremental. glop::GlopParameters parameters; - // TODO(user): Determine which simplex is better for this, dual or primal. + // Note(user): Primal simplex does better here since we have a limit on + // simplex iterations. So dual simplex sometimes fails to find a LP feasible + // solution. parameters.set_use_dual_simplex(false); parameters.set_max_number_of_iterations(2000); simplex_.SetParameters(parameters); @@ -461,10 +463,9 @@ bool FeasibilityPump::LockBasedRounding() { // We compute the number of locks based on variable coefficient in constraints // and constraint bounds. This doesn't change over time so we cache it. - if (var_has_more_up_locks_.empty()) { - std::vector up_locks(num_vars, 0); - std::vector down_locks(num_vars, 0); - var_has_more_up_locks_.resize(num_vars, false); + if (var_up_locks_.empty()) { + var_up_locks_.resize(num_vars, 0); + var_down_locks_.resize(num_vars, 0); for (int i = 0; i < num_vars; ++i) { for (const auto entry : lp_data_.GetSparseColumn(ColIndex(i))) { ColIndex slack = lp_data_.GetSlackVariable(entry.row()); @@ -475,23 +476,21 @@ bool FeasibilityPump::LockBasedRounding() { lp_data_.variable_upper_bounds()[slack] < glop::kInfinity; if (entry.coefficient() > 0) { - up_locks[i] += constraint_upper_bounded; - down_locks[i] += constraint_lower_bounded; + var_up_locks_[i] += constraint_upper_bounded; + var_down_locks_[i] += constraint_lower_bounded; } else { - up_locks[i] += constraint_lower_bounded; - down_locks[i] += constraint_upper_bounded; + var_up_locks_[i] += constraint_lower_bounded; + var_down_locks_[i] += constraint_upper_bounded; } } - var_has_more_up_locks_[i] = up_locks[i] > down_locks[i]; } } for (int i = 0; i < lp_solution_.size(); ++i) { - if (std::abs(lp_solution_[i] - std::round(lp_solution_[i])) < 0.1) { + if (std::abs(lp_solution_[i] - std::round(lp_solution_[i])) < 0.1 || + var_up_locks_[i] == var_down_locks_[i]) { integer_solution_[i] = static_cast(std::round(lp_solution_[i])); - } else if (var_has_more_up_locks_[i]) { - // TODO(user): When #up_locks and #down_locks are same, round to - // nearest integer. + } else if (var_up_locks_[i] > var_down_locks_[i]) { integer_solution_[i] = static_cast(std::floor(lp_solution_[i])); } else { integer_solution_[i] = static_cast(std::ceil(lp_solution_[i])); @@ -653,7 +652,10 @@ bool FeasibilityPump::PropagationRounding() { integer_encoder_->GetOrCreateLiteralAssociatedToEquality(var, value); } - if (!sat_solver_->FinishPropagation()) continue; + if (!sat_solver_->FinishPropagation()) { + model_is_unsat_ = true; + return false; + } sat_solver_->EnqueueDecisionAndBacktrackOnConflict(to_enqueue); if (sat_solver_->IsModelUnsat()) { @@ -692,10 +694,15 @@ void FeasibilityPump::FillIntegerSolutionStats() { if (activity > integer_lp_[i].ub || activity < integer_lp_[i].lb) { integer_solution_is_feasible_ = false; num_infeasible_constraints_++; + const int64 ub_infeasibility = activity > integer_lp_[i].ub.value() + ? activity - integer_lp_[i].ub.value() + : 0; + const int64 lb_infeasibility = activity < integer_lp_[i].lb.value() + ? integer_lp_[i].lb.value() - activity + : 0; integer_solution_infeasibility_ = std::max(integer_solution_infeasibility_, - std::max(activity - integer_lp_[i].ub.value(), - integer_lp_[i].lb.value() - activity)); + std::max(ub_infeasibility, lb_infeasibility)); } } } diff --git a/ortools/sat/feasibility_pump.h b/ortools/sat/feasibility_pump.h index 2fe7e056f9..8721d65461 100644 --- a/ortools/sat/feasibility_pump.h +++ b/ortools/sat/feasibility_pump.h @@ -179,10 +179,11 @@ class FeasibilityPump { // True if the variable was binary before we apply scaling. std::vector var_is_binary_; - // True if variable has more number of constraints restricting it to take - // higher values than number of constraints restricting it to take lower + // The following lock information is computed only once. + // Number of constraints restricting variable to take higher (resp. lower) // values. - std::vector var_has_more_up_locks_; + std::vector var_up_locks_; + std::vector var_down_locks_; // We need to remember what to optimize if an objective is given, because // then we will switch the objective between feasibility and optimization. diff --git a/ortools/sat/optimization.cc b/ortools/sat/optimization.cc index e011e22ea0..fb818552b5 100644 --- a/ortools/sat/optimization.cc +++ b/ortools/sat/optimization.cc @@ -29,13 +29,12 @@ #include "absl/strings/str_format.h" #include "ortools/base/cleanup.h" #include "ortools/base/int_type.h" -#include "ortools/base/integral_types.h" #include "ortools/base/logging.h" #include "ortools/base/macros.h" #include "ortools/base/map_util.h" #include "ortools/base/stl_util.h" #include "ortools/base/timer.h" -#if !defined(__PORTABLE_PLATFORM__) && (defined(USE_CBC) || defined(USE_SCIP)) +#if !defined(__PORTABLE_PLATFORM__) && defined(USE_SCIP) #include "ortools/linear_solver/linear_solver.h" #include "ortools/linear_solver/linear_solver.pb.h" #endif // __PORTABLE_PLATFORM__ @@ -1298,6 +1297,45 @@ SatSolver::Status FindCores(std::vector assumptions, return SatSolver::ASSUMPTIONS_UNSAT; } +// Slightly different algo than FindCores() which aim to extract more cores, but +// not necessarily non-overlaping ones. +SatSolver::Status FindMultipleCoresForMaxHs( + std::vector assumptions, Model* model, + std::vector>* cores) { + cores->clear(); + SatSolver* sat_solver = model->GetOrCreate(); + TimeLimit* limit = model->GetOrCreate(); + do { + if (limit->LimitReached()) return SatSolver::LIMIT_REACHED; + + const SatSolver::Status result = + ResetAndSolveIntegerProblem(assumptions, model); + if (result != SatSolver::ASSUMPTIONS_UNSAT) return result; + std::vector core = sat_solver->GetLastIncompatibleDecisions(); + if (sat_solver->parameters().minimize_core()) { + MinimizeCoreWithPropagation(sat_solver, &core); + } + CHECK(!core.empty()); + cores->push_back(core); + if (!sat_solver->parameters().find_multiple_cores()) break; + + // Pick a random literal from the core and remove it from the set of + // assumptions. + CHECK(!core.empty()); + auto* random = model->GetOrCreate(); + const Literal random_literal = + core[absl::Uniform(*random, 0, core.size())]; + for (int i = 0; i < assumptions.size(); ++i) { + if (assumptions[i] == random_literal) { + std::swap(assumptions[i], assumptions.back()); + assumptions.pop_back(); + break; + } + } + } while (!assumptions.empty()); + return SatSolver::ASSUMPTIONS_UNSAT; +} + } // namespace CoreBasedOptimizer::CoreBasedOptimizer( @@ -1741,10 +1779,14 @@ SatSolver::Status CoreBasedOptimizer::Optimize() { // // TODO(user): remove code duplication with MinimizeWithCoreAndLazyEncoding(); SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( - IntegerVariable objective_var, std::vector variables, - std::vector coefficients, + const ObjectiveDefinition& objective_definition, const std::function& feasible_solution_observer, Model* model) { -#if !defined(__PORTABLE_PLATFORM__) && (defined(USE_CBC) || defined(USE_SCIP)) +#if !defined(__PORTABLE_PLATFORM__) && defined(USE_SCIP) + + IntegerVariable objective_var = objective_definition.objective_var; + std::vector variables = objective_definition.vars; + std::vector coefficients = objective_definition.coeffs; + SatSolver* sat_solver = model->GetOrCreate(); IntegerTrail* integer_trail = model->GetOrCreate(); IntegerEncoder* integer_encoder = model->GetOrCreate(); @@ -1779,16 +1821,12 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( // This is the "generalized" hitting set problem we will solve. Each time // we find a core, a new constraint will be added to this problem. MPModelRequest request; -#if defined(USE_SCIP) request.set_solver_specific_parameters("limits/gap = 0"); request.set_solver_type(MPModelRequest::SCIP_MIXED_INTEGER_PROGRAMMING); -#else // USE_CBC - request.set_solver_type(MPModelRequest::CBC_MIXED_INTEGER_PROGRAMMING); -#endif // USE_CBC or USE_SCIP MPModelProto& hs_model = *request.mutable_model(); - const int num_variables = variables.size(); - for (int i = 0; i < num_variables; ++i) { + const int num_variables_in_objective = variables.size(); + for (int i = 0; i < num_variables_in_objective; ++i) { if (coefficients[i] < 0) { variables[i] = NegationOf(variables[i]); coefficients[i] = -coefficients[i]; @@ -1801,12 +1839,6 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( var_proto->set_is_integer(true); } -// The MIP solver. -#if defined(USE_SCIP) - MPSolver solver("HS solver", MPSolver::SCIP_MIXED_INTEGER_PROGRAMMING); -#else // USE_CBC - MPSolver solver("HS solver", MPSolver::CBC_MIXED_INTEGER_PROGRAMMING); -#endif // USE_CBC or USE_SCIP MPSolutionResponse response; // This is used by the "stratified" approach. We will only consider terms with @@ -1817,7 +1849,7 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( // TODO(user): The core is returned in the same order as the assumptions, // so we don't really need this map, we could just do a linear scan to // recover which node are part of the core. - std::map assumption_to_index; + std::map> assumption_to_indices; // New Booleans variable in the MIP model to represent X >= cte. std::map, int> created_var; @@ -1831,11 +1863,14 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( // not really done incrementally. It might be hard to improve though. // // TODO(user): deal with time limit. - solver.SolveWithProto(request, &response); + MPSolver::SolveWithProto(request, &response); CHECK_EQ(response.status(), MPSolverResponseStatus::MPSOLVER_OPTIMAL); + + const IntegerValue mip_objective( + static_cast(std::round(response.objective_value()))); VLOG(1) << "constraints: " << hs_model.constraint_size() - << " variables: " << hs_model.variable_size() - << " mip_lower_bound: " << response.objective_value() + << " variables: " << hs_model.variable_size() << " hs_lower_bound: " + << objective_definition.ScaleIntegerObjective(mip_objective) << " strat: " << stratified_threshold; // Update the objective lower bound with our current bound. @@ -1843,10 +1878,8 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( // Note(user): This is not needed for correctness, but it might cause // more propagation and is nice to have for reporting/logging purpose. if (!integer_trail->Enqueue( - IntegerLiteral::GreaterOrEqual( - objective_var, - IntegerValue(static_cast(response.objective_value()))), - {}, {})) { + IntegerLiteral::GreaterOrEqual(objective_var, mip_objective), {}, + {})) { result = SatSolver::INFEASIBLE; break; } @@ -1854,9 +1887,9 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( sat_solver->Backtrack(0); sat_solver->SetAssumptionLevel(0); std::vector assumptions; - assumption_to_index.clear(); + assumption_to_indices.clear(); IntegerValue next_stratified_threshold(0); - for (int i = 0; i < num_variables; ++i) { + for (int i = 0; i < num_variables_in_objective; ++i) { const IntegerValue hs_value( static_cast(response.variable_value(i))); if (hs_value == integer_trail->UpperBound(variables[i])) continue; @@ -1866,22 +1899,20 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( next_stratified_threshold = std::max(next_stratified_threshold, coefficients[i]); } else { + // It is possible that different variables have the same associated + // literal. So we do need to consider this case. assumptions.push_back(integer_encoder->GetOrCreateAssociatedLiteral( IntegerLiteral::LowerOrEqual(variables[i], hs_value))); - gtl::InsertOrDie(&assumption_to_index, assumptions.back().Index(), i); + assumption_to_indices[assumptions.back().Index()].push_back(i); } } // No assumptions with the current stratified_threshold? use the new one. - if (assumptions.empty()) { - if (next_stratified_threshold > 0) { - CHECK_LT(next_stratified_threshold, stratified_threshold); - stratified_threshold = next_stratified_threshold; - --iter; // "false" iteration, the lower bound does not increase. - continue; - } else { - return SatSolver::INFEASIBLE; - } + if (assumptions.empty() && next_stratified_threshold > 0) { + CHECK_LT(next_stratified_threshold, stratified_threshold); + stratified_threshold = next_stratified_threshold; + --iter; // "false" iteration, the lower bound does not increase. + continue; } // TODO(user): we could also randomly shuffle the assumptions to find more @@ -1889,10 +1920,7 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( // // TODO(user): Use the real weights and exploit the extra cores. std::vector> cores; - std::vector assumption_weights(assumptions.size(), - IntegerValue(1)); - result = FindCores(assumptions, assumption_weights, stratified_threshold, - model, &cores); + result = FindMultipleCoresForMaxHs(assumptions, model, &cores); if (result == SatSolver::FEASIBLE) { if (!process_solution()) return SatSolver::INFEASIBLE; if (parameters.stop_after_first_solution()) { @@ -1914,10 +1942,11 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( sat_solver->SetAssumptionLevel(0); for (const std::vector& core : cores) { if (core.size() == 1) { - const int index = - gtl::FindOrDie(assumption_to_index, core.front().Index()); - hs_model.mutable_variable(index)->set_lower_bound( - integer_trail->LowerBound(variables[index]).value()); + for (const int index : + gtl::FindOrDie(assumption_to_indices, core.front().Index())) { + hs_model.mutable_variable(index)->set_lower_bound( + integer_trail->LowerBound(variables[index]).value()); + } continue; } @@ -1925,44 +1954,46 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( MPConstraintProto* ct = hs_model.add_constraint(); ct->set_lower_bound(1.0); for (const Literal lit : core) { - const int index = gtl::FindOrDie(assumption_to_index, lit.Index()); - const double lb = integer_trail->LowerBound(variables[index]).value(); - const double hs_value = response.variable_value(index); - if (hs_value == lb) { - ct->add_var_index(index); - ct->add_coefficient(1.0); - ct->set_lower_bound(ct->lower_bound() + lb); - } else { - const std::pair key = {index, hs_value}; - if (!gtl::ContainsKey(created_var, key)) { - const int new_var_index = hs_model.variable_size(); - created_var[key] = new_var_index; + for (const int index : + gtl::FindOrDie(assumption_to_indices, lit.Index())) { + const double lb = integer_trail->LowerBound(variables[index]).value(); + const double hs_value = response.variable_value(index); + if (hs_value == lb) { + ct->add_var_index(index); + ct->add_coefficient(1.0); + ct->set_lower_bound(ct->lower_bound() + lb); + } else { + const std::pair key = {index, hs_value}; + if (!gtl::ContainsKey(created_var, key)) { + const int new_var_index = hs_model.variable_size(); + created_var[key] = new_var_index; - MPVariableProto* new_var = hs_model.add_variable(); - new_var->set_lower_bound(0); - new_var->set_upper_bound(1); - new_var->set_is_integer(true); + MPVariableProto* new_var = hs_model.add_variable(); + new_var->set_lower_bound(0); + new_var->set_upper_bound(1); + new_var->set_is_integer(true); - // (new_var == 1) => x > hs_value. - // (x - lb) - (hs_value - lb + 1) * new_var >= 0. - MPConstraintProto* implication = hs_model.add_constraint(); - implication->set_lower_bound(lb); - implication->add_var_index(index); - implication->add_coefficient(1.0); - implication->add_var_index(new_var_index); - implication->add_coefficient(lb - hs_value - 1); + // (new_var == 1) => x > hs_value. + // (x - lb) - (hs_value - lb + 1) * new_var >= 0. + MPConstraintProto* implication = hs_model.add_constraint(); + implication->set_lower_bound(lb); + implication->add_var_index(index); + implication->add_coefficient(1.0); + implication->add_var_index(new_var_index); + implication->add_coefficient(lb - hs_value - 1); + } + ct->add_var_index(gtl::FindOrDieNoPrint(created_var, key)); + ct->add_coefficient(1.0); } - ct->add_var_index(gtl::FindOrDieNoPrint(created_var, key)); - ct->add_coefficient(1.0); } } } } return result; -#else // !__PORTABLE_PLATFORM__ && (USE_CBC || USE_SCIP) +#else // !__PORTABLE_PLATFORM__ && USE_SCIP LOG(FATAL) << "Not supported."; -#endif // __PORTABLE_PLATFORM || (!USE_CBC && !USE_SCIP) +#endif // !__PORTABLE_PLATFORM__ && USE_SCIP } } // namespace sat diff --git a/ortools/sat/optimization.h b/ortools/sat/optimization.h index b1b6e4d8a7..bdaade5da5 100644 --- a/ortools/sat/optimization.h +++ b/ortools/sat/optimization.h @@ -11,11 +11,6 @@ // See the License for the specific language governing permissions and // limitations under the License. -// Optimization algorithms to solve a LinearBooleanProblem by using the SAT -// solver as a black-box. -// -// TODO(user): Currently, only the MINIMIZATION problem type is supported. - #ifndef OR_TOOLS_SAT_OPTIMIZATION_H_ #define OR_TOOLS_SAT_OPTIMIZATION_H_ @@ -23,6 +18,7 @@ #include #include "ortools/sat/boolean_problem.pb.h" +#include "ortools/sat/cp_model_loader.h" #include "ortools/sat/integer.h" #include "ortools/sat/integer_search.h" #include "ortools/sat/model.h" @@ -237,8 +233,7 @@ class CoreBasedOptimizer { // TODO(user): This function brings dependency to the SCIP MIP solver which is // quite big, maybe we should find a way not to do that. SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( - IntegerVariable objective_var, std::vector variables, - std::vector coefficients, + const ObjectiveDefinition& objective_definition, const std::function& feasible_solution_observer, Model* model); } // namespace sat diff --git a/ortools/sat/presolve_context.h b/ortools/sat/presolve_context.h index 5b13741001..c3fb5b660d 100644 --- a/ortools/sat/presolve_context.h +++ b/ortools/sat/presolve_context.h @@ -32,6 +32,7 @@ namespace sat { // We use some special constraint index in our variable <-> constraint graph. constexpr int kObjectiveConstraint = -1; constexpr int kAffineRelationConstraint = -2; +constexpr int kAssumptionsConstraint = -3; struct PresolveOptions { bool log_info = true; @@ -321,6 +322,14 @@ class PresolveContext { return interval_usage_[c]; } + // Make sure we never delete an "assumption" literal by using a special + // constraint for that. + void RegisterVariablesUsedInAssumptions() { + for (const int ref : working_model->assumptions()) { + var_to_constraints_[PositiveRef(ref)].insert(kAssumptionsConstraint); + } + } + // For each variables, list the constraints that just enforce a lower bound // (resp. upper bound) on that variable. If all the constraints in which a // variable appear are in the same direction, then we can usually fix a diff --git a/ortools/sat/synchronization.cc b/ortools/sat/synchronization.cc index d730d3f80f..4adc1f5e1d 100644 --- a/ortools/sat/synchronization.cc +++ b/ortools/sat/synchronization.cc @@ -275,6 +275,14 @@ void SharedResponseManager::NotifyThatImprovingProblemIsInfeasible( if (log_updates_) LogNewSatSolution("Done", wall_timer_.Get(), worker_info); } +void SharedResponseManager::AddUnsatCore(const std::vector& core) { + absl::MutexLock mutex_lock(&mutex_); + best_response_.clear_sufficient_assumptions_for_infeasibility(); + for (const int ref : core) { + best_response_.add_sufficient_assumptions_for_infeasibility(ref); + } +} + IntegerValue SharedResponseManager::GetInnerObjectiveLowerBound() { absl::MutexLock mutex_lock(&mutex_); return IntegerValue(inner_objective_lower_bound_); diff --git a/ortools/sat/synchronization.h b/ortools/sat/synchronization.h index 349b15f2db..4ae63d7c85 100644 --- a/ortools/sat/synchronization.h +++ b/ortools/sat/synchronization.h @@ -248,6 +248,10 @@ class SharedResponseManager { // reported. We check for this case in NewSolution(). void NotifyThatImprovingProblemIsInfeasible(const std::string& worker_info); + // Adds to the shared response a subset of assumptions that are enough to + // make the problem infeasible. + void AddUnsatCore(const std::vector& core); + // Sets the statistics in the response to the one of the solver inside the // given in-memory model. This does nothing if the model is nullptr. //