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

[CP-SAT] Support optional intervals in no_overlap_2d; add energetic cut to no_overlap_2d; disable cumulative relaxation of no_overlap_2d by default; fix #3012

Browse Source
This commit is contained in:
Laurent Perron
2021-12-13 10:47:57 +01:00
parent 93934e3660
commit dbe177a0fb
9 changed files with 360 additions and 20 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
ortools/sat/cp_model_loader.cc
View File

@@ -1213,7 +1213,9 @@ void LoadNoOverlap2dConstraint(const ConstraintProto& ct, Model* m) {
mapping->Intervals(ct.no_overlap_2d().y_intervals());
m->Add(NonOverlappingRectangles(
x_intervals, y_intervals,
!ct.no_overlap_2d().boxes_with_null_area_can_overlap()));
!ct.no_overlap_2d().boxes_with_null_area_can_overlap(),
m->GetOrCreate<SatParameters>()
->use_cumulative_constraint_in_no_overlap_2d_constraint()));
}
void LoadCumulativeConstraint(const ConstraintProto& ct, Model* m) {

10
ortools/sat/diffn.cc
View File

@@ -228,6 +228,8 @@ bool NonOverlappingRectanglesEnergyPropagator::FailWhenEnergyIsTooLarge(
const auto add_box_energy_in_rectangle_reason = [&](int b) {
x_.AddEnergyMinInIntervalReason(b, area.x_min, area.x_max);
y_.AddEnergyMinInIntervalReason(b, area.y_min, area.y_max);
x_.AddPresenceReason(b);
y_.AddPresenceReason(b);
};
for (int i = 0; i < neighbors_.size(); ++i) {
@@ -508,6 +510,8 @@ bool NonOverlappingRectanglesDisjunctivePropagator::Propagate() {
global_x_.AddReasonForBeingBefore(box1, box2);
global_x_.AddReasonForBeingBefore(box2, box1);
global_y_.ClearReason();
global_y_.AddPresenceReason(box1);
global_y_.AddPresenceReason(box2);
global_y_.AddReasonForBeingBefore(box1, box2);
global_y_.AddReasonForBeingBefore(box2, box1);
global_x_.ImportOtherReasons(global_y_);
@@ -531,6 +535,8 @@ bool NonOverlappingRectanglesDisjunctivePropagator::PropagateTwoBoxes() {
const IntegerValue left_end_min = x_.EndMin(left);
if (left_end_min > x_.StartMin(right)) {
x_.ClearReason();
x_.AddPresenceReason(left);
x_.AddPresenceReason(right);
x_.AddReasonForBeingBefore(left, right);
x_.AddEndMinReason(left, left_end_min);
RETURN_IF_FALSE(x_.IncreaseStartMin(right, left_end_min));
@@ -540,6 +546,8 @@ bool NonOverlappingRectanglesDisjunctivePropagator::PropagateTwoBoxes() {
const IntegerValue right_start_max = x_.StartMax(right);
if (right_start_max < x_.EndMax(left)) {
x_.ClearReason();
x_.AddPresenceReason(left);
x_.AddPresenceReason(right);
x_.AddReasonForBeingBefore(left, right);
x_.AddStartMaxReason(right, right_start_max);
RETURN_IF_FALSE(x_.DecreaseEndMax(left, right_start_max));
@@ -551,6 +559,8 @@ bool NonOverlappingRectanglesDisjunctivePropagator::PropagateTwoBoxes() {
switch (state) {
case 0: { // Conflict.
x_.ClearReason();
x_.AddPresenceReason(0);
x_.AddPresenceReason(1);
x_.AddReasonForBeingBefore(0, 1);
x_.AddReasonForBeingBefore(1, 0);
return x_.ReportConflict();

29
ortools/sat/diffn.h
View File

@@ -14,6 +14,7 @@
#ifndef OR_TOOLS_SAT_DIFFN_H_
#define OR_TOOLS_SAT_DIFFN_H_
#include <functional>
#include <vector>
#include "ortools/base/int_type.h"
@@ -167,8 +168,32 @@ inline std::function<void(Model*)> NonOverlappingRectangles(
model->TakeOwnership(constraint);
if (add_cumulative_relaxation) {
AddCumulativeRelaxation(x, x_helper, y_helper, model);
AddCumulativeRelaxation(y, y_helper, x_helper, model);
// We must first check if the cumulative relaxation is possible.
bool x_has_sole_optional_intervals = false;
bool y_has_sole_optional_intervals = false;
for (int i = 0; i < x.size(); ++i) {
if (x_helper->IsOptional(i) && y_helper->IsOptional(i) &&
x_helper->PresenceLiteral(i) != y_helper->PresenceLiteral(i)) {
// Abort as the task would be conditioned by two literals.
return;
}
if (x_helper->IsOptional(i) && !y_helper->IsOptional(i)) {
// We cannot use x_size as the demand of the cumulative based on
// the y_intervals.
x_has_sole_optional_intervals = true;
}
if (y_helper->IsOptional(i) && !x_helper->IsOptional(i)) {
// We cannot use y_size as the demand of the cumulative based on
// the y_intervals.
y_has_sole_optional_intervals = true;
}
}
if (!y_has_sole_optional_intervals) {
AddCumulativeRelaxation(x, x_helper, y_helper, model);
}
if (!x_has_sole_optional_intervals) {
AddCumulativeRelaxation(y, y_helper, x_helper, model);
}
}
};
}

3
ortools/sat/diffn_util.cc
View File

@@ -66,6 +66,9 @@ bool ReportEnergyConflict(Rectangle bounding_box, absl::Span<const int> boxes,
x->AddEnergyMinInIntervalReason(b, bounding_box.x_min, bounding_box.x_max);
y->AddEnergyMinInIntervalReason(b, bounding_box.y_min, bounding_box.y_max);
x->AddPresenceReason(b);
y->AddPresenceReason(b);
total_energy += x->SizeMin(b) * y->SizeMin(b);
// We abort early if a subset of boxes is enough.

69
ortools/sat/linear_relaxation.cc
View File

@@ -28,6 +28,7 @@
#include "ortools/sat/implied_bounds.h"
#include "ortools/sat/integer.h"
#include "ortools/sat/integer_expr.h"
#include "ortools/sat/intervals.h"
#include "ortools/sat/linear_constraint.h"
#include "ortools/sat/linear_programming_constraint.h"
#include "ortools/sat/sat_base.h"
@@ -802,14 +803,30 @@ void AppendNoOverlap2dRelaxation(const ConstraintProto& ct, Model* model,
LinearConstraintBuilder lc(model, IntegerValue(0), max_area);
for (int i = 0; i < ct.no_overlap_2d().x_intervals_size(); ++i) {
DCHECK(!intervals_repository->IsOptional(x_intervals[i]));
DCHECK(!intervals_repository->IsOptional(y_intervals[i]));
LinearConstraintBuilder linear_energy(model);
if (DetectLinearEncodingOfProducts(x_sizes[i], y_sizes[i], model,
&linear_energy)) {
lc.AddLinearExpression(linear_energy.BuildExpression());
} else {
lc.AddQuadraticLowerBound(x_sizes[i], y_sizes[i], integer_trail);
if (intervals_repository->IsPresent(x_intervals[i]) &&
intervals_repository->IsPresent(y_intervals[i])) {
LinearConstraintBuilder linear_energy(model);
if (DetectLinearEncodingOfProducts(x_sizes[i], y_sizes[i], model,
&linear_energy)) {
lc.AddLinearExpression(linear_energy.BuildExpression());
} else {
lc.AddQuadraticLowerBound(x_sizes[i], y_sizes[i], integer_trail);
}
} else if ((intervals_repository->PresenceLiteral(x_intervals[i]) ==
intervals_repository->PresenceLiteral(y_intervals[i])) ||
intervals_repository->IsPresent(x_intervals[i]) ||
intervals_repository->IsPresent(y_intervals[i])) {
// We have only one active literal.
const Literal presence_literal =
intervals_repository->IsPresent(x_intervals[i])
? intervals_repository->PresenceLiteral(y_intervals[i])
: intervals_repository->PresenceLiteral(x_intervals[i]);
const IntegerValue area_min =
integer_trail->LevelZeroLowerBound(x_sizes[i]) *
integer_trail->LevelZeroLowerBound(y_sizes[i]);
if (area_min != 0 && !lc.AddLiteralTerm(presence_literal, area_min)) {
return;
}
}
}
relaxation->linear_constraints.push_back(lc.Build());
@@ -1114,8 +1131,7 @@ void TryToLinearizeConstraint(const CpModelProto& model_proto,
// Adds an energetic relaxation (sum of areas fits in bounding box).
AppendNoOverlap2dRelaxation(ct, model, relaxation);
if (linearization_level > 1) {
// Adds a completion time cut generator.
// TODO(user): investigate adding an energetic cut generator.
// Adds a completion time cut generator and an energetic cut generator.
AddNoOverlap2dCutGenerator(ct, model, relaxation);
}
break;
@@ -1291,10 +1307,39 @@ void AddNoOverlap2dCutGenerator(const ConstraintProto& ct, Model* m,
mapping->Intervals(ct.no_overlap_2d().x_intervals());
std::vector<IntervalVariable> y_intervals =
mapping->Intervals(ct.no_overlap_2d().y_intervals());
// TODO(user): We can add CumulativeEnergyCuts for no_overlap_2d if boxes
// do not have a fixed size.
relaxation->cut_generators.push_back(
CreateNoOverlap2dCompletionTimeCutGenerator(x_intervals, y_intervals, m));
// Checks if at least one rectangle has a variable dimension or is optional.
IntervalsRepository* intervals_repository =
m->GetOrCreate<IntervalsRepository>();
bool has_variable_part = false;
for (int i = 0; i < x_intervals.size(); ++i) {
if (intervals_repository->IsOptional(x_intervals[i]) ||
intervals_repository->IsOptional(y_intervals[i])) {
has_variable_part = true;
break;
}
if (intervals_repository->MinSize(x_intervals[i]) !=
intervals_repository->MaxSize(x_intervals[i]) ||
intervals_repository->MinSize(y_intervals[i]) !=
intervals_repository->MaxSize(y_intervals[i])) {
has_variable_part = true;
break;
}
}
if (has_variable_part) {
std::vector<LinearExpression> energies;
std::vector<AffineExpression> x_sizes;
std::vector<AffineExpression> y_sizes;
for (int i = 0; i < x_intervals.size(); ++i) {
x_sizes.push_back(intervals_repository->Size(x_intervals[i]));
y_sizes.push_back(intervals_repository->Size(y_intervals[i]));
}
FillEnergies(x_sizes, y_sizes, m, &energies);
relaxation->cut_generators.push_back(CreateNoOverlap2dEnergyCutGenerator(
x_intervals, y_intervals, energies, m));
}
}
void AddLinMaxCutGenerator(const ConstraintProto& ct, Model* m,

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

@@ -1760,6 +1760,9 @@ class CpModel(object):
on a plane. Each rectangle is aligned with the X and Y axis, and is defined
by two intervals which represent its projection onto the X and Y axis.
Furthermore, one box is optional if at least one of the x or y interval is
optional.
Args:
x_intervals: The X coordinates of the rectangles.
y_intervals: The Y coordinates of the rectangles.
@@ -1782,8 +1785,8 @@ class CpModel(object):
for all t:
sum(demands[i]
if (start(intervals[t]) <= t < end(intervals[t])) and
(t is present)) <= capacity
if (start(intervals[i]) <= t < end(intervals[i])) and
(intervals[i] is present)) <= capacity
Args:
intervals: The list of intervals.

10
ortools/sat/sat_parameters.proto
View File

@@ -24,7 +24,7 @@ option csharp_namespace = "Google.OrTools.Sat";
// Contains the definitions for all the sat algorithm parameters and their
// default values.
//
// NEXT TAG: 200
// NEXT TAG: 201
message SatParameters {
// In some context, like in a portfolio of search, it makes sense to name a
// given parameters set for logging purpose.
@@ -615,6 +615,14 @@ message SatParameters {
optional bool use_disjunctive_constraint_in_cumulative_constraint = 80
[default = true];
// When this is true, the no_overlap_2d constraint is reinforced with
// propagators from the cumulative constraints. It consists of ignoring the
// position of rectangles in one position and projecting the no_overlap_2d on
// the other dimension to create a cumulative constraint. This is done on both
// axis.This additional level supplements the default level of reasoning.
optional bool use_cumulative_constraint_in_no_overlap_2d_constraint = 200
[default = false];
// A non-negative level indicating the type of constraints we consider in the
// LP relaxation. At level zero, no LP relaxation is used. At level 1, only
// the linear constraint and full encoding are added. At level 2, we also add

224
ortools/sat/scheduling_cuts.cc
View File

@@ -1074,5 +1074,229 @@ CutGenerator CreateNoOverlap2dCompletionTimeCutGenerator(
return result;
}
void GenerateNoOverlap2dEnergyCut(
const absl::StrongVector<IntegerVariable, double>& lp_values, Model* model,
IntegerTrail* integer_trail, IntegerEncoder* encoder,
LinearConstraintManager* manager,
const std::vector<LinearExpression>& energies, absl::Span<int> boxes,
const std::string& cut_name, SchedulingConstraintHelper* x_helper,
SchedulingConstraintHelper* y_helper) {
// Sort tasks by StartMin, tie breaked by EndMax.
std::sort(boxes.begin(), boxes.end(), [x_helper](int a, int b) {
return x_helper->StartMin(a) < x_helper->StartMin(b) ||
(x_helper->StartMin(a) == x_helper->StartMin(b) &&
x_helper->EndMax(a) < x_helper->EndMax(b));
});
for (int i1 = 0; i1 + 1 < boxes.size(); ++i1) {
// For each start time, we will keep the most violated cut
// generated while scanning the residual intervals.
int end_index_of_max_violation = -1;
double max_relative_violation = 1.01;
IntegerValue x_min_of_max_violation(0);
IntegerValue x_max_of_max_violation(0);
IntegerValue y_min_of_max_violation(0);
IntegerValue y_max_of_max_violation(0);
// Accumulates intervals and check for potential cuts.
double energy_lp = 0.0;
IntegerValue x_min = kMaxIntegerValue;
IntegerValue x_max = kMinIntegerValue;
IntegerValue y_min = kMaxIntegerValue;
IntegerValue y_max = kMinIntegerValue;
// We sort all tasks (start_min(task) >= start_min(start_index) by
// increasing end max.
std::vector<int> residual_intervals(boxes.begin() + i1, boxes.end());
std::sort(residual_intervals.begin(), residual_intervals.end(),
[&](int a, int b) {
return x_helper->EndMax(a) < x_helper->EndMax(b);
});
// Let's process residual tasks and evaluate the cut violation of
// the cut at each step. We follow the same structure as the cut
// creation code below.
for (int i2 = 0; i2 < residual_intervals.size(); ++i2) {
const int t = residual_intervals[i2];
if (x_helper->IsPresent(t) && y_helper->IsPresent(t)) {
if (EnergyIsDefined(energies[t])) {
energy_lp += energies[t].LpValue(lp_values);
} else { // demand and size are not fixed.
energy_lp += x_helper->Sizes()[t].LpValue(lp_values) *
ToDouble(y_helper->SizeMin(t));
energy_lp += y_helper->Sizes()[t].LpValue(lp_values) *
ToDouble(x_helper->SizeMin(t));
energy_lp -= ToDouble(x_helper->SizeMin(t) * y_helper->SizeMin(t));
}
} else {
const Literal lit = x_helper->IsOptional(t)
? x_helper->PresenceLiteral(t)
: y_helper->PresenceLiteral(t);
// TODO(user): Use the energy min if better than size_min *
// demand_min, here and when building the cut.
energy_lp += GetLiteralLpValue(lit, lp_values, encoder) *
ToDouble(x_helper->SizeMin(t) * y_helper->SizeMin(t));
}
// Update bounding box.
x_min = std::min(x_min, x_helper->StartMin(t));
x_max = std::max(x_max, x_helper->EndMax(t));
y_min = std::min(y_min, y_helper->StartMin(t));
y_max = std::max(y_max, y_helper->EndMax(t));
// Compute the violation of the potential cut.
const double relative_violation =
energy_lp / ToDouble((x_max - x_min) * (y_max - y_min));
if (relative_violation > max_relative_violation) {
end_index_of_max_violation = i2;
max_relative_violation = relative_violation;
x_min_of_max_violation = x_min;
x_max_of_max_violation = x_max;
y_min_of_max_violation = y_min;
y_max_of_max_violation = y_max;
}
}
if (end_index_of_max_violation == -1) return;
// A maximal violated cut has been found.
bool cut_generated = true;
bool has_opt_cuts = false;
bool has_quadratic_cuts = false;
bool use_energy = false;
LinearConstraintBuilder cut(
model, kMinIntegerValue,
(x_max_of_max_violation - x_min_of_max_violation) *
(y_max_of_max_violation - y_min_of_max_violation));
// Build the cut.
for (int i2 = 0; i2 <= end_index_of_max_violation; ++i2) {
const int t = residual_intervals[i2];
if (x_helper->IsPresent(t) && y_helper->IsPresent(t)) {
if (EnergyIsDefined(energies[t])) {
// We favor the energy info instead of the McCormick relaxation.
cut.AddLinearExpression(energies[t]);
use_energy = true;
} else {
// This will add linear term if the size is fixed.
cut.AddQuadraticLowerBound(x_helper->Sizes()[t], y_helper->Sizes()[t],
integer_trail);
has_quadratic_cuts = true;
}
} else {
// TODO(user): use the offset of the energy expression if better
// than size_min * demand_min.
has_opt_cuts = true;
if (!x_helper->SizeIsFixed(t) && !y_helper->SizeIsFixed(t)) {
has_quadratic_cuts = true;
}
const Literal lit = x_helper->IsOptional(t)
? x_helper->PresenceLiteral(t)
: y_helper->PresenceLiteral(t);
if (!cut.AddLiteralTerm(lit,
x_helper->SizeMin(t) * y_helper->SizeMin(t))) {
cut_generated = false;
break;
}
}
}
if (cut_generated) {
std::string full_name = cut_name;
if (has_opt_cuts) full_name.append("_opt");
if (has_quadratic_cuts) full_name.append("_quad");
if (use_energy) full_name.append("_energy");
manager->AddCut(cut.Build(), full_name, lp_values);
}
}
}
CutGenerator CreateNoOverlap2dEnergyCutGenerator(
const std::vector<IntervalVariable>& x_intervals,
const std::vector<IntervalVariable>& y_intervals,
const std::vector<LinearExpression>& energies, Model* model) {
CutGenerator result;
SchedulingConstraintHelper* x_helper =
new SchedulingConstraintHelper(x_intervals, model);
model->TakeOwnership(x_helper);
SchedulingConstraintHelper* y_helper =
new SchedulingConstraintHelper(y_intervals, model);
model->TakeOwnership(y_helper);
AddIntegerVariableFromIntervals(x_helper, model, &result.vars);
AddIntegerVariableFromIntervals(y_helper, model, &result.vars);
Trail* trail = model->GetOrCreate<Trail>();
IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
IntegerEncoder* encoder = model->GetOrCreate<IntegerEncoder>();
result.generate_cuts =
[integer_trail, trail, encoder, x_helper, y_helper, model, energies](
const absl::StrongVector<IntegerVariable, double>& lp_values,
LinearConstraintManager* manager) {
if (trail->CurrentDecisionLevel() > 0) return true;
if (!x_helper->SynchronizeAndSetTimeDirection(true)) return false;
if (!y_helper->SynchronizeAndSetTimeDirection(true)) return false;
const int num_boxes = x_helper->NumTasks();
std::vector<int> active_boxes;
std::vector<Rectangle> cached_rectangles(num_boxes);
for (int box = 0; box < num_boxes; ++box) {
if (y_helper->IsAbsent(box) || y_helper->IsAbsent(box)) continue;
// We cannot consider boxes controlled by 2 enforcement literals.
if (x_helper->IsOptional(box) && y_helper->IsOptional(box) &&
x_helper->PresenceLiteral(box) !=
y_helper->PresenceLiteral(box)) {
continue;
}
// TODO(user): It might be possible/better to use some shifted value
// here, but for now this code is not in the hot spot, so better be
// defensive and only do connected components on really disjoint
// boxes.
Rectangle& rectangle = cached_rectangles[box];
rectangle.x_min = x_helper->StartMin(box);
rectangle.x_max = x_helper->EndMax(box);
rectangle.y_min = y_helper->StartMin(box);
rectangle.y_max = y_helper->EndMax(box);
active_boxes.push_back(box);
}
if (active_boxes.size() <= 1) return true;
std::vector<absl::Span<int>> components =
GetOverlappingRectangleComponents(cached_rectangles,
absl::MakeSpan(active_boxes));
for (absl::Span<int> boxes : components) {
if (boxes.size() <= 1) continue;
if (!x_helper->SynchronizeAndSetTimeDirection(true)) return false;
if (!y_helper->SynchronizeAndSetTimeDirection(true)) return false;
GenerateNoOverlap2dEnergyCut(
lp_values, model, integer_trail, encoder, manager, energies,
boxes, "NoOverlap2dXEnergy", x_helper, y_helper);
GenerateNoOverlap2dEnergyCut(
lp_values, model, integer_trail, encoder, manager, energies,
boxes, "NoOverlap2dYEnergy", y_helper, x_helper);
if (!x_helper->SynchronizeAndSetTimeDirection(false)) return false;
if (!y_helper->SynchronizeAndSetTimeDirection(false)) return false;
GenerateNoOverlap2dEnergyCut(
lp_values, model, integer_trail, encoder, manager, energies,
boxes, "NoOverlap2dXEnergyMirror", x_helper, y_helper);
GenerateNoOverlap2dEnergyCut(
lp_values, model, integer_trail, encoder, manager, energies,
boxes, "NoOverlap2dYEnergyMirror", y_helper, x_helper);
}
return true;
};
return result;
}
} // namespace sat
} // namespace operations_research

24
ortools/sat/scheduling_cuts.h
View File

@@ -32,8 +32,8 @@
namespace operations_research {
namespace sat {
// For a given set of intervals and demands, we compute the maximum energy of
// each task and make sure it is less than the span of the intervals * its
// For a given set of intervals and demands, we compute the energy of
// each task and make sure their sum fits in the span of the intervals * its
// capacity.
//
// If an interval is optional, it contributes
@@ -89,6 +89,26 @@ CutGenerator CreateNoOverlap2dCompletionTimeCutGenerator(
const std::vector<IntervalVariable>& x_intervals,
const std::vector<IntervalVariable>& y_intervals, Model* model);
// Energetic cuts for the no_overlap_2d constraint.
//
// For a given set of rectangles, we compute the area of each rectangle
// and make sure their sum is less than the area of the bounding interval.
//
// If an interval is optional, it contributes
// min_size_x * min_size_y * presence_literal
// amount of total area.
//
// If an interval is performed, we use the linear area formulation (if
// defined, that is if different from a constant -1), or the McCormick
// relaxation of the size_x * size_y if not defined.
//
// The maximum area is the area of the bounding rectangle of each intervals
// at level 0.
CutGenerator CreateNoOverlap2dEnergyCutGenerator(
const std::vector<IntervalVariable>& x_intervals,
const std::vector<IntervalVariable>& y_intervals,
const std::vector<LinearExpression>& energies, Model* model);
// For a given set of intervals, we first compute the min and max of all
// intervals. Then we create a cut that indicates that all intervals must fit
// in that span.