diff --git a/makefiles/Makefile.gen.mk b/makefiles/Makefile.gen.mk index 59e181e33f..84ca30a589 100644 --- a/makefiles/Makefile.gen.mk +++ b/makefiles/Makefile.gen.mk @@ -1108,6 +1108,7 @@ SAT_DEPS = \ $(SRC_DIR)/ortools/sat/drat_proof_handler.h \ $(SRC_DIR)/ortools/sat/drat_writer.h \ $(SRC_DIR)/ortools/sat/encoding.h \ + $(SRC_DIR)/ortools/sat/feasibility_pump.h \ $(SRC_DIR)/ortools/sat/implied_bounds.h \ $(SRC_DIR)/ortools/sat/integer_expr.h \ $(SRC_DIR)/ortools/sat/integer.h \ @@ -1173,6 +1174,7 @@ SAT_LIB_OBJS = \ $(OBJ_DIR)/sat/drat_proof_handler.$O \ $(OBJ_DIR)/sat/drat_writer.$O \ $(OBJ_DIR)/sat/encoding.$O \ + $(OBJ_DIR)/sat/feasibility_pump.$O \ $(OBJ_DIR)/sat/implied_bounds.$O \ $(OBJ_DIR)/sat/integer.$O \ $(OBJ_DIR)/sat/integer_expr.$O \ @@ -1716,6 +1718,23 @@ objs/sat/encoding.$O: ortools/sat/encoding.cc ortools/sat/encoding.h \ ortools/util/integer_pq.h | $(OBJ_DIR)/sat $(CCC) $(CFLAGS) -c $(SRC_DIR)$Sortools$Ssat$Sencoding.cc $(OBJ_OUT)$(OBJ_DIR)$Ssat$Sencoding.$O +objs/sat/feasibility_pump.$O: ortools/sat/feasibility_pump.cc ortools/sat/feasibility_pump.h \ + ortools/base/int_type.h ortools/base/macros.h \ + ortools/base/integral_types.h ortools/base/logging.h \ + ortools/gen/ortools/sat/boolean_problem.pb.h ortools/sat/pb_constraint.h \ + ortools/base/int_type_indexed_vector.h ortools/sat/model.h \ + ortools/base/map_util.h ortools/base/typeid.h ortools/sat/sat_base.h \ + ortools/util/bitset.h ortools/gen/ortools/sat/sat_parameters.pb.h \ + ortools/util/stats.h ortools/base/timer.h ortools/base/basictypes.h \ + ortools/sat/sat_solver.h ortools/base/hash.h ortools/sat/clause.h \ + ortools/sat/drat_proof_handler.h ortools/sat/drat_checker.h \ + ortools/sat/drat_writer.h ortools/base/file.h \ + ortools/util/random_engine.h ortools/util/time_limit.h \ + ortools/base/commandlineflags.h ortools/util/running_stat.h \ + ortools/sat/restart.h ortools/sat/sat_decision.h \ + ortools/util/integer_pq.h | $(OBJ_DIR)/sat + $(CCC) $(CFLAGS) -c $(SRC_DIR)$Sortools$Ssat$Sfeasibility_pump.cc $(OBJ_OUT)$(OBJ_DIR)$Ssat$Sfeasibility_pump.$O + objs/sat/implied_bounds.$O: ortools/sat/implied_bounds.cc \ ortools/sat/implied_bounds.h ortools/base/int_type.h \ ortools/base/macros.h ortools/base/int_type_indexed_vector.h \ @@ -4100,4 +4119,3 @@ $(GEN_DIR)/ortools/constraint_solver/solver_parameters.pb.h: \ $(OBJ_DIR)/constraint_solver/solver_parameters.pb.$O: \ $(GEN_DIR)/ortools/constraint_solver/solver_parameters.pb.cc | $(OBJ_DIR)/constraint_solver $(CCC) $(CFLAGS) -c $(GEN_PATH)$Sortools$Sconstraint_solver$Ssolver_parameters.pb.cc $(OBJ_OUT)$(OBJ_DIR)$Sconstraint_solver$Ssolver_parameters.pb.$O - diff --git a/ortools/sat/BUILD b/ortools/sat/BUILD index 0e21bb46a4..83a81d6727 100644 --- a/ortools/sat/BUILD +++ b/ortools/sat/BUILD @@ -170,6 +170,7 @@ cc_library( ":cp_model_utils", ":drat_checker", ":drat_proof_handler", + ":feasibility_pump", ":integer", ":integer_expr", ":integer_search", @@ -1132,6 +1133,27 @@ cc_library( ], ) +cc_library( + name = "feasibility_pump", + srcs = ["feasibility_pump.cc"], + hdrs = ["feasibility_pump.h"], + visibility = ["//visibility:public"], + deps = [ + ":cp_model_loader", + ":cp_model_cc_proto", + ":integer", + ":linear_constraint", + ":synchronization", + ":util", + "//ortools/base", + "//ortools/glop:revised_simplex", + "//ortools/lp_data", + "//ortools/lp_data:base", + "//ortools/lp_data:lp_data_utils", + "//ortools/util:saturated_arithmetic", + ], +) + cc_library( name = "rins", srcs = ["rins.cc"], diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc index efbd4804ac..31dfba0286 100644 --- a/ortools/sat/cp_model_lns.cc +++ b/ortools/sat/cp_model_lns.cc @@ -13,6 +13,7 @@ #include "ortools/sat/cp_model_lns.h" +#include #include #include #include @@ -586,6 +587,10 @@ Neighborhood SchedulingTimeWindowNeighborhoodGenerator::Generate( } bool RelaxationInducedNeighborhoodGenerator::ReadyToGenerate() const { + if (incomplete_solutions_ != nullptr) { + return incomplete_solutions_->HasNewSolution(); + } + if (response_manager_ != nullptr) { if (response_manager_->SolutionsRepository().NumSolutions() == 0) { return false; @@ -618,6 +623,10 @@ Neighborhood RelaxationInducedNeighborhoodGenerator::Generate( (relaxation_solutions_ != nullptr && relaxation_solutions_->NumSolutions() > 0); + const bool incomplete_solution_available = + (incomplete_solutions_ != nullptr && + incomplete_solutions_->HasNewSolution()); + if (!lp_solution_available && !relaxation_solution_available) { return neighborhood; } @@ -632,14 +641,19 @@ Neighborhood RelaxationInducedNeighborhoodGenerator::Generate( ? random_bool(*random) : lp_solution_available; if (use_lp_relaxation) { - rins_neighborhood = GetRINSNeighborhood(response_manager_, - /*relaxation_solutions=*/nullptr, - lp_solutions_, random); + rins_neighborhood = + GetRINSNeighborhood(response_manager_, + /*relaxation_solutions=*/nullptr, lp_solutions_, + incomplete_solutions_, random); + neighborhood.source_info = + incomplete_solution_available ? "incomplete" : "lp"; } else { CHECK(relaxation_solution_available); - rins_neighborhood = - GetRINSNeighborhood(response_manager_, relaxation_solutions_, - /*lp_solutions=*/nullptr, random); + rins_neighborhood = GetRINSNeighborhood( + response_manager_, relaxation_solutions_, + /*lp_solutions=*/nullptr, incomplete_solutions_, random); + neighborhood.source_info = + incomplete_solution_available ? "incomplete" : "relaxation"; } if (rins_neighborhood.fixed_vars.empty() && diff --git a/ortools/sat/cp_model_lns.h b/ortools/sat/cp_model_lns.h index 40c418628c..14daf4c6de 100644 --- a/ortools/sat/cp_model_lns.h +++ b/ortools/sat/cp_model_lns.h @@ -49,6 +49,10 @@ struct Neighborhood { // TODO(user): Make sure that the id is unique for each generated // neighborhood for each generator. int64 id = 0; + + // Used for identifying the source of the neighborhood if it is generated + // using solution repositories. + std::string source_info = ""; }; // Contains pre-computed information about a given CpModelProto that is meant @@ -404,27 +408,36 @@ class SchedulingTimeWindowNeighborhoodGenerator : public NeighborhoodGenerator { double difficulty, random_engine_t* random) final; }; -// Generates a neighborhood by fixing the variables who have same solution value -// as their linear relaxation. This was published in "Exploring relaxation -// induced neighborhoods to improve MIP solutions" 2004 by E. Danna et. +// Generates a neighborhood by fixing the variables to solutions reported in +// various repositories. This is inspired from RINS published in "Exploring +// relaxation induced neighborhoods to improve MIP solutions" 2004 by E. Danna +// et. // -// If use_only_relaxation_values is true, this generates a neighborhood using -// only the linear/general relaxation values. The domain of the variables are -// reduced to the integer values around their lp solution/relaxation solution -// values. This was published in "RENS – The Relaxation Enforced Neighborhood" -// 2009 by Timo Berthold. +// If incomplete_solutions is provided, this generates a neighborhood by fixing +// the variable values to a solution in the SharedIncompleteSolutionManager and +// ignores the other repositories. +// +// Otherwise, if response_manager is not provided, this generates a neighborhood +// using only the linear/general relaxation values. The domain of the variables +// are reduced to the integer values around their lp solution/relaxation +// solution values. This was published in "RENS – The Relaxation Enforced +// Neighborhood" 2009 by Timo Berthold. class RelaxationInducedNeighborhoodGenerator : public NeighborhoodGenerator { public: explicit RelaxationInducedNeighborhoodGenerator( NeighborhoodGeneratorHelper const* helper, const SharedResponseManager* response_manager, const SharedRelaxationSolutionRepository* relaxation_solutions, - const SharedLPSolutionRepository* lp_solutions, const std::string& name) + const SharedLPSolutionRepository* lp_solutions, + SharedIncompleteSolutionManager* incomplete_solutions, + const std::string& name) : NeighborhoodGenerator(name, helper), response_manager_(response_manager), relaxation_solutions_(relaxation_solutions), - lp_solutions_(lp_solutions) { - CHECK(lp_solutions_ != nullptr || relaxation_solutions_ != nullptr); + lp_solutions_(lp_solutions), + incomplete_solutions_(incomplete_solutions) { + CHECK(lp_solutions_ != nullptr || relaxation_solutions_ != nullptr || + incomplete_solutions != nullptr); } // Both initial solution and difficulty values are ignored. @@ -434,9 +447,11 @@ class RelaxationInducedNeighborhoodGenerator : public NeighborhoodGenerator { // Returns true if the required solutions are available. bool ReadyToGenerate() const override; + private: const SharedResponseManager* response_manager_; const SharedRelaxationSolutionRepository* relaxation_solutions_; const SharedLPSolutionRepository* lp_solutions_; + SharedIncompleteSolutionManager* incomplete_solutions_; }; // Generates a relaxation of the original model by removing a consecutive span diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 8d078433fc..7e44e309f1 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -25,6 +25,7 @@ #include #include "ortools/base/cleanup.h" +#include "ortools/sat/feasibility_pump.h" #if !defined(__PORTABLE_PLATFORM__) #include "absl/synchronization/notification.h" @@ -551,9 +552,8 @@ void TryToAddCutGenerators(const CpModelProto& model_proto, } // namespace -// Adds one LinearProgrammingConstraint per connected component of the model. -IntegerVariable AddLPConstraints(const CpModelProto& model_proto, - int linearization_level, Model* m) { +LinearRelaxation ComputeLinearRelaxation(const CpModelProto& model_proto, + int linearization_level, Model* m) { LinearRelaxation relaxation; // Linearize the constraints. @@ -670,6 +670,14 @@ IntegerVariable AddLPConstraints(const CpModelProto& model_proto, VLOG(3) << relaxation.linear_constraints.size() << " constraints in the LP relaxation."; VLOG(3) << relaxation.cut_generators.size() << " cuts generators."; + return relaxation; +} + +// Adds one LinearProgrammingConstraint per connected component of the model. +IntegerVariable AddLPConstraints(const CpModelProto& model_proto, + int linearization_level, Model* m) { + const LinearRelaxation relaxation = + ComputeLinearRelaxation(model_proto, linearization_level, m); // The bipartite graph of LP constraints might be disconnected: // make a partition of the variables into connected components. @@ -719,6 +727,7 @@ IntegerVariable AddLPConstraints(const CpModelProto& model_proto, // as much as possible the objective bound by using any bounds the LP give // us on one of its components. This is critical on the zephyrus problems for // instance. + auto* mapping = m->GetOrCreate(); for (int i = 0; i < model_proto.objective().coeffs_size(); ++i) { const IntegerVariable var = mapping->Integer(model_proto.objective().vars(i)); @@ -1044,7 +1053,7 @@ void RegisterObjectiveBestBoundExport( // Registers a callback to import new objective bounds. It will be called each // time the search main loop is back to level zero. Note that it the presence of -// assumptions, this will not happend until the set of assumptions is changed. +// assumptions, this will not happen until the set of assumptions is changed. void RegisterObjectiveBoundsImport( SharedResponseManager* shared_response_manager, Model* model) { auto* solver = model->GetOrCreate(); @@ -1098,6 +1107,55 @@ void RegisterObjectiveBoundsImport( import_objective_bounds); } +void LoadFeasibilityPump(const CpModelProto& model_proto, + SharedResponseManager* shared_response_manager, + Model* model) { + CHECK(shared_response_manager != nullptr); + + auto* sat_solver = model->GetOrCreate(); + + // Simple function for the few places where we do "return unsat()". + const auto unsat = [shared_response_manager, sat_solver, model] { + sat_solver->NotifyThatModelIsUnsat(); + shared_response_manager->NotifyThatImprovingProblemIsInfeasible( + absl::StrCat(model->GetOrCreate()->worker_name, + " [loading]")); + }; + + auto* mapping = model->GetOrCreate(); + const SatParameters& parameters = *(model->GetOrCreate()); + const bool view_all_booleans_as_integers = + (parameters.linearization_level() >= 2) || + (parameters.search_branching() == SatParameters::FIXED_SEARCH && + model_proto.search_strategy().empty()); + mapping->CreateVariables(model_proto, view_all_booleans_as_integers, model); + + // Check the model is still feasible before continuing. + if (sat_solver->IsModelUnsat()) return unsat(); + + // We don't load any constraint here. + + if (parameters.linearization_level() == 0) return; + + // Add linear constraints to Feasibility Pump. + const LinearRelaxation relaxation = ComputeLinearRelaxation( + model_proto, parameters.linearization_level(), model); + const int num_lp_constraints = relaxation.linear_constraints.size(); + auto* feasibility_pump = model->GetOrCreate(); + for (int i = 0; i < num_lp_constraints; i++) { + feasibility_pump->AddLinearConstraint(relaxation.linear_constraints[i]); + } + + if (model_proto.has_objective()) { + for (int i = 0; i < model_proto.objective().coeffs_size(); ++i) { + const IntegerVariable var = + mapping->Integer(model_proto.objective().vars(i)); + const int64 coeff = model_proto.objective().coeffs(i); + feasibility_pump->SetObjectiveCoefficient(var, IntegerValue(coeff)); + } + } +} + // Loads a CpModelProto inside the given model. // This should only be called once on a given 'Model' class. // @@ -1836,6 +1894,7 @@ struct SharedClasses { SharedResponseManager* response; SharedRelaxationSolutionRepository* relaxation_solutions; SharedLPSolutionRepository* lp_solutions; + SharedIncompleteSolutionManager* incomplete_solutions; bool SearchIsDone() { if (response->ProblemIsSolved()) return true; @@ -1877,6 +1936,11 @@ class FullProblemSolver : public SubSolver { local_model_->Register(shared->lp_solutions); } + if (shared->incomplete_solutions != nullptr) { + local_model_->Register( + shared->incomplete_solutions); + } + // Level zero variable bounds sharing. if (shared_->bounds != nullptr) { RegisterVariableBoundsLevelZeroExport( @@ -1982,6 +2046,118 @@ class FullProblemSolver : public SubSolver { bool previous_task_is_completed_ ABSL_GUARDED_BY(mutex_) = true; }; +class FeasibilityPumpSolver : public SubSolver { + public: + FeasibilityPumpSolver(int id, const SatParameters& local_parameters, + SharedClasses* shared) + : SubSolver(id, "feasibility_pump"), + shared_(shared), + local_model_(absl::make_unique()) { + // Setup the local model parameters and time limit. + local_model_->Add(NewSatParameters(local_parameters)); + shared_->time_limit->UpdateLocalLimit( + local_model_->GetOrCreate()); + + // Stores info that will be used for logs in the local model. + WorkerInfo* worker_info = local_model_->GetOrCreate(); + worker_info->worker_name = "feasibility_pump"; + worker_info->worker_id = id; + + if (shared->response != nullptr) { + local_model_->Register(shared->response); + } + + if (shared->relaxation_solutions != nullptr) { + local_model_->Register( + shared->relaxation_solutions); + } + + if (shared->lp_solutions != nullptr) { + local_model_->Register(shared->lp_solutions); + } + + if (shared->incomplete_solutions != nullptr) { + local_model_->Register( + shared->incomplete_solutions); + } + + // Level zero variable bounds sharing. + if (shared_->bounds != nullptr) { + RegisterVariableBoundsLevelZeroImport( + *shared_->model_proto, shared_->bounds, local_model_.get()); + } + } + + bool TaskIsAvailable() override { + if (shared_->SearchIsDone()) return false; + absl::MutexLock mutex_lock(&mutex_); + return previous_task_is_completed_; + } + + std::function GenerateTask(int64 task_id) override { + return [this]() { + { + absl::MutexLock mutex_lock(&mutex_); + if (!previous_task_is_completed_) return; + previous_task_is_completed_ = false; + } + { + absl::MutexLock mutex_lock(&mutex_); + if (solving_first_chunk_) { + LoadFeasibilityPump(*shared_->model_proto, shared_->response, + local_model_.get()); + solving_first_chunk_ = false; + // Abort first chunk and allow to schedule the next. + previous_task_is_completed_ = true; + return; + } + } + + auto* time_limit = local_model_->GetOrCreate(); + const double saved_dtime = time_limit->GetElapsedDeterministicTime(); + auto* feasibility_pump = local_model_->Mutable(); + feasibility_pump->Solve(); + + { + absl::MutexLock mutex_lock(&mutex_); + deterministic_time_since_last_synchronize_ += + time_limit->GetElapsedDeterministicTime() - saved_dtime; + } + + // Abort if the problem is solved. + if (shared_->SearchIsDone()) { + shared_->time_limit->Stop(); + return; + } + + absl::MutexLock mutex_lock(&mutex_); + previous_task_is_completed_ = true; + }; + } + + void Synchronize() override { + absl::MutexLock mutex_lock(&mutex_); + deterministic_time_ += deterministic_time_since_last_synchronize_; + shared_->time_limit->AdvanceDeterministicTime( + deterministic_time_since_last_synchronize_); + deterministic_time_since_last_synchronize_ = 0.0; + } + + private: + SharedClasses* shared_; + std::unique_ptr local_model_; + + absl::Mutex mutex_; + + // The first chunk is special. It is the one in which we load the linear + // constraints. + bool solving_first_chunk_ ABSL_GUARDED_BY(mutex_) = true; + + double deterministic_time_since_last_synchronize_ ABSL_GUARDED_BY(mutex_) = + 0.0; + bool previous_task_is_completed_ ABSL_GUARDED_BY(mutex_) = true; +}; + // A Subsolver that generate LNS solve from a given neighborhood. class LnsSolver : public SubSolver { public: @@ -2058,9 +2234,13 @@ class LnsSolver : public SubSolver { const double fully_solved_proportion = static_cast(generator_->num_fully_solved_calls()) / std::max(int64{1}, generator_->num_calls()); + std::string source_info = name(); + if (!neighborhood.source_info.empty()) { + absl::StrAppend(&source_info, "_", neighborhood.source_info); + } const std::string solution_info = absl::StrFormat( - "%s(d=%0.2f s=%i t=%0.2f p=%0.2f)", name(), data.difficulty, task_id, - data.deterministic_limit, fully_solved_proportion); + "%s(d=%0.2f s=%i t=%0.2f p=%0.2f)", source_info, data.difficulty, + task_id, data.deterministic_limit, fully_solved_proportion); SatParameters local_params(parameters_); local_params.set_max_deterministic_time(data.deterministic_limit); @@ -2283,6 +2463,14 @@ void SolveCpModelParallel(const CpModelProto& model_proto, /*num_solutions_to_keep=*/10); global_model->Register(shared_lp_solutions.get()); + std::unique_ptr shared_incomplete_solutions; + if (global_model->GetOrCreate()->use_feasibility_pump()) { + shared_incomplete_solutions = + absl::make_unique(); + global_model->Register( + shared_incomplete_solutions.get()); + } + SharedClasses shared; shared.model_proto = &model_proto; shared.wall_timer = wall_timer; @@ -2291,6 +2479,7 @@ void SolveCpModelParallel(const CpModelProto& model_proto, shared.response = shared_response_manager; shared.relaxation_solutions = shared_relaxation_solutions.get(); shared.lp_solutions = shared_lp_solutions.get(); + shared.incomplete_solutions = shared_incomplete_solutions.get(); // The list of all the SubSolver that will be used in this parallel search. std::vector> subsolvers; @@ -2353,6 +2542,12 @@ void SolveCpModelParallel(const CpModelProto& model_proto, } } + // Add FeasibilityPumpSolver if enabled. + if (parameters.use_feasibility_pump() && !parameters.interleave_search()) { + subsolvers.push_back(absl::make_unique( + /*id=*/subsolvers.size(), parameters, &shared)); + } + // Add LNS SubSolver(s). // Add the NeighborhoodGeneratorHelper as a special subsolver so that its @@ -2418,7 +2613,8 @@ void SolveCpModelParallel(const CpModelProto& model_proto, /*id=*/subsolvers.size(), absl::make_unique( helper, shared.response, shared.relaxation_solutions, - shared.lp_solutions, absl::StrCat("rins_lns_", strategy_name)), + shared.lp_solutions, /*incomplete_solutions=*/nullptr, + absl::StrCat("rins_lns_", strategy_name)), local_params, helper, &shared)); // RENS. @@ -2426,9 +2622,11 @@ void SolveCpModelParallel(const CpModelProto& model_proto, /*id=*/subsolvers.size(), absl::make_unique( helper, /*respons_manager=*/nullptr, shared.relaxation_solutions, - shared.lp_solutions, absl::StrCat("rens_lns_", strategy_name)), + shared.lp_solutions, shared.incomplete_solutions, + absl::StrCat("rens_lns_", strategy_name)), local_params, helper, &shared)); } + if (parameters.use_relaxation_lns()) { subsolvers.push_back(absl::make_unique( /*id=*/subsolvers.size(), diff --git a/ortools/sat/feasibility_pump.cc b/ortools/sat/feasibility_pump.cc new file mode 100644 index 0000000000..67e6050a4b --- /dev/null +++ b/ortools/sat/feasibility_pump.cc @@ -0,0 +1,441 @@ +// Copyright 2010-2018 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "ortools/sat/feasibility_pump.h" + +#include + +#include "ortools/base/integral_types.h" +#include "ortools/lp_data/lp_types.h" +#include "ortools/sat/cp_model.pb.h" +#include "ortools/util/saturated_arithmetic.h" + +namespace operations_research { +namespace sat { + +using glop::ColIndex; +using glop::Fractional; +using glop::RowIndex; + +const double FeasibilityPump::kCpEpsilon = 1e-4; + +FeasibilityPump::FeasibilityPump(Model* model) + : sat_parameters_(*(model->GetOrCreate())), + time_limit_(model->GetOrCreate()), + integer_trail_(model->GetOrCreate()), + trail_(model->GetOrCreate()), + integer_encoder_(model->GetOrCreate()), + incomplete_solutions_(model->Mutable()), + 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. + parameters.set_use_dual_simplex(false); + parameters.set_max_number_of_iterations(2000); + simplex_.SetParameters(parameters); + lp_data_.Clear(); + integer_lp_.clear(); +} + +FeasibilityPump::~FeasibilityPump() { + VLOG(1) << "Feasibility Pump Total number of simplex iterations: " + << total_num_simplex_iterations_; +} + +void FeasibilityPump::AddLinearConstraint(const LinearConstraint& ct) { + // We still create the mirror variable right away though. + for (const IntegerVariable var : ct.vars) { + GetOrCreateMirrorVariable(PositiveVariable(var)); + } + + integer_lp_.push_back(LinearConstraintInternal()); + LinearConstraintInternal& new_ct = integer_lp_.back(); + new_ct.lb = ct.lb; + new_ct.ub = ct.ub; + const int size = ct.vars.size(); + CHECK_LE(ct.lb, ct.ub); + for (int i = 0; i < size; ++i) { + // We only use positive variable inside this class. + IntegerVariable var = ct.vars[i]; + IntegerValue coeff = ct.coeffs[i]; + if (!VariableIsPositive(var)) { + var = NegationOf(var); + coeff = -coeff; + } + new_ct.terms.push_back({GetOrCreateMirrorVariable(var), coeff}); + } + // Important to keep lp_data_ "clean". + std::sort(new_ct.terms.begin(), new_ct.terms.end()); +} + +void FeasibilityPump::SetObjectiveCoefficient(IntegerVariable ivar, + IntegerValue coeff) { + objective_is_defined_ = true; + const IntegerVariable pos_var = + VariableIsPositive(ivar) ? ivar : NegationOf(ivar); + if (ivar != pos_var) coeff = -coeff; + + const auto it = mirror_lp_variable_.find(pos_var); + if (it == mirror_lp_variable_.end()) return; + const glop::ColIndex col = it->second; + integer_objective_.push_back({col, coeff}); + objective_infinity_norm_ = + std::max(objective_infinity_norm_, IntTypeAbs(coeff)); +} + +glop::ColIndex FeasibilityPump::GetOrCreateMirrorVariable( + IntegerVariable positive_variable) { + DCHECK(VariableIsPositive(positive_variable)); + + const auto it = mirror_lp_variable_.find(positive_variable); + if (it == mirror_lp_variable_.end()) { + const int model_var = + mapping_->GetProtoVariableFromIntegerVariable(positive_variable); + model_vars_size_ = std::max(model_vars_size_, model_var + 1); + + const glop::ColIndex col(integer_variables_.size()); + mirror_lp_variable_[positive_variable] = col; + integer_variables_.push_back(positive_variable); + lp_solution_.push_back(std::numeric_limits::infinity()); + integer_solution_.push_back(kMaxIntegerValue.value()); + + return col; + } + return it->second; +} + +void FeasibilityPump::PrintStats() { + if (lp_solution_is_set_) { + VLOG(2) << "Fractionality: " << lp_solution_fractionality_; + } else { + VLOG(2) << "Fractionality: NA"; + VLOG(2) << "simplex status: " << simplex_.GetProblemStatus(); + } + + if (integer_solution_is_set_) { + VLOG(2) << "#Infeasible const: " << num_infeasible_constraints_; + VLOG(2) << "Infeasibility: " << integer_solution_infeasibility_; + } else { + VLOG(2) << "Infeasibility: NA"; + } +} + +void FeasibilityPump::Solve() { + if (lp_data_.num_variables() == 0) { + InitializeWorkingLP(); + } + UpdateBoundsOfLpVariables(); + lp_solution_is_set_ = false; + integer_solution_is_set_ = false; + + // Restore the original objective + for (ColIndex col(0); col < lp_data_.num_variables(); ++col) { + lp_data_.SetObjectiveCoefficient(col, 0.0); + } + for (const auto& term : integer_objective_) { + lp_data_.SetObjectiveCoefficient(term.first, ToDouble(term.second)); + } + + // TODO(user): Add cycle detection. + mixing_factor_ = 1.0; + for (int i = 0; i < max_fp_iterations_; ++i) { + L1DistanceMinimize(); + if (!SolveLp()) break; + if (lp_solution_is_integer_) break; + SimpleRounding(); + // We don't end this loop if the integer solutions is feasible in hope to + // get better solution. + if (integer_solution_is_feasible_) MaybePushToRepo(); + } + + PrintStats(); + MaybePushToRepo(); +} + +void FeasibilityPump::MaybePushToRepo() { + if (incomplete_solutions_ == nullptr) return; + + std::vector lp_solution(model_vars_size_, + std::numeric_limits::infinity()); + // TODO(user): Consider adding solutions that have low fractionality. + if (lp_solution_is_integer_) { + // Fill the solution using LP solution values. + for (const IntegerVariable positive_var : integer_variables_) { + const int model_var = + mapping_->GetProtoVariableFromIntegerVariable(positive_var); + if (model_var >= 0 && model_var < model_vars_size_) { + lp_solution[model_var] = GetLPSolutionValue(positive_var); + } + } + incomplete_solutions_->AddNewSolution(lp_solution); + } + + if (integer_solution_is_feasible_) { + // Fill the solution using Integer solution values. + for (const IntegerVariable positive_var : integer_variables_) { + const int model_var = + mapping_->GetProtoVariableFromIntegerVariable(positive_var); + if (model_var >= 0 && model_var < model_vars_size_) { + lp_solution[model_var] = GetIntegerSolutionValue(positive_var); + } + } + incomplete_solutions_->AddNewSolution(lp_solution); + } +} + +// ---------------------------------------------------------------- +// -------------------LPSolving------------------------------------ +// ---------------------------------------------------------------- + +void FeasibilityPump::InitializeWorkingLP() { + lp_data_.Clear(); + // Create variables. + for (int i = 0; i < integer_variables_.size(); ++i) { + CHECK_EQ(glop::ColIndex(i), lp_data_.CreateNewVariable()); + lp_data_.SetVariableType(glop::ColIndex(i), + glop::LinearProgram::VariableType::INTEGER); + } + + // Add constraints. + for (const LinearConstraintInternal& ct : integer_lp_) { + const ConstraintIndex row = lp_data_.CreateNewConstraint(); + lp_data_.SetConstraintBounds(row, ToDouble(ct.lb), ToDouble(ct.ub)); + for (const auto& term : ct.terms) { + lp_data_.SetCoefficient(row, term.first, ToDouble(term.second)); + } + } + + // Add objective. + for (const auto& term : integer_objective_) { + lp_data_.SetObjectiveCoefficient(term.first, ToDouble(term.second)); + } + + const int num_vars = integer_variables_.size(); + for (int i = 0; i < num_vars; i++) { + const IntegerVariable cp_var = integer_variables_[i]; + const double lb = ToDouble(integer_trail_->LevelZeroLowerBound(cp_var)); + const double ub = ToDouble(integer_trail_->LevelZeroUpperBound(cp_var)); + lp_data_.SetVariableBounds(glop::ColIndex(i), lb, ub); + } + + objective_normalization_factor_ = 0.0; + glop::ColIndexVector integer_variables; + const ColIndex num_cols = lp_data_.num_variables(); + for (ColIndex col : lp_data_.IntegerVariablesList()) { + if (!lp_data_.IsVariableBinary(col)) { + integer_variables.push_back(col); + } + + // The aim of this normalization value is to compute a coefficient of the + // d_i variables that should be minimized. + objective_normalization_factor_ += + std::abs(lp_data_.GetObjectiveCoefficientForMinimizationVersion(col)); + } + CHECK_GT(lp_data_.IntegerVariablesList().size(), 0); + objective_normalization_factor_ = + objective_normalization_factor_ / lp_data_.IntegerVariablesList().size(); + + if (!integer_variables.empty()) { + // Update the LpProblem with norm variables and constraints. + norm_variables_.assign(num_cols, ColIndex(-1)); + norm_lhs_constraints_.assign(num_cols, RowIndex(-1)); + norm_rhs_constraints_.assign(num_cols, RowIndex(-1)); + // For each integer non-binary variable x_i we introduce one new variable + // d_i subject to two new constraints: + // d_i - x_i >= -round(x'_i) + // d_i + x_i >= +round(x'_i) + // That's round(x'_i) - d_i <= x_i <= round(x'_i) + d_i, where d_i is an + // unbounded non-negative, and x'_i is the value of variable i from the + // previous solution obtained during the projection step. Consequently + // coefficients of the constraints are set here, but bounds of the + // constraints are updated at each iteration of the feasibility pump. Also + // coefficients of the objective are set here: x_i's are not present in the + // objective (i.e., coefficients set to 0.0), and d_i's are present in the + // objective with coefficients set to 1.0. + // Note that the treatment of integer non-binary variables is different + // from the treatment of binary variables. Binary variables do not impose + // any extra variables, nor extra constraints, but their objective + // coefficients are changed in the linear projection steps. + for (const ColIndex col : integer_variables) { + const ColIndex norm_variable = lp_data_.CreateNewVariable(); + norm_variables_[col] = norm_variable; + lp_data_.SetVariableBounds(norm_variable, 0.0, glop::kInfinity); + const RowIndex row_a = lp_data_.CreateNewConstraint(); + norm_lhs_constraints_[col] = row_a; + lp_data_.SetCoefficient(row_a, norm_variable, 1.0); + lp_data_.SetCoefficient(row_a, col, -1.0); + const RowIndex row_b = lp_data_.CreateNewConstraint(); + norm_rhs_constraints_[col] = row_b; + lp_data_.SetCoefficient(row_b, norm_variable, 1.0); + lp_data_.SetCoefficient(row_b, col, 1.0); + } + } + + // TODO(user): Try to add scaling. + lp_data_.AddSlackVariablesWhereNecessary( + /*detect_integer_constraints=*/false); +} + +void FeasibilityPump::L1DistanceMinimize() { + std::vector new_obj_coeffs(lp_data_.num_variables().value(), 0.0); + + // Set the original subobjective. The coefficients are scaled by mixing factor + // and the offset remains at 0 (because it does not affect the solution). + const ColIndex num_cols(lp_data_.objective_coefficients().size()); + for (ColIndex col(0); col < num_cols; ++col) { + new_obj_coeffs[col.value()] = + mixing_factor_ * lp_data_.objective_coefficients()[col]; + } + + // Set the norm subobjective. The coefficients are scaled by 1 - mixing factor + // and the offset remains at 0 (because it does not affect the solution). + for (const ColIndex col : lp_data_.IntegerVariablesList()) { + if (lp_data_.IsVariableBinary(col)) { + const Fractional objective_coefficient = + mixing_factor_ * lp_data_.objective_coefficients()[col] + + (1 - mixing_factor_) * objective_normalization_factor_ * + (1 - 2 * integer_solution_[col.value()]); + new_obj_coeffs[col.value()] = objective_coefficient; + } else { // The variable is integer. + // Update the bounds of the constraints added in + // InitializeIntegerVariables() (see there for more details): + // d_i - x_i >= -round(x'_i) + // d_i + x_i >= +round(x'_i) + + // TODO(user): We change both the objective and the bounds, thus + // breaking the incrementality. Handle integer variables differently, + // e.g., intensify rounding, or use soft fixing from: Fischetti, Lodi, + // "Local Branching", Math Program Ser B 98:23-47 (2003). + const Fractional objective_coefficient = + (1 - mixing_factor_) * objective_normalization_factor_; + new_obj_coeffs[norm_variables_[col].value()] = objective_coefficient; + // At this point, constraint bounds have already been transformed into + // bounds of slack variables. Instead of updating the constraints, we need + // to update the slack variables corresponding to them. + const ColIndex norm_a_slack_variable = + lp_data_.GetSlackVariable(norm_lhs_constraints_[col]); + lp_data_.SetVariableBounds(norm_a_slack_variable, -glop::kInfinity, + integer_solution_[col.value()]); + const ColIndex norm_b_slack_variable = + lp_data_.GetSlackVariable(norm_rhs_constraints_[col]); + lp_data_.SetVariableBounds(norm_b_slack_variable, -glop::kInfinity, + -integer_solution_[col.value()]); + } + } + for (ColIndex col(0); col < lp_data_.num_variables(); ++col) { + lp_data_.SetObjectiveCoefficient(col, new_obj_coeffs[col.value()]); + } + // TODO(user): Tune this or expose as parameter. + mixing_factor_ *= 0.8; +} + +bool FeasibilityPump::SolveLp() { + const int num_vars = integer_variables_.size(); + VLOG(3) << "LP relaxation: " << lp_data_.GetDimensionString() << "."; + + const auto status = simplex_.Solve(lp_data_, time_limit_); + total_num_simplex_iterations_ += simplex_.GetNumberOfIterations(); + if (!status.ok()) { + VLOG(1) << "The LP solver encountered an error: " << status.error_message(); + simplex_.ClearStateForNextSolve(); + return false; + } + + VLOG(3) << "simplex status: " << simplex_.GetProblemStatus(); + lp_solution_fractionality_ = 0.0; + if (simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL || + simplex_.GetProblemStatus() == glop::ProblemStatus::DUAL_FEASIBLE || + simplex_.GetProblemStatus() == glop::ProblemStatus::PRIMAL_FEASIBLE || + simplex_.GetProblemStatus() == glop::ProblemStatus::IMPRECISE) { + lp_solution_is_set_ = true; + for (int i = 0; i < num_vars; i++) { + const double value = simplex_.GetVariableValue(glop::ColIndex(i)); + lp_solution_[i] = value; + lp_solution_fractionality_ = std::max( + lp_solution_fractionality_, std::abs(value - std::round(value))); + } + + // Compute the objective value. + lp_objective_ = 0; + for (const auto& term : integer_objective_) { + lp_objective_ += lp_solution_[term.first.value()] * term.second.value(); + } + lp_solution_is_integer_ = lp_solution_fractionality_ < kCpEpsilon; + } + return true; +} + +void FeasibilityPump::UpdateBoundsOfLpVariables() { + const int num_vars = integer_variables_.size(); + for (int i = 0; i < num_vars; i++) { + const IntegerVariable cp_var = integer_variables_[i]; + const double lb = ToDouble(integer_trail_->LowerBound(cp_var)); + const double ub = ToDouble(integer_trail_->UpperBound(cp_var)); + lp_data_.SetVariableBounds(glop::ColIndex(i), lb, ub); + } +} + +double FeasibilityPump::GetLPSolutionValue(IntegerVariable variable) const { + return lp_solution_[gtl::FindOrDie(mirror_lp_variable_, variable).value()]; +} + +// ---------------------------------------------------------------- +// -------------------Rounding------------------------------------- +// ---------------------------------------------------------------- + +int64 FeasibilityPump::GetIntegerSolutionValue(IntegerVariable variable) const { + return integer_solution_[gtl::FindOrDie(mirror_lp_variable_, variable) + .value()]; +} + +void FeasibilityPump::SimpleRounding() { + if (!lp_solution_is_set_) return; + integer_solution_is_set_ = true; + for (int i = 0; i < lp_solution_.size(); ++i) { + integer_solution_[i] = static_cast(std::round(lp_solution_[i])); + } + // Compute the objective value. + integer_solution_objective_ = 0; + for (const auto& term : integer_objective_) { + integer_solution_objective_ += + integer_solution_[term.first.value()] * term.second.value(); + } + + integer_solution_is_feasible_ = true; + num_infeasible_constraints_ = 0; + integer_solution_infeasibility_ = 0; + for (glop::RowIndex i(0); i < integer_lp_.size(); ++i) { + int64 activity = 0; + for (const auto& term : integer_lp_[i].terms) { + const int64 prod = + CapProd(integer_solution_[term.first.value()], term.second.value()); + if (prod <= kint64min || prod >= kint64max) { + activity = prod; + break; + } + activity = CapAdd(activity, prod); + if (activity <= kint64min || activity >= kint64max) break; + } + if (activity > integer_lp_[i].ub || activity < integer_lp_[i].lb) { + integer_solution_is_feasible_ = false; + num_infeasible_constraints_++; + integer_solution_infeasibility_ = + std::max(integer_solution_infeasibility_, + std::max(activity - integer_lp_[i].ub.value(), + integer_lp_[i].lb.value() - activity)); + } + } +} + +} // namespace sat +} // namespace operations_research diff --git a/ortools/sat/feasibility_pump.h b/ortools/sat/feasibility_pump.h new file mode 100644 index 0000000000..a91685a9e3 --- /dev/null +++ b/ortools/sat/feasibility_pump.h @@ -0,0 +1,185 @@ +// Copyright 2010-2018 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#ifndef OR_TOOLS_SAT_FEASIBILITY_PUMP_H_ +#define OR_TOOLS_SAT_FEASIBILITY_PUMP_H_ + +#include "ortools/glop/revised_simplex.h" +#include "ortools/lp_data/lp_data.h" +#include "ortools/lp_data/lp_data_utils.h" +#include "ortools/lp_data/lp_types.h" +#include "ortools/sat/cp_model_loader.h" +#include "ortools/sat/integer.h" +#include "ortools/sat/linear_constraint.h" +#include "ortools/sat/synchronization.h" +#include "ortools/sat/util.h" + +namespace operations_research { +namespace sat { + +struct FeasibilityPumpLpSolution + : public gtl::ITIVector { + FeasibilityPumpLpSolution() {} +}; + +class FeasibilityPump { + public: + explicit FeasibilityPump(Model* model); + ~FeasibilityPump(); + + typedef glop::RowIndex ConstraintIndex; + + void SetMaxFPIterations(int max_iter) { + max_fp_iterations_ = std::max(1, max_iter); + } + + // Add a new linear constraint to this LP. + void AddLinearConstraint(const LinearConstraint& ct); + + // Set the coefficient of the variable in the objective. Calling it twice will + // overwrite the previous value. Note that this doesn't set the objective + // coefficient if the variable doesn't appear in any constraints. So this has + // to be called after all the constraints are added. + void SetObjectiveCoefficient(IntegerVariable ivar, IntegerValue coeff); + + // Returns the LP value of a variable in the current + // solution. These functions should only be called when HasSolution() is true. + bool HasLPSolution() const { return lp_solution_is_set_; } + double LPSolutionObjectiveValue() const { return lp_objective_; } + double GetLPSolutionValue(IntegerVariable variable) const; + bool LPSolutionIsInteger() const { return lp_solution_is_integer_; } + double LPSolutionFractionality() const { return lp_solution_fractionality_; } + + // Returns the Integer solution value of a variable in the current rounded + // solution. These functions should only be called when HasIntegerSolution() + // is true. + bool HasIntegerSolution() const { return integer_solution_is_set_; } + int64 IntegerSolutionObjectiveValue() const { + return integer_solution_objective_; + } + bool IntegerSolutionIsFeasible() const { + return integer_solution_is_feasible_; + } + int64 GetIntegerSolutionValue(IntegerVariable variable) const; + + void Solve(); + + private: + // Solve the LP, returns false if something went wrong in the LP solver. + bool SolveLp(); + + // Round the fractional LP solution values to nearest integer values. + void SimpleRounding(); + + // Loads the lp_data_. + void InitializeWorkingLP(); + + // Changes the LP objective and bounds of the norm constraints so the new + // objective also tries to minimize the distance to the rounded solution. + void L1DistanceMinimize(); + + // Stores the solutions in the shared repository. Stores LP solution if it is + // integer and stores the integer solution if it is feasible. + void MaybePushToRepo(); + + void PrintStats(); + + // Shortcut for an integer linear expression type. + using LinearExpression = std::vector>; + + // Gets or creates an LP variable that mirrors a model variable. + // The variable should be a positive reference. + glop::ColIndex GetOrCreateMirrorVariable(IntegerVariable positive_variable); + + // Updates the bounds of the LP variables from the CP bounds. + void UpdateBoundsOfLpVariables(); + + // This epsilon is related to the precision of the value returned by the LP + // once they have been scaled back into the CP domain. So for large domain or + // cost coefficient, we may have some issues. + static const double kCpEpsilon; + + // Initial problem in integer form. + // We always sort the inner vectors by increasing glop::ColIndex. + struct LinearConstraintInternal { + IntegerValue lb; + IntegerValue ub; + LinearExpression terms; + }; + LinearExpression integer_objective_; + IntegerValue objective_infinity_norm_ = IntegerValue(0); + double objective_normalization_factor_ = 0.0; + double mixing_factor_ = 1.0; + + gtl::ITIVector integer_lp_; + int model_vars_size_ = 0; + + // Underlying LP solver API. + glop::LinearProgram lp_data_; + glop::RevisedSimplex simplex_; + + glop::ColMapping norm_variables_; + glop::ColToRowMapping norm_lhs_constraints_; + glop::ColToRowMapping norm_rhs_constraints_; + + // Structures used for mirroring IntegerVariables inside the underlying LP + // solver: an integer variable var is mirrored by mirror_lp_variable_[var]. + // Note that these indices are dense in [0, mirror_lp_variable_.size()] so + // they can be used as vector indices. + std::vector integer_variables_; + absl::flat_hash_map mirror_lp_variable_; + + // We need to remember what to optimize if an objective is given, because + // then we will switch the objective between feasibility and optimization. + bool objective_is_defined_ = false; + + // Singletons from Model. + const SatParameters& sat_parameters_; + TimeLimit* time_limit_; + IntegerTrail* integer_trail_; + Trail* trail_; + IntegerEncoder* integer_encoder_; + SharedIncompleteSolutionManager* incomplete_solutions_; + const CpModelMapping* mapping_; + + // Last OPTIMAL/Feasible solution found by a call to the underlying LP solver. + bool lp_solution_is_set_ = false; + bool lp_solution_is_integer_ = false; + double lp_objective_; + std::vector lp_solution_; + std::vector best_lp_solution_; + // We use max fractionality of all variables. + double lp_solution_fractionality_; + + // Rounded Integer solution. This might not be feasible. + bool integer_solution_is_set_ = false; + bool integer_solution_is_feasible_ = false; + int64 integer_solution_objective_; + std::vector integer_solution_; + std::vector best_integer_solution_; + int num_infeasible_constraints_; + // We use max infeasibility of all constraints. + int64 integer_solution_infeasibility_; + + // Sum of all simplex iterations performed by this class. This is useful to + // test the incrementality and compare to other solvers. + int64 total_num_simplex_iterations_ = 0; + + // TODO(user): Tune default value. Expose as parameter. + int max_fp_iterations_ = 20; +}; + +} // namespace sat +} // namespace operations_research + +#endif // OR_TOOLS_SAT_FEASIBILITY_PUMP_H_ diff --git a/ortools/sat/rins.cc b/ortools/sat/rins.cc index 082aa5eab9..c601948a53 100644 --- a/ortools/sat/rins.cc +++ b/ortools/sat/rins.cc @@ -85,33 +85,50 @@ std::vector GetGeneralRelaxationValues( } return relaxation_values; } + +std::vector GetIncompleteSolutionValues( + SharedIncompleteSolutionManager* incomplete_solutions) { + std::vector empty_solution_values; + + if (incomplete_solutions == nullptr || + !incomplete_solutions->HasNewSolution()) { + return empty_solution_values; + } + + return incomplete_solutions->GetNewSolution(); +} } // namespace RINSNeighborhood GetRINSNeighborhood( const SharedResponseManager* response_manager, const SharedRelaxationSolutionRepository* relaxation_solutions, - const SharedLPSolutionRepository* lp_solutions, random_engine_t* random) { + const SharedLPSolutionRepository* lp_solutions, + SharedIncompleteSolutionManager* incomplete_solutions, + random_engine_t* random) { RINSNeighborhood rins_neighborhood; const bool use_only_relaxation_values = (response_manager == nullptr || response_manager->SolutionsRepository().NumSolutions() == 0); - if (use_only_relaxation_values && lp_solutions == nullptr) { + if (use_only_relaxation_values && lp_solutions == nullptr && + incomplete_solutions == nullptr) { // As of now RENS doesn't generate good neighborhoods from integer // relaxation solutions. return rins_neighborhood; } std::vector relaxation_values; - if (lp_solutions == nullptr) { + if (incomplete_solutions != nullptr) { + relaxation_values = GetIncompleteSolutionValues(incomplete_solutions); + } else if (lp_solutions != nullptr) { + relaxation_values = GetLPRelaxationValues(lp_solutions, random); + } else { CHECK(relaxation_solutions != nullptr) << "No relaxation solutions repository or lp solutions repository " "provided."; relaxation_values = GetGeneralRelaxationValues(relaxation_solutions, random); - } else { - relaxation_values = GetLPRelaxationValues(lp_solutions, random); } if (relaxation_values.empty()) return rins_neighborhood; diff --git a/ortools/sat/rins.h b/ortools/sat/rins.h index df0e4848ac..4131714ea2 100644 --- a/ortools/sat/rins.h +++ b/ortools/sat/rins.h @@ -60,17 +60,24 @@ struct RINSNeighborhood { }; // Helper method to create a RINS neighborhood by fixing variables with same -// values in relaxation solution and last solution. Uses lp values if -// 'use_lp_relaxation' is true, otherwise uses general relaxation value stored -// in the repository. +// values in relaxation solution and the current best solution in the +// response_manager. Prioritizes repositories in following order to get a +// relaxation solution. +// 1. incomplete_solutions +// 2. lp_solutions +// 3. relaxation_solutions // -// If use_only_relaxation_values is true, this Generates a RENS neighborhood by +// If response_manager is not provided, this generates a RENS neighborhood by // ignoring the solutions and using the relaxation values. The domain of the -// variables are reduced to integer values around relaxation values. +// variables are reduced to integer values around relaxation values. If the +// relaxation value is integer, then we fix the domain of the variable to that +// value. RINSNeighborhood GetRINSNeighborhood( const SharedResponseManager* response_manager, const SharedRelaxationSolutionRepository* relaxation_solutions, - const SharedLPSolutionRepository* lp_solutions, random_engine_t* random); + const SharedLPSolutionRepository* lp_solutions, + SharedIncompleteSolutionManager* incomplete_solutions, + random_engine_t* random); // Adds the current LP solution to the pool. void RecordLPRelaxationValues(Model* model); diff --git a/ortools/sat/samples/AssignmentSat.java b/ortools/sat/samples/AssignmentSat.java new file mode 100644 index 0000000000..0e52f3330b --- /dev/null +++ b/ortools/sat/samples/AssignmentSat.java @@ -0,0 +1,118 @@ +// Copyright 2010-2018 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// CP-SAT example that solves an assignment problem. +// [START program] +package com.google.ortools.sat.samples; +// [START import] +import com.google.ortools.sat.CpModel; +import com.google.ortools.sat.CpSolver; +import com.google.ortools.sat.CpSolverStatus; +import com.google.ortools.sat.IntVar; +import com.google.ortools.sat.LinearExpr; +// [END import] + +/** Assignment problem. */ +public class AssignmentSat { + static { + System.loadLibrary("jniortools"); + } + + public static void main(String[] args) { + // Data + // [START data_model] + int[][] costs = { + {90, 80, 75, 70}, + {35, 85, 55, 65}, + {125, 95, 90, 95}, + {45, 110, 95, 115}, + {50, 100, 90, 100}, + }; + final int numWorkers = costs.length; + final int numTasks = costs[0].length; + // [END data_model] + + // Model + // [START model] + CpModel model = new CpModel(); + // [END model] + + // Variables + // [START variables] + IntVar[][] x = new IntVar[numWorkers][numTasks]; + // Variables in a 1-dim array. + IntVar[] xFlat = new IntVar[numWorkers * numTasks]; + int[] costsFlat = new int[numWorkers * numTasks]; + for (int i = 0; i < numWorkers; ++i) { + for (int j = 0; j < numTasks; ++j) { + x[i][j] = model.newIntVar(0, 1, ""); + int k = i * numTasks + j; + xFlat[k] = x[i][j]; + costsFlat[k] = costs[i][j]; + } + } + // [END variables] + + // Constraints + // [START constraints] + // Each worker is assigned to at most one task. + for (int i = 0; i < numWorkers; ++i) { + IntVar[] vars = new IntVar[numTasks]; + for (int j = 0; j < numTasks; ++j) { + vars[j] = x[i][j]; + } + model.addLessOrEqual(LinearExpr.sum(vars), 1); + } + // Each task is assigned to exactly one worker. + for (int j = 0; j < numTasks; ++j) { + // LinearExpr taskSum; + IntVar[] vars = new IntVar[numWorkers]; + for (int i = 0; i < numWorkers; ++i) { + vars[i] = x[i][j]; + } + model.addEquality(LinearExpr.sum(vars), 1); + } + // [END constraints] + + // Objective + // [START objective] + model.minimize(LinearExpr.scalProd(xFlat, costsFlat)); + // [END objective] + + // Solve + // [START solve] + CpSolver solver = new CpSolver(); + CpSolverStatus status = solver.solve(model); + // [END solve] + + // Print solution. + // [START print_solution] + // Check that the problem has a feasible solution. + if (status == CpSolverStatus.OPTIMAL || status == CpSolverStatus.FEASIBLE) { + System.out.println("Total cost: " + solver.objectiveValue() + "\n"); + for (int i = 0; i < numWorkers; ++i) { + for (int j = 0; j < numTasks; ++j) { + if (solver.value(x[i][j]) == 1) { + System.out.println( + "Worker " + i + " assigned to task " + j + ". Cost: " + costs[i][j]); + } + } + } + } else { + System.err.println("No solution found."); + } + // [END print_solution] + } + + private AssignmentSat() {} +} diff --git a/ortools/sat/samples/assignment_sat.cc b/ortools/sat/samples/assignment_sat.cc new file mode 100644 index 0000000000..e5c640b8ec --- /dev/null +++ b/ortools/sat/samples/assignment_sat.cc @@ -0,0 +1,110 @@ +// Copyright 2010-2018 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// [START program] +// [START import] +#include "ortools/sat/cp_model.h" +// [END import] +namespace operations_research { +namespace sat { + +void IntegerProgrammingExample() { + // Data + // [START data_model] + const std::vector> costs{ + {90, 80, 75, 70}, {35, 85, 55, 65}, {125, 95, 90, 95}, + {45, 110, 95, 115}, {50, 100, 90, 100}, + }; + const int num_workers = costs.size(); + const int num_tasks = costs[0].size(); + // [END data_model] + + // Model + // [START model] + CpModelBuilder cp_model; + // [END model] + + // Variables + // [START variables] + // x[i][j] is an array of Boolean variables. x[i][j] is true + // if worker i is assigned to task j. + std::vector> x(num_workers, + std::vector(num_tasks)); + for (int i = 0; i < num_workers; ++i) { + for (int j = 0; j < num_tasks; ++j) { + x[i][j] = cp_model.NewBoolVar(); + } + } + // [END variables] + + // Constraints + // [START constraints] + // Each worker is assigned to at most one task. + for (int i = 0; i < num_workers; ++i) { + LinearExpr worker_sum; + for (int j = 0; j < num_tasks; ++j) { + worker_sum.AddTerm(x[i][j], 1); + } + cp_model.AddLessOrEqual(worker_sum, 1); + } + // Each task is assigned to exactly one worker. + for (int j = 0; j < num_tasks; ++j) { + LinearExpr task_sum; + for (int i = 0; i < num_workers; ++i) { + task_sum.AddTerm(x[i][j], 1); + } + cp_model.AddEquality(task_sum, 1); + } + // [END constraints] + + // Objective + // [START objective] + LinearExpr total_cost; + for (int i = 0; i < num_workers; ++i) { + for (int j = 0; j < num_tasks; ++j) { + total_cost.AddTerm(x[i][j], costs[i][j]); + } + } + cp_model.Minimize(total_cost); + // [END objective] + + // Solve + // [START solve] + const CpSolverResponse response = Solve(cp_model.Build()); + // [END solve] + + // Print solution. + // [START print_solution] + if (response.status() == CpSolverStatus::INFEASIBLE) { + LOG(FATAL) << "No solution found."; + } + + LOG(INFO) << "Total cost: " << response.objective_value(); + LOG(INFO); + for (int i = 0; i < num_workers; ++i) { + for (int j = 0; j < num_tasks; ++j) { + if (SolutionBooleanValue(response, x[i][j])) { + LOG(INFO) << "Task " << i << " assigned to worker " << j + << ". Cost: " << costs[i][j]; + } + } + } + // [END print_solution] +} +} // namespace sat +} // namespace operations_research + +int main(int argc, char** argv) { + operations_research::sat::IntegerProgrammingExample(); + return EXIT_SUCCESS; +} diff --git a/ortools/sat/samples/assignment_sat.py b/ortools/sat/samples/assignment_sat.py new file mode 100644 index 0000000000..5b3290d3f2 --- /dev/null +++ b/ortools/sat/samples/assignment_sat.py @@ -0,0 +1,92 @@ +# Copyright 2010-2018 Google LLC +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +"""Solve a simple assignment problem.""" +# [START program] +# [START import] +from ortools.sat.python import cp_model +# [END import] + + +def main(): + # Data + # [START data_model] + costs = [ + [90, 80, 75, 70], + [35, 85, 55, 65], + [125, 95, 90, 95], + [45, 110, 95, 115], + [50, 100, 90, 100], + ] + num_workers = len(costs) + num_tasks = len(costs[0]) + # [END data_model] + + # Model + # [START model] + model = cp_model.CpModel() + # [END model] + + # Variables + # [START variables] + x = [] + for i in range(num_workers): + t = [] + for j in range(num_tasks): + t.append(model.NewBoolVar('x[%i,%i]' % (i, j))) + x.append(t) + # [END variables] + + # Constraints + # [START constraints] + # Each worker is assigned to at most one task. + for i in range(num_workers): + model.Add(sum(x[i][j] for j in range(num_tasks)) <= 1) + + # Each task is assigned to exactly one worker. + for j in range(num_tasks): + model.Add(sum(x[i][j] for i in range(num_workers)) == 1) + # [END constraints] + + # Objective + # [START objective] + objective_terms = [] + for i in range(num_workers): + for j in range(num_tasks): + objective_terms.append(costs[i][j] * x[i][j]) + model.Minimize(sum(objective_terms)) + # [END objective] + + # Solve + # [START solve] + solver = cp_model.CpSolver() + status = solver.Solve(model) + # [END solve] + + # Print solution. + # [START print_solution] + if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE: + print('Total cost = %i' % solver.ObjectiveValue()) + print() + for i in range(num_workers): + for j in range(num_tasks): + if solver.BooleanValue(x[i][j]): + print('Worker ', i, ' assigned to task ', j, ' Cost = ', + costs[i][j]) + else: + print('No solution found.') + # [END print_solution] + + +if __name__ == '__main__': + main() +# [END program] diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index ee5cdf1e44..e5c6e585ca 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -21,7 +21,7 @@ option java_multiple_files = true; // Contains the definitions for all the sat algorithm parameters and their // default values. // -// NEXT TAG: 164 +// NEXT TAG: 165 message SatParameters { // ========================================================================== // Branching and polarity @@ -808,6 +808,9 @@ message SatParameters { // Turns on relaxation induced neighborhood generator. optional bool use_rins_lns = 129 [default = true]; + // Adds a feasibility pump subsolver along with lns subsolvers. + optional bool use_feasibility_pump = 164 [default = false]; + // Turns on a lns worker which solves relaxed version of the original problem // by removing constraints from the problem in order to get better bounds. optional bool use_relaxation_lns = 150 [default = false]; diff --git a/ortools/sat/synchronization.cc b/ortools/sat/synchronization.cc index 6f06161b47..5912168127 100644 --- a/ortools/sat/synchronization.cc +++ b/ortools/sat/synchronization.cc @@ -63,7 +63,7 @@ void SharedRelaxationSolutionRepository::NewRelaxationSolution( } void SharedLPSolutionRepository::NewLPSolution( - std::vector lp_solution) { + const std::vector& lp_solution) { // Note that the Add() method already applies mutex lock. So we don't need it // here. @@ -80,6 +80,27 @@ void SharedLPSolutionRepository::NewLPSolution( AddInternal(solution); } +bool SharedIncompleteSolutionManager::HasNewSolution() const { + absl::MutexLock mutex_lock(&mutex_); + return !solutions_.empty(); +} + +std::vector SharedIncompleteSolutionManager::GetNewSolution() { + absl::MutexLock mutex_lock(&mutex_); + std::vector solution; + if (solutions_.empty()) return solution; + + solution = std::move(solutions_.back()); + solutions_.pop_back(); + return solution; +} + +void SharedIncompleteSolutionManager::AddNewSolution( + const std::vector& lp_solution) { + absl::MutexLock mutex_lock(&mutex_); + solutions_.push_back(lp_solution); +} + // TODO(user): Experiments and play with the num_solutions_to_keep parameter. SharedResponseManager::SharedResponseManager(bool log_updates, bool enumerate_all_solutions, diff --git a/ortools/sat/synchronization.h b/ortools/sat/synchronization.h index acf44cda09..a5643630f8 100644 --- a/ortools/sat/synchronization.h +++ b/ortools/sat/synchronization.h @@ -136,7 +136,28 @@ class SharedLPSolutionRepository : public SharedSolutionRepository { explicit SharedLPSolutionRepository(int num_solutions_to_keep) : SharedSolutionRepository(num_solutions_to_keep) {} - void NewLPSolution(std::vector lp_solution); + void NewLPSolution(const std::vector& lp_solution); +}; + +// Set of partly filled solutions. They are meant to be finished by some lns +// worker. +// +// The solutions are stored as a vector of doubles. The value at index i +// represents the solution value of model variable indexed i. Note that some +// values can be infinity which should be interpreted as 'unknown' solution +// value for that variable. These solutions can not necessarily be completed to +// complete feasible solutions. +class SharedIncompleteSolutionManager { + public: + bool HasNewSolution() const; + std::vector GetNewSolution(); + + void AddNewSolution(const std::vector& lp_solution); + + private: + // New solutions are added and removed from the back. + std::vector> solutions_; + mutable absl::Mutex mutex_; }; // Manages the global best response kept by the solver.