always-cautious/ortools-clone
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

[CP-SAT] add pandas support for intervals; new sample that uses them; improve the meta-heuristics in feasibility_jump

Browse Source
This commit is contained in:
Laurent Perron
2023-10-09 18:07:11 +02:00
parent 5ace9d7c4a
commit 6d1a1b455c
19 changed files with 897 additions and 579 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
ortools/sat/BUILD.bazel
View File

@@ -205,7 +205,8 @@ cc_library(
":util",
"//ortools/algorithms:binary_search",
"//ortools/util:sorted_interval_list",
"@com_google_absl//absl/log",
"@com_google_absl//absl/functional:bind_front",
"@com_google_absl//absl/log",
"@com_google_absl//absl/random",
"@com_google_absl//absl/strings",
"@com_google_absl//absl/time",

38
ortools/sat/constraint_violation.cc
View File

@@ -1284,11 +1284,13 @@ void AddCircuitFlowConstraints(LinearIncrementalEvaluator& linear_evaluator,
LsEvaluator::LsEvaluator(const CpModelProto& cp_model) : cp_model_(cp_model) {
var_to_constraints_.resize(cp_model_.variables_size());
jump_value_optimal_.resize(cp_model_.variables_size(), true);
num_violated_constraint_per_var_.assign(cp_model_.variables_size(), 0);
std::vector<bool> ignored_constraints(cp_model_.constraints_size(), false);
std::vector<ConstraintProto> additional_constraints;
CompileConstraintsAndObjective(ignored_constraints, additional_constraints);
BuildVarConstraintGraph();
pos_in_violated_constraints_.assign(NumEvaluatorConstraints(), -1);
}
LsEvaluator::LsEvaluator(
@@ -1297,8 +1299,10 @@ LsEvaluator::LsEvaluator(
: cp_model_(cp_model) {
var_to_constraints_.resize(cp_model_.variables_size());
jump_value_optimal_.resize(cp_model_.variables_size(), true);
num_violated_constraint_per_var_.assign(cp_model_.variables_size(), 0);
CompileConstraintsAndObjective(ignored_constraints, additional_constraints);
BuildVarConstraintGraph();
pos_in_violated_constraints_.assign(NumEvaluatorConstraints(), -1);
}
void LsEvaluator::BuildVarConstraintGraph() {
@@ -1497,7 +1501,11 @@ void LsEvaluator::CompileConstraintsAndObjective(
bool LsEvaluator::ReduceObjectiveBounds(int64_t lb, int64_t ub) {
if (!cp_model_.has_objective()) return false;
return linear_evaluator_.ReduceBounds(/*c=*/0, lb, ub);
if (linear_evaluator_.ReduceBounds(/*c=*/0, lb, ub)) {
UpdateViolatedList(0);
return true;
}
return false;
}
void LsEvaluator::OverwriteCurrentSolution(absl::Span<const int64_t> solution) {
@@ -1667,27 +1675,6 @@ double LsEvaluator::WeightedViolationDelta(absl::Span<const double> weights,
return linear_evaluator_.WeightedViolationDelta(weights, var, delta) +
WeightedNonLinearViolationDelta(weights, var, delta);
}
// TODO(user): Speed-up:
// - Use a row oriented representation of the model (could reuse the Apply
// methods on proto).
// - Maintain the list of violated constraints ?
std::vector<int> LsEvaluator::VariablesInViolatedConstraints() const {
std::vector<int> variables;
for (int var = 0; var < cp_model_.variables_size(); ++var) {
if (linear_evaluator_.AppearsInViolatedConstraints(var)) {
variables.push_back(var);
} else {
for (const int ct_index : var_to_constraints_[var]) {
if (constraints_[ct_index]->violation() > 0) {
variables.push_back(var);
break;
}
}
}
}
return variables;
}
bool LsEvaluator::VariableOnlyInLinearConstraintWithConvexViolationChange(
int var) const {
@@ -1695,6 +1682,7 @@ bool LsEvaluator::VariableOnlyInLinearConstraintWithConvexViolationChange(
}
void LsEvaluator::RecomputeViolatedList(bool linear_only) {
num_violated_constraint_per_var_.assign(cp_model_.variables_size(), 0);
violated_constraints_.clear();
pos_in_violated_constraints_.assign(NumEvaluatorConstraints(), -1);
@@ -1713,6 +1701,9 @@ void LsEvaluator::UpdateViolatedList(const int c) {
if (pos != -1) return;
pos_in_violated_constraints_[c] = violated_constraints_.size();
violated_constraints_.push_back(c);
for (const int v : ConstraintToVars(c)) {
num_violated_constraint_per_var_[v] += 1;
}
return;
}
@@ -1723,6 +1714,9 @@ void LsEvaluator::UpdateViolatedList(const int c) {
violated_constraints_[pos] = last;
pos_in_violated_constraints_[c] = -1;
violated_constraints_.pop_back();
for (const int v : ConstraintToVars(c)) {
num_violated_constraint_per_var_[v] -= 1;
}
}
} // namespace sat

11
ortools/sat/constraint_violation.h
View File

@@ -359,19 +359,19 @@ class LsEvaluator {
return &linear_evaluator_;
}
// Returns the set of variables appearing in a violated constraint.
std::vector<int> VariablesInViolatedConstraints() const;
// List of the currently violated constraints.
// - It is initialized by RecomputeViolatedList()
// - And incrementally maintained by UpdateVariableValue()
//
// The order depends on the algorithm used and shouldn't be
// The order depends on the algorithm used and shouldn't be relied on.
void RecomputeViolatedList(bool linear_only);
const std::vector<int>& ViolatedConstraints() const {
return violated_constraints_;
}
// Returns the number of constraints in ViolatedConstraints containing `var`.
int NumViolatedConstraintsForVar(int var) const {
return num_violated_constraint_per_var_[var];
}
// Indicates if the computed jump value is always the best choice.
bool VariableOnlyInLinearConstraintWithConvexViolationChange(int var) const;
@@ -435,6 +435,7 @@ class LsEvaluator {
// - violated_constraints_[pos_in_violated_constraints_[c]] = c
std::vector<int> pos_in_violated_constraints_;
std::vector<int> violated_constraints_;
std::vector<int> num_violated_constraint_per_var_;
// Constraint index and violation delta for the last update.
std::vector<std::pair<int, int64_t>> last_update_violation_changes_;

8
ortools/sat/docs/README.md
View File

@@ -1,13 +1,9 @@
| [home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) |
| ----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | -------------------------------- |
[home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Troubleshooting](troubleshooting.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html)
----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | ------------------------------------- | ---------------------------------------------------------------------------------------
# Using the CP-SAT solver
https://developers.google.com/optimization/cp/cp_solver
## Documentation structure
This document presents modeling recipes for the CP-SAT solver.

8
ortools/sat/docs/boolean_logic.md
View File

@@ -1,13 +1,9 @@
| [home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) |
| ----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | -------------------------------- |
[home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Troubleshooting](troubleshooting.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html)
----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | ------------------------------------- | ---------------------------------------------------------------------------------------
# Boolean logic recipes for the CP-SAT solver.
https://developers.google.com/optimization/
## Introduction
The CP-SAT solver can express Boolean variables and constraints. A **Boolean

7
ortools/sat/docs/channeling.md
View File

@@ -1,12 +1,9 @@
| [home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) |
| ----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | -------------------------------- |
[home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Troubleshooting](troubleshooting.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html)
----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | ------------------------------------- | ---------------------------------------------------------------------------------------
# Channeling constraints
https://developers.google.com/optimization/
A *channeling constraint* links variables inside a model. They're used when you
want to express a complicated relationship between variables, such as "if this
variable satisfies a condition, force another variable to a particular value".

7
ortools/sat/docs/integer_arithmetic.md
View File

@@ -1,12 +1,9 @@
| [home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) |
| ----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | -------------------------------- |
[home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Troubleshooting](troubleshooting.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html)
----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | ------------------------------------- | ---------------------------------------------------------------------------------------
# Integer arithmetic recipes for the CP-SAT solver.
https://developers.google.com/optimization/
## Introduction
The CP-SAT solver can express integer variables and constraints.

7
ortools/sat/docs/model.md
View File

@@ -1,12 +1,9 @@
| [home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) |
| ----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | -------------------------------- |
[home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Troubleshooting](troubleshooting.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html)
----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | ------------------------------------- | ---------------------------------------------------------------------------------------
# Model manipulation
https://developers.google.com/optimization/
## Introduction
In all languages, the CpModel class is a thin wrapper around a

27
ortools/sat/docs/scheduling.md
View File

@@ -1,12 +1,9 @@
| [home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) |
| ----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | -------------------------------- |
[home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Troubleshooting](troubleshooting.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html)
----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | ------------------------------------- | ---------------------------------------------------------------------------------------
# Scheduling recipes for the CP-SAT solver.
https://developers.google.com/optimization/
## Introduction
Scheduling in Operations Research involves problems of tasks, resources and
@@ -754,7 +751,7 @@ def main():
index=tasks_df.index,
starts=starts,
sizes=tasks_df.duration,
performed_literals=performed,
are_present=performed,
)
# Set up the profile. We use fixed (intervals, demands) to fill in the space
@@ -795,7 +792,7 @@ def main():
else:
print(f"task {task} is not performed")
elif status == cp_model.INFEASIBLE:
print("The problem is infeasible")
print("No solution found")
else:
print("Something is wrong, check the status and the log of the solve")
@@ -808,7 +805,7 @@ if __name__ == "__main__":
## Ranking tasks in a disjunctive resource
To rank intervals in a NoOverlap constraint, we will count the number of
To rank intervals in a no_overlap constraint, we will count the number of
performed intervals that precede each interval.
This is slightly complicated if some interval variables are optional. To
@@ -1455,8 +1452,8 @@ public class RankingSampleSat
## Ranking tasks in a disjunctive resource with a circuit constraint.
To rank intervals in a NoOverlap constraint, we will use a circuit constraint to
perform the transitive reduction from precedences to successors.
To rank intervals in a no_overlap constraint, we will use a circuit constraint
to perform the transitive reduction from precedences to successors.
This is slightly complicated if some interval variables are optional, and we
need to take into account the case where no task is performed.
@@ -1728,10 +1725,10 @@ SchedulingWithCalendarSampleSat()
## Detecting if two intervals overlap.
We want a Boolean variable to be true iff two intervals overlap. To enforce
this, we will create 3 Boolean variables, link two of them to the relative
positions of the two intervals, and define the third one using the other two
Boolean variables.
We want a Boolean variable to be true if and only if two intervals overlap. To
enforce this, we will create 3 Boolean variables, link two of them to the
relative positions of the two intervals, and define the third one using the
other two Boolean variables.
There are two ways of linking the three Boolean variables. The first version
uses one clause and two implications. Propagation will be faster using this
@@ -1855,5 +1852,3 @@ OverlappingIntervals()
## Precedences between intervals
## Convex hull of a set of intervals
## Reservoir constraint

7
ortools/sat/docs/solver.md
View File

@@ -1,12 +1,9 @@
| [home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) |
| ----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | -------------------------------- |
[home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Troubleshooting](troubleshooting.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html)
----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | ------------------------------------- | ---------------------------------------------------------------------------------------
# Solving a CP-SAT model
https://developers.google.com/optimization/
## Changing the parameters of the solver
The SatParameters protobuf encapsulates the set of parameters of a CP-SAT

268
ortools/sat/docs/troubleshooting.md Normal file
View File

@@ -0,0 +1,268 @@
# Troubleshooting
## Enable logging
CP-SAT supports extensive logging of its internals.
To enable it, just set the parameter `log_search_progress` to true.
A good explanation of the log can be found in the
[cpsat-primer](https://github.com/d-krupke/cpsat-primer/blob/main/understanding_the_log.md).
## Improving performance with multiple workers
CP-SAT is built with parallelism in mind. While single worker is a good
solution, the solver is tuned to reached its best performance with multiple
workers.
We distinguish multiple thresholds:
- **[8 workers]** This is the minimal setting for parallel search. It blends
workers with different linear relaxations (none, default, maximal), core
based search if applicable, a quick_restart subsolver, a feasibility_jump
first solution subsolver, and dedicated Large Neighborhood Search subsolvers
to find improving solutions.
- **[16 workers]** Bumping to 16 workers adds a continuous probing subsolver,
more first solution subsolvers (random search and feasibility_jump), two
dual subsolvers dedicated to improving the lower bound of the objective
(when minimizing), and more LNS subsolvers.
- **[24 workers]** Adds more dual subsolvers, no_lp and max_lp variations of
previous subsolvers, more feasibility_jump first solution subsolvers, and
more LNS subsolvers.
- **[32 workers and more]** Adds more LNS workers, and more first solution
subsolvers.
## Solve status
Solving a problem yields the following possible status (CpSolverStatus):
- **[UNKNOWN]** The status of the model is still unknown. A search limit has
been reached before any of the statuses below could be determined.
- **[MODEL_INVALID]** The given CpModelProto didn't pass the validation step.
You can get a detailed error by calling `ValidateCpModel(model_proto)`in
C++, or `model.Validate()` in other languages.
- **[FEASIBLE]** A feasible solution has been found. But the search was
stopped before we could prove optimality or before we enumerated all
solutions of a feasibility problem (if asked).
- **[INFEASIBLE]** The problem has been proven infeasible.
- **[OPTIMAL]** An optimal feasible solution has been found. More generally,
this status represent a success. So we also return OPTIMAL if we find a
solution for a pure feasibility problem or if a gap limit has been specified
and we return a solution within this limit. In the case where we need to
return all the feasible solution, this status will only be returned if
solve() enumerated all of them; If we stopped before, we will return
FEASIBLE.
## Debugging infeasible models
Oftentimes, solving a model yields an infeasibility status. While this can be a
valid answer, it can happen because of either a data issue, or a model issue.
Here a simple common sense processes to help diagnose the source of the
infeasibility.
- reduce the model size if this is a parameter
- remove all constraints while keeping the problem infeasible
- play with the domain of the variables to enlarge them
- inject a known feasible solution and try to find where it breaks
- add assumptions to check the soundness of your data (capacity is >= 0, no
crazy ranges of values, ...)
## Using assumptions to explain infeasibility
CP-SAT implements a system of assumptions to find the root cause of an
infeasible problem. To use it, one must instrument constraints with literals and
add these literals to the set of assumptions.
Note that this set is minimized but not guaranteed to be minimal. To find a
minimal unsatisfiable set (MUS), you must minimize the (weighted) sum of these
assumption literals instead of using the assumptions mechanism.
Finally, be aware that solving with assumptions is not compatible with
parallelism. Therefore, the number of workers must be set to 1.
### Python code sample
```python
#!/usr/bin/env python3
"""Code sample that solves a model and gets the infeasibility assumptions."""
from ortools.sat.python import cp_model
def main():
"""Showcases assumptions."""
# Creates the model.
model = cp_model.CpModel()
# Creates the variables.
x = model.NewIntVar(0, 10, "x")
y = model.NewIntVar(0, 10, "y")
z = model.NewIntVar(0, 10, "z")
a = model.NewBoolVar("a")
b = model.NewBoolVar("b")
c = model.NewBoolVar("c")
# Creates the constraints.
model.Add(x > y).OnlyEnforceIf(a)
model.Add(y > z).OnlyEnforceIf(b)
model.Add(z > x).OnlyEnforceIf(c)
# Add assumptions
model.AddAssumptions([a, b, c])
# Creates a solver and solves.
solver = cp_model.CpSolver()
status = solver.Solve(model)
# Print solution.
print(f"Status = {solver.StatusName(status)}")
if status == cp_model.INFEASIBLE:
print(
"SufficientAssumptionsForInfeasibility = "
f"{solver.SufficientAssumptionsForInfeasibility()}"
)
if __name__ == "__main__":
main()
```
### C++ code samples
```cpp
#include <stdlib.h>
#include "absl/types/span.h"
#include "ortools/base/logging.h"
#include "ortools/sat/cp_model.h"
#include "ortools/sat/cp_model.pb.h"
#include "ortools/sat/cp_model_solver.h"
#include "ortools/util/sorted_interval_list.h"
namespace operations_research {
namespace sat {
void AssumptionsSampleSat() {
CpModelBuilder cp_model;
const Domain domain(0, 10);
const IntVar x = cp_model.NewIntVar(domain).WithName("x");
const IntVar y = cp_model.NewIntVar(domain).WithName("y");
const IntVar z = cp_model.NewIntVar(domain).WithName("z");
const BoolVar a = cp_model.NewBoolVar().WithName("a");
const BoolVar b = cp_model.NewBoolVar().WithName("b");
const BoolVar c = cp_model.NewBoolVar().WithName("c");
cp_model.AddGreaterThan(x, y).OnlyEnforceIf(a);
cp_model.AddGreaterThan(y, z).OnlyEnforceIf(b);
cp_model.AddGreaterThan(z, x).OnlyEnforceIf(c);
// Add assumptions
cp_model.AddAssumptions({a, b, c});
// Solving part.
const CpSolverResponse response = Solve(cp_model.Build());
// Print solution.
LOG(INFO) << CpSolverResponseStats(response);
if (response.status() == CpSolverStatus::INFEASIBLE) {
for (const int index :
response.sufficient_assumptions_for_infeasibility()) {
LOG(INFO) << index;
}
}
}
} // namespace sat
} // namespace operations_research
int main(int argc, char** argv) {
operations_research::sat::AssumptionsSampleSat();
return EXIT_SUCCESS;
}
```
### Java code samples
```java
package com.google.ortools.sat.samples;
import com.google.ortools.Loader;
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.Literal;
/** Minimal CP-SAT example to showcase assumptions. */
public class AssumptionsSampleSat {
public static void main(String[] args) {
Loader.loadNativeLibraries();
// Create the model.
CpModel model = new CpModel();
// Create the variables.
IntVar x = model.newIntVar(0, 10, "x");
IntVar y = model.newIntVar(0, 10, "y");
IntVar z = model.newIntVar(0, 10, "z");
Literal a = model.newBoolVar("a");
Literal b = model.newBoolVar("b");
Literal c = model.newBoolVar("c");
// Creates the constraints.
model.addGreaterThan(x, y).onlyEnforceIf(a);
model.addGreaterThan(y, z).onlyEnforceIf(b);
model.addGreaterThan(z, x).onlyEnforceIf(c);
// Add assumptions
model.addAssumptions(new Literal[] {a, b, c});
// Create a solver and solve the model.
CpSolver solver = new CpSolver();
CpSolverStatus status = solver.solve(model);
// Print solution.
// Check that the problem is infeasible.
if (status == CpSolverStatus.INFEASIBLE) {
System.out.println(solver.sufficientAssumptionsForInfeasibility());
}
}
private AssumptionsSampleSat() {}
}
```
### C\# code samples
```cs
using System;
using Google.OrTools.Sat;
public class AssumptionsSampleSat
{
static void Main()
{
// Creates the model.
CpModel model = new CpModel();
// Creates the variables.
IntVar x = model.NewIntVar(0, 10, "x");
IntVar y = model.NewIntVar(0, 10, "y");
IntVar z = model.NewIntVar(0, 10, "z");
ILiteral a = model.NewBoolVar("a");
ILiteral b = model.NewBoolVar("b");
ILiteral c = model.NewBoolVar("c");
// Creates the constraints.
model.Add(x > y).OnlyEnforceIf(a);
model.Add(y > z).OnlyEnforceIf(b);
model.Add(z > x).OnlyEnforceIf(c);
// Add assumptions
model.AddAssumptions(new ILiteral[] { a, b, c });
// Creates a solver and solves the model.
CpSolver solver = new CpSolver();
CpSolverStatus status = solver.Solve(model);
Console.WriteLine(solver.SufficientAssumptionsForInfeasibility());
}
}
```

793
ortools/sat/feasibility_jump.cc
View File

File diff suppressed because it is too large Load Diff

173
ortools/sat/feasibility_jump.h
View File

@@ -20,9 +20,11 @@
#include <limits>
#include <memory>
#include <string>
#include <utility>
#include <vector>
#include "absl/time/time.h"
#include "absl/functional/any_invocable.h"
#include "absl/functional/bind_front.h"
#include "absl/types/span.h"
#include "ortools/sat/constraint_violation.h"
#include "ortools/sat/linear_model.h"
@@ -36,6 +38,55 @@ namespace operations_research::sat {
class CompoundMoveBuilder;
// This class lazily caches the results of `compute_jump(var)` which returns a
// <delta, score> pair.
// Variables' scores can be manually modified using MutableScores (if the
// optimal jump is known not to change), or marked for recomputation on the next
// call to GetJump(var) by calling Recompute.
class JumpTable {
public:
explicit JumpTable(
absl::AnyInvocable<std::pair<int64_t, double>(int)> compute_jump);
void RecomputeAll(int num_variables);
// Gets the current jump delta and score, recomputing if necessary.
std::pair<int64_t, double> GetJump(int var);
// If the new optimum value and score is known, users can update it directly.
// e.g. after weight rescaling, or after changing a binary variable.
void SetJump(int var, int64_t delta, double score);
// Recompute the jump for `var` when `GetJump(var)` is next called.
void Recompute(int var);
// Returns true if the jump score for `var` might be negative.
bool PossiblyGood(int var) const;
// Advanced usage, allows users to read possibly stale deltas for incremental
// score updates.
absl::Span<const int64_t> Deltas() const {
return absl::MakeConstSpan(deltas_);
}
absl::Span<const double> Scores() const {
return absl::MakeConstSpan(scores_);
}
absl::Span<double> MutableScores() { return absl::MakeSpan(scores_); }
// Note if you have very high weights (e.g. when using decay), the tolerances
// in this function are likely too tight.
bool JumpIsUpToDate(int var); // For debugging and testing.
private:
absl::AnyInvocable<std::pair<int64_t, double>(int)> compute_jump_;
// For each variable, we store:
// - A jump delta which represents a change in variable value:
// new_value = old + delta, which is non-zero except for fixed variables.
// - The associated weighted feasibility violation change if we take this
// jump.
std::vector<int64_t> deltas_;
std::vector<double> scores_;
std::vector<bool> needs_recomputation_;
};
// Implements and heuristic similar to the one described in the paper:
// "Feasibility Jump: an LP-free Lagrangian MIP heuristic", Bjørnar
// Luteberget, Giorgio Sartor, 2023, Mathematical Programming Computation.
@@ -62,7 +113,11 @@ class FeasibilityJumpSolver : public SubSolver {
shared_response_(shared_response),
shared_bounds_(shared_bounds),
shared_stats_(shared_stats),
random_(params_) {}
random_(params_),
linear_jumps_(
absl::bind_front(&FeasibilityJumpSolver::ComputeLinearJump, this)),
general_jumps_(absl::bind_front(
&FeasibilityJumpSolver::ComputeGeneralJump, this)) {}
// If VLOG_IS_ON(1), it will export a bunch of statistics.
~FeasibilityJumpSolver() override;
@@ -103,51 +158,60 @@ class FeasibilityJumpSolver : public SubSolver {
void PerturbateCurrentSolution();
std::string OneLineStats() const;
// Linear only.
bool JumpIsUpToDate(int var); // Debug.
void RecomputeJump(int var);
void RecomputeAllJumps();
void MarkJumpsThatNeedsToBeRecomputed(int changed_var);
// Returns the weighted violation delta plus epsilon * the objective delta.
double ComputeScore(absl::Span<const double> scan_weights, int var,
int64_t delta, bool linear_only) const;
// As above, but uses appropriate scan weights based on the current value of
// `use_compound_moves_`.
double ComputeScore(int var, int64_t delta, bool linear_only) const {
return ComputeScore(use_compound_moves_ ? compound_weights_ : weights_, var,
delta, linear_only);
}
// Computes the optimal value for variable v, considering only the violation
// of linear constraints.
std::pair<int64_t, double> ComputeLinearJump(int var);
// Computes the optimal value for variable v, considering all constraints
// (assuming violation functions are convex).
std::pair<int64_t, double> ComputeGeneralJump(int var);
// Marks all variables whose jump value may have changed due to the last
// update, except for `changed var`.
void MarkJumpsThatNeedToBeRecomputed(int changed_var, JumpTable& jumps);
// Moves.
bool DoSomeLinearIterations();
bool DoSomeGeneralIterations();
// Returns true if an improving move was found.
bool ScanRelevantVariables(int* improving_var, int64_t* improving_value,
double* improving_score, bool* time_limit_crossed);
bool ScanRelevantVariables(int num_to_scan, JumpTable& jumps, int* var,
int64_t* value, double* score);
// Increases the weight of the currently infeasible constraints.
void UpdateViolatedConstraintWeights();
// Ensures jumps remains consistent.
void UpdateViolatedConstraintWeights(JumpTable& jumps);
// Returns true if changing this variable to its jump value reduces weighted
// violation, or has no impact on weighted violation but reduces the
// objective.
bool IsGood(int var) const;
void UpdateNumViolatedConstraintsPerVar();
void RecomputeVarsToScan();
void RecomputeVarsToScan(JumpTable&);
// Returns true if it is possible that `var` may have value that reduces
// weighted violation or improve the objective.
// Note that this is independent of the actual weights used.
bool ShouldScan(int var) const;
void AddVarToScan(int var);
bool MustRecomputeJumpOnGeneralUpdate(int var) const;
bool ShouldScan(const JumpTable& jumps, int var) const;
void AddVarToScan(const JumpTable&, int var);
// Resets the weights used to find compound moves. Only modifies weights that
// have changed, but the invariant below holds for all constraints:
// Resets the weights used to find compound moves.
// Ensures the following invariant holds afterwards:
// compound_weights[c] = weights_[c] if c is violated, and epsilon *
// weights_[c] otherwise.
void ResetChangedCompoundWeights();
// Returns true if we should accept a single variable move.
bool ShouldAcceptUnitMove(int var, int64_t delta, double score);
// Returns true if we should push this change to move_.
// `novelty` is the total discount applied to the score due to using
// `cumulative_weights_`, should always be positive (modulo floating-point
// errors).
bool ShouldExtendCompoundMove(int var, int64_t delta, double score,
double novelty);
bool ShouldExtendCompoundMove(double score, double novelty);
// Validates each element in num_violated_constraints_per_var_ against
// evaluator_->ViolatedConstraints.
@@ -168,19 +232,15 @@ class FeasibilityJumpSolver : public SubSolver {
bool is_initialized_ = false;
std::atomic<bool> model_is_supported_ = true;
std::atomic<bool> task_generated_ = false;
bool time_limit_crossed_ = false;
std::unique_ptr<LsEvaluator> evaluator_;
std::vector<Domain> var_domains_;
std::vector<bool> var_has_two_values_;
std::vector<bool> var_occurs_in_non_linear_constraint_;
// For each variable, we store:
// - A jump delta which represent a change in variable value:
// new_value = old + delta, which is non-zero except for fixed variable.
// - The associated weighted feasibility violation change if we take this
// jump.
std::vector<bool> jump_need_recomputation_;
std::vector<int64_t> jump_deltas_;
std::vector<double> jump_scores_;
JumpTable linear_jumps_;
JumpTable general_jumps_;
std::vector<double> for_weight_update_;
// The score of a solution is just the sum of infeasibility of each
@@ -199,11 +259,9 @@ class FeasibilityJumpSolver : public SubSolver {
// weight like for SAT activities.
double bump_value_ = 1.0;
// Sparse list of jumps with a potential delta < 0.0.
// If jump_need_recomputation_[var] is true, we lazily recompute the exact
// delta as we randomly pick variables from here.
std::vector<bool> in_good_jumps_;
std::vector<int> good_jumps_;
// A list of variables that might be relevant to check for improving jumps.
std::vector<bool> in_vars_to_scan_;
std::vector<int> vars_to_scan_;
// We restart each time our local deterministic time crosses this.
double dtime_restart_threshold_ = 0.0;
@@ -219,9 +277,14 @@ class FeasibilityJumpSolver : public SubSolver {
// Each time we reset the weights, randomly decide if we will use compound
// moves or not.
bool use_compound_moves_ = false;
// We randomize the beam width each batch, choosing either 1 or 2, biased
// towards 1.
int compound_move_beam_search_width_ = 1;
// Limit the discrepancy in compound move search (i.e. limit the number of
// backtracks to any ancestor of the current leaf). This is set to 0 whenever
// a new incumbent is found or weights are updated, and increased at fixed
// point.
// Weights are only increased if no moves are found with discrepancy 2.
// Empirically we have seen very few moves applied with discrepancy > 2.
int compound_move_max_discrepancy_ = 0;
// Statistics
int64_t num_batches_ = 0;
@@ -235,14 +298,9 @@ class FeasibilityJumpSolver : public SubSolver {
int64_t num_solutions_imported_ = 0;
int64_t num_weight_updates_ = 0;
// A list of variables that might be relevant to check in general iterations.
std::vector<bool> in_vars_to_scan_;
std::vector<int> vars_to_scan_;
std::unique_ptr<CompoundMoveBuilder> move_;
// Counts the number of violated constraints each var is in.
// Only maintained in general iterations as jump_score_ is used for filtering
// moves in linear iterations.
std::vector<int> num_violated_constraints_per_var_;
// Info on the last solution loaded.
@@ -281,29 +339,15 @@ class CompoundMoveBuilder {
return stack_.empty() ? 0.0 : stack_.back().cumulative_score;
}
// Returns the sum of objective deltas from the root.
double ObjectiveDelta() const {
return stack_.empty() ? 0.0 : stack_.back().cumulative_objective_delta;
}
// Returns the minimum score of any child of the current node in the stack,
// with a maximum of 0.
// NB: We assume all improving moves will immediately be committed, so this
// always returns 0 when the stack is empty.
double BestChildScore() const {
return stack_.empty() ? 0.0 : stack_.back().best_child_score;
}
// Returns the number of times a child has been pushed to the current node.
// NB: Always returns 0 when the stack is empty.
int NumChildrenExplored() const {
return stack_.empty() ? 0 : stack_.back().num_children;
// Returns the number of backtracks to any ancestor of the current leaf.
int Discrepancy() const {
return stack_.empty() ? 0 : stack_.back().discrepancy;
}
// Returns true if this move reduces weighted violation, or improves the
// objective without increasing violation.
bool IsImproving() const;
// Returns the number of backtracking moves that have been applied.
int NumBacktracks() const { return num_backtracks_; }
@@ -316,13 +360,12 @@ class CompoundMoveBuilder {
// Instead of tracking this on the compound move, we store the partial sums
// to avoid numerical issues causing negative scores after backtracking.
double cumulative_score;
double cumulative_objective_delta;
// Used to avoid infinite loops, this tracks the best score of any immediate
// children (and not deeper descendants) to avoid re-exploring the same
// prefix.
double best_child_score = 0.0;
int num_children = 0;
int discrepancy = 0;
};
LsEvaluator* evaluator_;
std::vector<bool> var_on_stack_;

43
ortools/sat/python/cp_model.py
View File

@@ -2146,10 +2146,11 @@ class CpModel:
"""Creates an optional interval var from start, size, end, and is_present.
An optional interval variable is a constraint, that is itself used in other
constraints like NoOverlap. This constraint is protected by an is_present
constraints like NoOverlap. This constraint is protected by a presence
literal that indicates if it is active or not.
Internally, it ensures that `is_present` implies `start + size == end`.
Internally, it ensures that `is_present` implies `start + size ==
end`.
Args:
start: The start of the interval. It can be an integer value, or an
@@ -2193,7 +2194,7 @@ class CpModel:
starts: Union[LinearExprT, pd.Series],
sizes: Union[LinearExprT, pd.Series],
ends: Union[LinearExprT, pd.Series],
performed_literals: Union[LiteralT, pd.Series],
are_present: Union[LiteralT, pd.Series],
) -> pd.Series:
"""Creates a series of interval variables with the given name.
@@ -2209,6 +2210,9 @@ class CpModel:
ends (Union[LinearExprT, pd.Series]): The ends of each interval in the
set. If a `pd.Series` is passed in, it will be based on the
corresponding values of the pd.Series.
are_present (Union[LiteralT, pd.Series]): The performed literal of each
interval in the set. If a `pd.Series` is passed in, it will be based on
the corresponding values of the pd.Series.
Returns:
pd.Series: The interval variable set indexed by its corresponding
@@ -2227,9 +2231,7 @@ class CpModel:
starts = _ConvertToLinearExprSeriesAndValidateIndex(starts, index)
sizes = _ConvertToLinearExprSeriesAndValidateIndex(sizes, index)
ends = _ConvertToLinearExprSeriesAndValidateIndex(ends, index)
performed_literals = _ConvertToLiteralSeriesAndValidateIndex(
performed_literals, index
)
are_present = _ConvertToLiteralSeriesAndValidateIndex(are_present, index)
interval_array = []
for i in index:
@@ -2238,14 +2240,18 @@ class CpModel:
start=starts[i],
size=sizes[i],
end=ends[i],
is_present=performed_literals[i],
is_present=are_present[i],
name=f"{name}[{i}]",
)
)
return pd.Series(index=index, data=interval_array)
def NewOptionalFixedSizeIntervalVar(
self, start: LinearExprT, size: IntegralT, is_present: LiteralT, name: str
self,
start: LinearExprT,
size: IntegralT,
is_present: LiteralT,
name: str,
) -> IntervalVar:
"""Creates an interval variable from start, and a fixed size.
@@ -2273,7 +2279,12 @@ class CpModel:
)
is_present_index = self.GetOrMakeBooleanIndex(is_present)
return IntervalVar(
self.__model, start_expr, size_expr, end_expr, is_present_index, name
self.__model,
start_expr,
size_expr,
end_expr,
is_present_index,
name,
)
def NewOptionalFixedSizeIntervalVarSeries(
@@ -2282,7 +2293,7 @@ class CpModel:
index: pd.Index,
starts: Union[LinearExprT, pd.Series],
sizes: Union[IntegralT, pd.Series],
performed_literals: Union[LiteralT, pd.Series],
are_present: Union[LiteralT, pd.Series],
) -> pd.Series:
"""Creates a series of interval variables with the given name.
@@ -2295,9 +2306,9 @@ class CpModel:
sizes (Union[IntegralT, pd.Series]): The fixed size of each interval in
the set. If a `pd.Series` is passed in, it will be based on the
corresponding values of the pd.Series.
performed_literals (Union[LiteralT, pd.Series]): The performed literal of
each interval in the set. If a `pd.Series` is passed in, it will be
based on the corresponding values of the pd.Series.
are_present (Union[LiteralT, pd.Series]): The performed literal of each
interval in the set. If a `pd.Series` is passed in, it will be based on
the corresponding values of the pd.Series.
Returns:
pd.Series: The interval variable set indexed by its corresponding
@@ -2315,16 +2326,14 @@ class CpModel:
starts = _ConvertToLinearExprSeriesAndValidateIndex(starts, index)
sizes = _ConvertToIntegralSeriesAndValidateIndex(sizes, index)
performed_literals = _ConvertToLiteralSeriesAndValidateIndex(
performed_literals, index
)
are_present = _ConvertToLiteralSeriesAndValidateIndex(are_present, index)
interval_array = []
for i in index:
interval_array.append(
self.NewOptionalFixedSizeIntervalVar(
start=starts[i],
size=sizes[i],
is_present=performed_literals[i],
is_present=are_present[i],
name=f"{name}[{i}]",
)
)

58
ortools/sat/python/cp_model_test.py
View File

@@ -1444,6 +1444,64 @@ class CpModelTest(absltest.TestCase):
solution = solver.BooleanValues(x)
self.assertTrue((solution.values == [False, True, False]).all())
def testFixedSizeIntervalVarSeries(self):
print("testFixedSizeIntervalVarSeries")
df = pd.DataFrame([2, 4, 6], columns=["size"])
model = cp_model.CpModel()
starts = model.NewIntVarSeries(
name="starts", index=df.index, lower_bounds=0, upper_bounds=5
)
presences = model.NewBoolVarSeries(name="rresences", index=df.index)
fixed_size_intervals = model.NewFixedSizeIntervalVarSeries(
name="fixed_size_intervals",
index=df.index,
starts=starts,
sizes=df.size,
)
opt_fixed_size_intervals = model.NewOptionalFixedSizeIntervalVarSeries(
name="fixed_size_intervals",
index=df.index,
starts=starts,
sizes=df.size,
are_present=presences,
)
model.AddNoOverlap(
fixed_size_intervals.to_list() + opt_fixed_size_intervals.to_list()
)
self.assertLen(model.Proto().constraints, 7)
def testIntervalVarSeries(self):
print("testIntervalVarSeries")
df = pd.DataFrame([2, 4, 6], columns=["size"])
model = cp_model.CpModel()
starts = model.NewIntVarSeries(
name="starts", index=df.index, lower_bounds=0, upper_bounds=5
)
sizes = model.NewIntVarSeries(
name="sizes", index=df.index, lower_bounds=2, upper_bounds=4
)
ends = model.NewIntVarSeries(
name="ends", index=df.index, lower_bounds=0, upper_bounds=10
)
presences = model.NewBoolVarSeries(name="rresences", index=df.index)
intervals = model.NewIntervalVarSeries(
name="fixed_size_intervals",
index=df.index,
starts=starts,
sizes=sizes,
ends=ends,
)
opt_intervals = model.NewOptionalIntervalVarSeries(
name="fixed_size_intervals",
index=df.index,
starts=starts,
sizes=sizes,
ends=ends,
are_present=presences,
)
model.AddNoOverlap(intervals.to_list() + opt_intervals.to_list())
self.assertLen(model.Proto().constraints, 13)
if __name__ == "__main__":
absltest.main()

4
ortools/sat/samples/assignment_sat.py
View File

@@ -92,8 +92,10 @@ def main():
selected = data.loc[solver.BooleanValues(x).loc[lambda x: x].index]
for unused_index, row in selected.iterrows():
print(f"{row.task} assigned to {row.worker} with a cost of {row.cost}")
elif status == cp_model.INFEASIBLE:
print("No solution found")
else:
print("No solution found.")
print("Something is wrong, check the status and the log of the solve")
# [END print_solution]

8
ortools/sat/samples/bin_packing_sat.py
View File

@@ -105,13 +105,13 @@ def main():
# [END objective]
# [START solve]
# Create the solver with the CP-SAT backend, and solve the model.
# Create the solver and solve the model.
solver = cp_model.CpSolver()
status = solver.Solve(model)
# [END solve]
# [START print_solution]
if status == cp_model.OPTIMAL:
if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
print(f"Number of bins used = {solver.ObjectiveValue()}")
x_values = solver.BooleanValues(x)
@@ -129,8 +129,10 @@ def main():
print(f"Total packed weight: {items.weight.sum()}")
print()
print(f"Time = {solver.WallTime()} seconds")
elif status == cp_model.INFEASIBLE:
print("No solution found")
else:
print("The problem does not have an optimal solution.")
print("Something is wrong, check the status and the log of the solve")
# [END print_solution]

4
ortools/sat/samples/cumulative_variable_profile_sample_sat.py
View File

@@ -117,7 +117,7 @@ def main():
index=tasks_df.index,
starts=starts,
sizes=tasks_df.duration,
performed_literals=performed,
are_present=performed,
)
# [END variables]
@@ -166,7 +166,7 @@ def main():
else:
print(f"task {task} is not performed")
elif status == cp_model.INFEASIBLE:
print("The problem is infeasible")
print("No solution found")
else:
print("Something is wrong, check the status and the log of the solve")
# [END print_solution]

2
ortools/sat/synchronization.h
View File

@@ -400,7 +400,7 @@ class SharedResponseManager {
CpSolverStatus synchronized_best_status_ ABSL_GUARDED_BY(mutex_) =
CpSolverStatus::UNKNOWN;
std::vector<int> unsat_cores_ ABSL_GUARDED_BY(mutex_);
SharedSolutionRepository<int64_t> solutions_ ABSL_GUARDED_BY(mutex_);
SharedSolutionRepository<int64_t> solutions_; // Thread-safe.
int num_solutions_ ABSL_GUARDED_BY(mutex_) = 0;
int64_t inner_objective_lower_bound_ ABSL_GUARDED_BY(mutex_) =