OR-Tools  8.2
optimization.cc
Go to the documentation of this file.
1// Copyright 2010-2018 Google LLC
2// Licensed under the Apache License, Version 2.0 (the "License");
3// you may not use this file except in compliance with the License.
4// You may obtain a copy of the License at
5//
6// http://www.apache.org/licenses/LICENSE-2.0
7//
8// Unless required by applicable law or agreed to in writing, software
9// distributed under the License is distributed on an "AS IS" BASIS,
10// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11// See the License for the specific language governing permissions and
12// limitations under the License.
13
15
16#include <algorithm>
17#include <cstdio>
18#include <cstdlib>
19#include <deque>
20#include <limits>
21#include <map>
22#include <memory>
23#include <set>
24#include <string>
25#include <utility>
26
27#include "absl/random/random.h"
28#include "absl/strings/str_cat.h"
29#include "absl/strings/str_format.h"
33#include "ortools/base/macros.h"
36#include "ortools/base/timer.h"
37#if !defined(__PORTABLE_PLATFORM__) && defined(USE_SCIP)
39#include "ortools/linear_solver/linear_solver.pb.h"
40#endif // __PORTABLE_PLATFORM__
41#include "ortools/base/random.h"
47#include "ortools/sat/sat_parameters.pb.h"
48#include "ortools/sat/util.h"
50
51namespace operations_research {
52namespace sat {
53
54namespace {
55
56// Used to log messages to stdout or to the normal logging framework according
57// to the given LogBehavior value.
58class Logger {
59 public:
60 explicit Logger(LogBehavior v) : use_stdout_(v == STDOUT_LOG) {}
61 void Log(const std::string& message) {
62 if (use_stdout_) {
63 absl::PrintF("%s\n", message);
64 } else {
65 LOG(INFO) << message;
66 }
67 }
68
69 private:
70 bool use_stdout_;
71};
72
73// Outputs the current objective value in the cnf output format.
74// Note that this function scale the given objective.
75std::string CnfObjectiveLine(const LinearBooleanProblem& problem,
76 Coefficient objective) {
77 const double scaled_objective =
78 AddOffsetAndScaleObjectiveValue(problem, objective);
79 return absl::StrFormat("o %d", static_cast<int64>(scaled_objective));
80}
81
82struct LiteralWithCoreIndex {
83 LiteralWithCoreIndex(Literal l, int i) : literal(l), core_index(i) {}
84 Literal literal;
86};
87
88// Deletes the given indices from a vector. The given indices must be sorted in
89// increasing order. The order of the non-deleted entries in the vector is
90// preserved.
91template <typename Vector>
92void DeleteVectorIndices(const std::vector<int>& indices, Vector* v) {
93 int new_size = 0;
94 int indices_index = 0;
95 for (int i = 0; i < v->size(); ++i) {
96 if (indices_index < indices.size() && i == indices[indices_index]) {
97 ++indices_index;
98 } else {
99 (*v)[new_size] = (*v)[i];
100 ++new_size;
101 }
102 }
103 v->resize(new_size);
104}
105
106// In the Fu & Malik algorithm (or in WPM1), when two cores overlap, we
107// artifically introduce symmetries. More precisely:
108//
109// The picture below shows two cores with index 0 and 1, with one blocking
110// variable per '-' and with the variables ordered from left to right (by their
111// assumptions index). The blocking variables will be the one added to "relax"
112// the core for the next iteration.
113//
114// 1: -------------------------------
115// 0: ------------------------------------
116//
117// The 2 following assignment of the blocking variables are equivalent.
118// Remember that exactly one blocking variable per core must be assigned to 1.
119//
120// 1: ----------------------1--------
121// 0: --------1---------------------------
122//
123// and
124//
125// 1: ---------------------------1---
126// 0: ---1--------------------------------
127//
128// This class allows to add binary constraints excluding the second possibility.
129// Basically, each time a new core is added, if two of its blocking variables
130// (b1, b2) have the same assumption index of two blocking variables from
131// another core (c1, c2), then we forbid the assignment c1 true and b2 true.
132//
133// Reference: C Ansótegui, ML Bonet, J Levy, "Sat-based maxsat algorithms",
134// Artificial Intelligence, 2013 - Elsevier.
135class FuMalikSymmetryBreaker {
136 public:
137 FuMalikSymmetryBreaker() {}
138
139 // Must be called before a new core is processed.
140 void StartResolvingNewCore(int new_core_index) {
141 literal_by_core_.resize(new_core_index);
142 for (int i = 0; i < new_core_index; ++i) {
143 literal_by_core_[i].clear();
144 }
145 }
146
147 // This should be called for each blocking literal b of the new core. The
148 // assumption_index identify the soft clause associated to the given blocking
149 // literal. Note that between two StartResolvingNewCore() calls,
150 // ProcessLiteral() is assumed to be called with different assumption_index.
151 //
152 // Changing the order of the calls will not change the correctness, but will
153 // change the symmetry-breaking clause produced.
154 //
155 // Returns a set of literals which can't be true at the same time as b (under
156 // symmetry breaking).
157 std::vector<Literal> ProcessLiteral(int assumption_index, Literal b) {
158 if (assumption_index >= info_by_assumption_index_.size()) {
159 info_by_assumption_index_.resize(assumption_index + 1);
160 }
161
162 // Compute the function result.
163 // info_by_assumption_index_[assumption_index] will contain all the pairs
164 // (blocking_literal, core) of the previous resolved cores at the same
165 // assumption index as b.
166 std::vector<Literal> result;
167 for (LiteralWithCoreIndex data :
168 info_by_assumption_index_[assumption_index]) {
169 // literal_by_core_ will contain all the blocking literal of a given core
170 // with an assumption_index that was used in one of the ProcessLiteral()
171 // calls since the last StartResolvingNewCore().
172 //
173 // Note that there can be only one such literal by core, so we will not
174 // add duplicates.
175 result.insert(result.end(), literal_by_core_[data.core_index].begin(),
176 literal_by_core_[data.core_index].end());
177 }
178
179 // Update the internal data structure.
180 for (LiteralWithCoreIndex data :
181 info_by_assumption_index_[assumption_index]) {
182 literal_by_core_[data.core_index].push_back(data.literal);
183 }
184 info_by_assumption_index_[assumption_index].push_back(
185 LiteralWithCoreIndex(b, literal_by_core_.size()));
186 return result;
187 }
188
189 // Deletes the given assumption indices.
190 void DeleteIndices(const std::vector<int>& indices) {
191 DeleteVectorIndices(indices, &info_by_assumption_index_);
192 }
193
194 // This is only used in WPM1 to forget all the information related to a given
195 // assumption_index.
196 void ClearInfo(int assumption_index) {
197 CHECK_LE(assumption_index, info_by_assumption_index_.size());
198 info_by_assumption_index_[assumption_index].clear();
199 }
200
201 // This is only used in WPM1 when a new assumption_index is created.
202 void AddInfo(int assumption_index, Literal b) {
203 CHECK_GE(assumption_index, info_by_assumption_index_.size());
204 info_by_assumption_index_.resize(assumption_index + 1);
205 info_by_assumption_index_[assumption_index].push_back(
206 LiteralWithCoreIndex(b, literal_by_core_.size()));
207 }
208
209 private:
210 std::vector<std::vector<LiteralWithCoreIndex>> info_by_assumption_index_;
211 std::vector<std::vector<Literal>> literal_by_core_;
212
213 DISALLOW_COPY_AND_ASSIGN(FuMalikSymmetryBreaker);
214};
215
216} // namespace
217
219 std::vector<Literal>* core) {
220 if (solver->IsModelUnsat()) return;
221 std::set<LiteralIndex> moved_last;
222 std::vector<Literal> candidate(core->begin(), core->end());
223
224 solver->Backtrack(0);
225 solver->SetAssumptionLevel(0);
226 while (!limit->LimitReached()) {
227 // We want each literal in candidate to appear last once in our propagation
228 // order. We want to do that while maximizing the reutilization of the
229 // current assignment prefix, that is minimizing the number of
230 // decision/progagation we need to perform.
231 const int target_level = MoveOneUnprocessedLiteralLast(
232 moved_last, solver->CurrentDecisionLevel(), &candidate);
233 if (target_level == -1) break;
234 solver->Backtrack(target_level);
235 while (!solver->IsModelUnsat() && !limit->LimitReached() &&
236 solver->CurrentDecisionLevel() < candidate.size()) {
237 const Literal decision = candidate[solver->CurrentDecisionLevel()];
238 if (solver->Assignment().LiteralIsTrue(decision)) {
239 candidate.erase(candidate.begin() + solver->CurrentDecisionLevel());
240 continue;
241 } else if (solver->Assignment().LiteralIsFalse(decision)) {
242 // This is a "weird" API to get the subset of decisions that caused
243 // this literal to be false with reason analysis.
245 candidate = solver->GetLastIncompatibleDecisions();
246 break;
247 } else {
248 solver->EnqueueDecisionAndBackjumpOnConflict(decision);
249 }
250 }
251 if (candidate.empty() || solver->IsModelUnsat()) return;
252 moved_last.insert(candidate.back().Index());
253 }
254
255 solver->Backtrack(0);
256 solver->SetAssumptionLevel(0);
257 if (candidate.size() < core->size()) {
258 VLOG(1) << "minimization " << core->size() << " -> " << candidate.size();
259 core->assign(candidate.begin(), candidate.end());
260 }
261}
262
263// This algorithm works by exploiting the unsat core returned by the SAT solver
264// when the problem is UNSAT. It starts by trying to solve the decision problem
265// where all the objective variables are set to their value with minimal cost,
266// and relax in each step some of these fixed variables until the problem
267// becomes satisfiable.
269 const LinearBooleanProblem& problem,
270 SatSolver* solver,
271 std::vector<bool>* solution) {
272 Logger logger(log);
273 FuMalikSymmetryBreaker symmetry;
274
275 // blocking_clauses will contains a set of clauses that are currently added to
276 // the initial problem.
277 //
278 // Initially, each clause just contains a literal associated to an objective
279 // variable with non-zero cost. Setting all these literals to true will lead
280 // to the lowest possible objective.
281 //
282 // During the algorithm, "blocking" literals will be added to each clause.
283 // Moreover each clause will contain an extra "assumption" literal stored in
284 // the separate assumptions vector (in its negated form).
285 //
286 // The meaning of a given clause will always be:
287 // If the assumption literal and all blocking literals are false, then the
288 // "objective" literal (which is the first one in the clause) must be true.
289 // When the "objective" literal is true, its variable (which have a non-zero
290 // cost) is set to the value that minimize the objective cost.
291 //
292 // ex: If a variable "x" as a cost of 3, its cost contribution is smaller when
293 // it is set to false (since it will contribute to zero instead of 3).
294 std::vector<std::vector<Literal>> blocking_clauses;
295 std::vector<Literal> assumptions;
296
297 // Initialize blocking_clauses and assumptions.
298 const LinearObjective& objective = problem.objective();
299 CHECK_GT(objective.coefficients_size(), 0);
300 const Coefficient unique_objective_coeff(std::abs(objective.coefficients(0)));
301 for (int i = 0; i < objective.literals_size(); ++i) {
302 CHECK_EQ(std::abs(objective.coefficients(i)), unique_objective_coeff)
303 << "The basic Fu & Malik algorithm needs constant objective coeffs.";
304 const Literal literal(objective.literals(i));
305
306 // We want to minimize the cost when this literal is true.
307 const Literal min_literal =
308 objective.coefficients(i) > 0 ? literal.Negated() : literal;
309 blocking_clauses.push_back(std::vector<Literal>(1, min_literal));
310
311 // Note that initialy, we do not create any extra variables.
312 assumptions.push_back(min_literal);
313 }
314
315 // Print the number of variable with a non-zero cost.
316 logger.Log(absl::StrFormat("c #weights:%u #vars:%d #constraints:%d",
317 assumptions.size(), problem.num_variables(),
318 problem.constraints_size()));
319
320 // Starts the algorithm. Each loop will solve the problem under the given
321 // assumptions, and if unsat, will relax exactly one of the objective
322 // variables (from the unsat core) to be in its "costly" state. When the
323 // algorithm terminates, the number of iterations is exactly the minimal
324 // objective value.
325 for (int iter = 0;; ++iter) {
326 const SatSolver::Status result =
327 solver->ResetAndSolveWithGivenAssumptions(assumptions);
328 if (result == SatSolver::FEASIBLE) {
329 ExtractAssignment(problem, *solver, solution);
330 Coefficient objective = ComputeObjectiveValue(problem, *solution);
331 logger.Log(CnfObjectiveLine(problem, objective));
332 return SatSolver::FEASIBLE;
333 }
334 if (result != SatSolver::ASSUMPTIONS_UNSAT) return result;
335
336 // The interesting case: we have an unsat core.
337 //
338 // We need to add new "blocking" variables b_i for all the objective
339 // variable appearing in the core. Moreover, we will only relax as little
340 // as possible (to not miss the optimal), so we will enforce that the sum
341 // of the b_i is exactly one.
342 std::vector<Literal> core = solver->GetLastIncompatibleDecisions();
343 MinimizeCore(solver, &core);
344 solver->Backtrack(0);
345
346 // Print the search progress.
347 logger.Log(absl::StrFormat("c iter:%d core:%u", iter, core.size()));
348
349 // Special case for a singleton core.
350 if (core.size() == 1) {
351 // Find the index of the "objective" variable that need to be fixed in
352 // its "costly" state.
353 const int index =
354 std::find(assumptions.begin(), assumptions.end(), core[0]) -
355 assumptions.begin();
356 CHECK_LT(index, assumptions.size());
357
358 // Fix it. We also fix all the associated blocking variables if any.
359 if (!solver->AddUnitClause(core[0].Negated())) {
361 }
362 for (Literal b : blocking_clauses[index]) {
363 if (!solver->AddUnitClause(b.Negated())) return SatSolver::INFEASIBLE;
364 }
365
366 // Erase this entry from the current "objective"
367 std::vector<int> to_delete(1, index);
368 DeleteVectorIndices(to_delete, &assumptions);
369 DeleteVectorIndices(to_delete, &blocking_clauses);
370 symmetry.DeleteIndices(to_delete);
371 } else {
372 symmetry.StartResolvingNewCore(iter);
373
374 // We will add 2 * |core.size()| variables.
375 const int old_num_variables = solver->NumVariables();
376 if (core.size() == 2) {
377 // Special case. If core.size() == 2, we can use only one blocking
378 // variable (the other one beeing its negation). This actually do happen
379 // quite often in practice, so it is worth it.
380 solver->SetNumVariables(old_num_variables + 3);
381 } else {
382 solver->SetNumVariables(old_num_variables + 2 * core.size());
383 }
384
385 // Temporary vectors for the constraint (sum new blocking variable == 1).
386 std::vector<LiteralWithCoeff> at_most_one_constraint;
387 std::vector<Literal> at_least_one_constraint;
388
389 // This will be set to false if the problem becomes unsat while adding a
390 // new clause. This is unlikely, but may be possible.
391 bool ok = true;
392
393 // Loop over the core.
394 int index = 0;
395 for (int i = 0; i < core.size(); ++i) {
396 // Since the assumptions appear in order in the core, we can find the
397 // relevant "objective" variable efficiently with a simple linear scan
398 // in the assumptions vector (done with index).
399 index =
400 std::find(assumptions.begin() + index, assumptions.end(), core[i]) -
401 assumptions.begin();
402 CHECK_LT(index, assumptions.size());
403
404 // The new blocking and assumption variables for this core entry.
405 const Literal a(BooleanVariable(old_num_variables + i), true);
406 Literal b(BooleanVariable(old_num_variables + core.size() + i), true);
407 if (core.size() == 2) {
408 b = Literal(BooleanVariable(old_num_variables + 2), true);
409 if (i == 1) b = b.Negated();
410 }
411
412 // Symmetry breaking clauses.
413 for (Literal l : symmetry.ProcessLiteral(index, b)) {
414 ok &= solver->AddBinaryClause(l.Negated(), b.Negated());
415 }
416
417 // Note(user): There is more than one way to encode the algorithm in
418 // SAT. Here we "delete" the old blocking clause and add a new one. In
419 // the WPM1 algorithm below, the blocking clause is decomposed into
420 // 3-SAT and we don't need to delete anything.
421
422 // First, fix the old "assumption" variable to false, which has the
423 // effect of deleting the old clause from the solver.
424 if (assumptions[index].Variable() >= problem.num_variables()) {
425 CHECK(solver->AddUnitClause(assumptions[index].Negated()));
426 }
427
428 // Add the new blocking variable.
429 blocking_clauses[index].push_back(b);
430
431 // Add the new clause to the solver. Temporary including the
432 // assumption, but removing it right afterwards.
433 blocking_clauses[index].push_back(a);
434 ok &= solver->AddProblemClause(blocking_clauses[index]);
435 blocking_clauses[index].pop_back();
436
437 // For the "== 1" constraint on the blocking literals.
438 at_most_one_constraint.push_back(LiteralWithCoeff(b, 1.0));
439 at_least_one_constraint.push_back(b);
440
441 // The new assumption variable replace the old one.
442 assumptions[index] = a.Negated();
443 }
444
445 // Add the "<= 1" side of the "== 1" constraint.
446 ok &= solver->AddLinearConstraint(false, Coefficient(0), true,
447 Coefficient(1.0),
448 &at_most_one_constraint);
449
450 // TODO(user): The algorithm does not really need the >= 1 side of this
451 // constraint. Initial investigation shows that it doesn't really help,
452 // but investigate more.
453 if (/* DISABLES CODE */ (false)) {
454 ok &= solver->AddProblemClause(at_least_one_constraint);
455 }
456
457 if (!ok) {
458 LOG(INFO) << "Infeasible while adding a clause.";
460 }
461 }
462 }
463}
464
466 const LinearBooleanProblem& problem,
467 SatSolver* solver,
468 std::vector<bool>* solution) {
469 Logger logger(log);
470 FuMalikSymmetryBreaker symmetry;
471
472 // The current lower_bound on the cost.
473 // It will be correct after the initialization.
474 Coefficient lower_bound(static_cast<int64>(problem.objective().offset()));
475 Coefficient upper_bound(kint64max);
476
477 // The assumption literals and their associated cost.
478 std::vector<Literal> assumptions;
479 std::vector<Coefficient> costs;
480 std::vector<Literal> reference;
481
482 // Initialization.
483 const LinearObjective& objective = problem.objective();
484 CHECK_GT(objective.coefficients_size(), 0);
485 for (int i = 0; i < objective.literals_size(); ++i) {
486 const Literal literal(objective.literals(i));
487 const Coefficient coeff(objective.coefficients(i));
488
489 // We want to minimize the cost when the assumption is true.
490 // Note that initially, we do not create any extra variables.
491 if (coeff > 0) {
492 assumptions.push_back(literal.Negated());
493 costs.push_back(coeff);
494 } else {
495 assumptions.push_back(literal);
496 costs.push_back(-coeff);
497 lower_bound += coeff;
498 }
499 }
500 reference = assumptions;
501
502 // This is used by the "stratified" approach.
503 Coefficient stratified_lower_bound =
504 *std::max_element(costs.begin(), costs.end());
505
506 // Print the number of variables with a non-zero cost.
507 logger.Log(absl::StrFormat("c #weights:%u #vars:%d #constraints:%d",
508 assumptions.size(), problem.num_variables(),
509 problem.constraints_size()));
510
511 for (int iter = 0;; ++iter) {
512 // This is called "hardening" in the literature.
513 // Basically, we know that there is only hardening_threshold weight left
514 // to distribute, so any assumption with a greater cost than this can never
515 // be false. We fix it instead of treating it as an assumption.
516 solver->Backtrack(0);
517 const Coefficient hardening_threshold = upper_bound - lower_bound;
518 CHECK_GE(hardening_threshold, 0);
519 std::vector<int> to_delete;
520 int num_above_threshold = 0;
521 for (int i = 0; i < assumptions.size(); ++i) {
522 if (costs[i] > hardening_threshold) {
523 if (!solver->AddUnitClause(assumptions[i])) {
525 }
526 to_delete.push_back(i);
527 ++num_above_threshold;
528 } else {
529 // This impact the stratification heuristic.
530 if (solver->Assignment().LiteralIsTrue(assumptions[i])) {
531 to_delete.push_back(i);
532 }
533 }
534 }
535 if (!to_delete.empty()) {
536 logger.Log(absl::StrFormat("c fixed %u assumptions, %d with cost > %d",
537 to_delete.size(), num_above_threshold,
538 hardening_threshold.value()));
539 DeleteVectorIndices(to_delete, &assumptions);
540 DeleteVectorIndices(to_delete, &costs);
541 DeleteVectorIndices(to_delete, &reference);
542 symmetry.DeleteIndices(to_delete);
543 }
544
545 // This is the "stratification" part.
546 // Extract the assumptions with a cost >= stratified_lower_bound.
547 std::vector<Literal> assumptions_subset;
548 for (int i = 0; i < assumptions.size(); ++i) {
549 if (costs[i] >= stratified_lower_bound) {
550 assumptions_subset.push_back(assumptions[i]);
551 }
552 }
553
554 const SatSolver::Status result =
555 solver->ResetAndSolveWithGivenAssumptions(assumptions_subset);
556 if (result == SatSolver::FEASIBLE) {
557 // If not all assumptions were taken, continue with a lower stratified
558 // bound. Otherwise we have an optimal solution!
559 //
560 // TODO(user): Try more advanced variant where the bound is lowered by
561 // more than this minimal amount.
562 const Coefficient old_lower_bound = stratified_lower_bound;
563 for (Coefficient cost : costs) {
564 if (cost < old_lower_bound) {
565 if (stratified_lower_bound == old_lower_bound ||
566 cost > stratified_lower_bound) {
567 stratified_lower_bound = cost;
568 }
569 }
570 }
571
572 ExtractAssignment(problem, *solver, solution);
573 DCHECK(IsAssignmentValid(problem, *solution));
574 const Coefficient objective_offset(
575 static_cast<int64>(problem.objective().offset()));
576 const Coefficient objective = ComputeObjectiveValue(problem, *solution);
577 if (objective + objective_offset < upper_bound) {
578 logger.Log(CnfObjectiveLine(problem, objective));
579 upper_bound = objective + objective_offset;
580 }
581
582 if (stratified_lower_bound < old_lower_bound) continue;
583 return SatSolver::FEASIBLE;
584 }
585 if (result != SatSolver::ASSUMPTIONS_UNSAT) return result;
586
587 // The interesting case: we have an unsat core.
588 //
589 // We need to add new "blocking" variables b_i for all the objective
590 // variables appearing in the core. Moreover, we will only relax as little
591 // as possible (to not miss the optimal), so we will enforce that the sum
592 // of the b_i is exactly one.
593 std::vector<Literal> core = solver->GetLastIncompatibleDecisions();
594 MinimizeCore(solver, &core);
595 solver->Backtrack(0);
596
597 // Compute the min cost of all the assertions in the core.
598 // The lower bound will be updated by that much.
599 Coefficient min_cost = kCoefficientMax;
600 {
601 int index = 0;
602 for (int i = 0; i < core.size(); ++i) {
603 index =
604 std::find(assumptions.begin() + index, assumptions.end(), core[i]) -
605 assumptions.begin();
606 CHECK_LT(index, assumptions.size());
607 min_cost = std::min(min_cost, costs[index]);
608 }
609 }
610 lower_bound += min_cost;
611
612 // Print the search progress.
613 logger.Log(absl::StrFormat(
614 "c iter:%d core:%u lb:%d min_cost:%d strat:%d", iter, core.size(),
615 lower_bound.value(), min_cost.value(), stratified_lower_bound.value()));
616
617 // This simple line helps a lot on the packup-wpms instances!
618 //
619 // TODO(user): That was because of a bug before in the way
620 // stratified_lower_bound was decremented, not sure it helps that much now.
621 if (min_cost > stratified_lower_bound) {
622 stratified_lower_bound = min_cost;
623 }
624
625 // Special case for a singleton core.
626 if (core.size() == 1) {
627 // Find the index of the "objective" variable that need to be fixed in
628 // its "costly" state.
629 const int index =
630 std::find(assumptions.begin(), assumptions.end(), core[0]) -
631 assumptions.begin();
632 CHECK_LT(index, assumptions.size());
633
634 // Fix it.
635 if (!solver->AddUnitClause(core[0].Negated())) {
637 }
638
639 // Erase this entry from the current "objective".
640 std::vector<int> to_delete(1, index);
641 DeleteVectorIndices(to_delete, &assumptions);
642 DeleteVectorIndices(to_delete, &costs);
643 DeleteVectorIndices(to_delete, &reference);
644 symmetry.DeleteIndices(to_delete);
645 } else {
646 symmetry.StartResolvingNewCore(iter);
647
648 // We will add 2 * |core.size()| variables.
649 const int old_num_variables = solver->NumVariables();
650 if (core.size() == 2) {
651 // Special case. If core.size() == 2, we can use only one blocking
652 // variable (the other one beeing its negation). This actually do happen
653 // quite often in practice, so it is worth it.
654 solver->SetNumVariables(old_num_variables + 3);
655 } else {
656 solver->SetNumVariables(old_num_variables + 2 * core.size());
657 }
658
659 // Temporary vectors for the constraint (sum new blocking variable == 1).
660 std::vector<LiteralWithCoeff> at_most_one_constraint;
661 std::vector<Literal> at_least_one_constraint;
662
663 // This will be set to false if the problem becomes unsat while adding a
664 // new clause. This is unlikely, but may be possible.
665 bool ok = true;
666
667 // Loop over the core.
668 int index = 0;
669 for (int i = 0; i < core.size(); ++i) {
670 // Since the assumptions appear in order in the core, we can find the
671 // relevant "objective" variable efficiently with a simple linear scan
672 // in the assumptions vector (done with index).
673 index =
674 std::find(assumptions.begin() + index, assumptions.end(), core[i]) -
675 assumptions.begin();
676 CHECK_LT(index, assumptions.size());
677
678 // The new blocking and assumption variables for this core entry.
679 const Literal a(BooleanVariable(old_num_variables + i), true);
680 Literal b(BooleanVariable(old_num_variables + core.size() + i), true);
681 if (core.size() == 2) {
682 b = Literal(BooleanVariable(old_num_variables + 2), true);
683 if (i == 1) b = b.Negated();
684 }
685
686 // a false & b false => previous assumptions (which was false).
687 const Literal old_a = assumptions[index];
688 ok &= solver->AddTernaryClause(a, b, old_a);
689
690 // Optional. Also add the two implications a => x and b => x where x is
691 // the negation of the previous assumption variable.
692 ok &= solver->AddBinaryClause(a.Negated(), old_a.Negated());
693 ok &= solver->AddBinaryClause(b.Negated(), old_a.Negated());
694
695 // Optional. Also add the implication a => not(b).
696 ok &= solver->AddBinaryClause(a.Negated(), b.Negated());
697
698 // This is the difference with the Fu & Malik algorithm.
699 // If the soft clause protected by old_a has a cost greater than
700 // min_cost then:
701 // - its cost is disminished by min_cost.
702 // - an identical clause with cost min_cost is artifically added to
703 // the problem.
704 CHECK_GE(costs[index], min_cost);
705 if (costs[index] == min_cost) {
706 // The new assumption variable replaces the old one.
707 assumptions[index] = a.Negated();
708
709 // Symmetry breaking clauses.
710 for (Literal l : symmetry.ProcessLiteral(index, b)) {
711 ok &= solver->AddBinaryClause(l.Negated(), b.Negated());
712 }
713 } else {
714 // Since the cost of the given index changes, we need to start a new
715 // "equivalence" class for the symmetry breaking algo and clear the
716 // old one.
717 symmetry.AddInfo(assumptions.size(), b);
718 symmetry.ClearInfo(index);
719
720 // Reduce the cost of the old assumption.
721 costs[index] -= min_cost;
722
723 // We add the new assumption with a cost of min_cost.
724 //
725 // Note(user): I think it is nice that these are added after old_a
726 // because assuming old_a will implies all the derived assumptions to
727 // true, and thus they will never appear in a core until old_a is not
728 // an assumption anymore.
729 assumptions.push_back(a.Negated());
730 costs.push_back(min_cost);
731 reference.push_back(reference[index]);
732 }
733
734 // For the "<= 1" constraint on the blocking literals.
735 // Note(user): we don't add the ">= 1" side because it is not needed for
736 // the correctness and it doesn't seems to help.
737 at_most_one_constraint.push_back(LiteralWithCoeff(b, 1.0));
738
739 // Because we have a core, we know that at least one of the initial
740 // problem variables must be true. This seems to help a bit.
741 //
742 // TODO(user): Experiment more.
743 at_least_one_constraint.push_back(reference[index].Negated());
744 }
745
746 // Add the "<= 1" side of the "== 1" constraint.
747 ok &= solver->AddLinearConstraint(false, Coefficient(0), true,
748 Coefficient(1.0),
749 &at_most_one_constraint);
750
751 // Optional. Add the ">= 1" constraint on the initial problem variables.
752 ok &= solver->AddProblemClause(at_least_one_constraint);
753
754 if (!ok) {
755 LOG(INFO) << "Unsat while adding a clause.";
757 }
758 }
759 }
760}
761
763 const LinearBooleanProblem& problem,
764 int num_times, SatSolver* solver,
765 std::vector<bool>* solution) {
766 Logger logger(log);
767 const SatParameters initial_parameters = solver->parameters();
768
769 MTRandom random("A random seed.");
770 SatParameters parameters = initial_parameters;
771 TimeLimit time_limit(parameters.max_time_in_seconds());
772
773 // We start with a low conflict limit and increase it until we are able to
774 // solve the problem at least once. After this, the limit stays the same.
775 int max_number_of_conflicts = 5;
776 parameters.set_log_search_progress(false);
777
778 Coefficient min_seen(std::numeric_limits<int64>::max());
779 Coefficient max_seen(std::numeric_limits<int64>::min());
780 Coefficient best(min_seen);
781 for (int i = 0; i < num_times; ++i) {
782 solver->Backtrack(0);
784
785 parameters.set_max_number_of_conflicts(max_number_of_conflicts);
786 parameters.set_max_time_in_seconds(time_limit.GetTimeLeft());
787 parameters.set_random_seed(i);
788 solver->SetParameters(parameters);
789 solver->ResetDecisionHeuristic();
790
791 const bool use_obj = absl::Bernoulli(random, 1.0 / 4);
792 if (use_obj) UseObjectiveForSatAssignmentPreference(problem, solver);
793
794 const SatSolver::Status result = solver->Solve();
795 if (result == SatSolver::INFEASIBLE) {
796 // If the problem is INFEASIBLE after we over-constrained the objective,
797 // then we found an optimal solution, otherwise, even the decision problem
798 // is INFEASIBLE.
799 if (best == kCoefficientMax) return SatSolver::INFEASIBLE;
800 return SatSolver::FEASIBLE;
801 }
802 if (result == SatSolver::LIMIT_REACHED) {
803 // We augment the number of conflict until we have one feasible solution.
804 if (best == kCoefficientMax) ++max_number_of_conflicts;
805 if (time_limit.LimitReached()) return SatSolver::LIMIT_REACHED;
806 continue;
807 }
808
810 std::vector<bool> candidate;
811 ExtractAssignment(problem, *solver, &candidate);
812 CHECK(IsAssignmentValid(problem, candidate));
813 const Coefficient objective = ComputeObjectiveValue(problem, candidate);
814 if (objective < best) {
815 *solution = candidate;
816 best = objective;
817 logger.Log(CnfObjectiveLine(problem, objective));
818
819 // Overconstrain the objective.
820 solver->Backtrack(0);
821 if (!AddObjectiveConstraint(problem, false, Coefficient(0), true,
822 objective - 1, solver)) {
823 return SatSolver::FEASIBLE;
824 }
825 }
826 min_seen = std::min(min_seen, objective);
827 max_seen = std::max(max_seen, objective);
828
829 logger.Log(absl::StrCat(
830 "c ", objective.value(), " [", min_seen.value(), ", ", max_seen.value(),
831 "] objective_preference: ", use_obj ? "true" : "false", " ",
833 }
834
835 // Retore the initial parameter (with an updated time limit).
836 parameters = initial_parameters;
837 parameters.set_max_time_in_seconds(time_limit.GetTimeLeft());
838 solver->SetParameters(parameters);
840}
841
843 const LinearBooleanProblem& problem,
844 SatSolver* solver,
845 std::vector<bool>* solution) {
846 Logger logger(log);
847
848 // This has a big positive impact on most problems.
850
851 Coefficient objective = kCoefficientMax;
852 if (!solution->empty()) {
853 CHECK(IsAssignmentValid(problem, *solution));
854 objective = ComputeObjectiveValue(problem, *solution);
855 }
856 while (true) {
857 if (objective != kCoefficientMax) {
858 // Over constrain the objective.
859 solver->Backtrack(0);
860 if (!AddObjectiveConstraint(problem, false, Coefficient(0), true,
861 objective - 1, solver)) {
862 return SatSolver::FEASIBLE;
863 }
864 }
865
866 // Solve the problem.
867 const SatSolver::Status result = solver->Solve();
869 if (result == SatSolver::INFEASIBLE) {
870 if (objective == kCoefficientMax) return SatSolver::INFEASIBLE;
871 return SatSolver::FEASIBLE;
872 }
873 if (result == SatSolver::LIMIT_REACHED) {
875 }
876
877 // Extract the new best solution.
879 ExtractAssignment(problem, *solver, solution);
880 CHECK(IsAssignmentValid(problem, *solution));
881 const Coefficient old_objective = objective;
882 objective = ComputeObjectiveValue(problem, *solution);
883 CHECK_LT(objective, old_objective);
884 logger.Log(CnfObjectiveLine(problem, objective));
885 }
886}
887
889 LogBehavior log, const LinearBooleanProblem& problem, SatSolver* solver,
890 std::vector<bool>* solution) {
891 Logger logger(log);
892 std::deque<EncodingNode> repository;
893
894 // Create one initial node per variables with cost.
895 Coefficient offset(0);
896 std::vector<EncodingNode*> nodes =
897 CreateInitialEncodingNodes(problem.objective(), &offset, &repository);
898
899 // This algorithm only work with weights of the same magnitude.
900 CHECK(!nodes.empty());
901 const Coefficient reference = nodes.front()->weight();
902 for (const EncodingNode* n : nodes) CHECK_EQ(n->weight(), reference);
903
904 // Initialize the current objective.
905 Coefficient objective = kCoefficientMax;
906 Coefficient upper_bound = kCoefficientMax;
907 if (!solution->empty()) {
908 CHECK(IsAssignmentValid(problem, *solution));
909 objective = ComputeObjectiveValue(problem, *solution);
910 upper_bound = objective + offset;
911 }
912
913 // Print the number of variables with a non-zero cost.
914 logger.Log(absl::StrFormat("c #weights:%u #vars:%d #constraints:%d",
915 nodes.size(), problem.num_variables(),
916 problem.constraints_size()));
917
918 // Create the sorter network.
919 solver->Backtrack(0);
920 EncodingNode* root =
921 MergeAllNodesWithDeque(upper_bound, nodes, solver, &repository);
922 logger.Log(absl::StrFormat("c encoding depth:%d", root->depth()));
923
924 while (true) {
925 if (objective != kCoefficientMax) {
926 // Over constrain the objective by fixing the variable index - 1 of the
927 // root node to 0.
928 const int index = offset.value() + objective.value();
929 if (index == 0) return SatSolver::FEASIBLE;
930 solver->Backtrack(0);
931 if (!solver->AddUnitClause(root->literal(index - 1).Negated())) {
932 return SatSolver::FEASIBLE;
933 }
934 }
935
936 // Solve the problem.
937 const SatSolver::Status result = solver->Solve();
939 if (result == SatSolver::INFEASIBLE) {
940 if (objective == kCoefficientMax) return SatSolver::INFEASIBLE;
941 return SatSolver::FEASIBLE;
942 }
944
945 // Extract the new best solution.
947 ExtractAssignment(problem, *solver, solution);
948 CHECK(IsAssignmentValid(problem, *solution));
949 const Coefficient old_objective = objective;
950 objective = ComputeObjectiveValue(problem, *solution);
951 CHECK_LT(objective, old_objective);
952 logger.Log(CnfObjectiveLine(problem, objective));
953 }
954}
955
957 LogBehavior log, const LinearBooleanProblem& problem, SatSolver* solver,
958 std::vector<bool>* solution) {
959 Logger logger(log);
960 SatParameters parameters = solver->parameters();
961
962 // Create one initial nodes per variables with cost.
963 Coefficient offset(0);
964 std::deque<EncodingNode> repository;
965 std::vector<EncodingNode*> nodes =
966 CreateInitialEncodingNodes(problem.objective(), &offset, &repository);
967
968 // Initialize the bounds.
969 // This is in term of number of variables not at their minimal value.
970 Coefficient lower_bound(0);
971 Coefficient upper_bound(kint64max);
972 if (!solution->empty()) {
973 CHECK(IsAssignmentValid(problem, *solution));
974 upper_bound = ComputeObjectiveValue(problem, *solution) + offset;
975 }
976
977 // Print the number of variables with a non-zero cost.
978 logger.Log(absl::StrFormat("c #weights:%u #vars:%d #constraints:%d",
979 nodes.size(), problem.num_variables(),
980 problem.constraints_size()));
981
982 // This is used by the "stratified" approach.
983 Coefficient stratified_lower_bound(0);
984 if (parameters.max_sat_stratification() ==
985 SatParameters::STRATIFICATION_DESCENT) {
986 // In this case, we initialize it to the maximum assumption weights.
987 for (EncodingNode* n : nodes) {
988 stratified_lower_bound = std::max(stratified_lower_bound, n->weight());
989 }
990 }
991
992 // Start the algorithm.
993 int max_depth = 0;
994 std::string previous_core_info = "";
995 for (int iter = 0;; ++iter) {
996 const std::vector<Literal> assumptions = ReduceNodesAndExtractAssumptions(
997 upper_bound, stratified_lower_bound, &lower_bound, &nodes, solver);
998 if (assumptions.empty()) return SatSolver::FEASIBLE;
999
1000 // Display the progress.
1001 const std::string gap_string =
1002 (upper_bound == kCoefficientMax)
1003 ? ""
1004 : absl::StrFormat(" gap:%d", (upper_bound - lower_bound).value());
1005 logger.Log(
1006 absl::StrFormat("c iter:%d [%s] lb:%d%s assumptions:%u depth:%d", iter,
1007 previous_core_info,
1008 lower_bound.value() - offset.value() +
1009 static_cast<int64>(problem.objective().offset()),
1010 gap_string, nodes.size(), max_depth));
1011
1012 // Solve under the assumptions.
1013 const SatSolver::Status result =
1014 solver->ResetAndSolveWithGivenAssumptions(assumptions);
1015 if (result == SatSolver::FEASIBLE) {
1016 // Extract the new solution and save it if it is the best found so far.
1017 std::vector<bool> temp_solution;
1018 ExtractAssignment(problem, *solver, &temp_solution);
1019 CHECK(IsAssignmentValid(problem, temp_solution));
1020 const Coefficient obj = ComputeObjectiveValue(problem, temp_solution);
1021 if (obj + offset < upper_bound) {
1022 *solution = temp_solution;
1023 logger.Log(CnfObjectiveLine(problem, obj));
1024 upper_bound = obj + offset;
1025 }
1026
1027 // If not all assumptions were taken, continue with a lower stratified
1028 // bound. Otherwise we have an optimal solution.
1029 stratified_lower_bound =
1030 MaxNodeWeightSmallerThan(nodes, stratified_lower_bound);
1031 if (stratified_lower_bound > 0) continue;
1032 return SatSolver::FEASIBLE;
1033 }
1034 if (result != SatSolver::ASSUMPTIONS_UNSAT) return result;
1035
1036 // We have a new core.
1037 std::vector<Literal> core = solver->GetLastIncompatibleDecisions();
1038 if (parameters.minimize_core()) MinimizeCore(solver, &core);
1039
1040 // Compute the min weight of all the nodes in the core.
1041 // The lower bound will be increased by that much.
1042 const Coefficient min_weight = ComputeCoreMinWeight(nodes, core);
1043 previous_core_info =
1044 absl::StrFormat("core:%u mw:%d", core.size(), min_weight.value());
1045
1046 // Increase stratified_lower_bound according to the parameters.
1047 if (stratified_lower_bound < min_weight &&
1048 parameters.max_sat_stratification() ==
1049 SatParameters::STRATIFICATION_ASCENT) {
1050 stratified_lower_bound = min_weight;
1051 }
1052
1053 ProcessCore(core, min_weight, &repository, &nodes, solver);
1054 max_depth = std::max(max_depth, nodes.back()->depth());
1055 }
1056}
1057
1059 IntegerVariable objective_var,
1060 const std::function<void()>& feasible_solution_observer, Model* model) {
1061 SatSolver* sat_solver = model->GetOrCreate<SatSolver>();
1062 IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
1063 const SatParameters& parameters = *(model->GetOrCreate<SatParameters>());
1064
1065 // Simple linear scan algorithm to find the optimal.
1066 while (true) {
1068 if (result != SatSolver::FEASIBLE) return result;
1069
1070 // The objective is the current lower bound of the objective_var.
1071 const IntegerValue objective = integer_trail->LowerBound(objective_var);
1072
1073 // We have a solution!
1074 if (feasible_solution_observer != nullptr) {
1075 feasible_solution_observer();
1076 }
1077 if (parameters.stop_after_first_solution()) {
1079 }
1080
1081 // Restrict the objective.
1082 sat_solver->Backtrack(0);
1083 if (!integer_trail->Enqueue(
1084 IntegerLiteral::LowerOrEqual(objective_var, objective - 1), {},
1085 {})) {
1086 return SatSolver::INFEASIBLE;
1087 }
1088 }
1089}
1090
1092 IntegerVariable objective_var,
1093 const std::function<void()>& feasible_solution_observer, Model* model) {
1094 const SatParameters old_params = *model->GetOrCreate<SatParameters>();
1095 SatSolver* sat_solver = model->GetOrCreate<SatSolver>();
1096 IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
1097 IntegerEncoder* integer_encoder = model->GetOrCreate<IntegerEncoder>();
1098
1099 // Set the requested conflict limit.
1100 {
1101 SatParameters new_params = old_params;
1102 new_params.set_max_number_of_conflicts(
1103 old_params.binary_search_num_conflicts());
1104 *model->GetOrCreate<SatParameters>() = new_params;
1105 }
1106
1107 // The assumption (objective <= value) for values in
1108 // [unknown_min, unknown_max] reached the conflict limit.
1109 bool loop = true;
1110 IntegerValue unknown_min = integer_trail->UpperBound(objective_var);
1111 IntegerValue unknown_max = integer_trail->LowerBound(objective_var);
1112 while (loop) {
1113 sat_solver->Backtrack(0);
1114 const IntegerValue lb = integer_trail->LowerBound(objective_var);
1115 const IntegerValue ub = integer_trail->UpperBound(objective_var);
1116 unknown_min = std::min(unknown_min, ub);
1117 unknown_max = std::max(unknown_max, lb);
1118
1119 // We first refine the lower bound and then the upper bound.
1120 IntegerValue target;
1121 if (lb < unknown_min) {
1122 target = lb + (unknown_min - lb) / 2;
1123 } else if (unknown_max < ub) {
1124 target = ub - (ub - unknown_max) / 2;
1125 } else {
1126 VLOG(1) << "Binary-search, done.";
1127 break;
1128 }
1129 VLOG(1) << "Binary-search, objective: [" << lb << "," << ub << "]"
1130 << " tried: [" << unknown_min << "," << unknown_max << "]"
1131 << " target: obj<=" << target;
1132 SatSolver::Status result;
1133 if (target < ub) {
1134 const Literal assumption = integer_encoder->GetOrCreateAssociatedLiteral(
1135 IntegerLiteral::LowerOrEqual(objective_var, target));
1136 result = ResetAndSolveIntegerProblem({assumption}, model);
1137 } else {
1138 result = ResetAndSolveIntegerProblem({}, model);
1139 }
1140
1141 switch (result) {
1142 case SatSolver::INFEASIBLE: {
1143 loop = false;
1144 break;
1145 }
1147 // Update the objective lower bound.
1148 sat_solver->Backtrack(0);
1149 if (!integer_trail->Enqueue(
1150 IntegerLiteral::GreaterOrEqual(objective_var, target + 1), {},
1151 {})) {
1152 loop = false;
1153 }
1154 break;
1155 }
1156 case SatSolver::FEASIBLE: {
1157 // The objective is the current lower bound of the objective_var.
1158 const IntegerValue objective = integer_trail->LowerBound(objective_var);
1159 if (feasible_solution_observer != nullptr) {
1160 feasible_solution_observer();
1161 }
1162
1163 // We have a solution, restrict the objective upper bound to only look
1164 // for better ones now.
1165 sat_solver->Backtrack(0);
1166 if (!integer_trail->Enqueue(
1167 IntegerLiteral::LowerOrEqual(objective_var, objective - 1), {},
1168 {})) {
1169 loop = false;
1170 }
1171 break;
1172 }
1174 unknown_min = std::min(target, unknown_min);
1175 unknown_max = std::max(target, unknown_max);
1176 break;
1177 }
1178 }
1179 }
1180
1181 sat_solver->Backtrack(0);
1182 *model->GetOrCreate<SatParameters>() = old_params;
1183}
1184
1185namespace {
1186
1187// If the given model is unsat under the given assumptions, returns one or more
1188// non-overlapping set of assumptions, each set making the problem infeasible on
1189// its own (the cores).
1190//
1191// In presence of weights, we "generalize" the notions of disjoints core using
1192// the WCE idea describe in "Weight-Aware Core Extraction in SAT-Based MaxSAT
1193// solving" Jeremias Berg And Matti Jarvisalo.
1194//
1195// The returned status can be either:
1196// - ASSUMPTIONS_UNSAT if the set of returned core perfectly cover the given
1197// assumptions, in this case, we don't bother trying to find a SAT solution
1198// with no assumptions.
1199// - FEASIBLE if after finding zero or more core we have a solution.
1200// - LIMIT_REACHED if we reached the time-limit before one of the two status
1201// above could be decided.
1202//
1203// TODO(user): There is many way to combine the WCE and stratification
1204// heuristics. I didn't had time to properly compare the different approach. See
1205// the WCE papers for some ideas, but there is many more ways to try to find a
1206// lot of core at once and try to minimize the minimum weight of each of the
1207// cores.
1208SatSolver::Status FindCores(std::vector<Literal> assumptions,
1209 std::vector<IntegerValue> assumption_weights,
1210 IntegerValue stratified_threshold, Model* model,
1211 std::vector<std::vector<Literal>>* cores) {
1212 cores->clear();
1213 SatSolver* sat_solver = model->GetOrCreate<SatSolver>();
1214 TimeLimit* limit = model->GetOrCreate<TimeLimit>();
1215 do {
1216 if (limit->LimitReached()) return SatSolver::LIMIT_REACHED;
1217
1218 const SatSolver::Status result =
1219 ResetAndSolveIntegerProblem(assumptions, model);
1220 if (result != SatSolver::ASSUMPTIONS_UNSAT) return result;
1221 std::vector<Literal> core = sat_solver->GetLastIncompatibleDecisions();
1222 if (sat_solver->parameters().minimize_core()) {
1223 MinimizeCoreWithPropagation(limit, sat_solver, &core);
1224 }
1225 CHECK(!core.empty());
1226 cores->push_back(core);
1227 if (!sat_solver->parameters().find_multiple_cores()) break;
1228
1229 // Recover the original indices of the assumptions that are part of the
1230 // core.
1231 std::vector<int> indices;
1232 {
1233 std::set<Literal> temp(core.begin(), core.end());
1234 for (int i = 0; i < assumptions.size(); ++i) {
1235 if (gtl::ContainsKey(temp, assumptions[i])) {
1236 indices.push_back(i);
1237 }
1238 }
1239 }
1240
1241 // Remove min_weight from the weights of all the assumptions in the core.
1242 //
1243 // TODO(user): push right away the objective bound by that much? This should
1244 // be better in a multi-threading context as we can share more quickly the
1245 // better bound.
1246 IntegerValue min_weight = assumption_weights[indices.front()];
1247 for (const int i : indices) {
1248 min_weight = std::min(min_weight, assumption_weights[i]);
1249 }
1250 for (const int i : indices) {
1251 assumption_weights[i] -= min_weight;
1252 }
1253
1254 // Remove from assumptions all the one with a new weight smaller than the
1255 // current stratification threshold and see if we can find another core.
1256 int new_size = 0;
1257 for (int i = 0; i < assumptions.size(); ++i) {
1258 if (assumption_weights[i] < stratified_threshold) continue;
1259 assumptions[new_size] = assumptions[i];
1260 assumption_weights[new_size] = assumption_weights[i];
1261 ++new_size;
1262 }
1263 assumptions.resize(new_size);
1264 assumption_weights.resize(new_size);
1265 } while (!assumptions.empty());
1267}
1268
1269// Slightly different algo than FindCores() which aim to extract more cores, but
1270// not necessarily non-overlaping ones.
1271SatSolver::Status FindMultipleCoresForMaxHs(
1272 std::vector<Literal> assumptions, Model* model,
1273 std::vector<std::vector<Literal>>* cores) {
1274 cores->clear();
1275 SatSolver* sat_solver = model->GetOrCreate<SatSolver>();
1276 TimeLimit* limit = model->GetOrCreate<TimeLimit>();
1277 do {
1278 if (limit->LimitReached()) return SatSolver::LIMIT_REACHED;
1279
1280 const SatSolver::Status result =
1281 ResetAndSolveIntegerProblem(assumptions, model);
1282 if (result != SatSolver::ASSUMPTIONS_UNSAT) return result;
1283 std::vector<Literal> core = sat_solver->GetLastIncompatibleDecisions();
1284 if (sat_solver->parameters().minimize_core()) {
1285 MinimizeCoreWithPropagation(limit, sat_solver, &core);
1286 }
1287 CHECK(!core.empty());
1288 cores->push_back(core);
1289 if (!sat_solver->parameters().find_multiple_cores()) break;
1290
1291 // Pick a random literal from the core and remove it from the set of
1292 // assumptions.
1293 CHECK(!core.empty());
1294 auto* random = model->GetOrCreate<ModelRandomGenerator>();
1295 const Literal random_literal =
1296 core[absl::Uniform<int>(*random, 0, core.size())];
1297 for (int i = 0; i < assumptions.size(); ++i) {
1298 if (assumptions[i] == random_literal) {
1299 std::swap(assumptions[i], assumptions.back());
1300 assumptions.pop_back();
1301 break;
1302 }
1303 }
1304 } while (!assumptions.empty());
1306}
1307
1308} // namespace
1309
1311 IntegerVariable objective_var,
1312 const std::vector<IntegerVariable>& variables,
1313 const std::vector<IntegerValue>& coefficients,
1314 std::function<void()> feasible_solution_observer, Model* model)
1315 : parameters_(model->GetOrCreate<SatParameters>()),
1316 sat_solver_(model->GetOrCreate<SatSolver>()),
1317 time_limit_(model->GetOrCreate<TimeLimit>()),
1318 integer_trail_(model->GetOrCreate<IntegerTrail>()),
1319 integer_encoder_(model->GetOrCreate<IntegerEncoder>()),
1320 model_(model),
1321 objective_var_(objective_var),
1322 feasible_solution_observer_(std::move(feasible_solution_observer)) {
1323 CHECK_EQ(variables.size(), coefficients.size());
1324 for (int i = 0; i < variables.size(); ++i) {
1325 if (coefficients[i] > 0) {
1326 terms_.push_back({variables[i], coefficients[i]});
1327 } else if (coefficients[i] < 0) {
1328 terms_.push_back({NegationOf(variables[i]), -coefficients[i]});
1329 } else {
1330 continue; // coefficients[i] == 0
1331 }
1332 terms_.back().depth = 0;
1333 }
1334
1335 // This is used by the "stratified" approach. We will only consider terms with
1336 // a weight not lower than this threshold. The threshold will decrease as the
1337 // algorithm progress.
1338 stratification_threshold_ = parameters_->max_sat_stratification() ==
1339 SatParameters::STRATIFICATION_NONE
1340 ? IntegerValue(1)
1342}
1343
1344bool CoreBasedOptimizer::ProcessSolution() {
1345 // We don't assume that objective_var is linked with its linear term, so
1346 // we recompute the objective here.
1347 IntegerValue objective(0);
1348 for (ObjectiveTerm& term : terms_) {
1349 const IntegerValue value = integer_trail_->LowerBound(term.var);
1350 objective += term.weight * value;
1351
1352 // Also keep in term.cover_ub the minimum value for term.var that we have
1353 // seens amongst all the feasible solutions found so far.
1354 term.cover_ub = std::min(term.cover_ub, value);
1355 }
1356
1357 // We use the level zero upper bound of the objective to indicate an upper
1358 // limit for the solution objective we are looking for. Again, because the
1359 // objective_var is not assumed to be linked, it could take any value in the
1360 // current solution.
1361 if (objective > integer_trail_->LevelZeroUpperBound(objective_var_)) {
1362 return true;
1363 }
1364
1365 if (feasible_solution_observer_ != nullptr) {
1366 feasible_solution_observer_();
1367 }
1368 if (parameters_->stop_after_first_solution()) {
1369 stop_ = true;
1370 }
1371
1372 // Constrain objective_var. This has a better result when objective_var is
1373 // used in an LP relaxation for instance.
1374 sat_solver_->Backtrack(0);
1375 sat_solver_->SetAssumptionLevel(0);
1376 return integer_trail_->Enqueue(
1377 IntegerLiteral::LowerOrEqual(objective_var_, objective - 1), {}, {});
1378}
1379
1380bool CoreBasedOptimizer::PropagateObjectiveBounds() {
1381 // We assumes all terms (modulo stratification) at their lower-bound.
1382 bool some_bound_were_tightened = true;
1383 while (some_bound_were_tightened) {
1384 some_bound_were_tightened = false;
1385 if (!sat_solver_->ResetToLevelZero()) return false;
1386
1387 // Compute implied lb.
1388 IntegerValue implied_objective_lb(0);
1389 for (ObjectiveTerm& term : terms_) {
1390 const IntegerValue var_lb = integer_trail_->LowerBound(term.var);
1391 term.old_var_lb = var_lb;
1392 implied_objective_lb += term.weight * var_lb.value();
1393 }
1394
1395 // Update the objective lower bound with our current bound.
1396 if (implied_objective_lb > integer_trail_->LowerBound(objective_var_)) {
1397 if (!integer_trail_->Enqueue(IntegerLiteral::GreaterOrEqual(
1398 objective_var_, implied_objective_lb),
1399 {}, {})) {
1400 return false;
1401 }
1402
1403 some_bound_were_tightened = true;
1404 }
1405
1406 // The gap is used to propagate the upper-bound of all variable that are
1407 // in the current objective (Exactly like done in the propagation of a
1408 // linear constraint with the slack). When this fix a variable to its
1409 // lower bound, it is called "hardening" in the max-sat literature. This
1410 // has a really beneficial effect on some weighted max-sat problems like
1411 // the haplotyping-pedigrees ones.
1412 const IntegerValue gap =
1413 integer_trail_->UpperBound(objective_var_) - implied_objective_lb;
1414
1415 for (const ObjectiveTerm& term : terms_) {
1416 if (term.weight == 0) continue;
1417 const IntegerValue var_lb = integer_trail_->LowerBound(term.var);
1418 const IntegerValue var_ub = integer_trail_->UpperBound(term.var);
1419 if (var_lb == var_ub) continue;
1420
1421 // Hardening. This basically just propagate the implied upper bound on
1422 // term.var from the current best solution. Note that the gap is
1423 // non-negative and the weight positive here. The test is done in order
1424 // to avoid any integer overflow provided (ub - lb) do not overflow, but
1425 // this is a precondition in our cp-model.
1426 if (gap / term.weight < var_ub - var_lb) {
1427 some_bound_were_tightened = true;
1428 const IntegerValue new_ub = var_lb + gap / term.weight;
1429 CHECK_LT(new_ub, var_ub);
1430 CHECK(!integer_trail_->IsCurrentlyIgnored(term.var));
1431 if (!integer_trail_->Enqueue(
1432 IntegerLiteral::LowerOrEqual(term.var, new_ub), {}, {})) {
1433 return false;
1434 }
1435 }
1436 }
1437 }
1438 return true;
1439}
1440
1441// A basic algorithm is to take the next one, or at least the next one
1442// that invalidate the current solution. But to avoid corner cases for
1443// problem with a lot of terms all with different objective weights (in
1444// which case we will kind of introduce only one assumption per loop
1445// which is little), we use an heuristic and take the 90% percentile of
1446// the unique weights not yet included.
1447//
1448// TODO(user): There is many other possible heuristics here, and I
1449// didn't have the time to properly compare them.
1450void CoreBasedOptimizer::ComputeNextStratificationThreshold() {
1451 std::vector<IntegerValue> weights;
1452 for (ObjectiveTerm& term : terms_) {
1453 if (term.weight >= stratification_threshold_) continue;
1454 if (term.weight == 0) continue;
1455
1456 const IntegerValue var_lb = integer_trail_->LevelZeroLowerBound(term.var);
1457 const IntegerValue var_ub = integer_trail_->LevelZeroUpperBound(term.var);
1458 if (var_lb == var_ub) continue;
1459
1460 weights.push_back(term.weight);
1461 }
1462 if (weights.empty()) {
1463 stratification_threshold_ = IntegerValue(0);
1464 return;
1465 }
1466
1468 stratification_threshold_ =
1469 weights[static_cast<int>(std::floor(0.9 * weights.size()))];
1470}
1471
1472bool CoreBasedOptimizer::CoverOptimization() {
1473 // We set a fix deterministic time limit per all sub-solve and skip to the
1474 // next core if the sum of the subsolve is also over this limit.
1475 constexpr double max_dtime_per_core = 0.5;
1476 const double old_time_limit = parameters_->max_deterministic_time();
1477 parameters_->set_max_deterministic_time(max_dtime_per_core);
1478 auto cleanup = ::absl::MakeCleanup([old_time_limit, this]() {
1479 parameters_->set_max_deterministic_time(old_time_limit);
1480 });
1481
1482 for (const ObjectiveTerm& term : terms_) {
1483 // We currently skip the initial objective terms as there could be many
1484 // of them. TODO(user): provide an option to cover-optimize them? I
1485 // fear that this will slow down the solver too much though.
1486 if (term.depth == 0) continue;
1487
1488 // Find out the true lower bound of var. This is called "cover
1489 // optimization" in some of the max-SAT literature. It can helps on some
1490 // problem families and hurt on others, but the overall impact is
1491 // positive.
1492 const IntegerVariable var = term.var;
1493 IntegerValue best =
1494 std::min(term.cover_ub, integer_trail_->UpperBound(var));
1495
1496 // Note(user): this can happen in some corner case because each time we
1497 // find a solution, we constrain the objective to be smaller than it, so
1498 // it is possible that a previous best is now infeasible.
1499 if (best <= integer_trail_->LowerBound(var)) continue;
1500
1501 // Compute the global deterministic time for this core cover
1502 // optimization.
1503 const double deterministic_limit =
1504 time_limit_->GetElapsedDeterministicTime() + max_dtime_per_core;
1505
1506 // Simple linear scan algorithm to find the optimal of var.
1507 SatSolver::Status result;
1508 while (best > integer_trail_->LowerBound(var)) {
1509 const Literal assumption = integer_encoder_->GetOrCreateAssociatedLiteral(
1511 result = ResetAndSolveIntegerProblem({assumption}, model_);
1512 if (result != SatSolver::FEASIBLE) break;
1513
1514 best = integer_trail_->LowerBound(var);
1515 VLOG(1) << "cover_opt var:" << var << " domain:["
1516 << integer_trail_->LevelZeroLowerBound(var) << "," << best << "]";
1517 if (!ProcessSolution()) return false;
1518 if (!sat_solver_->ResetToLevelZero()) return false;
1519 if (stop_ ||
1520 time_limit_->GetElapsedDeterministicTime() > deterministic_limit) {
1521 break;
1522 }
1523 }
1524 if (result == SatSolver::INFEASIBLE) return false;
1525 if (result == SatSolver::ASSUMPTIONS_UNSAT) {
1526 // TODO(user): If we improve the lower bound of var, we should check
1527 // if our global lower bound reached our current best solution in
1528 // order to abort early if the optimal is proved.
1529 if (!integer_trail_->Enqueue(IntegerLiteral::GreaterOrEqual(var, best),
1530 {}, {})) {
1531 return false;
1532 }
1533 }
1534 }
1535
1536 if (!PropagateObjectiveBounds()) return false;
1537 return true;
1538}
1539
1541 // TODO(user): The core is returned in the same order as the assumptions,
1542 // so we don't really need this map, we could just do a linear scan to
1543 // recover which node are part of the core. This however needs to be properly
1544 // unit tested before usage.
1545 std::map<LiteralIndex, int> literal_to_term_index;
1546
1547 // Start the algorithm.
1548 stop_ = false;
1549 while (true) {
1550 // TODO(user): This always resets the solver to level zero.
1551 // Because of that we don't resume a solve in "chunk" perfectly. Fix.
1552 if (!PropagateObjectiveBounds()) return SatSolver::INFEASIBLE;
1553
1554 // Bulk cover optimization.
1555 //
1556 // TODO(user): If the search is aborted during this phase and we solve in
1557 // "chunk", we don't resume perfectly from where it was. Fix.
1558 if (parameters_->cover_optimization()) {
1559 if (!CoverOptimization()) return SatSolver::INFEASIBLE;
1560 if (stop_) return SatSolver::LIMIT_REACHED;
1561 }
1562
1563 // We assumes all terms (modulo stratification) at their lower-bound.
1564 std::vector<int> term_indices;
1565 std::vector<IntegerLiteral> integer_assumptions;
1566 std::vector<IntegerValue> assumption_weights;
1567 IntegerValue objective_offset(0);
1568 bool some_assumptions_were_skipped = false;
1569 for (int i = 0; i < terms_.size(); ++i) {
1570 const ObjectiveTerm term = terms_[i];
1571
1572 // TODO(user): These can be simply removed from the list.
1573 if (term.weight == 0) continue;
1574
1575 // Skip fixed terms.
1576 // We still keep them around for a proper lower-bound computation.
1577 //
1578 // TODO(user): we could keep an objective offset instead.
1579 const IntegerValue var_lb = integer_trail_->LowerBound(term.var);
1580 const IntegerValue var_ub = integer_trail_->UpperBound(term.var);
1581 if (var_lb == var_ub) {
1582 objective_offset += term.weight * var_lb.value();
1583 continue;
1584 }
1585
1586 // Only consider the terms above the threshold.
1587 if (term.weight >= stratification_threshold_) {
1588 integer_assumptions.push_back(
1589 IntegerLiteral::LowerOrEqual(term.var, var_lb));
1590 assumption_weights.push_back(term.weight);
1591 term_indices.push_back(i);
1592 } else {
1593 some_assumptions_were_skipped = true;
1594 }
1595 }
1596
1597 // No assumptions with the current stratification? use the next one.
1598 if (term_indices.empty() && some_assumptions_were_skipped) {
1599 ComputeNextStratificationThreshold();
1600 continue;
1601 }
1602
1603 // If there is only one or two assumptions left, we switch the algorithm.
1604 if (term_indices.size() <= 2 && !some_assumptions_were_skipped) {
1605 VLOG(1) << "Switching to linear scan...";
1606 if (!already_switched_to_linear_scan_) {
1607 already_switched_to_linear_scan_ = true;
1608 std::vector<IntegerVariable> constraint_vars;
1609 std::vector<int64> constraint_coeffs;
1610 for (const int index : term_indices) {
1611 constraint_vars.push_back(terms_[index].var);
1612 constraint_coeffs.push_back(terms_[index].weight.value());
1613 }
1614 constraint_vars.push_back(objective_var_);
1615 constraint_coeffs.push_back(-1);
1616 model_->Add(WeightedSumLowerOrEqual(constraint_vars, constraint_coeffs,
1617 -objective_offset.value()));
1618 }
1619
1621 objective_var_, feasible_solution_observer_, model_);
1622 }
1623
1624 // Display the progress.
1625 if (VLOG_IS_ON(1)) {
1626 int max_depth = 0;
1627 for (const ObjectiveTerm& term : terms_) {
1628 max_depth = std::max(max_depth, term.depth);
1629 }
1630 const int64 lb = integer_trail_->LowerBound(objective_var_).value();
1631 const int64 ub = integer_trail_->UpperBound(objective_var_).value();
1632 const int gap =
1633 lb == ub
1634 ? 0
1635 : static_cast<int>(std::ceil(
1636 100.0 * (ub - lb) / std::max(std::abs(ub), std::abs(lb))));
1637 VLOG(1) << absl::StrCat("unscaled_next_obj_range:[", lb, ",", ub,
1638 "]"
1639 " gap:",
1640 gap, "%", " assumptions:", term_indices.size(),
1641 " strat:", stratification_threshold_.value(),
1642 " depth:", max_depth);
1643 }
1644
1645 // Convert integer_assumptions to Literals.
1646 std::vector<Literal> assumptions;
1647 literal_to_term_index.clear();
1648 for (int i = 0; i < integer_assumptions.size(); ++i) {
1649 assumptions.push_back(integer_encoder_->GetOrCreateAssociatedLiteral(
1650 integer_assumptions[i]));
1651
1652 // Tricky: In some rare case, it is possible that the same literal
1653 // correspond to more that one assumptions. In this case, we can just
1654 // pick one of them when converting back a core to term indices.
1655 //
1656 // TODO(user): We can probably be smarter about the cost of the
1657 // assumptions though.
1658 literal_to_term_index[assumptions.back().Index()] = term_indices[i];
1659 }
1660
1661 // Solve under the assumptions.
1662 //
1663 // TODO(user): If the "search" is interupted while computing cores, we
1664 // currently do not resume it flawlessly. We however add any cores we found
1665 // before aborting.
1666 std::vector<std::vector<Literal>> cores;
1667 const SatSolver::Status result =
1668 FindCores(assumptions, assumption_weights, stratification_threshold_,
1669 model_, &cores);
1670 if (result == SatSolver::INFEASIBLE) return SatSolver::INFEASIBLE;
1671 if (result == SatSolver::FEASIBLE) {
1672 if (!ProcessSolution()) return SatSolver::INFEASIBLE;
1673 if (stop_) return SatSolver::LIMIT_REACHED;
1674 if (cores.empty()) {
1675 ComputeNextStratificationThreshold();
1676 if (stratification_threshold_ == 0) return SatSolver::INFEASIBLE;
1677 continue;
1678 }
1679 }
1680
1681 // Process the cores by creating new variables and transferring the minimum
1682 // weight of each core to it.
1683 if (!sat_solver_->ResetToLevelZero()) return SatSolver::INFEASIBLE;
1684 for (const std::vector<Literal>& core : cores) {
1685 // This just increase the lower-bound of the corresponding node, which
1686 // should already be done by the solver.
1687 if (core.size() == 1) continue;
1688
1689 // Compute the min weight of all the terms in the core. The lower bound
1690 // will be increased by that much because at least one assumption in the
1691 // core must be true. This is also why we can start at 1 for new_var_lb.
1692 bool ignore_this_core = false;
1693 IntegerValue min_weight = kMaxIntegerValue;
1694 IntegerValue max_weight(0);
1695 IntegerValue new_var_lb(1);
1696 IntegerValue new_var_ub(0);
1697 int new_depth = 0;
1698 for (const Literal lit : core) {
1699 const int index = gtl::FindOrDie(literal_to_term_index, lit.Index());
1700
1701 // When this happen, the core is now trivially "minimized" by the new
1702 // bound on this variable, so there is no point in adding it.
1703 if (terms_[index].old_var_lb <
1704 integer_trail_->LowerBound(terms_[index].var)) {
1705 ignore_this_core = true;
1706 break;
1707 }
1708
1709 const IntegerValue weight = terms_[index].weight;
1710 min_weight = std::min(min_weight, weight);
1711 max_weight = std::max(max_weight, weight);
1712 new_depth = std::max(new_depth, terms_[index].depth + 1);
1713 new_var_lb += integer_trail_->LowerBound(terms_[index].var);
1714 new_var_ub += integer_trail_->UpperBound(terms_[index].var);
1715 }
1716 if (ignore_this_core) continue;
1717
1718 VLOG(1) << absl::StrFormat(
1719 "core:%u weight:[%d,%d] domain:[%d,%d] depth:%d", core.size(),
1720 min_weight.value(), max_weight.value(), new_var_lb.value(),
1721 new_var_ub.value(), new_depth);
1722
1723 // We will "transfer" min_weight from all the variables of the core
1724 // to a new variable.
1725 const IntegerVariable new_var =
1726 integer_trail_->AddIntegerVariable(new_var_lb, new_var_ub);
1727 terms_.push_back({new_var, min_weight, new_depth});
1728 terms_.back().cover_ub = new_var_ub;
1729
1730 // Sum variables in the core <= new_var.
1731 {
1732 std::vector<IntegerVariable> constraint_vars;
1733 std::vector<int64> constraint_coeffs;
1734 for (const Literal lit : core) {
1735 const int index = gtl::FindOrDie(literal_to_term_index, lit.Index());
1736 terms_[index].weight -= min_weight;
1737 constraint_vars.push_back(terms_[index].var);
1738 constraint_coeffs.push_back(1);
1739 }
1740 constraint_vars.push_back(new_var);
1741 constraint_coeffs.push_back(-1);
1742 model_->Add(
1743 WeightedSumLowerOrEqual(constraint_vars, constraint_coeffs, 0));
1744 }
1745 }
1746
1747 // Abort if we reached the time limit. Note that we still add any cores we
1748 // found in case the solve is splitted in "chunk".
1749 if (result == SatSolver::LIMIT_REACHED) return result;
1750 }
1751}
1752
1753// TODO(user): take the MPModelRequest or MPModelProto directly, so that we can
1754// have initial constraints!
1755//
1756// TODO(user): remove code duplication with MinimizeWithCoreAndLazyEncoding();
1758 const ObjectiveDefinition& objective_definition,
1759 const std::function<void()>& feasible_solution_observer, Model* model) {
1760#if !defined(__PORTABLE_PLATFORM__) && defined(USE_SCIP)
1761
1762 IntegerVariable objective_var = objective_definition.objective_var;
1763 std::vector<IntegerVariable> variables = objective_definition.vars;
1764 std::vector<IntegerValue> coefficients = objective_definition.coeffs;
1765
1766 SatSolver* sat_solver = model->GetOrCreate<SatSolver>();
1767 IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
1768 IntegerEncoder* integer_encoder = model->GetOrCreate<IntegerEncoder>();
1769
1770 // This will be called each time a feasible solution is found.
1771 const auto process_solution = [&]() {
1772 // We don't assume that objective_var is linked with its linear term, so
1773 // we recompute the objective here.
1774 IntegerValue objective(0);
1775 for (int i = 0; i < variables.size(); ++i) {
1776 objective +=
1777 coefficients[i] * IntegerValue(model->Get(Value(variables[i])));
1778 }
1779 if (objective > integer_trail->UpperBound(objective_var)) return true;
1780
1781 if (feasible_solution_observer != nullptr) {
1782 feasible_solution_observer();
1783 }
1784
1785 // Constrain objective_var. This has a better result when objective_var is
1786 // used in an LP relaxation for instance.
1787 sat_solver->Backtrack(0);
1788 sat_solver->SetAssumptionLevel(0);
1789 if (!integer_trail->Enqueue(
1790 IntegerLiteral::LowerOrEqual(objective_var, objective - 1), {},
1791 {})) {
1792 return false;
1793 }
1794 return true;
1795 };
1796
1797 // This is the "generalized" hitting set problem we will solve. Each time
1798 // we find a core, a new constraint will be added to this problem.
1799 MPModelRequest request;
1800 request.set_solver_specific_parameters("limits/gap = 0");
1801 request.set_solver_type(MPModelRequest::SCIP_MIXED_INTEGER_PROGRAMMING);
1802
1803 MPModelProto& hs_model = *request.mutable_model();
1804 const int num_variables_in_objective = variables.size();
1805 for (int i = 0; i < num_variables_in_objective; ++i) {
1806 if (coefficients[i] < 0) {
1807 variables[i] = NegationOf(variables[i]);
1808 coefficients[i] = -coefficients[i];
1809 }
1810 const IntegerVariable var = variables[i];
1811 MPVariableProto* var_proto = hs_model.add_variable();
1812 var_proto->set_lower_bound(integer_trail->LowerBound(var).value());
1813 var_proto->set_upper_bound(integer_trail->UpperBound(var).value());
1814 var_proto->set_objective_coefficient(coefficients[i].value());
1815 var_proto->set_is_integer(true);
1816 }
1817
1818 MPSolutionResponse response;
1819
1820 // This is used by the "stratified" approach. We will only consider terms with
1821 // a weight not lower than this threshold. The threshold will decrease as the
1822 // algorithm progress.
1823 IntegerValue stratified_threshold = kMaxIntegerValue;
1824
1825 // TODO(user): The core is returned in the same order as the assumptions,
1826 // so we don't really need this map, we could just do a linear scan to
1827 // recover which node are part of the core.
1828 std::map<LiteralIndex, std::vector<int>> assumption_to_indices;
1829
1830 // New Booleans variable in the MIP model to represent X >= cte.
1831 std::map<std::pair<int, double>, int> created_var;
1832
1833 const SatParameters& parameters = *(model->GetOrCreate<SatParameters>());
1834
1835 // Start the algorithm.
1836 SatSolver::Status result;
1837 for (int iter = 0;; ++iter) {
1838 // TODO(user): Even though we keep the same solver, currently the solve is
1839 // not really done incrementally. It might be hard to improve though.
1840 //
1841 // TODO(user): deal with time limit.
1843 CHECK_EQ(response.status(), MPSolverResponseStatus::MPSOLVER_OPTIMAL);
1844
1845 const IntegerValue mip_objective(
1846 static_cast<int64>(std::round(response.objective_value())));
1847 VLOG(1) << "constraints: " << hs_model.constraint_size()
1848 << " variables: " << hs_model.variable_size() << " hs_lower_bound: "
1849 << objective_definition.ScaleIntegerObjective(mip_objective)
1850 << " strat: " << stratified_threshold;
1851
1852 // Update the objective lower bound with our current bound.
1853 //
1854 // Note(user): This is not needed for correctness, but it might cause
1855 // more propagation and is nice to have for reporting/logging purpose.
1856 if (!integer_trail->Enqueue(
1857 IntegerLiteral::GreaterOrEqual(objective_var, mip_objective), {},
1858 {})) {
1859 result = SatSolver::INFEASIBLE;
1860 break;
1861 }
1862
1863 sat_solver->Backtrack(0);
1864 sat_solver->SetAssumptionLevel(0);
1865 std::vector<Literal> assumptions;
1866 assumption_to_indices.clear();
1867 IntegerValue next_stratified_threshold(0);
1868 for (int i = 0; i < num_variables_in_objective; ++i) {
1869 const IntegerValue hs_value(
1870 static_cast<int64>(response.variable_value(i)));
1871 if (hs_value == integer_trail->UpperBound(variables[i])) continue;
1872
1873 // Only consider the terms above the threshold.
1874 if (coefficients[i] < stratified_threshold) {
1875 next_stratified_threshold =
1876 std::max(next_stratified_threshold, coefficients[i]);
1877 } else {
1878 // It is possible that different variables have the same associated
1879 // literal. So we do need to consider this case.
1880 assumptions.push_back(integer_encoder->GetOrCreateAssociatedLiteral(
1881 IntegerLiteral::LowerOrEqual(variables[i], hs_value)));
1882 assumption_to_indices[assumptions.back().Index()].push_back(i);
1883 }
1884 }
1885
1886 // No assumptions with the current stratified_threshold? use the new one.
1887 if (assumptions.empty() && next_stratified_threshold > 0) {
1888 CHECK_LT(next_stratified_threshold, stratified_threshold);
1889 stratified_threshold = next_stratified_threshold;
1890 --iter; // "false" iteration, the lower bound does not increase.
1891 continue;
1892 }
1893
1894 // TODO(user): we could also randomly shuffle the assumptions to find more
1895 // cores for only one MIP solve.
1896 //
1897 // TODO(user): Use the real weights and exploit the extra cores.
1898 std::vector<std::vector<Literal>> cores;
1899 result = FindMultipleCoresForMaxHs(assumptions, model, &cores);
1900 if (result == SatSolver::FEASIBLE) {
1901 if (!process_solution()) return SatSolver::INFEASIBLE;
1902 if (parameters.stop_after_first_solution()) {
1904 }
1905 if (cores.empty()) {
1906 // If not all assumptions were taken, continue with a lower stratified
1907 // bound. Otherwise we have an optimal solution.
1908 stratified_threshold = next_stratified_threshold;
1909 if (stratified_threshold == 0) break;
1910 --iter; // "false" iteration, the lower bound does not increase.
1911 continue;
1912 }
1913 } else if (result != SatSolver::ASSUMPTIONS_UNSAT) {
1914 break;
1915 }
1916
1917 sat_solver->Backtrack(0);
1918 sat_solver->SetAssumptionLevel(0);
1919 for (const std::vector<Literal>& core : cores) {
1920 if (core.size() == 1) {
1921 for (const int index :
1922 gtl::FindOrDie(assumption_to_indices, core.front().Index())) {
1923 hs_model.mutable_variable(index)->set_lower_bound(
1924 integer_trail->LowerBound(variables[index]).value());
1925 }
1926 continue;
1927 }
1928
1929 // Add the corresponding constraint to hs_model.
1930 MPConstraintProto* ct = hs_model.add_constraint();
1931 ct->set_lower_bound(1.0);
1932 for (const Literal lit : core) {
1933 for (const int index :
1934 gtl::FindOrDie(assumption_to_indices, lit.Index())) {
1935 const double lb = integer_trail->LowerBound(variables[index]).value();
1936 const double hs_value = response.variable_value(index);
1937 if (hs_value == lb) {
1938 ct->add_var_index(index);
1939 ct->add_coefficient(1.0);
1940 ct->set_lower_bound(ct->lower_bound() + lb);
1941 } else {
1942 const std::pair<int, double> key = {index, hs_value};
1943 if (!gtl::ContainsKey(created_var, key)) {
1944 const int new_var_index = hs_model.variable_size();
1945 created_var[key] = new_var_index;
1946
1947 MPVariableProto* new_var = hs_model.add_variable();
1948 new_var->set_lower_bound(0);
1949 new_var->set_upper_bound(1);
1950 new_var->set_is_integer(true);
1951
1952 // (new_var == 1) => x > hs_value.
1953 // (x - lb) - (hs_value - lb + 1) * new_var >= 0.
1954 MPConstraintProto* implication = hs_model.add_constraint();
1955 implication->set_lower_bound(lb);
1956 implication->add_var_index(index);
1957 implication->add_coefficient(1.0);
1958 implication->add_var_index(new_var_index);
1959 implication->add_coefficient(lb - hs_value - 1);
1960 }
1961 ct->add_var_index(gtl::FindOrDieNoPrint(created_var, key));
1962 ct->add_coefficient(1.0);
1963 }
1964 }
1965 }
1966 }
1967 }
1968
1969 return result;
1970#else // !__PORTABLE_PLATFORM__ && USE_SCIP
1971 LOG(FATAL) << "Not supported.";
1972#endif // !__PORTABLE_PLATFORM__ && USE_SCIP
1973}
1974
1975} // namespace sat
1976} // namespace operations_research
int64 min
Definition: alldiff_cst.cc:138
int64 max
Definition: alldiff_cst.cc:139
#define CHECK(condition)
Definition: base/logging.h:495
#define CHECK_LT(val1, val2)
Definition: base/logging.h:700
#define CHECK_EQ(val1, val2)
Definition: base/logging.h:697
#define CHECK_GE(val1, val2)
Definition: base/logging.h:701
#define CHECK_GT(val1, val2)
Definition: base/logging.h:702
#define CHECK_NE(val1, val2)
Definition: base/logging.h:698
#define LOG(severity)
Definition: base/logging.h:420
#define DCHECK(condition)
Definition: base/logging.h:884
#define CHECK_LE(val1, val2)
Definition: base/logging.h:699
#define VLOG(verboselevel)
Definition: base/logging.h:978
static void SolveWithProto(const MPModelRequest &model_request, MPSolutionResponse *response)
Solves the model encoded by a MPModelRequest protocol buffer and fills the solution encoded as a MPSo...
A simple class to enforce both an elapsed time limit and a deterministic time limit in the same threa...
Definition: time_limit.h:105
bool LimitReached()
Returns true when the external limit is true, or the deterministic time is over the deterministic lim...
Definition: time_limit.h:532
double GetElapsedDeterministicTime() const
Returns the elapsed deterministic time since the construction of this object.
Definition: time_limit.h:260
CoreBasedOptimizer(IntegerVariable objective_var, const std::vector< IntegerVariable > &variables, const std::vector< IntegerValue > &coefficients, std::function< void()> feasible_solution_observer, Model *model)
Literal literal(int i) const
Definition: encoding.h:78
Literal GetOrCreateAssociatedLiteral(IntegerLiteral i_lit)
Definition: integer.cc:202
ABSL_MUST_USE_RESULT bool Enqueue(IntegerLiteral i_lit, absl::Span< const Literal > literal_reason, absl::Span< const IntegerLiteral > integer_reason)
Definition: integer.cc:989
bool IsCurrentlyIgnored(IntegerVariable i) const
Definition: integer.h:625
IntegerValue UpperBound(IntegerVariable i) const
Definition: integer.h:1304
IntegerValue LevelZeroUpperBound(IntegerVariable var) const
Definition: integer.h:1355
IntegerVariable AddIntegerVariable(IntegerValue lower_bound, IntegerValue upper_bound)
Definition: integer.cc:603
IntegerValue LevelZeroLowerBound(IntegerVariable var) const
Definition: integer.h:1350
IntegerValue LowerBound(IntegerVariable i) const
Definition: integer.h:1300
Class that owns everything related to a particular optimization model.
Definition: sat/model.h:38
T Add(std::function< T(Model *)> f)
This makes it possible to have a nicer API on the client side, and it allows both of these forms:
Definition: sat/model.h:81
bool AddLinearConstraint(bool use_lower_bound, Coefficient lower_bound, bool use_upper_bound, Coefficient upper_bound, std::vector< LiteralWithCoeff > *cst)
Definition: sat_solver.cc:299
void SetNumVariables(int num_variables)
Definition: sat_solver.cc:64
bool AddTernaryClause(Literal a, Literal b, Literal c)
Definition: sat_solver.cc:191
const SatParameters & parameters() const
Definition: sat_solver.cc:110
Status ResetAndSolveWithGivenAssumptions(const std::vector< Literal > &assumptions)
Definition: sat_solver.cc:947
const VariablesAssignment & Assignment() const
Definition: sat_solver.h:362
void SetAssumptionLevel(int assumption_level)
Definition: sat_solver.cc:962
int EnqueueDecisionAndBackjumpOnConflict(Literal true_literal)
Definition: sat_solver.cc:499
void SetParameters(const SatParameters &parameters)
Definition: sat_solver.cc:115
bool AddBinaryClause(Literal a, Literal b)
Definition: sat_solver.cc:180
int EnqueueDecisionAndBacktrackOnConflict(Literal true_literal)
Definition: sat_solver.cc:861
void Backtrack(int target_level)
Definition: sat_solver.cc:888
std::vector< Literal > GetLastIncompatibleDecisions()
Definition: sat_solver.cc:1272
bool AddProblemClause(absl::Span< const Literal > literals)
Definition: sat_solver.cc:203
bool AddUnitClause(Literal true_literal)
Definition: sat_solver.cc:164
bool LiteralIsTrue(Literal literal) const
Definition: sat_base.h:150
bool LiteralIsFalse(Literal literal) const
Definition: sat_base.h:147
SatParameters parameters
SharedResponseManager * response
SharedTimeLimit * time_limit
const Constraint * ct
int64 value
IntVar * var
Definition: expr_array.cc:1858
GRBmodel * model
static const int64 kint64max
int64_t int64
A C++ wrapper that provides a simple and unified interface to several linear programming and mixed in...
const int INFO
Definition: log_severity.h:31
const int FATAL
Definition: log_severity.h:32
#define DISALLOW_COPY_AND_ASSIGN(TypeName)
Definition: macros.h:29
absl::Cleanup< absl::decay_t< Callback > > MakeCleanup(Callback &&callback)
Definition: cleanup.h:120
const Collection::value_type::second_type & FindOrDie(const Collection &collection, const typename Collection::value_type::first_type &key)
Definition: map_util.h:176
bool ContainsKey(const Collection &collection, const Key &key)
Definition: map_util.h:170
const Collection::value_type::second_type & FindOrDieNoPrint(const Collection &collection, const typename Collection::value_type::first_type &key)
Definition: map_util.h:186
void STLSortAndRemoveDuplicates(T *v, const LessFunc &less_func)
Definition: stl_util.h:58
std::function< void(Model *)> WeightedSumLowerOrEqual(const std::vector< IntegerVariable > &vars, const VectorInt &coefficients, int64 upper_bound)
Definition: integer_expr.h:299
bool AddObjectiveConstraint(const LinearBooleanProblem &problem, bool use_lower_bound, Coefficient lower_bound, bool use_upper_bound, Coefficient upper_bound, SatSolver *solver)
void RestrictObjectiveDomainWithBinarySearch(IntegerVariable objective_var, const std::function< void()> &feasible_solution_observer, Model *model)
SatSolver::Status ResetAndSolveIntegerProblem(const std::vector< Literal > &assumptions, Model *model)
SatSolver::Status SolveWithCardinalityEncodingAndCore(LogBehavior log, const LinearBooleanProblem &problem, SatSolver *solver, std::vector< bool > *solution)
Coefficient ComputeCoreMinWeight(const std::vector< EncodingNode * > &nodes, const std::vector< Literal > &core)
Definition: encoding.cc:418
EncodingNode * MergeAllNodesWithDeque(Coefficient upper_bound, const std::vector< EncodingNode * > &nodes, SatSolver *solver, std::deque< EncodingNode > *repository)
Definition: encoding.cc:263
int MoveOneUnprocessedLiteralLast(const std::set< LiteralIndex > &processed, int relevant_prefix_size, std::vector< Literal > *literals)
Definition: sat/util.cc:24
std::vector< Literal > ReduceNodesAndExtractAssumptions(Coefficient upper_bound, Coefficient stratified_lower_bound, Coefficient *lower_bound, std::vector< EncodingNode * > *nodes, SatSolver *solver)
Definition: encoding.cc:366
void UseObjectiveForSatAssignmentPreference(const LinearBooleanProblem &problem, SatSolver *solver)
std::function< int64(const Model &)> Value(IntegerVariable v)
Definition: integer.h:1487
std::function< int64(const Model &)> LowerBound(IntegerVariable v)
Definition: integer.h:1467
SatSolver::Status SolveWithLinearScan(LogBehavior log, const LinearBooleanProblem &problem, SatSolver *solver, std::vector< bool > *solution)
double AddOffsetAndScaleObjectiveValue(const LinearBooleanProblem &problem, Coefficient v)
void MinimizeCore(SatSolver *solver, std::vector< Literal > *core)
Definition: sat_solver.cc:2547
SatSolver::Status MinimizeWithHittingSetAndLazyEncoding(const ObjectiveDefinition &objective_definition, const std::function< void()> &feasible_solution_observer, Model *model)
SatSolver::Status SolveIntegerProblem(Model *model)
SatSolver::Status SolveWithWPM1(LogBehavior log, const LinearBooleanProblem &problem, SatSolver *solver, std::vector< bool > *solution)
bool IsAssignmentValid(const LinearBooleanProblem &problem, const std::vector< bool > &assignment)
std::vector< IntegerVariable > NegationOf(const std::vector< IntegerVariable > &vars)
Definition: integer.cc:27
void MinimizeCoreWithPropagation(TimeLimit *limit, SatSolver *solver, std::vector< Literal > *core)
void ProcessCore(const std::vector< Literal > &core, Coefficient min_weight, std::deque< EncodingNode > *repository, std::vector< EncodingNode * > *nodes, SatSolver *solver)
Definition: encoding.cc:445
Coefficient ComputeObjectiveValue(const LinearBooleanProblem &problem, const std::vector< bool > &assignment)
const Coefficient kCoefficientMax(std::numeric_limits< Coefficient::ValueType >::max())
SatSolver::Status SolveWithFuMalik(LogBehavior log, const LinearBooleanProblem &problem, SatSolver *solver, std::vector< bool > *solution)
Coefficient MaxNodeWeightSmallerThan(const std::vector< EncodingNode * > &nodes, Coefficient upper_bound)
Definition: encoding.cc:433
SatSolver::Status SolveWithRandomParameters(LogBehavior log, const LinearBooleanProblem &problem, int num_times, SatSolver *solver, std::vector< bool > *solution)
SatSolver::Status SolveWithCardinalityEncoding(LogBehavior log, const LinearBooleanProblem &problem, SatSolver *solver, std::vector< bool > *solution)
void ExtractAssignment(const LinearBooleanProblem &problem, const SatSolver &solver, std::vector< bool > *assignment)
std::vector< EncodingNode * > CreateInitialEncodingNodes(const std::vector< Literal > &literals, const std::vector< Coefficient > &coeffs, Coefficient *offset, std::deque< EncodingNode > *repository)
Definition: encoding.cc:302
void RandomizeDecisionHeuristic(URBG *random, SatParameters *parameters)
Definition: sat/util.h:100
constexpr IntegerValue kMaxIntegerValue(std::numeric_limits< IntegerValue::ValueType >::max() - 1)
SatSolver::Status MinimizeIntegerVariableWithLinearScanAndLazyEncoding(IntegerVariable objective_var, const std::function< void()> &feasible_solution_observer, Model *model)
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
std::string ProtobufShortDebugString(const P &message)
STL namespace.
int core_index
Definition: optimization.cc:85
Literal literal
Definition: optimization.cc:84
int index
Definition: pack.cc:508
int64 weight
Definition: pack.cc:509
int64 cost
std::vector< double > coefficients
static IntegerLiteral LowerOrEqual(IntegerVariable i, IntegerValue bound)
Definition: integer.h:1270
static IntegerLiteral GreaterOrEqual(IntegerVariable i, IntegerValue bound)
Definition: integer.h:1264
double ScaleIntegerObjective(IntegerValue value) const
std::string message
Definition: trace.cc:395
#define VLOG_IS_ON(verboselevel)
Definition: vlog_is_on.h:41