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

[CP-SAT] Change semantics/API of the reservoir constraint; added presolve on the reservoir constraint; experimental propagator for the reservoir constraint; fix #2293

Browse Source
This commit is contained in:
Laurent Perron
2020-12-24 06:41:44 +01:00
parent edb6c0b630
commit c80eec373b
18 changed files with 743 additions and 145 deletions
Download Patch File Download Diff File Expand all files Collapse all files

70
ortools/java/com/google/ortools/sat/CpModel.java
View File

@@ -499,31 +499,37 @@ public final class CpModel {
/**
* Adds {@code Reservoir(times, demands, minLevel, maxLevel)}.
*
* <p>Maintains a reservoir level within bounds. The water level starts at 0, and at any non
* negative time , it must be between minLevel and maxLevel. Furthermore, this constraints expect
* all times variables to be non negative. If the variable times[i] is assigned a value t, then
* the current level changes by demands[i] (which is constant) at the time t.
* <p>Maintains a reservoir level within bounds. The water level starts at 0, and at any times, it
* must be between minLevel and maxLevel. If the variable {@code times[i]} is assigned a value t,
* and if {@code actives[i]} is true, then the current level changes by {@code demands[i]} (which
* is constant) at the time t.
*
* <p>Note that {@code minLevel} can be greater than 0, or {@code maxLevel} can be less than 0. It
* just forces some demands to be executed at time 0 to make sure that we are within those bounds
* with the executed demands. Therefore, {@code forall t >= 0: minLevel <= sum(demands[i] if
* times[i] <= t) <= maxLevel}.
* <p>Note that {@code minLevel} must be less than 0, and {@code maxLevel} must be greater than 0.
* Therefore, {@code forall t : minLevel <= sum(demands[i] if times[i] <= t) <= maxLevel}.
*
* @param times a list of positive integer variables which specify the time of the filling or
* emptying the reservoir
* @param times a list of integer variables which specify the time of the filling or emptying the
* reservoir
* @param demands a list of integer values that specifies the amount of the emptying or feeling
* @param minLevel at any non negative time, the level of the reservoir must be greater of equal
* than the min level
* @param maxLevel at any non negative time, the level of the reservoir must be less or equal than
* the max level
* @param minLevel at any time, the level of the reservoir must be greater of equal than the min
* level. minLevel must me <= 0.
* @param maxLevel at any time, the level of the reservoir must be less or equal than the max
* level. maxLevel must be >= 0.
* @return an instance of the Constraint class
* @throws MismatchedArrayLengths if times and demands have different length
* @throws IllegalArgumentException if minLevel > 0
* @throws IllegalArgumentException if maxLevel < 0
*/
public Constraint addReservoirConstraint(
IntVar[] times, long[] demands, long minLevel, long maxLevel) throws MismatchedArrayLengths {
if (times.length != demands.length) {
throw new MismatchedArrayLengths("addReservoirConstraint", "times", "demands");
}
if (minLevel > 0) {
throw new IllegalArgumentException("addReservoirConstraint: minLevel must be <= 0");
}
if (maxLevel < 0) {
throw new IllegalArgumentException("addReservoirConstraint: maxLevel must be >= 0");
}
Constraint ct = new Constraint(modelBuilder);
ReservoirConstraintProto.Builder reservoir = ct.getBuilder().getReservoirBuilder();
for (IntVar var : times) {
@@ -549,27 +555,25 @@ public final class CpModel {
/**
* Adds {@code Reservoir(times, demands, actives, minLevel, maxLevel)}.
*
* <p>Maintains a reservoir level within bounds. The water level starts at 0, and at any non
* negative time , it must be between minLevel and maxLevel. Furthermore, this constraints expect
* all times variables to be non negative. If actives[i] is true, and if times[i] is assigned a
* value t, then the current level changes by demands[i] (which is constant) at the time t.
* <p>Maintains a reservoir level within bounds. The water level starts at 0, and at any time, it
* must be between minLevel and maxLevel. If the variable {@code times[i]} is assigned a value t,
* then the current level changes by {@code demands[i]} (which is constant) at the time t.
*
* <p>Note that {@code minLevel} can be greater than 0, or {@code maxLevel} can be less than 0. It
* just forces some demands to be executed at time 0 to make sure that we are within those bounds
* with the executed demands. Therefore, {@code forall t >= 0: minLevel <= sum(demands[i] *
* actives[i] if times[i] <= t) <= maxLevel}.
* <p>Note that {@code minLevel} must be less than 0, and {@code maxLevel} must be greater than 0.
* Therefore, {@code forall t : minLevel <= sum(demands[i] * actives[i] if times[i] <= t) <=
* maxLevel}.
*
* @param times a list of positive integer variables which specify the time of the filling or
* emptying the reservoir
* @param times a list of integer variables which specify the time of the filling or emptying the
* reservoir
* @param demands a list of integer values that specifies the amount of the emptying or feeling
* @param actives a list of integer variables that specifies if the operation actually takes
* place.
* @param minLevel at any non negative time, the level of the reservoir must be greater of equal
* than the min level
* @param maxLevel at any non negative time, the level of the reservoir must be less or equal than
* the max level
* @param minLevel at any time, the level of the reservoir must be greater of equal than the min
* level. minLevel must me <= 0.
* @param maxLevel at any time, the level of the reservoir must be less or equal than the max
* level. maxLevel must be >= 0.
* @return an instance of the Constraint class
* @throws MismatchedArrayLengths if times, demands, or actives have different length
* @throws IllegalArgumentException if minLevel > 0
* @throws IllegalArgumentException if maxLevel < 0
*/
public Constraint addReservoirConstraintWithActive(IntVar[] times, long[] demands,
IntVar[] actives, long minLevel, long maxLevel) throws MismatchedArrayLengths {
@@ -579,6 +583,12 @@ public final class CpModel {
if (times.length != actives.length) {
throw new MismatchedArrayLengths("addReservoirConstraint", "times", "actives");
}
if (minLevel > 0) {
throw new IllegalArgumentException("addReservoirConstraint: minLevel must be <= 0");
}
if (maxLevel < 0) {
throw new IllegalArgumentException("addReservoirConstraint: maxLevel must be >= 0");
}
Constraint ct = new Constraint(modelBuilder);
ReservoirConstraintProto.Builder reservoir = ct.getBuilder().getReservoirBuilder();

14
ortools/sat/cp_model.proto
View File

@@ -171,16 +171,16 @@ message CumulativeConstraintProto {
}
// Maintain a reservoir level within bounds. The water level starts at 0, and at
// any time >= 0, it must be within min_level, and max_level. Furthermore, this
// constraints expect all times variables to be >= 0.
// any time, it must be within [min_level, max_level].
//
// If the variable actives[i] is true, and if the variable times[i] is assigned
// a value t, then the current level changes by demands[i] (which is constant)
// at the time t.
//
// Note that level min can be > 0, or level max can be < 0. It just forces
// some demands to be executed at time 0 to make sure that we are within those
// bounds with the executed demands. Therefore, at any time t >= 0:
// at the time t. Therefore, at any time t:
// sum(demands[i] * actives[i] if times[i] <= t) in [min_level, max_level]
//
// Note that min level must be <= 0, and the max level must be >= 0. Please use
// fixed demands to simulate initial state.
//
// The array of boolean variables 'actives', if defined, indicates which actions
// are actually performed. If this array is not defined, then it is assumed that
// all actions will be performed.

20
ortools/sat/cp_model_checker.cc
View File

@@ -312,15 +312,19 @@ std::string ValidateReservoirConstraint(const CpModelProto& model,
return absl::StrCat("Times and demands fields must be of the same size: ",
ProtobufShortDebugString(ct));
}
for (const int t : ct.reservoir().times()) {
const IntegerVariableProto& time = model.variables(t);
for (const int64 bound : time.domain()) {
if (bound < 0) {
return absl::StrCat("Time variables must be >= 0 in constraint ",
ProtobufShortDebugString(ct));
}
}
if (ct.reservoir().min_level() > 0) {
return absl::StrCat(
"The min level of a reservoir must be <= 0. Please use fixed events to "
"setup initial state: ",
ProtobufShortDebugString(ct));
}
if (ct.reservoir().max_level() < 0) {
return absl::StrCat(
"The max level of a reservoir must be >= 0. Please use fixed events to "
"setup initial state: ",
ProtobufShortDebugString(ct));
}
int64 sum_abs = 0;
for (const int64 demand : ct.reservoir().demands()) {
sum_abs = CapAdd(sum_abs, std::abs(demand));

49
ortools/sat/cp_model_expand.cc
View File

@@ -126,49 +126,6 @@ void ExpandReservoir(ConstraintProto* ct, PresolveContext* context) {
sum->add_domain(reservoir.max_level());
}
// Constrains the reservoir level to be consistent at time 0.
// We need to do it only if 0 is not in [min_level..max_level].
// Otherwise, the regular propagation will already check it.
if (reservoir.min_level() > 0 || reservoir.max_level() < 0) {
auto* const sum_at_zero =
context->working_model->add_constraints()->mutable_linear();
for (int i = 0; i < num_events; ++i) {
const int active_i = is_active_literal(i);
if (context->LiteralIsFalse(active_i)) continue;
const int time_i = reservoir.times(i);
const int lesseq_0 = context->NewBoolVar();
// lesseq_0 => (time_i <= 0) && active_i is true
ConstraintProto* const lesseq = context->working_model->add_constraints();
lesseq->add_enforcement_literal(lesseq_0);
lesseq->mutable_linear()->add_vars(time_i);
lesseq->mutable_linear()->add_coeffs(1);
lesseq->mutable_linear()->add_domain(kint64min);
lesseq->mutable_linear()->add_domain(0);
if (!context->LiteralIsTrue(active_i)) {
context->AddImplication(lesseq_0, active_i);
}
// Not(lesseq_0) && active_i => (time_i >= 1)
ConstraintProto* const greater =
context->working_model->add_constraints();
greater->add_enforcement_literal(NegatedRef(lesseq_0));
greater->add_enforcement_literal(active_i);
greater->mutable_linear()->add_vars(time_i);
greater->mutable_linear()->add_coeffs(1);
greater->mutable_linear()->add_domain(1);
greater->mutable_linear()->add_domain(kint64max);
// Accumulate in the sum_at_zero constraint.
sum_at_zero->add_vars(lesseq_0);
sum_at_zero->add_coeffs(reservoir.demands(i));
}
sum_at_zero->add_domain(reservoir.min_level());
sum_at_zero->add_domain(reservoir.max_level());
}
ct->Clear();
context->UpdateRuleStats("reservoir: expanded");
}
@@ -1441,6 +1398,8 @@ void ExpandAllDiff(bool expand_non_permutations, ConstraintProto* ct,
} // namespace
void ExpandCpModel(PresolveContext* context) {
if (context->params().disable_constraint_expansion()) return;
if (context->ModelIsUnsat()) return;
// Make sure all domains are initialized.
@@ -1455,7 +1414,9 @@ void ExpandCpModel(PresolveContext* context) {
bool skip = false;
switch (ct->constraint_case()) {
case ConstraintProto::ConstraintCase::kReservoir:
ExpandReservoir(ct, context);
if (context->params().expand_reservoir_constraints()) {
ExpandReservoir(ct, context);
}
break;
case ConstraintProto::ConstraintCase::kIntMod:
ExpandIntMod(ct, context);

24
ortools/sat/cp_model_loader.cc
View File

@@ -46,6 +46,7 @@
#include "ortools/sat/sat_parameters.pb.h"
#include "ortools/sat/sat_solver.h"
#include "ortools/sat/table.h"
#include "ortools/sat/timetable.h"
#include "ortools/util/saturated_arithmetic.h"
#include "ortools/util/sorted_interval_list.h"
@@ -1345,6 +1346,26 @@ void LoadCumulativeConstraint(const ConstraintProto& ct, Model* m) {
m->Add(Cumulative(intervals, demands, capacity));
}
void LoadReservoirConstraint(const ConstraintProto& ct, Model* m) {
auto* mapping = m->GetOrCreate<CpModelMapping>();
auto* encoder = m->GetOrCreate<IntegerEncoder>();
std::vector<AffineExpression> times;
std::vector<IntegerValue> deltas;
std::vector<Literal> presences;
const int size = ct.reservoir().times().size();
for (int i = 0; i < size; ++i) {
times.push_back(mapping->Integer(ct.reservoir().times(i)));
deltas.push_back(IntegerValue(ct.reservoir().demands(i)));
if (!ct.reservoir().actives().empty()) {
presences.push_back(mapping->Literal(ct.reservoir().actives(i)));
} else {
presences.push_back(encoder->GetTrueLiteral());
}
}
AddReservoirConstraint(times, deltas, presences, ct.reservoir().min_level(),
ct.reservoir().max_level(), m);
}
// If a variable is constant and its value appear in no other variable domains,
// then the literal encoding the index and the one encoding the target at this
// value are equivalent.
@@ -1804,6 +1825,9 @@ bool LoadConstraint(const ConstraintProto& ct, Model* m) {
case ConstraintProto::ConstraintProto::kCumulative:
LoadCumulativeConstraint(ct, m);
return true;
case ConstraintProto::ConstraintProto::kReservoir:
LoadReservoirConstraint(ct, m);
return true;
case ConstraintProto::ConstraintProto::kElement:
LoadElementConstraint(ct, m);
return true;

1
ortools/sat/cp_model_loader.h
View File

@@ -268,6 +268,7 @@ void LoadIntMaxConstraint(const ConstraintProto& ct, Model* m);
void LoadNoOverlapConstraint(const ConstraintProto& ct, Model* m);
void LoadNoOverlap2dConstraint(const ConstraintProto& ct, Model* m);
void LoadCumulativeConstraint(const ConstraintProto& ct, Model* m);
void LoadReservoirConstraint(const ConstraintProto& ct, Model* m);
void LoadElementConstraintBounds(const ConstraintProto& ct, Model* m);
void LoadElementConstraintAC(const ConstraintProto& ct, Model* m);
void LoadElementConstraint(const ConstraintProto& ct, Model* m);

15
ortools/sat/cp_model_postsolve.cc
View File

@@ -85,13 +85,16 @@ void PostsolveLinear(const ConstraintProto& ct,
// The postsolve code is a bit involved if there is more than one free
// variable, we have to postsolve them one by one.
//
// Note that if there are some not free variable that are not assigned, the
// presolve should have made it sure that whatever the value of those
// variables the first variable can always be assigned to satisfy the
// constraint. In that case, the MultiplicationBy() below might be inexact,
// but that should be okay.
std::vector<Domain> to_add;
to_add.push_back(Domain(0));
for (int i = 0; i + 1 < free_vars.size(); ++i) {
bool exact = false;
Domain term =
(*domains)[free_vars[i]].MultiplicationBy(-free_coeffs[i], &exact);
CHECK(exact);
Domain term = (*domains)[free_vars[i]].MultiplicationBy(-free_coeffs[i]);
to_add.push_back(term.AdditionWith(to_add.back()));
}
for (int i = free_vars.size() - 1; i >= 0; --i) {
@@ -102,6 +105,10 @@ void PostsolveLinear(const ConstraintProto& ct,
const Domain domain = rhs.AdditionWith(to_add[i])
.InverseMultiplicationBy(coeff)
.IntersectionWith((*domains)[var]);
// TODO(user): I am not 100% that the algo here might cover all the presolve
// case, so if this fail, it might indicate an issue here and not in the
// presolve/solver code.
CHECK(!domain.IsEmpty()) << ct.ShortDebugString();
const int64 value = prefer_lower_value[var] ? domain.Min() : domain.Max();
(*domains)[var] = Domain(value);

174
ortools/sat/cp_model_presolve.cc
View File

@@ -1794,9 +1794,22 @@ bool CpModelPresolver::PropagateDomainsInLinear(int c, ConstraintProto* ct) {
// The variable now only appear in its definition and we can remove it
// because it was implied free.
//
// Tricky: If the linear constraint contains other variables that are only
// used here, then the postsolve needs more info. We do need to indicate
// that whatever the value of those other variables, we will have a way to
// assign var. We do that by putting it fist.
CHECK_EQ(context_->VarToConstraints(var).size(), 1);
context_->MarkVariableAsRemoved(var);
const int ct_index = context_->mapping_model->constraints().size();
*context_->mapping_model->add_constraints() = *ct;
LinearConstraintProto* mapping_linear_ct =
context_->mapping_model->mutable_constraints(ct_index)
->mutable_linear();
std::swap(mapping_linear_ct->mutable_vars()->at(0),
mapping_linear_ct->mutable_vars()->at(i));
std::swap(mapping_linear_ct->mutable_coeffs()->at(0),
mapping_linear_ct->mutable_coeffs()->at(i));
return RemoveConstraint(ct);
}
if (new_bounds) {
@@ -3459,9 +3472,164 @@ bool CpModelPresolver::PresolveAutomaton(ConstraintProto* ct) {
return false;
}
bool CpModelPresolver::PresolveReservoir(ConstraintProto* ct) {
if (context_->ModelIsUnsat()) return false;
if (HasEnforcementLiteral(*ct)) return false;
bool changed = false;
ReservoirConstraintProto& mutable_reservoir = *ct->mutable_reservoir();
if (mutable_reservoir.actives().empty()) {
const int true_literal = context_->GetOrCreateConstantVar(1);
for (int i = 0; i < mutable_reservoir.times_size(); ++i) {
mutable_reservoir.add_actives(true_literal);
}
changed = true;
}
const auto& demand_is_null = [&](int i) {
return mutable_reservoir.demands(i) == 0 ||
context_->LiteralIsFalse(mutable_reservoir.actives(i));
};
// Remove zero demands, and inactive events.
int num_zeros = 0;
for (int i = 0; i < mutable_reservoir.demands_size(); ++i) {
if (demand_is_null(i)) num_zeros++;
}
if (num_zeros > 0) { // Remove null events
changed = true;
int new_size = 0;
for (int i = 0; i < mutable_reservoir.demands_size(); ++i) {
if (demand_is_null(i)) continue;
mutable_reservoir.set_demands(new_size, mutable_reservoir.demands(i));
mutable_reservoir.set_times(new_size, mutable_reservoir.times(i));
mutable_reservoir.set_actives(new_size, mutable_reservoir.actives(i));
new_size++;
}
mutable_reservoir.mutable_demands()->Truncate(new_size);
mutable_reservoir.mutable_times()->Truncate(new_size);
mutable_reservoir.mutable_actives()->Truncate(new_size);
context_->UpdateRuleStats(
"reservoir: remove zero demands or inactive events.");
}
const int num_events = mutable_reservoir.demands_size();
int64 gcd = mutable_reservoir.demands().empty()
? 0
: std::abs(mutable_reservoir.demands(0));
int num_positives = 0;
int num_negatives = 0;
int64 sum_of_demands = 0;
int64 max_sum_of_positive_demands = 0;
int64 min_sum_of_negative_demands = 0;
for (int i = 0; i < num_events; ++i) {
const int64 demand = mutable_reservoir.demands(i);
sum_of_demands += demand;
gcd = MathUtil::GCD64(gcd, std::abs(demand));
if (demand > 0) {
num_positives++;
max_sum_of_positive_demands += demand;
} else {
DCHECK_LT(demand, 0);
num_negatives++;
min_sum_of_negative_demands += demand;
}
}
if (min_sum_of_negative_demands >= mutable_reservoir.min_level() &&
max_sum_of_positive_demands <= mutable_reservoir.max_level()) {
context_->UpdateRuleStats("reservoir: always feasible");
return RemoveConstraint(ct);
}
if (min_sum_of_negative_demands > mutable_reservoir.max_level() ||
max_sum_of_positive_demands < mutable_reservoir.min_level()) {
context_->UpdateRuleStats("reservoir: trivially infeasible");
return context_->NotifyThatModelIsUnsat();
}
if (min_sum_of_negative_demands > mutable_reservoir.min_level()) {
mutable_reservoir.set_min_level(min_sum_of_negative_demands);
context_->UpdateRuleStats(
"reservoir: increase min_level to reachable value");
}
if (max_sum_of_positive_demands < mutable_reservoir.max_level()) {
mutable_reservoir.set_max_level(max_sum_of_positive_demands);
context_->UpdateRuleStats("reservoir: reduce max_level to reachable value");
}
if (mutable_reservoir.min_level() <= 0 &&
mutable_reservoir.max_level() >= 0 &&
(num_positives == 0 || num_negatives == 0)) {
// If all demands have the same sign, and if the initial state is
// always feasible, we do not care about the order, just the sum.
auto* const sum =
context_->working_model->add_constraints()->mutable_linear();
int64 fixed_contrib = 0;
for (int i = 0; i < mutable_reservoir.demands_size(); ++i) {
const int64 demand = mutable_reservoir.demands(i);
DCHECK_NE(demand, 0);
const int active = mutable_reservoir.actives(i);
if (RefIsPositive(active)) {
sum->add_vars(active);
sum->add_coeffs(demand);
} else {
sum->add_vars(PositiveRef(active));
sum->add_coeffs(-demand);
fixed_contrib += demand;
}
}
sum->add_domain(mutable_reservoir.min_level() - fixed_contrib);
sum->add_domain(mutable_reservoir.max_level() - fixed_contrib);
context_->UpdateRuleStats("reservoir: converted to linear");
return RemoveConstraint(ct);
}
if (gcd > 1) {
for (int i = 0; i < mutable_reservoir.demands_size(); ++i) {
mutable_reservoir.set_demands(i, mutable_reservoir.demands(i) / gcd);
}
// Adjust min and max levels.
// max level is always rounded down.
// min level is always rounded up.
const Domain reduced_domain =
Domain({mutable_reservoir.min_level(), mutable_reservoir.max_level()})
.InverseMultiplicationBy(gcd);
mutable_reservoir.set_min_level(reduced_domain.Min());
mutable_reservoir.set_max_level(reduced_domain.Max());
context_->UpdateRuleStats("reservoir: simplify demands and levels by gcd.");
}
if (num_positives == 1 && num_negatives > 0) {
context_->UpdateRuleStats(
"TODO reservoir: one producer, multiple consumers.");
}
absl::flat_hash_set<std::pair<int, int>> time_active_set;
for (int i = 0; i < mutable_reservoir.demands_size(); ++i) {
const std::pair<int, int> key = std::make_pair(
mutable_reservoir.times(i), mutable_reservoir.actives(i));
if (time_active_set.contains(key)) {
context_->UpdateRuleStats("TODO reservoir: merge synchronized events.");
break;
} else {
time_active_set.insert(key);
}
}
return changed;
}
// TODO(user): It is probably more efficient to keep all the bool_and in a
// global place during all the presolve, and just output them at the end rather
// than modifying more than once the proto.
// global place during all the presolve, and just output them at the end
// rather than modifying more than once the proto.
void CpModelPresolver::ExtractBoolAnd() {
absl::flat_hash_map<int, int> ref_to_bool_and;
const int num_constraints = context_->working_model->constraints_size();
@@ -4255,6 +4423,8 @@ bool CpModelPresolver::PresolveOneConstraint(int c) {
return PresolveRoutes(ct);
case ConstraintProto::ConstraintCase::kAutomaton:
return PresolveAutomaton(ct);
case ConstraintProto::ConstraintCase::kReservoir:
return PresolveReservoir(ct);
default:
return false;
}

1
ortools/sat/cp_model_presolve.h
View File

@@ -99,6 +99,7 @@ class CpModelPresolver {
bool PresolveRoutes(ConstraintProto* ct);
bool PresolveCumulative(ConstraintProto* ct);
bool PresolveNoOverlap(ConstraintProto* ct);
bool PresolveReservoir(ConstraintProto* ct);
bool PresolveAllDiff(ConstraintProto* ct);
bool PresolveTable(ConstraintProto* ct);
bool PresolveElement(ConstraintProto* ct);

39
ortools/sat/cumulative.cc
View File

@@ -231,5 +231,44 @@ std::function<void(Model*)> CumulativeTimeDecomposition(
};
}
std::function<void(Model*)> CumulativeUsingReservoir(
const std::vector<IntervalVariable>& vars,
const std::vector<AffineExpression>& demands, AffineExpression capacity,
SchedulingConstraintHelper* helper) {
return [=](Model* model) {
if (vars.empty()) return;
auto* integer_trail = model->GetOrCreate<IntegerTrail>();
auto* encoder = model->GetOrCreate<IntegerEncoder>();
auto* intervals = model->GetOrCreate<IntervalsRepository>();
CHECK(integer_trail->IsFixed(capacity));
const IntegerValue fixed_capacity(
integer_trail->UpperBound(capacity).value());
std::vector<AffineExpression> times;
std::vector<IntegerValue> deltas;
std::vector<Literal> presences;
const int num_tasks = vars.size();
for (int t = 0; t < num_tasks; ++t) {
CHECK(integer_trail->IsFixed(demands[t]));
times.push_back(intervals->StartVar(vars[t]));
deltas.push_back(integer_trail->LowerBound(demands[t]));
times.push_back(intervals->EndVar(vars[t]));
deltas.push_back(-integer_trail->LowerBound(demands[t]));
if (intervals->IsOptional(vars[t])) {
presences.push_back(intervals->PresenceLiteral(vars[t]));
presences.push_back(intervals->PresenceLiteral(vars[t]));
} else {
presences.push_back(encoder->GetTrueLiteral());
presences.push_back(encoder->GetTrueLiteral());
}
}
AddReservoirConstraint(times, deltas, presences, 0, fixed_capacity.value(),
model);
};
}
} // namespace sat
} // namespace operations_research

6
ortools/sat/cumulative.h
View File

@@ -57,6 +57,12 @@ std::function<void(Model*)> CumulativeTimeDecomposition(
const std::vector<AffineExpression>& demands, AffineExpression capacity,
SchedulingConstraintHelper* helper = nullptr);
// Another testing code, same assumptions as the CumulativeTimeDecomposition().
std::function<void(Model*)> CumulativeUsingReservoir(
const std::vector<IntervalVariable>& vars,
const std::vector<AffineExpression>& demands, AffineExpression capacity,
SchedulingConstraintHelper* helper = nullptr);
} // namespace sat
} // namespace operations_research

29
ortools/sat/integer.cc
View File

@@ -993,6 +993,35 @@ bool IntegerTrail::Enqueue(IntegerLiteral i_lit,
integer_trail_.size());
}
bool IntegerTrail::ConditionalEnqueue(
Literal lit, IntegerLiteral i_lit, std::vector<Literal>* literal_reason,
std::vector<IntegerLiteral>* integer_reason) {
const VariablesAssignment& assignment = trail_->Assignment();
if (assignment.LiteralIsFalse(lit)) return true;
// We can always push var if the optional literal is the same.
//
// TODO(user): we can also push lit.var if its presence implies lit.
if (lit.Index() == OptionalLiteralIndex(i_lit.var)) {
return Enqueue(i_lit, *literal_reason, *integer_reason);
}
if (assignment.LiteralIsTrue(lit)) {
literal_reason->push_back(lit.Negated());
return Enqueue(i_lit, *literal_reason, *integer_reason);
}
if (IntegerLiteralIsFalse(i_lit)) {
integer_reason->push_back(
IntegerLiteral::LowerOrEqual(i_lit.var, i_lit.bound - 1));
EnqueueLiteral(lit.Negated(), *literal_reason, *integer_reason);
return true;
}
// We can't push anything in this case.
return true;
}
bool IntegerTrail::Enqueue(IntegerLiteral i_lit,
absl::Span<const Literal> literal_reason,
absl::Span<const IntegerLiteral> integer_reason,

10
ortools/sat/integer.h
View File

@@ -724,6 +724,16 @@ class IntegerTrail : public SatPropagator {
IntegerLiteral i_lit, absl::Span<const Literal> literal_reason,
absl::Span<const IntegerLiteral> integer_reason);
// Pushes the given integer literal assuming that the Boolean literal is true.
// This can do a few things:
// - If lit it true, add it to the reason and push the integer bound.
// - If the bound is infeasible, push lit to false.
// - If the underlying variable is optional and also controlled by lit, push
// the bound even if lit is not assigned.
ABSL_MUST_USE_RESULT bool ConditionalEnqueue(
Literal lit, IntegerLiteral i_lit, std::vector<Literal>* literal_reason,
std::vector<IntegerLiteral>* integer_reason);
// Same as Enqueue(), but takes an extra argument which if smaller than
// integer_trail_.size() is interpreted as the trail index of an old Enqueue()
// that had the same reason as this one. Note that the given Span must still

27
ortools/sat/intervals.cc
View File

@@ -384,31 +384,12 @@ bool SchedulingConstraintHelper::PushIntegerLiteralIfTaskPresent(
int t, IntegerLiteral lit) {
if (IsAbsent(t)) return true;
AddOtherReason(t);
// TODO(user): we can also push lit.var if its presence implies the interval
// presence.
if (IsOptional(t) && integer_trail_->OptionalLiteralIndex(lit.var) !=
reason_for_presence_[t]) {
if (IsPresent(t)) {
// We can still push, but we do need the presence reason.
AddPresenceReason(t);
} else {
// In this case we cannot push lit.var, but we may detect the interval
// absence.
if (lit.bound > integer_trail_->UpperBound(lit.var)) {
integer_reason_.push_back(
IntegerLiteral::LowerOrEqual(lit.var, lit.bound - 1));
if (!PushTaskAbsence(t)) return false;
}
return true;
}
}
ImportOtherReasons();
if (!integer_trail_->Enqueue(lit, literal_reason_, integer_reason_)) {
return false;
if (IsOptional(t)) {
return integer_trail_->ConditionalEnqueue(
PresenceLiteral(t), lit, &literal_reason_, &integer_reason_);
}
return true;
return integer_trail_->Enqueue(lit, literal_reason_, integer_reason_);
}
// We also run directly the precedence propagator for this variable so that when

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

@@ -1037,38 +1037,48 @@ class CpModel(object):
"""Adds Reservoir(times, demands, min_level, max_level).
Maintains a reservoir level within bounds. The water level starts at 0, and
at any time >= 0, it must be between min_level and max_level. Furthermore,
this constraint expects all times variables to be >= 0.
at any time, it must be between min_level and max_level.
If the variable `times[i]` is assigned a value t, then the current level
changes by `demands[i]`, which is constant, at time t.
Note that level min can be > 0, or level max can be < 0. It just forces
some demands to be executed at time 0 to make sure that we are within those
bounds with the executed demands. Therefore, at any time t >= 0:
Note that min level must be <= 0, and the max level must be >= 0. Please
use fixed demands to simulate initial state.
sum(demands[i] if times[i] <= t) in [min_level, max_level]
Therefore, at any time:
sum(demands[i] if times[i] <= t) in [min_level, max_level]
Args:
times: A list of positive integer variables which specify the time of the
times: A list of integer variables which specify the time of the
filling or emptying the reservoir.
demands: A list of integer values that specifies the amount of the
emptying or filling.
min_level: At any time >= 0, the level of the reservoir must be greater of
min_level: At any time, the level of the reservoir must be greater or
equal than the min level.
max_level: At any time >= 0, the level of the reservoir must be less or
equal than the max level.
max_level: At any time, the level of the reservoir must be less or equal
than the max level.
Returns:
An instance of the `Constraint` class.
Raises:
ValueError: if max_level < min_level.
ValueError: if max_level < 0.
ValueError: if min_level > 0
"""
if max_level < min_level:
return ValueError(
'Reservoir constraint must have a max_level >= min_level')
if max_level < 0:
return ValueError('Reservoir constraint must have a max_level >= 0')
if min_level > 0:
return ValueError('Reservoir constraint must have a min_level <= 0')
ct = Constraint(self.__model.constraints)
model_ct = self.__model.constraints[ct.Index()]
model_ct.reservoir.times.extend([self.GetOrMakeIndex(x) for x in times])
@@ -1081,46 +1091,56 @@ class CpModel(object):
min_level, max_level):
"""Adds Reservoir(times, demands, actives, min_level, max_level).
Maintain a reservoir level within bounds. The water level starts at 0, and
at
any time >= 0, it must be within min_level, and max_level. Furthermore, this
constraints expect all times variables to be >= 0.
If `actives[i]` is true, and if `times[i]` is assigned a value t, then the
level of the reservoir changes by `demands[i]`, which is constant, at
time t.
Maintains a reservoir level within bounds. The water level starts at 0, and
at any time, it must be between min_level and max_level.
Note that level_min can be > 0, or level_max can be < 0. It just forces
some demands to be executed at time 0 to make sure that we are within those
bounds with the executed demands. Therefore, at any time t >= 0:
If the variable `times[i]` is assigned a value t, and `actives[i]` is
`True`, then the current level changes by `demands[i]`, which is constant,
at time t.
Note that min level must be <= 0, and the max level must be >= 0. Please
use fixed demands to simulate initial state.
Therefore, at any time:
sum(demands[i] * actives[i] if times[i] <= t) in [min_level, max_level]
sum(demands[i] * actives[i] if times[i] <= t) in [min_level, max_level]
The array of boolean variables 'actives', if defined, indicates which
actions are actually performed.
Args:
times: A list of positive integer variables which specify the time of the
times: A list of integer variables which specify the time of the
filling or emptying the reservoir.
demands: A list of integer values that specifies the amount of the
emptying or filling.
actives: a list of boolean variables. They indicates if the
emptying/refilling events actually take place.
min_level: At any time >= 0, the level of the reservoir must be greater of
min_level: At any time, the level of the reservoir must be greater or
equal than the min level.
max_level: At any time >= 0, the level of the reservoir must be less or
equal than the max level.
max_level: At any time, the level of the reservoir must be less or equal
than the max level.
Returns:
An instance of the `Constraint` class.
Raises:
ValueError: if max_level < min_level.
ValueError: if max_level < 0.
ValueError: if min_level > 0
"""
if max_level < min_level:
return ValueError(
'Reservoir constraint must have a max_level >= min_level')
if max_level < 0:
return ValueError('Reservoir constraint must have a max_level >= 0')
if min_level > 0:
return ValueError('Reservoir constraint must have a min_level <= 0')
ct = Constraint(self.__model.constraints)
model_ct = self.__model.constraints[ct.Index()]
model_ct.reservoir.times.extend([self.GetOrMakeIndex(x) for x in times])

10
ortools/sat/sat_parameters.proto
View File

@@ -21,7 +21,7 @@ option java_multiple_files = true;
// Contains the definitions for all the sat algorithm parameters and their
// default values.
//
// NEXT TAG: 181
// NEXT TAG: 183
message SatParameters {
// In some context, like in a portfolio of search, it makes sense to name a
// given parameters set for logging purpose.
@@ -463,6 +463,14 @@ message SatParameters {
// Permutations (#Variables = #Values) are always expanded.
optional bool expand_alldiff_constraints = 170 [default = false];
// If true, expand the reservoir constraints by creating booleans for all
// possible precedences between event and encoding the constraint.
optional bool expand_reservoir_constraints = 182 [default = true];
// If true, it disable all constraint expansion.
// This should only be used to test the presolve of expanded constraints.
optional bool disable_constraint_expansion = 181 [default = false];
// During presolve, we use a maximum clique heuristic to merge together
// no-overlap constraints or at most one constraints. This code can be slow,
// so we have a limit in place on the number of explored nodes in the

252
ortools/sat/timetable.cc
View File

@@ -24,6 +24,258 @@
namespace operations_research {
namespace sat {
void AddReservoirConstraint(std::vector<AffineExpression> times,
std::vector<IntegerValue> deltas,
std::vector<Literal> presences, int64 min_level,
int64 max_level, Model* model) {
// We only create a side if it can fail.
IntegerValue min_possible(0);
IntegerValue max_possible(0);
for (const IntegerValue d : deltas) {
if (d > 0) {
max_possible += d;
} else {
min_possible += d;
}
}
if (max_possible > max_level) {
model->TakeOwnership(new ReservoirTimeTabling(
times, deltas, presences, IntegerValue(max_level), model));
}
if (min_possible < min_level) {
for (IntegerValue& ref : deltas) ref = -ref;
model->TakeOwnership(new ReservoirTimeTabling(
times, deltas, presences, IntegerValue(-min_level), model));
}
}
ReservoirTimeTabling::ReservoirTimeTabling(
const std::vector<AffineExpression>& times,
const std::vector<IntegerValue>& deltas,
const std::vector<Literal>& presences, IntegerValue capacity, Model* model)
: times_(times),
deltas_(deltas),
presences_(presences),
capacity_(capacity),
assignment_(model->GetOrCreate<Trail>()->Assignment()),
integer_trail_(model->GetOrCreate<IntegerTrail>()) {
auto* watcher = model->GetOrCreate<GenericLiteralWatcher>();
const int id = watcher->Register(this);
const int num_events = times.size();
for (int e = 0; e < num_events; e++) {
if (deltas_[e] > 0) {
watcher->WatchUpperBound(times_[e].var, id);
watcher->WatchLiteral(presences_[e], id);
}
if (deltas_[e] < 0) {
watcher->WatchLowerBound(times_[e].var, id);
watcher->WatchLiteral(presences_[e].Negated(), id);
}
}
watcher->NotifyThatPropagatorMayNotReachFixedPointInOnePass(id);
}
bool ReservoirTimeTabling::Propagate() {
const int num_events = times_.size();
if (!BuildProfile()) return false;
for (int e = 0; e < num_events; e++) {
if (assignment_.LiteralIsFalse(presences_[e])) continue;
// For positive delta, we can maybe increase the min.
if (deltas_[e] > 0 && !TryToIncreaseMin(e)) return false;
// For negative delta, we can maybe decrease the max.
if (deltas_[e] < 0 && !TryToDecreaseMax(e)) return false;
}
return true;
}
// We compute the lowest possible profile at time t.
//
// TODO(user): If we have precedences between events, we should be able to do
// more.
bool ReservoirTimeTabling::BuildProfile() {
// Starts by copying the "events" in the profile and sort them by time.
profile_.clear();
const int num_events = times_.size();
profile_.emplace_back(kMinIntegerValue, IntegerValue(0)); // Sentinel.
for (int e = 0; e < num_events; e++) {
if (deltas_[e] > 0) {
// Only consider present event for positive delta.
if (!assignment_.LiteralIsTrue(presences_[e])) continue;
const IntegerValue ub = integer_trail_->UpperBound(times_[e]);
profile_.push_back({ub, deltas_[e]});
} else if (deltas_[e] < 0) {
// Only consider non-absent event for negative delta.
if (assignment_.LiteralIsFalse(presences_[e])) continue;
profile_.push_back({integer_trail_->LowerBound(times_[e]), deltas_[e]});
}
}
profile_.emplace_back(kMaxIntegerValue, IntegerValue(0)); // Sentinel.
std::sort(profile_.begin(), profile_.end());
// Accumulate delta and collapse entries.
int last = 0;
for (const ProfileRectangle& rect : profile_) {
if (rect.start == profile_[last].start) {
profile_[last].height += rect.height;
} else {
++last;
profile_[last].start = rect.start;
profile_[last].height = rect.height + profile_[last - 1].height;
}
}
profile_.resize(last + 1);
// Conflict?
for (const ProfileRectangle& rect : profile_) {
if (rect.height <= capacity_) continue;
FillReasonForProfileAtGivenTime(rect.start);
return integer_trail_->ReportConflict(literal_reason_, integer_reason_);
}
return true;
}
// TODO(user): Minimize with how high the profile needs to be. We can also
// remove from the reason the absence of a negative event provided that the
// level zero min of the event is greater than t anyway.
//
// TODO(user): Make sure the code work with fixed time since pushing always
// true/false literal to the reason is not completely supported.
void ReservoirTimeTabling::FillReasonForProfileAtGivenTime(
IntegerValue t, int event_to_ignore) {
integer_reason_.clear();
literal_reason_.clear();
const int num_events = times_.size();
for (int e = 0; e < num_events; e++) {
if (e == event_to_ignore) continue;
if (deltas_[e] > 0) {
if (!assignment_.LiteralIsTrue(presences_[e])) continue;
if (integer_trail_->UpperBound(times_[e]) > t) continue;
integer_reason_.push_back(times_[e].LowerOrEqual(t));
literal_reason_.push_back(presences_[e].Negated());
} else if (deltas_[e] < 0) {
if (assignment_.LiteralIsFalse(presences_[e])) {
literal_reason_.push_back(presences_[e]);
} else if (integer_trail_->LowerBound(times_[e]) > t) {
integer_reason_.push_back(times_[e].GreaterOrEqual(t + 1));
}
}
}
}
// Note that a negative event will always be in the profile, even if its
// presence is still not settled.
bool ReservoirTimeTabling::TryToDecreaseMax(int event) {
CHECK_LT(deltas_[event], 0);
const IntegerValue start = integer_trail_->LowerBound(times_[event]);
const IntegerValue end = integer_trail_->UpperBound(times_[event]);
// We already tested for conflict in BuildProfile().
if (start == end) return true;
// Find the profile rectangle that overlaps the start of the given event.
// The sentinel prevents out of bound exceptions.
DCHECK(std::is_sorted(profile_.begin(), profile_.end()));
int rec_id =
std::upper_bound(profile_.begin(), profile_.end(), start,
[&](IntegerValue value, const ProfileRectangle& rect) {
return value < rect.start;
}) -
profile_.begin();
--rec_id;
bool push = false;
IntegerValue new_end = end;
for (; profile_[rec_id].start < end; ++rec_id) {
if (profile_[rec_id].height - deltas_[event] > capacity_) {
new_end = profile_[rec_id].start;
push = true;
break;
}
}
if (!push) return true;
// The reason is simply why the capacity at new_end (without the event)
// would overflow.
FillReasonForProfileAtGivenTime(new_end, event);
// Note(user): I don't think this is possible since it would have been
// detected at profile construction, but then, since the bound might have been
// updated, better be defensive.
if (new_end < start) {
integer_reason_.push_back(times_[event].GreaterOrEqual(new_end + 1));
return integer_trail_->ReportConflict(literal_reason_, integer_reason_);
}
// First, the task MUST be present, otherwise we have a conflict.
//
// TODO(user): We actually need to look after 'end' to potentially push the
// presence in more situation.
if (!assignment_.LiteralIsTrue(presences_[event])) {
integer_trail_->EnqueueLiteral(presences_[event], literal_reason_,
integer_reason_);
}
// Push new_end too. Note that we don't need the presence reason.
return integer_trail_->Enqueue(times_[event].LowerOrEqual(new_end),
literal_reason_, integer_reason_);
}
bool ReservoirTimeTabling::TryToIncreaseMin(int event) {
CHECK_GT(deltas_[event], 0);
const IntegerValue start = integer_trail_->LowerBound(times_[event]);
const IntegerValue end = integer_trail_->UpperBound(times_[event]);
// We already tested for conflict in BuildProfile().
if (start == end) return true;
// Find the profile rectangle containing the end of the given event.
// The sentinel prevents out of bound exceptions.
//
// TODO(user): If the task is no present, we should actually look at the
// maximum profile after end to maybe push its absence.
DCHECK(std::is_sorted(profile_.begin(), profile_.end()));
int rec_id =
std::upper_bound(profile_.begin(), profile_.end(), end,
[&](IntegerValue value, const ProfileRectangle& rect) {
return value < rect.start;
}) -
profile_.begin();
--rec_id;
bool push = false;
IntegerValue new_start = start;
if (profile_[rec_id].height + deltas_[event] > capacity_) {
if (!assignment_.LiteralIsTrue(presences_[event])) {
// Push to false since it wasn't part of the profile and cannot fit.
push = true;
new_start = end + 1;
} else if (profile_[rec_id].start < end) {
// It must be at end in this case.
push = true;
new_start = end;
}
}
if (!push) {
for (; profile_[rec_id].start > start; --rec_id) {
if (profile_[rec_id - 1].height + deltas_[event] > capacity_) {
push = true;
new_start = profile_[rec_id].start;
break;
}
}
}
if (!push) return true;
// The reason is simply the capacity at new_start - 1;
FillReasonForProfileAtGivenTime(new_start - 1, event);
return integer_trail_->ConditionalEnqueue(
presences_[event], times_[event].GreaterOrEqual(new_start),
&literal_reason_, &integer_reason_);
}
TimeTablingPerTask::TimeTablingPerTask(
const std::vector<AffineExpression>& demands, AffineExpression capacity,
IntegerTrail* integer_trail, SchedulingConstraintHelper* helper)

77
ortools/sat/timetable.h
View File

@@ -24,6 +24,80 @@
namespace operations_research {
namespace sat {
// Adds a reservoir constraint to the model. Note that to account for level not
// containing zero at time zero, we might needs to create an artificial fixed
// event.
//
// This instantiate one or more ReservoirTimeTabling class to perform the
// propagation.
void AddReservoirConstraint(std::vector<AffineExpression> times,
std::vector<IntegerValue> deltas,
std::vector<Literal> presences, int64 min_level,
int64 max_level, Model* model);
// The piecewise constant function must be below the given capacity. The initial
// function value is zero. Note that a negative capacity will thus be trivially
// infeasible.
//
// Note that we take for the definition of the function at time t to be the sum
// of all delta with time <= t. But because we check for the capacity over the
// full horizon, we could have taken < t with no behavior change.
class ReservoirTimeTabling : public PropagatorInterface {
public:
ReservoirTimeTabling(const std::vector<AffineExpression>& times,
const std::vector<IntegerValue>& deltas,
const std::vector<Literal>& presences,
IntegerValue capacity, Model* model);
bool Propagate() final;
private:
// The rectangle will be ordered by start, and the end of each rectangle
// will be equal to the start of the next one. The height correspond to the
// one from start (inclusive) until the next one (exclusive).
struct ProfileRectangle {
ProfileRectangle() {}
ProfileRectangle(IntegerValue start, IntegerValue height)
: start(start), height(height) {}
bool operator<(const ProfileRectangle& other) const {
return start < other.start;
}
/* const */ IntegerValue start = IntegerValue(0);
/* const */ IntegerValue height = IntegerValue(0);
};
// Builds the profile and increases the lower bound of the capacity
// variable accordingly.
bool BuildProfile();
// Explanation of the profile minimum value at time t, eventually ignoring the
// given event.
void FillReasonForProfileAtGivenTime(IntegerValue t,
int event_to_ignore = -1);
// Tries to tighten the min/max time of the given event depending on the sign
// of the delta associated with this event.
bool TryToIncreaseMin(int event);
bool TryToDecreaseMax(int event);
// Input.
std::vector<AffineExpression> times_;
std::vector<IntegerValue> deltas_;
std::vector<Literal> presences_;
IntegerValue capacity_;
// Model class.
const VariablesAssignment& assignment_;
IntegerTrail* integer_trail_;
// Temporary data.
std::vector<Literal> literal_reason_;
std::vector<IntegerLiteral> integer_reason_;
std::vector<ProfileRectangle> profile_;
};
// A strongly quadratic version of Time Tabling filtering. This propagator
// is similar to the CumulativeTimeTable propagator of the constraint solver.
class TimeTablingPerTask : public PropagatorInterface {
@@ -38,7 +112,8 @@ class TimeTablingPerTask : public PropagatorInterface {
private:
// The rectangle will be ordered by start, and the end of each rectangle
// will be equal to the start of the next one.
// will be equal to the start of the next one. The height correspond to the
// one from start (inclusive) until the next one (exclusive).
struct ProfileRectangle {
/* const */ IntegerValue start;
/* const */ IntegerValue height;