From 6a99deb2a37da0048de27ced5da501075ec462c0 Mon Sep 17 00:00:00 2001 From: Corentin Le Molgat Date: Mon, 23 Aug 2021 16:08:30 +0200 Subject: [PATCH] Export constraint_solver to github --- .../constraint_solver/constraint_solveri.h | 3 +- ortools/constraint_solver/local_search.cc | 15 +- ortools/constraint_solver/routing.cc | 138 ++- ortools/constraint_solver/routing.h | 142 ++- ortools/constraint_solver/routing_breaks.cc | 2 +- ortools/constraint_solver/routing_filters.cc | 69 +- ortools/constraint_solver/routing_filters.h | 3 +- ortools/constraint_solver/routing_flow.cc | 4 +- .../routing_lp_scheduling.cc | 219 +++-- .../constraint_solver/routing_lp_scheduling.h | 111 ++- .../routing_neighborhoods.cc | 109 +-- .../constraint_solver/routing_neighborhoods.h | 8 +- .../constraint_solver/routing_parameters.cc | 14 + .../routing_parameters.proto | 11 +- ortools/constraint_solver/routing_sat.cc | 811 +++++++++++++++--- 15 files changed, 1340 insertions(+), 319 deletions(-) diff --git a/ortools/constraint_solver/constraint_solveri.h b/ortools/constraint_solver/constraint_solveri.h index 559f27a0a6..c1c8a745be 100644 --- a/ortools/constraint_solver/constraint_solveri.h +++ b/ortools/constraint_solver/constraint_solveri.h @@ -3479,7 +3479,8 @@ LocalSearchFilter* MakePathStateFilter(Solver* solver, // - Synchronize() must be called before. // - Revert() must be called before. LocalSearchFilter* MakeUnaryDimensionFilter( - Solver* solver, std::unique_ptr checker); + Solver* solver, std::unique_ptr checker, + const std::string& dimension_name); #endif // !defined(SWIG) diff --git a/ortools/constraint_solver/local_search.cc b/ortools/constraint_solver/local_search.cc index e688367728..2f5a3ed365 100644 --- a/ortools/constraint_solver/local_search.cc +++ b/ortools/constraint_solver/local_search.cc @@ -3182,9 +3182,11 @@ namespace { class UnaryDimensionFilter : public LocalSearchFilter { public: - std::string DebugString() const override { return "UnaryDimensionFilter"; } - explicit UnaryDimensionFilter(std::unique_ptr checker) - : checker_(std::move(checker)) {} + std::string DebugString() const override { return name_; } + UnaryDimensionFilter(std::unique_ptr checker, + const std::string& dimension_name) + : checker_(std::move(checker)), + name_(absl::StrCat("UnaryDimensionFilter(", dimension_name, ")")) {} bool Accept(const Assignment* delta, const Assignment* deltadelta, int64_t objective_min, int64_t objective_max) override { @@ -3198,13 +3200,16 @@ class UnaryDimensionFilter : public LocalSearchFilter { private: std::unique_ptr checker_; + const std::string name_; }; } // namespace LocalSearchFilter* MakeUnaryDimensionFilter( - Solver* solver, std::unique_ptr checker) { - UnaryDimensionFilter* filter = new UnaryDimensionFilter(std::move(checker)); + Solver* solver, std::unique_ptr checker, + const std::string& dimension_name) { + UnaryDimensionFilter* filter = + new UnaryDimensionFilter(std::move(checker), dimension_name); return solver->RevAlloc(filter); } diff --git a/ortools/constraint_solver/routing.cc b/ortools/constraint_solver/routing.cc index b53818d472..9da174c44b 100644 --- a/ortools/constraint_solver/routing.cc +++ b/ortools/constraint_solver/routing.cc @@ -100,6 +100,8 @@ class TwoOpt; namespace operations_research { namespace { +using ResourceGroup = RoutingModel::ResourceGroup; + // A decision builder which tries to assign values to variables as close as // possible to target values first. // TODO(user): Move to CP solver. @@ -338,8 +340,11 @@ class SetCumulsFromGlobalDimensionCosts : public DecisionBuilder { SetCumulsFromGlobalDimensionCosts( const std::vector>* global_optimizers, + const std::vector>* + global_mp_optimizers, SearchMonitor* monitor, bool optimize_and_pack = false) : global_optimizers_(*global_optimizers), + global_mp_optimizers_(*global_mp_optimizers), monitor_(monitor), optimize_and_pack_(optimize_and_pack) {} Decision* Next(Solver* const solver) override { @@ -347,7 +352,8 @@ class SetCumulsFromGlobalDimensionCosts : public DecisionBuilder { // order to postpone the Fail() call until after the for loop, so there are // no memory leaks related to the cumul_values vector. bool should_fail = false; - for (const auto& global_optimizer : global_optimizers_) { + for (int i = 0; i < global_optimizers_.size(); ++i) { + const auto& global_optimizer = global_optimizers_[i]; const RoutingDimension* dimension = global_optimizer->dimension(); RoutingModel* const model = dimension->model(); @@ -359,16 +365,34 @@ class SetCumulsFromGlobalDimensionCosts : public DecisionBuilder { std::vector cumul_values; std::vector break_start_end_values; - const bool cumuls_optimized = + const DimensionSchedulingStatus status = optimize_and_pack_ ? global_optimizer->ComputePackedCumuls(next, &cumul_values, &break_start_end_values) : global_optimizer->ComputeCumuls(next, &cumul_values, &break_start_end_values); - if (!cumuls_optimized) { + if (status == DimensionSchedulingStatus::INFEASIBLE) { should_fail = true; break; } + // If relaxation is not feasible, try the MILP optimizer. + if (status == DimensionSchedulingStatus::RELAXED_OPTIMAL_ONLY) { + cumul_values.clear(); + break_start_end_values.clear(); + DCHECK(global_mp_optimizers_[i] != nullptr); + const DimensionSchedulingStatus mp_status = + optimize_and_pack_ + ? global_mp_optimizers_[i]->ComputePackedCumuls( + next, &cumul_values, &break_start_end_values) + : global_mp_optimizers_[i]->ComputeCumuls( + next, &cumul_values, &break_start_end_values); + if (mp_status != DimensionSchedulingStatus::OPTIMAL) { + should_fail = true; + break; + } + } else { + DCHECK(status == DimensionSchedulingStatus::OPTIMAL); + } // Concatenate cumul_values and break_start_end_values into cp_values, // generate corresponding cp_variables vector. std::vector cp_variables = dimension->cumuls(); @@ -409,6 +433,8 @@ class SetCumulsFromGlobalDimensionCosts : public DecisionBuilder { private: const std::vector>& global_optimizers_; + const std::vector>& + global_mp_optimizers_; SearchMonitor* const monitor_; const bool optimize_and_pack_; }; @@ -447,7 +473,7 @@ const Assignment* RoutingModel::PackCumulsOfOptimizerDimensionsFromAssignment( /*optimize_and_pack=*/true))); decision_builders.push_back( solver_->RevAlloc(new SetCumulsFromGlobalDimensionCosts( - &global_dimension_optimizers_, + &global_dimension_optimizers_, &global_dimension_mp_optimizers_, GetOrCreateLargeNeighborhoodSearchLimit(), /*optimize_and_pack=*/true))); decision_builders.push_back( @@ -1181,6 +1207,18 @@ GlobalDimensionCumulOptimizer* RoutingModel::GetMutableGlobalCumulOptimizer( return global_dimension_optimizers_[optimizer_index].get(); } +GlobalDimensionCumulOptimizer* RoutingModel::GetMutableGlobalCumulMPOptimizer( + const RoutingDimension& dimension) const { + const DimensionIndex dim_index = GetDimensionIndex(dimension.name()); + if (dim_index < 0 || dim_index >= global_optimizer_index_.size() || + global_optimizer_index_[dim_index] < 0) { + return nullptr; + } + const int optimizer_index = global_optimizer_index_[dim_index]; + DCHECK_LT(optimizer_index, global_dimension_mp_optimizers_.size()); + return global_dimension_mp_optimizers_[optimizer_index].get(); +} + LocalDimensionCumulOptimizer* RoutingModel::GetMutableLocalCumulOptimizer( const RoutingDimension& dimension) const { const DimensionIndex dim_index = GetDimensionIndex(dimension.name()); @@ -1229,6 +1267,72 @@ RoutingDimension* RoutingModel::GetMutableDimension( return nullptr; } +// ResourceGroup +ResourceGroup::Attributes::Attributes() + : start_domain_(Domain::AllValues()), end_domain_(Domain::AllValues()) { + /// The default attributes have unconstrained start/end domains. +} + +RoutingModel::ResourceGroup::Attributes::Attributes(Domain start_domain, + Domain end_domain) + : start_domain_(std::move(start_domain)), + end_domain_(std::move(end_domain)) {} + +const ResourceGroup::Attributes& +ResourceGroup::Resource::GetDimensionAttributes( + const RoutingDimension* dimension) const { + DimensionIndex dimension_index = model_->GetDimensionIndex(dimension->name()); + DCHECK_NE(dimension_index, kNoDimension); + return gtl::FindWithDefault(dimension_attributes_, dimension_index, + GetDefaultAttributes()); +} + +void ResourceGroup::Resource::SetDimensionAttributes( + Attributes attributes, const RoutingDimension* dimension) { + DCHECK(dimension_attributes_.empty()) + << "As of 2021/07, each resource can only constrain a single dimension."; + + const DimensionIndex dimension_index = + model_->GetDimensionIndex(dimension->name()); + DCHECK_NE(dimension_index, kNoDimension); + DCHECK(!dimension_attributes_.contains(dimension_index)); + dimension_attributes_[dimension_index] = std::move(attributes); +} + +const ResourceGroup::Attributes& ResourceGroup::Resource::GetDefaultAttributes() + const { + static const Attributes* const kAttributes = new Attributes(); + return *kAttributes; +} + +RoutingModel::ResourceGroup* const RoutingModel::AddResourceGroup() { + resource_groups_.push_back(absl::make_unique(this)); + return resource_groups_.back().get(); +} + +void RoutingModel::ResourceGroup::AddResource( + Attributes attributes, const RoutingDimension* dimension) { + resources_.push_back(Resource(model_)); + resources_.back().SetDimensionAttributes(std::move(attributes), dimension); + + const DimensionIndex dimension_index = + model_->GetDimensionIndex(dimension->name()); + DCHECK_NE(dimension_index, kNoDimension); + affected_dimension_indices_.insert(dimension_index); + + DCHECK_EQ(affected_dimension_indices_.size(), 1) + << "As of 2021/07, each ResourceGroup can only affect a single " + "RoutingDimension at a time."; +} + +const std::vector& RoutingModel::GetDimensionResourceGroupIndices( + const RoutingDimension* dimension) const { + DCHECK(closed_); + const DimensionIndex dim = GetDimensionIndex(dimension->name()); + DCHECK_NE(dim, kNoDimension); + return dimension_resource_group_indices_[dim]; +} + void RoutingModel::SetArcCostEvaluatorOfAllVehicles(int evaluator_index) { CHECK_LT(0, vehicles_); for (int i = 0; i < vehicles_; ++i) { @@ -1719,8 +1823,9 @@ void RoutingModel::AddPickupAndDelivery(int64_t pickup, int64_t delivery) { void RoutingModel::AddPickupAndDeliverySets( DisjunctionIndex pickup_disjunction, DisjunctionIndex delivery_disjunction) { - AddPickupAndDeliverySetsInternal(GetDisjunctionIndices(pickup_disjunction), - GetDisjunctionIndices(delivery_disjunction)); + AddPickupAndDeliverySetsInternal( + GetDisjunctionNodeIndices(pickup_disjunction), + GetDisjunctionNodeIndices(delivery_disjunction)); pickup_delivery_disjunctions_.push_back( {pickup_disjunction, delivery_disjunction}); } @@ -2124,6 +2229,15 @@ void RoutingModel::CloseModelWithParameters( for (RoutingDimension* const dimension : dimensions_) { dimension->CloseModel(UsesLightPropagation(parameters)); } + + dimension_resource_group_indices_.resize(dimensions_.size()); + for (int rg_index = 0; rg_index < resource_groups_.size(); rg_index++) { + for (DimensionIndex dim_index : + resource_groups_[rg_index]->GetAffectedDimensionIndices()) { + dimension_resource_group_indices_[dim_index].push_back(rg_index); + } + } + ComputeCostClasses(parameters); ComputeVehicleClasses(); ComputeVehicleTypes(); @@ -2705,7 +2819,8 @@ const Assignment* RoutingModel::SolveFromAssignmentsWithParameters( } } - if (parameters.use_cp_sat() == BOOL_TRUE) { + if (parameters.use_cp_sat() == BOOL_TRUE || + parameters.use_generalized_cp_sat() == BOOL_TRUE) { const int solution_count = collect_assignments_->solution_count(); Assignment* const cp_solution = solution_count >= 1 ? collect_assignments_->solution(solution_count - 1) @@ -4354,7 +4469,11 @@ void RoutingModel::StoreDimensionCumulOptimizers( // Use global optimizer. global_optimizer_index_[dim] = global_dimension_optimizers_.size(); global_dimension_optimizers_.push_back( - absl::make_unique(dimension)); + absl::make_unique( + dimension, parameters.continuous_scheduling_solver())); + global_dimension_mp_optimizers_.push_back( + absl::make_unique( + dimension, parameters.mixed_integer_scheduling_solver())); packed_dimensions_collector_assignment->Add(dimension->cumuls()); if (!AllTransitsPositive(*dimension)) { dimension->SetOffsetForGlobalOptimizer(0); @@ -4514,7 +4633,8 @@ DecisionBuilder* RoutingModel::CreateSolutionFinalizer(SearchLimit* lns_limit) { if (!global_dimension_optimizers_.empty()) { decision_builders.push_back( solver_->RevAlloc(new SetCumulsFromGlobalDimensionCosts( - &global_dimension_optimizers_, lns_limit))); + &global_dimension_optimizers_, &global_dimension_mp_optimizers_, + lns_limit))); } decision_builders.push_back( CreateFinalizerForMinimizedAndMaximizedVariables()); diff --git a/ortools/constraint_solver/routing.h b/ortools/constraint_solver/routing.h index 2637de8dc6..4577b4c003 100644 --- a/ortools/constraint_solver/routing.h +++ b/ortools/constraint_solver/routing.h @@ -176,6 +176,7 @@ #include "ortools/base/integral_types.h" #include "ortools/base/logging.h" #include "ortools/base/macros.h" +#include "ortools/base/map_util.h" #include "ortools/base/strong_vector.h" #include "ortools/constraint_solver/constraint_solver.h" #include "ortools/constraint_solver/constraint_solveri.h" @@ -375,6 +376,75 @@ class RoutingModel { // clang-format on }; + /// A ResourceGroup defines a set of available Resources with attributes on + /// one or multiple dimensions. + /// For every ResourceGroup in the model, each (used) vehicle in the solution + /// must be assigned to exactly 1 resource, and each resource can in turn be + /// assigned to at most 1 vehicle. This vehicle-to-resource assignment will + /// apply the corresponding Attributes to the dimensions affected by the + /// resource group. + /// NOTE: As of 2021/07, each ResourceGroup can only affect a single + /// RoutingDimension at a time, i.e. all Resources in a group must apply + /// attributes to the same single dimension. + class ResourceGroup { + public: + /// Attributes for a dimension. + class Attributes { + public: + Attributes(); + Attributes(Domain start_domain, Domain end_domain); + + const Domain& start_domain() const { return start_domain_; } + const Domain& end_domain() const { return end_domain_; } + + private: + /// The following domains constrain the dimension start/end cumul of the + /// the vehicle assigned to this resource: + /// start_domain_.Min() <= cumul[Start(v)] <= start_domain_.Max() + Domain start_domain_; + /// end_domain_.Min() <= cumul[End(v)] <= end_domain_.Max() + Domain end_domain_; + }; + + /// A Resource sets attributes (costs/constraints) for a set of dimensions. + class Resource { + public: + const Attributes& GetDimensionAttributes( + const RoutingDimension* dimension) const; + + private: + explicit Resource(const RoutingModel* model) : model_(model) {} + + void SetDimensionAttributes(Attributes attributes, + const RoutingDimension* dimension); + const Attributes& GetDefaultAttributes() const; + + const RoutingModel* const model_; + absl::flat_hash_map dimension_attributes_; + + friend class ResourceGroup; + }; + + explicit ResourceGroup(const RoutingModel* model) : model_(model) {} + + /// Add a Resource with the given attributes for the corresponding + /// dimension. + void AddResource(Attributes attributes, const RoutingDimension* dimension); + + const std::vector& GetResources() const { return resources_; } + const absl::flat_hash_set& GetAffectedDimensionIndices() + const { + return affected_dimension_indices_; + } + int Size() const { return resources_.size(); } + + private: + const RoutingModel* const model_; + std::vector resources_; + /// All indices of dimensions affected by this resource group. + absl::flat_hash_set affected_dimension_indices_; + }; + /// Constant used to express a hard constraint instead of a soft penalty. static const int64_t kNoPenalty; @@ -563,6 +633,10 @@ class RoutingModel { GetGlobalDimensionCumulOptimizers() const { return global_dimension_optimizers_; } + const std::vector >& + GetGlobalDimensionCumulMPOptimizers() const { + return global_dimension_mp_optimizers_; + } const std::vector >& GetLocalDimensionCumulOptimizers() const { return local_dimension_optimizers_; @@ -577,6 +651,8 @@ class RoutingModel { /// or nullptr if there is none. GlobalDimensionCumulOptimizer* GetMutableGlobalCumulOptimizer( const RoutingDimension& dimension) const; + GlobalDimensionCumulOptimizer* GetMutableGlobalCumulMPOptimizer( + const RoutingDimension& dimension) const; LocalDimensionCumulOptimizer* GetMutableLocalCumulOptimizer( const RoutingDimension& dimension) const; LocalDimensionCumulOptimizer* GetMutableLocalCumulMPOptimizer( @@ -603,6 +679,19 @@ class RoutingModel { const std::string& GetPrimaryConstrainedDimension() const { return primary_constrained_dimension_; } + + ResourceGroup* const AddResourceGroup(); + // clang-format off + const std::vector >& GetResourceGroups() + const { + return resource_groups_; + } + // clang-format on + /// Returns the indices of resource groups for this dimension. This method can + /// only be called after the model has been closed. + const std::vector& GetDimensionResourceGroupIndices( + const RoutingDimension* dimension) const; + /// Adds a disjunction constraint on the indices: exactly 'max_cardinality' of /// the indices are active. Start and end indices of any vehicle cannot be /// part of a disjunction. @@ -644,7 +733,7 @@ class RoutingModel { #if !defined(SWIGPYTHON) /// Returns the variable indices of the nodes in the disjunction of index /// 'index'. - const std::vector& GetDisjunctionIndices( + const std::vector& GetDisjunctionNodeIndices( DisjunctionIndex index) const { return disjunctions_[index].indices; } @@ -875,9 +964,9 @@ class RoutingModel { } /// Get the "unperformed" penalty of a node. This is only well defined if the - /// node is only part of a single Disjunction involving only itself, and that - /// disjunction has a penalty. In all other cases, including forced active - /// nodes, this returns 0. + /// node is only part of a single Disjunction, and that disjunction has a + /// penalty. For forced active nodes returns max int64_t. In all other cases, + /// this returns 0. int64_t UnperformedPenalty(int64_t var_index) const; /// Same as above except that it returns default_value instead of 0 when /// penalty is not well defined (default value is passed as first argument to @@ -1279,6 +1368,21 @@ class RoutingModel { DCHECK(closed_); return vehicle_class_index_of_vehicle_[vehicle]; } + /// Returns a vehicle of the given vehicle class, and -1 if there are no + /// vehicles for this class. + int GetVehicleOfClass(VehicleClassIndex vehicle_class) const { + DCHECK(closed_); + const RoutingModel::VehicleTypeContainer& vehicle_type_container = + GetVehicleTypeContainer(); + if (vehicle_class.value() >= GetVehicleClassesCount() || + vehicle_type_container.vehicles_per_vehicle_class[vehicle_class.value()] + .empty()) { + return -1; + } + return vehicle_type_container + .vehicles_per_vehicle_class[vehicle_class.value()] + .front(); + } /// Returns the number of different vehicle classes in the model. int GetVehicleClassesCount() const { return vehicle_classes_.size(); } /// Returns variable indices of nodes constrained to be on the same route. @@ -1744,12 +1848,23 @@ class RoutingModel { /// Dimensions absl::flat_hash_map dimension_name_to_index_; absl::StrongVector dimensions_; + /// Resource Groups. + /// If resource_groups_ is not empty, then for each group of resources, each + /// (used) vehicle must be assigned to exactly 1 resource, and each resource + /// can in turn be assigned to at most 1 vehicle. // clang-format off + std::vector > resource_groups_; + /// Stores the set of resource groups related to each dimension. + absl::StrongVector > + dimension_resource_group_indices_; + /// TODO(user): Define a new Dimension[Global|Local]OptimizerIndex type /// and use it to define ITIVectors and for the dimension to optimizer index /// mappings below. std::vector > global_dimension_optimizers_; + std::vector > + global_dimension_mp_optimizers_; absl::StrongVector global_optimizer_index_; std::vector > local_dimension_optimizers_; @@ -1941,6 +2056,7 @@ class RoutingModel { friend class RoutingDimension; friend class RoutingModelInspector; + friend class ResourceGroup::Resource; DISALLOW_COPY_AND_ASSIGN(RoutingModel); }; @@ -2458,6 +2574,16 @@ class RoutingDimension { return model_->TransitCallback( class_evaluators_[vehicle_to_class_[vehicle]]); } + + /// Returns the callback evaluating the transit value between two node indices + /// for a given vehicle class. + const RoutingModel::TransitCallback2& class_transit_evaluator( + RoutingVehicleClassIndex vehicle_class) const { + const int vehicle = model_->GetVehicleOfClass(vehicle_class); + DCHECK_NE(vehicle, -1); + return transit_evaluator(vehicle); + } + /// Returns the unary callback evaluating the transit value between two node /// indices for a given vehicle. If the corresponding callback is not unary, /// returns a null callback. @@ -2683,6 +2809,14 @@ class RoutingDimension { int64_t GetSpanCostCoefficientForVehicle(int vehicle) const { return vehicle_span_cost_coefficients_[vehicle]; } +#ifndef SWIG + int64_t GetSpanCostCoefficientForVehicleClass( + RoutingVehicleClassIndex vehicle_class) const { + const int vehicle = model_->GetVehicleOfClass(vehicle_class); + DCHECK_NE(vehicle, -1); + return GetSpanCostCoefficientForVehicle(vehicle); + } +#endif // SWIG #ifndef SWIG const std::vector& vehicle_span_cost_coefficients() const { return vehicle_span_cost_coefficients_; diff --git a/ortools/constraint_solver/routing_breaks.cc b/ortools/constraint_solver/routing_breaks.cc index ddbf6968b1..06cae183c8 100644 --- a/ortools/constraint_solver/routing_breaks.cc +++ b/ortools/constraint_solver/routing_breaks.cc @@ -334,7 +334,7 @@ bool DisjunctivePropagator::DistanceDuration(Tasks* tasks) { // Skip breaks that cannot be performed after start. int index_break_by_emax = tasks->num_chain_tasks; while (index_break_by_emax < num_tasks && - tasks->end_max[index_break_by_emax] <= tasks->end_max[route_start]) { + tasks->end_max[index_break_by_emax] <= tasks->end_min[route_start]) { ++index_break_by_emax; } // Special case: no breaks after start. diff --git a/ortools/constraint_solver/routing_filters.cc b/ortools/constraint_solver/routing_filters.cc index 8a84d4ea37..480d4ea30f 100644 --- a/ortools/constraint_solver/routing_filters.cc +++ b/ortools/constraint_solver/routing_filters.cc @@ -212,7 +212,7 @@ class NodeDisjunctionFilter : public IntVarLocalSearchFilter { const int inactive_nodes = current_inactive_nodes + disjunction_inactive_delta.second; const int max_inactive_cardinality = - routing_model_.GetDisjunctionIndices(disjunction_index).size() - + routing_model_.GetDisjunctionNodeIndices(disjunction_index).size() - routing_model_.GetDisjunctionMaxCardinality(disjunction_index); // Too many inactive nodes. if (inactive_nodes > max_inactive_cardinality) { @@ -258,7 +258,7 @@ class NodeDisjunctionFilter : public IntVarLocalSearchFilter { active_per_disjunction_[i] = 0; inactive_per_disjunction_[i] = 0; const std::vector& disjunction_indices = - routing_model_.GetDisjunctionIndices(i); + routing_model_.GetDisjunctionNodeIndices(i); for (const int64_t index : disjunction_indices) { const bool index_synced = IsVarSynced(index); if (index_synced) { @@ -2199,7 +2199,7 @@ void AppendLightWeightDimensionFilters( std::move(demands), std::move(node_capacity)); const auto kAccept = LocalSearchFilterManager::FilterEventType::kAccept; LocalSearchFilter* filter = MakeUnaryDimensionFilter( - dimension->model()->solver(), std::move(checker)); + dimension->model()->solver(), std::move(checker), dimension->name()); filters->push_back({filter, kAccept}); } } @@ -2275,6 +2275,7 @@ void AppendDimensionCumulFilters( DCHECK(model.GetMutableGlobalCumulOptimizer(dimension) != nullptr); filters->push_back({MakeGlobalLPCumulFilter( model.GetMutableGlobalCumulOptimizer(dimension), + model.GetMutableGlobalCumulMPOptimizer(dimension), filter_objective_cost), kAccept}); } else if (use_cumul_bounds_propagator_filter[d]) { @@ -2609,6 +2610,7 @@ class LPCumulFilter : public IntVarLocalSearchFilter { public: LPCumulFilter(const std::vector& nexts, GlobalDimensionCumulOptimizer* optimizer, + GlobalDimensionCumulOptimizer* mp_optimizer, bool filter_objective_cost); bool Accept(const Assignment* delta, const Assignment* deltadelta, int64_t objective_min, int64_t objective_max) override; @@ -2621,6 +2623,7 @@ class LPCumulFilter : public IntVarLocalSearchFilter { private: GlobalDimensionCumulOptimizer& optimizer_; + GlobalDimensionCumulOptimizer& mp_optimizer_; const bool filter_objective_cost_; int64_t synchronized_cost_without_transit_; int64_t delta_cost_without_transit_; @@ -2630,9 +2633,11 @@ class LPCumulFilter : public IntVarLocalSearchFilter { LPCumulFilter::LPCumulFilter(const std::vector& nexts, GlobalDimensionCumulOptimizer* optimizer, + GlobalDimensionCumulOptimizer* mp_optimizer, bool filter_objective_cost) : IntVarLocalSearchFilter(nexts), optimizer_(*optimizer), + mp_optimizer_(*mp_optimizer), filter_objective_cost_(filter_objective_cost), synchronized_cost_without_transit_(-1), delta_cost_without_transit_(-1), @@ -2662,12 +2667,30 @@ bool LPCumulFilter::Accept(const Assignment* delta, if (!filter_objective_cost_) { // No need to compute the cost of the LP, only verify its feasibility. delta_cost_without_transit_ = 0; - return optimizer_.IsFeasible(next_accessor); + const DimensionSchedulingStatus status = + optimizer_.ComputeCumuls(next_accessor, nullptr, nullptr); + if (status == DimensionSchedulingStatus::OPTIMAL) return true; + if (status == DimensionSchedulingStatus::RELAXED_OPTIMAL_ONLY && + mp_optimizer_.ComputeCumuls(next_accessor, nullptr, nullptr) == + DimensionSchedulingStatus::OPTIMAL) { + return true; + } + return false; } - if (!optimizer_.ComputeCumulCostWithoutFixedTransits( - next_accessor, &delta_cost_without_transit_)) { - // Infeasible. + const DimensionSchedulingStatus status = + optimizer_.ComputeCumulCostWithoutFixedTransits( + next_accessor, &delta_cost_without_transit_); + if (status == DimensionSchedulingStatus::INFEASIBLE) { + delta_cost_without_transit_ = std::numeric_limits::max(); + return false; + } + if (delta_cost_without_transit_ > objective_max) return false; + + if (status == DimensionSchedulingStatus::RELAXED_OPTIMAL_ONLY && + mp_optimizer_.ComputeCumulCostWithoutFixedTransits( + next_accessor, &delta_cost_without_transit_) != + DimensionSchedulingStatus::OPTIMAL) { delta_cost_without_transit_ = std::numeric_limits::max(); return false; } @@ -2682,17 +2705,28 @@ void LPCumulFilter::OnSynchronize(const Assignment* delta) { // TODO(user): Try to optimize this so the LP is not called when the last // computed delta cost corresponds to the solution being synchronized. const RoutingModel& model = *optimizer_.dimension()->model(); - if (!optimizer_.ComputeCumulCostWithoutFixedTransits( - [this, &model](int64_t index) { - return IsVarSynced(index) ? Value(index) - : model.IsStart(index) ? model.End(model.VehicleIndex(index)) - : index; - }, - &synchronized_cost_without_transit_)) { + const auto& next_accessor = [this, &model](int64_t index) { + return IsVarSynced(index) ? Value(index) + : model.IsStart(index) ? model.End(model.VehicleIndex(index)) + : index; + }; + + const DimensionSchedulingStatus status = + optimizer_.ComputeCumulCostWithoutFixedTransits( + next_accessor, &synchronized_cost_without_transit_); + if (status == DimensionSchedulingStatus::INFEASIBLE) { // TODO(user): This should only happen if the LP solver times out. // DCHECK the fail wasn't due to an infeasible model. synchronized_cost_without_transit_ = 0; } + if (status == DimensionSchedulingStatus::RELAXED_OPTIMAL_ONLY && + mp_optimizer_.ComputeCumulCostWithoutFixedTransits( + next_accessor, &synchronized_cost_without_transit_) != + DimensionSchedulingStatus::OPTIMAL) { + // TODO(user): This should only happen if the MP solver times out. + // DCHECK the fail wasn't due to an infeasible model. + synchronized_cost_without_transit_ = 0; + } } int64_t LPCumulFilter::GetSynchronizedObjectiveValue() const { @@ -2702,10 +2736,11 @@ int64_t LPCumulFilter::GetSynchronizedObjectiveValue() const { } // namespace IntVarLocalSearchFilter* MakeGlobalLPCumulFilter( - GlobalDimensionCumulOptimizer* optimizer, bool filter_objective_cost) { + GlobalDimensionCumulOptimizer* optimizer, + GlobalDimensionCumulOptimizer* mp_optimizer, bool filter_objective_cost) { const RoutingModel& model = *optimizer->dimension()->model(); - return model.solver()->RevAlloc( - new LPCumulFilter(model.Nexts(), optimizer, filter_objective_cost)); + return model.solver()->RevAlloc(new LPCumulFilter( + model.Nexts(), optimizer, mp_optimizer, filter_objective_cost)); } namespace { diff --git a/ortools/constraint_solver/routing_filters.h b/ortools/constraint_solver/routing_filters.h index 6c63213d96..eeba1e508e 100644 --- a/ortools/constraint_solver/routing_filters.h +++ b/ortools/constraint_solver/routing_filters.h @@ -66,7 +66,8 @@ IntVarLocalSearchFilter* MakeCumulBoundsPropagatorFilter( /// Returns a filter checking global linear constraints and costs. IntVarLocalSearchFilter* MakeGlobalLPCumulFilter( - GlobalDimensionCumulOptimizer* optimizer, bool filter_objective_cost); + GlobalDimensionCumulOptimizer* optimizer, + GlobalDimensionCumulOptimizer* mp_optimizer, bool filter_objective_cost); /// Returns a filter checking the current solution using CP propagation. IntVarLocalSearchFilter* MakeCPFeasibilityFilter(RoutingModel* routing_model); diff --git a/ortools/constraint_solver/routing_flow.cc b/ortools/constraint_solver/routing_flow.cc index 18ff386a47..4a3e86183a 100644 --- a/ortools/constraint_solver/routing_flow.cc +++ b/ortools/constraint_solver/routing_flow.cc @@ -58,7 +58,7 @@ bool RoutingModel::IsMatchingModel() const { absl::flat_hash_set disjunction_nodes; for (DisjunctionIndex i(0); i < GetNumberOfDisjunctions(); ++i) { if (GetDisjunctionMaxCardinality(i) > 1) return false; - for (int64_t node : GetDisjunctionIndices(i)) { + for (int64_t node : GetDisjunctionNodeIndices(i)) { if (!disjunction_nodes.insert(node).second) return false; } } @@ -225,7 +225,7 @@ bool RoutingModel::SolveMatchingModel( disjunction_to_flow_nodes.back().push_back(num_flow_nodes); num_flow_nodes++; } else { - for (int n : GetDisjunctionIndices(disjunctions.back())) { + for (int n : GetDisjunctionNodeIndices(disjunctions.back())) { in_disjunction[n] = true; flow_to_non_pd[num_flow_nodes] = n; disjunction_to_flow_nodes.back().push_back(num_flow_nodes); diff --git a/ortools/constraint_solver/routing_lp_scheduling.cc b/ortools/constraint_solver/routing_lp_scheduling.cc index e6ad6d04a6..000f32ed98 100644 --- a/ortools/constraint_solver/routing_lp_scheduling.cc +++ b/ortools/constraint_solver/routing_lp_scheduling.cc @@ -171,7 +171,10 @@ LocalDimensionCumulOptimizer::LocalDimensionCumulOptimizer( case RoutingSearchParameters::GLOP: { const glop::GlopParameters parameters = GetGlopParametersForLocalLP(); for (int vehicle = 0; vehicle < vehicles; ++vehicle) { - solver_[vehicle] = absl::make_unique(parameters); + // TODO(user): Instead of passing false, detect if the relaxation + // will always violate the MIPL constraints. + solver_[vehicle] = + absl::make_unique(false, parameters); } break; } @@ -516,7 +519,7 @@ DimensionSchedulingStatus DimensionCumulOptimizerCore::OptimizeSingleRoute( return status; } -bool DimensionCumulOptimizerCore::Optimize( +DimensionSchedulingStatus DimensionCumulOptimizerCore::Optimize( const std::function& next_accessor, RoutingLinearSolverWrapper* solver, std::vector* cumul_values, std::vector* break_values, int64_t* cost, int64_t* transit_cost, @@ -532,23 +535,34 @@ bool DimensionCumulOptimizerCore::Optimize( if (propagator_ != nullptr && !propagator_->PropagateCumulBounds(next_accessor, cumul_offset)) { - return false; + return DimensionSchedulingStatus::INFEASIBLE; } int64_t total_transit_cost = 0; int64_t total_cost_offset = 0; const RoutingModel* model = dimension()->model(); + int num_needed_resources = 0; + int num_available_resources = model->vehicles(); + const auto& resource_groups = model->GetResourceGroups(); + for (int rg_index : model->GetDimensionResourceGroupIndices(dimension_)) { + num_available_resources = + std::min(num_available_resources, resource_groups[rg_index]->Size()); + } + DCHECK_GT(num_available_resources, 0); for (int vehicle = 0; vehicle < model->vehicles(); vehicle++) { int64_t route_transit_cost = 0; int64_t route_cost_offset = 0; - const bool optimize_vehicle_costs = - optimize_costs && - (!model->IsEnd(next_accessor(model->Start(vehicle))) || - model->AreEmptyRouteCostsConsideredForVehicle(vehicle)); + const bool vehicle_is_used = + !model->IsEnd(next_accessor(model->Start(vehicle))) || + model->AreEmptyRouteCostsConsideredForVehicle(vehicle); + if (vehicle_is_used && ++num_needed_resources > num_available_resources) { + return DimensionSchedulingStatus::INFEASIBLE; + } + const bool optimize_vehicle_costs = optimize_costs && vehicle_is_used; if (!SetRouteCumulConstraints(vehicle, next_accessor, cumul_offset, optimize_vehicle_costs, solver, &route_transit_cost, &route_cost_offset)) { - return false; + return DimensionSchedulingStatus::INFEASIBLE; } total_transit_cost = CapAdd(total_transit_cost, route_transit_cost); total_cost_offset = CapAdd(total_cost_offset, route_cost_offset); @@ -558,13 +572,17 @@ bool DimensionCumulOptimizerCore::Optimize( *transit_cost = total_transit_cost; } - SetGlobalConstraints(has_vehicles_being_optimized, solver); + SetGlobalConstraints(next_accessor, cumul_offset, + has_vehicles_being_optimized, solver); - if (solver->Solve(model->RemainingTime()) == - DimensionSchedulingStatus::INFEASIBLE) { - return false; + const DimensionSchedulingStatus status = + solver->Solve(model->RemainingTime()); + if (status == DimensionSchedulingStatus::INFEASIBLE) { + return status; } + // TODO(user): In case the status is RELAXED_OPTIMAL_ONLY, check we can + // safely avoid filling variable and cost values. SetValuesFromLP(index_to_cumul_variable_, cumul_offset, solver, cumul_values); SetValuesFromLP(all_break_variables_, cumul_offset, solver, break_values); @@ -575,34 +593,38 @@ bool DimensionCumulOptimizerCore::Optimize( if (clear_lp) { solver->Clear(); } - return true; + return status; } -bool DimensionCumulOptimizerCore::OptimizeAndPack( +DimensionSchedulingStatus DimensionCumulOptimizerCore::OptimizeAndPack( const std::function& next_accessor, RoutingLinearSolverWrapper* solver, std::vector* cumul_values, std::vector* break_values) { - // Note: We pass a non-nullptr cost to the Optimize() method so the costs are - // optimized by the LP. + // Note: We pass a non-nullptr cost to the Optimize() method so the costs + // are optimized by the solver. int64_t cost = 0; - if (!Optimize(next_accessor, solver, - /*cumul_values=*/nullptr, /*break_values=*/nullptr, &cost, - /*transit_cost=*/nullptr, /*clear_lp=*/false)) { - return false; + if (Optimize(next_accessor, solver, + /*cumul_values=*/nullptr, /*break_values=*/nullptr, &cost, + /*transit_cost=*/nullptr, + /*clear_lp=*/false) == DimensionSchedulingStatus::INFEASIBLE) { + return DimensionSchedulingStatus::INFEASIBLE; } - std::vector vehicles(dimension()->model()->vehicles()); std::iota(vehicles.begin(), vehicles.end(), 0); - if (PackRoutes(std::move(vehicles), solver) == - DimensionSchedulingStatus::INFEASIBLE) { - return false; + // Subtle: Even if the status was RELAXED_OPTIMAL_ONLY we try to pack just in + // case packing manages to make the solution completely feasible. + DimensionSchedulingStatus status = PackRoutes(vehicles, solver); + if (status == DimensionSchedulingStatus::INFEASIBLE) { + return status; } + // TODO(user): In case the status is RELAXED_OPTIMAL_ONLY, check we can + // safely avoid filling variable values. const int64_t global_offset = dimension_->GetGlobalOptimizerOffset(); SetValuesFromLP(index_to_cumul_variable_, global_offset, solver, cumul_values); SetValuesFromLP(all_break_variables_, global_offset, solver, break_values); solver->Clear(); - return true; + return status; } DimensionSchedulingStatus @@ -650,15 +672,7 @@ DimensionSchedulingStatus DimensionCumulOptimizerCore::PackRoutes( DCHECK(solver->SolutionIsInteger()); // Minimize the route end times without increasing the cost. - const int objective_ct = - solver->CreateNewConstraint(0, solver->GetObjectiveValue()); - - for (int variable = 0; variable < solver->NumVariables(); variable++) { - const double coefficient = solver->GetObjectiveCoefficient(variable); - if (coefficient != 0) { - solver->SetCoefficient(objective_ct, variable, coefficient); - } - } + solver->AddObjectiveConstraint(); solver->ClearObjective(); for (int vehicle : vehicles) { solver->SetObjectiveCoefficient( @@ -1286,6 +1300,7 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints( } void DimensionCumulOptimizerCore::SetGlobalConstraints( + const std::function& next_accessor, int64_t cumul_offset, bool optimize_costs, RoutingLinearSolverWrapper* solver) { // Global span cost = // global_span_cost_coefficient * (max_end_cumul - min_start_cumul). @@ -1316,6 +1331,95 @@ void DimensionCumulOptimizerCore::SetGlobalConstraints( solver->SetCoefficient(ct, second_cumul_var, 1); solver->SetCoefficient(ct, first_cumul_var, -1); } + + const RoutingModel& model = *dimension_->model(); + if (!solver->IsCPSATSolver()) { + // The resource attributes conditional constraints can only be added with + // the CP-SAT MIP solver. + return; + } + + using ResourceGroup = RoutingModel::ResourceGroup; + const auto& resource_groups = model.GetResourceGroups(); + for (int rg_index : model.GetDimensionResourceGroupIndices(dimension_)) { + // Resource domain constraints: + // Every (used) vehicle must be assigned to exactly one resource in this + // group, and each resource must be assigned to at most 1 vehicle. + // For every resource r with Attributes A = resources[r].attributes(dim) + // and every vehicle v, assign(r, v) == 1 --> + // A.start_domain.Min() <= cumul[Start(v)] <= A.start_domain.Max(), + // and + // A.end_domain.Min() <= cumul[End(v)] <= A.end_domain.Max(). + const std::vector& resources = + resource_groups[rg_index]->GetResources(); + + static const int kNoConstraint = -1; + // Assignment constraints for vehicles: each (used) vehicle must have + // exactly one resource assigned to it. + std::vector vehicle_constraints(model.vehicles(), kNoConstraint); + for (int v = 0; v < model.vehicles(); v++) { + if (model.IsEnd(next_accessor(model.Start(v))) && + !model.AreEmptyRouteCostsConsideredForVehicle(v)) { + // We don't assign a driver to unused vehicles. + continue; + } + vehicle_constraints[v] = solver->CreateNewConstraint(1, 1); + } + // Assignment constraints for resources: each resource must be assigned to + // at most one (used) vehicle. + std::vector resource_constraints(resources.size(), kNoConstraint); + for (int r = 0; r < resources.size(); r++) { + const ResourceGroup::Attributes& attributes = + resources[r].GetDimensionAttributes(dimension_); + if (attributes.start_domain().Max() < cumul_offset || + attributes.end_domain().Max() < cumul_offset) { + // This resource's domain has a cumul max lower than the offset, so it's + // not possible to restrict any vehicle start/end to this domain; skip + // it. + continue; + } + resource_constraints[r] = solver->CreateNewConstraint(0, 1); + } + // Create assignment variables, add them to the corresponding constraints, + // and create the reified constraints assign(r, v) == 1 --> + // A(r).start_domain.Min() <= cumul[Start(v)] <= A(r).start_domain.Max(), + // and + // A(r).end_domain.Min() <= cumul[End(v)] <= A(r).end_domain.Max(). + for (int r = 0; r < resources.size(); r++) { + if (resource_constraints[r] == kNoConstraint) continue; + const ResourceGroup::Attributes& attributes = + resources[r].GetDimensionAttributes(dimension_); + for (int v = 0; v < dimension_->model()->vehicles(); v++) { + if (vehicle_constraints[v] == kNoConstraint) continue; + + const int assign_r_to_v = solver->AddVariable(0, 1); + solver->SetCoefficient(vehicle_constraints[v], assign_r_to_v, 1); + solver->SetCoefficient(resource_constraints[r], assign_r_to_v, 1); + + const auto& add_domain_constraint = + [&solver, cumul_offset, assign_r_to_v](const Domain& domain, + int cumul_variable) { + if (domain == Domain::AllValues()) { + return; + } + const int64_t cumul_min = + std::max(CapSub(domain.Min(), cumul_offset), 0); + const int64_t cumul_max = + domain.Max() == std::numeric_limits::max() + ? std::numeric_limits::max() + : CapSub(domain.Max(), cumul_offset); + DCHECK_GE(cumul_max, 0); + const int cumul_constraint = solver->AddLinearConstraint( + cumul_min, cumul_max, {{cumul_variable, 1}}); + solver->SetEnforcementLiteral(cumul_constraint, assign_r_to_v); + }; + add_domain_constraint(attributes.start_domain(), + index_to_cumul_variable_[model.Start(v)]); + add_domain_constraint(attributes.end_domain(), + index_to_cumul_variable_[model.End(v)]); + } + } + } } void DimensionCumulOptimizerCore::SetValuesFromLP( @@ -1336,30 +1440,45 @@ void DimensionCumulOptimizerCore::SetValuesFromLP( } GlobalDimensionCumulOptimizer::GlobalDimensionCumulOptimizer( - const RoutingDimension* dimension) + const RoutingDimension* dimension, + RoutingSearchParameters::SchedulingSolver solver_type) : optimizer_core_(dimension, /*use_precedence_propagator=*/ !dimension->GetNodePrecedences().empty()) { - solver_ = - absl::make_unique(GetGlopParametersForGlobalLP()); + switch (solver_type) { + case RoutingSearchParameters::GLOP: { + solver_ = absl::make_unique( + /*is_relaxation=*/!dimension->model() + ->GetDimensionResourceGroupIndices(dimension) + .empty(), + GetGlopParametersForGlobalLP()); + break; + } + case RoutingSearchParameters::CP_SAT: { + solver_ = absl::make_unique(); + break; + } + default: + LOG(DFATAL) << "Unrecognized solver type: " << solver_type; + } } -bool GlobalDimensionCumulOptimizer::ComputeCumulCostWithoutFixedTransits( +DimensionSchedulingStatus +GlobalDimensionCumulOptimizer::ComputeCumulCostWithoutFixedTransits( const std::function& next_accessor, int64_t* optimal_cost_without_transits) { int64_t cost = 0; int64_t transit_cost = 0; - if (optimizer_core_.Optimize(next_accessor, solver_.get(), nullptr, nullptr, - &cost, &transit_cost)) { - if (optimal_cost_without_transits != nullptr) { - *optimal_cost_without_transits = CapSub(cost, transit_cost); - } - return true; + DimensionSchedulingStatus status = optimizer_core_.Optimize( + next_accessor, solver_.get(), nullptr, nullptr, &cost, &transit_cost); + if (status != DimensionSchedulingStatus::INFEASIBLE && + optimal_cost_without_transits != nullptr) { + *optimal_cost_without_transits = CapSub(cost, transit_cost); } - return false; + return status; } -bool GlobalDimensionCumulOptimizer::ComputeCumuls( +DimensionSchedulingStatus GlobalDimensionCumulOptimizer::ComputeCumuls( const std::function& next_accessor, std::vector* optimal_cumuls, std::vector* optimal_breaks) { @@ -1367,13 +1486,7 @@ bool GlobalDimensionCumulOptimizer::ComputeCumuls( optimal_breaks, nullptr, nullptr); } -bool GlobalDimensionCumulOptimizer::IsFeasible( - const std::function& next_accessor) { - return optimizer_core_.Optimize(next_accessor, solver_.get(), nullptr, - nullptr, nullptr, nullptr); -} - -bool GlobalDimensionCumulOptimizer::ComputePackedCumuls( +DimensionSchedulingStatus GlobalDimensionCumulOptimizer::ComputePackedCumuls( const std::function& next_accessor, std::vector* packed_cumuls, std::vector* packed_breaks) { return optimizer_core_.OptimizeAndPack(next_accessor, solver_.get(), diff --git a/ortools/constraint_solver/routing_lp_scheduling.h b/ortools/constraint_solver/routing_lp_scheduling.h index 2f62a5a5ff..8586654d1e 100644 --- a/ortools/constraint_solver/routing_lp_scheduling.h +++ b/ortools/constraint_solver/routing_lp_scheduling.h @@ -174,6 +174,7 @@ class RoutingLinearSolverWrapper { virtual int CreateNewConstraint(int64_t lower_bound, int64_t upper_bound) = 0; virtual void SetCoefficient(int ct, int index, double coefficient) = 0; virtual bool IsCPSATSolver() = 0; + virtual void AddObjectiveConstraint() = 0; virtual void AddMaximumConstraint(int max_var, std::vector vars) = 0; virtual void AddProductConstraint(int product_var, std::vector vars) = 0; virtual void SetEnforcementLiteral(int ct, int condition) = 0; @@ -236,7 +237,8 @@ class RoutingLinearSolverWrapper { class RoutingGlopWrapper : public RoutingLinearSolverWrapper { public: - explicit RoutingGlopWrapper(const glop::GlopParameters& parameters) { + RoutingGlopWrapper(bool is_relaxation, const glop::GlopParameters& parameters) + : is_relaxation_(is_relaxation) { lp_solver_.SetParameters(parameters); linear_program_.SetMaximizationProblem(false); } @@ -307,8 +309,17 @@ class RoutingGlopWrapper : public RoutingLinearSolverWrapper { coefficient); } bool IsCPSATSolver() override { return false; } - void AddMaximumConstraint(int max_var, std::vector vars) override{}; - void AddProductConstraint(int product_var, std::vector vars) override{}; + void AddObjectiveConstraint() override { + const int ct = CreateNewConstraint(0, GetObjectiveValue()); + for (int variable = 0; variable < NumVariables(); variable++) { + const double coefficient = GetObjectiveCoefficient(variable); + if (coefficient != 0) { + SetCoefficient(ct, variable, coefficient); + } + } + } + void AddMaximumConstraint(int max_var, std::vector vars) override {} + void AddProductConstraint(int product_var, std::vector vars) override {} void SetEnforcementLiteral(int ct, int condition) override{}; DimensionSchedulingStatus Solve(absl::Duration duration_limit) override { lp_solver_.GetMutableParameters()->set_max_time_in_seconds( @@ -327,6 +338,9 @@ class RoutingGlopWrapper : public RoutingLinearSolverWrapper { linear_program_.Clear(); return DimensionSchedulingStatus::INFEASIBLE; } + if (is_relaxation_) { + return DimensionSchedulingStatus::RELAXED_OPTIMAL_ONLY; + } for (const auto& allowed_interval : allowed_intervals_) { const double value_double = GetValue(allowed_interval.first); const int64_t value = @@ -354,6 +368,7 @@ class RoutingGlopWrapper : public RoutingLinearSolverWrapper { } private: + const bool is_relaxation_; glop::LinearProgram linear_program_; glop::LPSolver lp_solver_; absl::flat_hash_map> @@ -362,7 +377,7 @@ class RoutingGlopWrapper : public RoutingLinearSolverWrapper { class RoutingCPSatWrapper : public RoutingLinearSolverWrapper { public: - RoutingCPSatWrapper() : objective_offset_(0), first_constraint_to_offset_(0) { + RoutingCPSatWrapper() : first_constraint_to_offset_(0) { parameters_.set_num_search_workers(1); // Keeping presolve but with 0 iterations; as of 11/2019 it is // significantly faster than both full presolve and no presolve. @@ -377,7 +392,6 @@ class RoutingCPSatWrapper : public RoutingLinearSolverWrapper { model_.Clear(); response_.Clear(); objective_coefficients_.clear(); - objective_offset_ = 0; variable_offset_.clear(); constraint_offset_.clear(); first_constraint_to_offset_ = 0; @@ -436,17 +450,15 @@ class RoutingCPSatWrapper : public RoutingLinearSolverWrapper { sat::CpObjectiveProto* const objective = model_.mutable_objective(); objective->add_vars(index); objective->add_coeffs(coefficient); - objective_offset_ += coefficient * variable_offset_[index]; + objective->set_offset(objective->offset() + + coefficient * variable_offset_[index]); } double GetObjectiveCoefficient(int index) const override { return (index < objective_coefficients_.size()) ? objective_coefficients_[index] : 0; } - void ClearObjective() override { - model_.mutable_objective()->Clear(); - objective_offset_ = 0; - } + void ClearObjective() override { model_.mutable_objective()->Clear(); } int NumVariables() const override { return model_.variables_size(); } int CreateNewConstraint(int64_t lower_bound, int64_t upper_bound) override { const int ct_index = model_.constraints_size(); @@ -471,6 +483,19 @@ class RoutingCPSatWrapper : public RoutingLinearSolverWrapper { CapProd(variable_offset_[index], coefficient)); } bool IsCPSATSolver() override { return true; } + void AddObjectiveConstraint() override { + const double scaling_factor = + model_.has_objective() && model_.objective().scaling_factor() != 0 + ? model_.objective().scaling_factor() + : 1; + const int ct = + CreateNewConstraint(0, response_.objective_value() / scaling_factor - + model_.objective().offset()); + const sat::CpObjectiveProto& objective = model_.objective(); + for (int i = 0; i < objective.vars_size(); ++i) { + SetCoefficient(ct, objective.vars(i), objective.coeffs(i)); + } + } void AddMaximumConstraint(int max_var, std::vector vars) override { sat::LinearArgumentProto* const ct = model_.add_constraints()->mutable_lin_max(); @@ -505,6 +530,21 @@ class RoutingCPSatWrapper : public RoutingLinearSolverWrapper { ct->set_domain(1, CapSub(ct->domain(1), constraint_offset_[ct_index])); } first_constraint_to_offset_ = constraint_offset_.size(); + + // Scale the objective + sat::CpObjectiveProto* const objective = model_.mutable_objective(); + int64_t gcd(0); + for (int64_t coeff : objective->coeffs()) { + gcd = MathUtil::GCD64(gcd, std::abs(coeff)); + } + if (gcd > 1) { + for (int i = 0; i < objective->coeffs_size(); ++i) { + objective->set_coeffs(i, objective->coeffs(i) / gcd); + } + objective->set_offset(objective->offset() / gcd); + objective->set_scaling_factor(gcd); + } + parameters_.set_max_time_in_seconds(absl::ToDoubleSeconds(duration_limit)); VLOG(2) << model_.DebugString(); if (hint_.vars_size() == model_.variables_size()) { @@ -527,8 +567,7 @@ class RoutingCPSatWrapper : public RoutingLinearSolverWrapper { return DimensionSchedulingStatus::INFEASIBLE; } int64_t GetObjectiveValue() const override { - return MathUtil::FastInt64Round(response_.objective_value() + - objective_offset_); + return MathUtil::FastInt64Round(response_.objective_value()); } double GetValue(int index) const override { return response_.solution(index) + variable_offset_[index]; @@ -540,7 +579,6 @@ class RoutingCPSatWrapper : public RoutingLinearSolverWrapper { sat::CpSolverResponse response_; sat::SatParameters parameters_; std::vector objective_coefficients_; - double objective_offset_; std::vector variable_offset_; std::vector constraint_offset_; int first_constraint_to_offset_; @@ -564,16 +602,16 @@ class DimensionCumulOptimizerCore { std::vector* break_values, int64_t* cost, int64_t* transit_cost, bool clear_lp = true); - bool Optimize(const std::function& next_accessor, - RoutingLinearSolverWrapper* solver, - std::vector* cumul_values, - std::vector* break_values, int64_t* cost, - int64_t* transit_cost, bool clear_lp = true); + DimensionSchedulingStatus Optimize( + const std::function& next_accessor, + RoutingLinearSolverWrapper* solver, std::vector* cumul_values, + std::vector* break_values, int64_t* cost, int64_t* transit_cost, + bool clear_lp = true); - bool OptimizeAndPack(const std::function& next_accessor, - RoutingLinearSolverWrapper* solver, - std::vector* cumul_values, - std::vector* break_values); + DimensionSchedulingStatus OptimizeAndPack( + const std::function& next_accessor, + RoutingLinearSolverWrapper* solver, std::vector* cumul_values, + std::vector* break_values); DimensionSchedulingStatus OptimizeAndPackSingleRoute( int vehicle, const std::function& next_accessor, @@ -613,8 +651,10 @@ class DimensionCumulOptimizerCore { // SetRouteCumulConstraints() has been called on all routes, so that // index_to_cumul_variable_ and min_start/max_end_cumul_ are correctly // initialized. - void SetGlobalConstraints(bool optimize_costs, - RoutingLinearSolverWrapper* solver); + void SetGlobalConstraints( + const std::function& next_accessor, + int64_t cumul_offset, bool optimize_costs, + RoutingLinearSolverWrapper* solver); void SetValuesFromLP(const std::vector& lp_variables, int64_t offset, RoutingLinearSolverWrapper* solver, @@ -709,13 +749,15 @@ class LocalDimensionCumulOptimizer { class GlobalDimensionCumulOptimizer { public: - explicit GlobalDimensionCumulOptimizer(const RoutingDimension* dimension); + GlobalDimensionCumulOptimizer( + const RoutingDimension* dimension, + RoutingSearchParameters::SchedulingSolver solver_type); // If feasible, computes the optimal cost of the entire model with regards to // the optimizer_core_'s dimension costs, minimizing cumul soft lower/upper // bound costs and vehicle/global span costs, and stores it in "optimal_cost" // (if not null). // Returns true iff all the constraints can be respected. - bool ComputeCumulCostWithoutFixedTransits( + DimensionSchedulingStatus ComputeCumulCostWithoutFixedTransits( const std::function& next_accessor, int64_t* optimal_cost_without_transits); // If feasible, computes the optimal values for cumul and break variables, @@ -723,20 +765,17 @@ class GlobalDimensionCumulOptimizer { // costs, stores them in "optimal_cumuls" (if not null) and optimal breaks, // and returns true. // Returns false if the routes are not feasible. - bool ComputeCumuls(const std::function& next_accessor, - std::vector* optimal_cumuls, - std::vector* optimal_breaks); - - // Returns true iff the routes resulting from the next_accessor are feasible - // wrt the constraints on the optimizer_core_.dimension()'s cumuls. - bool IsFeasible(const std::function& next_accessor); + DimensionSchedulingStatus ComputeCumuls( + const std::function& next_accessor, + std::vector* optimal_cumuls, + std::vector* optimal_breaks); // Similar to ComputeCumuls, but also tries to pack the cumul values on all // routes, such that the cost remains the same, the cumuls of route ends are // minimized, and then the cumuls of the starts of the routes are maximized. - bool ComputePackedCumuls(const std::function& next_accessor, - std::vector* packed_cumuls, - std::vector* packed_breaks); + DimensionSchedulingStatus ComputePackedCumuls( + const std::function& next_accessor, + std::vector* packed_cumuls, std::vector* packed_breaks); const RoutingDimension* dimension() const { return optimizer_core_.dimension(); diff --git a/ortools/constraint_solver/routing_neighborhoods.cc b/ortools/constraint_solver/routing_neighborhoods.cc index f4a617dcf6..defb1d9980 100644 --- a/ortools/constraint_solver/routing_neighborhoods.cc +++ b/ortools/constraint_solver/routing_neighborhoods.cc @@ -677,17 +677,24 @@ FilteredHeuristicLocalSearchOperator::FilteredHeuristicLocalSearchOperator( bool keep_inverse_values) : IntVarLocalSearchOperator(heuristic->model()->Nexts(), keep_inverse_values), - model_(*heuristic->model()), - removed_nodes_(model_.Size()), + model_(heuristic->model()), + removed_nodes_(model_->Size()), heuristic_(std::move(heuristic)), - consider_vehicle_vars_(!model_.CostsAreHomogeneousAcrossVehicles()) { + consider_vehicle_vars_(!model_->CostsAreHomogeneousAcrossVehicles()) { if (consider_vehicle_vars_) { - AddVars(model_.VehicleVars()); + AddVars(model_->VehicleVars()); } } bool FilteredHeuristicLocalSearchOperator::MakeOneNeighbor() { while (IncrementPosition()) { + if (model_->CheckLimit()) { + // NOTE: Even though the limit is checked in the BuildSolutionFromRoutes() + // method of the heuristics, we still check it here to avoid calling + // IncrementPosition() and building a solution for every possible position + // if the time limit is reached. + return false; + } // NOTE: No need to call RevertChanges() here as MakeChangeAndInsertNodes() // will always return true if any change was made. if (MakeChangesAndInsertNodes()) { @@ -715,13 +722,13 @@ bool FilteredHeuristicLocalSearchOperator::MakeChangesAndInsertNodes() { bool has_change = false; const std::vector& elements = result_assignment->IntVarContainer().elements(); - for (int vehicle = 0; vehicle < model_.vehicles(); vehicle++) { - int64_t node_index = model_.Start(vehicle); - while (!model_.IsEnd(node_index)) { + for (int vehicle = 0; vehicle < model_->vehicles(); vehicle++) { + int64_t node_index = model_->Start(vehicle); + while (!model_->IsEnd(node_index)) { // NOTE: When building the solution in the heuristic, Next vars are added // to the assignment at the position corresponding to their index. const IntVarElement& node_element = elements[node_index]; - DCHECK_EQ(node_element.Var(), model_.NextVar(node_index)); + DCHECK_EQ(node_element.Var(), model_->NextVar(node_index)); const int64_t new_node_value = node_element.Value(); DCHECK_NE(new_node_value, node_index); @@ -742,7 +749,7 @@ bool FilteredHeuristicLocalSearchOperator::MakeChangesAndInsertNodes() { // the heuristic. for (int64_t node : removed_nodes_.PositionsSetAtLeastOnce()) { const IntVarElement& node_element = elements[node]; - DCHECK_EQ(node_element.Var(), model_.NextVar(node)); + DCHECK_EQ(node_element.Var(), model_->NextVar(node)); if (node_element.Value() == node) { DCHECK_NE(OldValue(node), node); has_change = true; @@ -786,11 +793,11 @@ bool FilteredHeuristicPathLNSOperator::IncrementPosition() { } bool FilteredHeuristicPathLNSOperator::CurrentRouteIsEmpty() const { - return model_.IsEnd(OldValue(model_.Start(current_route_))); + return model_->IsEnd(OldValue(model_->Start(current_route_))); } void FilteredHeuristicPathLNSOperator::IncrementCurrentRouteToNextNonEmpty() { - const int num_routes = model_.vehicles(); + const int num_routes = model_->vehicles(); do { ++current_route_ %= num_routes; if (current_route_ == last_route_) { @@ -802,8 +809,8 @@ void FilteredHeuristicPathLNSOperator::IncrementCurrentRouteToNextNonEmpty() { std::function FilteredHeuristicPathLNSOperator::SetupNextAccessorForNeighbor() { - const int64_t start_node = model_.Start(current_route_); - const int64_t end_node = model_.End(current_route_); + const int64_t start_node = model_->Start(current_route_); + const int64_t end_node = model_->End(current_route_); int64_t node = Value(start_node); while (node != end_node) { @@ -829,30 +836,30 @@ RelocatePathAndHeuristicInsertUnperformedOperator:: void RelocatePathAndHeuristicInsertUnperformedOperator::OnStart() { has_unperformed_nodes_ = false; - last_node_on_route_.resize(model_.vehicles()); + last_node_on_route_.resize(model_->vehicles()); routes_to_relocate_.clear(); empty_routes_.clear(); std::vector empty_vehicle_of_vehicle_class_added( - model_.GetVehicleClassesCount(), false); - for (int64_t node = 0; node < model_.Size(); node++) { + model_->GetVehicleClassesCount(), false); + for (int64_t node = 0; node < model_->Size(); node++) { const int64_t next = OldValue(node); if (next == node) { has_unperformed_nodes_ = true; continue; } - if (model_.IsEnd(next)) { - last_node_on_route_[model_.VehicleIndex(next)] = node; + if (model_->IsEnd(next)) { + last_node_on_route_[model_->VehicleIndex(next)] = node; } } - for (int vehicle = 0; vehicle < model_.vehicles(); vehicle++) { - const int64_t next = OldValue(model_.Start(vehicle)); - if (!model_.IsEnd(next)) { + for (int vehicle = 0; vehicle < model_->vehicles(); vehicle++) { + const int64_t next = OldValue(model_->Start(vehicle)); + if (!model_->IsEnd(next)) { routes_to_relocate_.push_back(vehicle); continue; } const int vehicle_class = - model_.GetVehicleClassIndexOfVehicle(vehicle).value(); + model_->GetVehicleClassIndexOfVehicle(vehicle).value(); if (!empty_vehicle_of_vehicle_class_added[vehicle_class]) { empty_routes_.push_back(vehicle); empty_vehicle_of_vehicle_class_added[vehicle_class] = true; @@ -897,19 +904,19 @@ RelocatePathAndHeuristicInsertUnperformedOperator:: SetupNextAccessorForNeighbor() { const int empty_route = empty_routes_[empty_route_index_]; const int relocated_route = routes_to_relocate_[route_to_relocate_index_]; - if (model_.GetVehicleClassIndexOfVehicle(empty_route) == - model_.GetVehicleClassIndexOfVehicle(relocated_route)) { + if (model_->GetVehicleClassIndexOfVehicle(empty_route) == + model_->GetVehicleClassIndexOfVehicle(relocated_route)) { // Don't try to relocate the route to an empty vehicle of the same class. return nullptr; } - const int64_t empty_start_node = model_.Start(empty_route); - const int64_t empty_end_node = model_.End(empty_route); + const int64_t empty_start_node = model_->Start(empty_route); + const int64_t empty_end_node = model_->End(empty_route); - const int64_t relocated_route_start = model_.Start(relocated_route); + const int64_t relocated_route_start = model_->Start(relocated_route); const int64_t first_relocated_node = OldValue(relocated_route_start); const int64_t last_relocated_node = last_node_on_route_[relocated_route]; - const int64_t relocated_route_end = model_.End(relocated_route); + const int64_t relocated_route_end = model_->End(relocated_route); return [this, empty_start_node, empty_end_node, first_relocated_node, last_relocated_node, relocated_route_start, @@ -927,33 +934,33 @@ FilteredHeuristicCloseNodesLNSOperator::FilteredHeuristicCloseNodesLNSOperator( std::unique_ptr heuristic, int num_close_nodes) : FilteredHeuristicLocalSearchOperator(std::move(heuristic), /*keep_inverse_values*/ true), - pickup_delivery_pairs_(model_.GetPickupAndDeliveryPairs()), + pickup_delivery_pairs_(model_->GetPickupAndDeliveryPairs()), current_node_(0), last_node_(0), - close_nodes_(model_.Size()), - new_nexts_(model_.Size()), - changed_nexts_(model_.Size()), - new_prevs_(model_.Size()), - changed_prevs_(model_.Size()) { - const int64_t size = model_.Size(); + close_nodes_(model_->Size()), + new_nexts_(model_->Size()), + changed_nexts_(model_->Size()), + new_prevs_(model_->Size()), + changed_prevs_(model_->Size()) { + const int64_t size = model_->Size(); const int64_t max_num_neighbors = - std::max(0, size - 1 - model_.vehicles()); + std::max(0, size - 1 - model_->vehicles()); const int64_t num_closest_neighbors = std::min(num_close_nodes, max_num_neighbors); DCHECK_GE(num_closest_neighbors, 0); if (num_closest_neighbors == 0) return; - const int64_t num_cost_classes = model_.GetCostClassesCount(); + const int64_t num_cost_classes = model_->GetCostClassesCount(); for (int64_t node = 0; node < size; node++) { - if (model_.IsStart(node) || model_.IsEnd(node)) continue; + if (model_->IsStart(node) || model_->IsEnd(node)) continue; std::vector> costed_after_nodes; costed_after_nodes.reserve(size); for (int64_t after_node = 0; after_node < size; after_node++) { - if (model_.IsStart(after_node) || model_.IsEnd(after_node) || + if (model_->IsStart(after_node) || model_->IsEnd(after_node) || after_node == node) { continue; } @@ -961,7 +968,7 @@ FilteredHeuristicCloseNodesLNSOperator::FilteredHeuristicCloseNodesLNSOperator( // NOTE: We don't consider the 'always-zero' cost class when searching for // closest neighbors. for (int cost_class = 1; cost_class < num_cost_classes; cost_class++) { - total_cost += model_.GetArcCostForClass(node, after_node, cost_class); + total_cost += model_->GetArcCostForClass(node, after_node, cost_class); } costed_after_nodes.emplace_back(total_cost, after_node); } @@ -987,12 +994,12 @@ bool FilteredHeuristicCloseNodesLNSOperator::IncrementPosition() { just_started_ = false; return true; } - ++current_node_ %= model_.Size(); + ++current_node_ %= model_->Size(); return current_node_ != last_node_; } void FilteredHeuristicCloseNodesLNSOperator::RemoveNode(int64_t node) { - DCHECK(!model_.IsEnd(node) && !model_.IsStart(node)); + DCHECK(!model_->IsEnd(node) && !model_->IsStart(node)); DCHECK_NE(Value(node), node); DCHECK(IsActive(node)); @@ -1001,7 +1008,7 @@ void FilteredHeuristicCloseNodesLNSOperator::RemoveNode(int64_t node) { const int64_t next = Next(node); changed_nexts_.Set(prev); new_nexts_[prev] = next; - if (next < model_.Size()) { + if (next < model_->Size()) { changed_prevs_.Set(next); new_prevs_[next] = prev; } @@ -1013,7 +1020,7 @@ void FilteredHeuristicCloseNodesLNSOperator::RemoveNodeAndActiveSibling( RemoveNode(node); for (int64_t sibling_node : GetActiveSiblings(node)) { - if (!model_.IsStart(sibling_node) && !model_.IsEnd(sibling_node)) { + if (!model_->IsStart(sibling_node) && !model_->IsEnd(sibling_node)) { RemoveNode(sibling_node); } } @@ -1026,7 +1033,7 @@ std::vector FilteredHeuristicCloseNodesLNSOperator::GetActiveSiblings( // alternative deliveries or pickups for this index pair. std::vector active_siblings; for (std::pair index_pair : - model_.GetPickupIndexPairs(node)) { + model_->GetPickupIndexPairs(node)) { for (int64_t sibling_delivery : pickup_delivery_pairs_[index_pair.first].second) { if (IsActive(sibling_delivery)) { @@ -1036,7 +1043,7 @@ std::vector FilteredHeuristicCloseNodesLNSOperator::GetActiveSiblings( } } for (std::pair index_pair : - model_.GetDeliveryIndexPairs(node)) { + model_->GetDeliveryIndexPairs(node)) { for (int64_t sibling_pickup : pickup_delivery_pairs_[index_pair.first].first) { if (IsActive(sibling_pickup)) { @@ -1050,10 +1057,10 @@ std::vector FilteredHeuristicCloseNodesLNSOperator::GetActiveSiblings( std::function FilteredHeuristicCloseNodesLNSOperator::SetupNextAccessorForNeighbor() { - if (model_.IsStart(current_node_)) { + if (model_->IsStart(current_node_)) { return nullptr; } - DCHECK(!model_.IsEnd(current_node_)); + DCHECK(!model_->IsEnd(current_node_)); changed_nexts_.SparseClearAll(); changed_prevs_.SparseClearAll(); @@ -1135,7 +1142,7 @@ FilteredHeuristicExpensiveChainLNSOperator::SetupNextAccessorForNeighbor() { } bool FilteredHeuristicExpensiveChainLNSOperator::IncrementRoute() { - ++current_route_ %= model_.vehicles(); + ++current_route_ %= model_->vehicles(); return current_route_ != last_route_; } @@ -1219,9 +1226,9 @@ bool FilteredHeuristicExpensiveChainLNSOperator:: FindMostExpensiveChainsOnRemainingRoutes() { do { if (FindMostExpensiveArcsOnRoute( - num_arcs_to_consider_, model_.Start(current_route_), + num_arcs_to_consider_, model_->Start(current_route_), [this](int64_t i) { return OldValue(i); }, - [this](int64_t node) { return model_.IsEnd(node); }, + [this](int64_t node) { return model_->IsEnd(node); }, arc_cost_for_route_start_, &most_expensive_arc_starts_and_ranks_, ¤t_expensive_arc_indices_)) { return true; diff --git a/ortools/constraint_solver/routing_neighborhoods.h b/ortools/constraint_solver/routing_neighborhoods.h index 085012ef78..948c2d0a27 100644 --- a/ortools/constraint_solver/routing_neighborhoods.h +++ b/ortools/constraint_solver/routing_neighborhoods.h @@ -385,7 +385,7 @@ class FilteredHeuristicLocalSearchOperator : public IntVarLocalSearchOperator { // TODO(user): Remove the dependency from RoutingModel by storing an // IntVarFilteredHeuristic here instead and storing information on path // start/ends like PathOperator does (instead of relying on the model). - const RoutingModel& model_; + RoutingModel* const model_; /// Keeps track of removed nodes when making a neighbor. SparseBitset<> removed_nodes_; @@ -393,7 +393,7 @@ class FilteredHeuristicLocalSearchOperator : public IntVarLocalSearchOperator { bool MakeOneNeighbor() override; bool MakeChangesAndInsertNodes(); - int64_t VehicleVarIndex(int64_t node) const { return model_.Size() + node; } + int64_t VehicleVarIndex(int64_t node) const { return model_->Size() + node; } const std::unique_ptr heuristic_; const bool consider_vehicle_vars_; @@ -528,7 +528,7 @@ class FilteredHeuristicCloseNodesLNSOperator void RemoveNodeAndActiveSibling(int64_t node); bool IsActive(int64_t node) const { - DCHECK_LT(node, model_.Size()); + DCHECK_LT(node, model_->Size()); return Value(node) != node && !removed_nodes_[node]; } @@ -538,7 +538,7 @@ class FilteredHeuristicCloseNodesLNSOperator return changed_prevs_[node] ? new_prevs_[node] : InverseValue(node); } int64_t Next(int64_t node) const { - DCHECK(!model_.IsEnd(node)); + DCHECK(!model_->IsEnd(node)); return changed_nexts_[node] ? new_nexts_[node] : Value(node); } diff --git a/ortools/constraint_solver/routing_parameters.cc b/ortools/constraint_solver/routing_parameters.cc index 7bebced587..8cec12c974 100644 --- a/ortools/constraint_solver/routing_parameters.cc +++ b/ortools/constraint_solver/routing_parameters.cc @@ -26,6 +26,7 @@ #include "ortools/constraint_solver/constraint_solver.h" #include "ortools/constraint_solver/routing_enums.pb.h" #include "ortools/constraint_solver/solver_parameters.pb.h" +#include "ortools/sat/sat_parameters.pb.h" #include "ortools/util/optional_boolean.pb.h" namespace operations_research { @@ -107,6 +108,8 @@ RoutingSearchParameters DefaultRoutingSearchParameters() { "use_depth_first_search: false " "use_cp: BOOL_TRUE " "use_cp_sat: BOOL_FALSE " + "use_generalized_cp_sat: BOOL_FALSE " + "sat_parameters { linearization_level: 2 num_search_workers: 1 } " "continuous_scheduling_solver: GLOP " "mixed_integer_scheduling_solver: CP_SAT " "optimization_step: 0.0 " @@ -370,6 +373,17 @@ std::string FindErrorInRoutingSearchParameters( } } + { + const sat::SatParameters& sat_parameters = + search_parameters.sat_parameters(); + if (sat_parameters.enumerate_all_solutions() && + (sat_parameters.num_search_workers() > 1 || + sat_parameters.interleave_search())) { + return "sat_parameters.enumerate_all_solutions cannot be true in parallel" + " search"; + } + } + return ""; // = Valid (No error). } diff --git a/ortools/constraint_solver/routing_parameters.proto b/ortools/constraint_solver/routing_parameters.proto index 33c20cb803..2fbefad860 100644 --- a/ortools/constraint_solver/routing_parameters.proto +++ b/ortools/constraint_solver/routing_parameters.proto @@ -24,6 +24,7 @@ option csharp_namespace = "Google.OrTools.ConstraintSolver"; import "google/protobuf/duration.proto"; import "ortools/constraint_solver/routing_enums.proto"; import "ortools/constraint_solver/solver_parameters.proto"; +import "ortools/sat/sat_parameters.proto"; import "ortools/util/optional_boolean.proto"; package operations_research; @@ -34,7 +35,7 @@ package operations_research; // then the routing library will pick its preferred value for that parameter // automatically: this should be the case for most parameters. // To see those "default" parameters, call GetDefaultRoutingSearchParameters(). -// Next ID: 47 +// Next ID: 49 message RoutingSearchParameters { // First solution strategies, used as starting point of local search. FirstSolutionStrategy.Value first_solution_strategy = 1; @@ -423,6 +424,14 @@ message RoutingSearchParameters { // remaining and will use the CP solution as a hint for the CP-SAT search. // As of 5/2019, only TSP models can be solved. OptionalBoolean use_cp_sat = 27; + // If true, use the CP-SAT solver to find a solution on generalized routing + // model. If use_cp is also true, the CP-SAT solver will be run after the CP + // solver if there is time remaining and will use the CP solution as a hint + // for the CP-SAT search. + OptionalBoolean use_generalized_cp_sat = 47; + // If use_cp_sat or use_generalized_cp_sat is true, contains the SAT algorithm + // parameters which will be used. + sat.SatParameters sat_parameters = 48; // Underlying solver to use in dimension scheduling, respectively for // continuous and mixed models. enum SchedulingSolver { diff --git a/ortools/constraint_solver/routing_sat.cc b/ortools/constraint_solver/routing_sat.cc index be43fb4e33..e7c30c270e 100644 --- a/ortools/constraint_solver/routing_sat.cc +++ b/ortools/constraint_solver/routing_sat.cc @@ -20,6 +20,7 @@ #include #include +#include "absl/container/flat_hash_map.h" #include "absl/time/time.h" #include "ortools/base/integral_types.h" #include "ortools/base/logging.h" @@ -27,11 +28,13 @@ #include "ortools/constraint_solver/constraint_solver.h" #include "ortools/constraint_solver/routing.h" #include "ortools/constraint_solver/routing_parameters.pb.h" +#include "ortools/constraint_solver/routing_types.h" #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/cp_model_solver.h" #include "ortools/sat/integer.h" #include "ortools/sat/model.h" #include "ortools/sat/sat_parameters.pb.h" +#include "ortools/util/optional_boolean.pb.h" #include "ortools/util/saturated_arithmetic.h" namespace operations_research { @@ -55,6 +58,34 @@ int AddVariable(CpModelProto* cp_model, int64_t lb, int64_t ub) { return index; } +// Adds a linear constraint, enforcing +// enforcement_literals -> lower_bound <= sum variable * coeff <= upper_bound. +void AddLinearConstraint( + CpModelProto* cp_model, int64_t lower_bound, int64_t upper_bound, + const std::vector>& variable_coeffs, + const std::vector& enforcement_literals) { + CHECK_LE(lower_bound, upper_bound); + ConstraintProto* ct = cp_model->add_constraints(); + for (const int enforcement_literal : enforcement_literals) { + ct->add_enforcement_literal(enforcement_literal); + } + LinearConstraintProto* arg = ct->mutable_linear(); + arg->add_domain(lower_bound); + arg->add_domain(upper_bound); + for (const auto [var, coeff] : variable_coeffs) { + arg->add_vars(var); + arg->add_coeffs(coeff); + } +} + +// Adds a linear constraint, enforcing +// lower_bound <= sum variable * coeff <= upper_bound. +void AddLinearConstraint( + CpModelProto* cp_model, int64_t lower_bound, int64_t upper_bound, + const std::vector>& variable_coeffs) { + AddLinearConstraint(cp_model, lower_bound, upper_bound, variable_coeffs, {}); +} + // Returns the unique depot node used in the CP-SAT models (as of 01/2020). int64_t GetDepotFromModel(const RoutingModel& model) { return model.Start(0); } @@ -112,15 +143,9 @@ void AddDimensions(const RoutingModel& model, const ArcVarMap& arc_vars, if (tail == head || model.IsStart(tail) || model.IsStart(head)) continue; // arc[tail][head] -> cumuls[head] >= cumuls[tail] + transit. // This is a relaxation of the model as it does not consider slack max. - ConstraintProto* ct = cp_model->add_constraints(); - ct->add_enforcement_literal(arc_var.second); - LinearConstraintProto* arg = ct->mutable_linear(); - arg->add_domain(transit(tail, head)); - arg->add_domain(std::numeric_limits::max()); - arg->add_vars(cumuls[tail]); - arg->add_coeffs(-1); - arg->add_vars(cumuls[head]); - arg->add_coeffs(1); + AddLinearConstraint( + cp_model, transit(tail, head), std::numeric_limits::max(), + {{cumuls[head], 1}, {cumuls[tail], -1}}, {arc_var.second}); } } } @@ -142,15 +167,8 @@ std::vector CreateRanks(const RoutingModel& model, const int head = arc_var.first.head; if (tail == head || head == depot) continue; // arc[tail][head] -> ranks[head] == ranks[tail] + 1. - ConstraintProto* ct = cp_model->add_constraints(); - ct->add_enforcement_literal(arc_var.second); - LinearConstraintProto* arg = ct->mutable_linear(); - arg->add_domain(1); - arg->add_domain(1); - arg->add_vars(ranks[tail]); - arg->add_coeffs(-1); - arg->add_vars(ranks[head]); - arg->add_coeffs(1); + AddLinearConstraint(cp_model, 1, 1, {{ranks[head], 1}, {ranks[tail], -1}}, + {arc_var.second}); } return ranks; } @@ -175,25 +193,14 @@ std::vector CreateVehicleVars(const RoutingModel& model, if (tail == head || head == depot) continue; if (tail == depot) { // arc[depot][head] -> vehicles[head] == head. - ConstraintProto* ct = cp_model->add_constraints(); - ct->add_enforcement_literal(arc_var.second); - LinearConstraintProto* arg = ct->mutable_linear(); - arg->add_domain(head); - arg->add_domain(head); - arg->add_vars(vehicles[head]); - arg->add_coeffs(1); + AddLinearConstraint(cp_model, head, head, {{vehicles[head], 1}}, + {arc_var.second}); continue; } // arc[tail][head] -> vehicles[head] == vehicles[tail]. - ConstraintProto* ct = cp_model->add_constraints(); - ct->add_enforcement_literal(arc_var.second); - LinearConstraintProto* arg = ct->mutable_linear(); - arg->add_domain(0); - arg->add_domain(0); - arg->add_vars(vehicles[tail]); - arg->add_coeffs(-1); - arg->add_vars(vehicles[head]); - arg->add_coeffs(1); + AddLinearConstraint(cp_model, 0, 0, + {{vehicles[head], 1}, {vehicles[tail], -1}}, + {arc_var.second}); } return vehicles; } @@ -208,28 +215,12 @@ void AddPickupDeliveryConstraints(const RoutingModel& model, for (const auto& pairs : model.GetPickupAndDeliveryPairs()) { const int64_t pickup = pairs.first[0]; const int64_t delivery = pairs.second[0]; - { - // ranks[pickup] + 1 <= ranks[delivery]. - ConstraintProto* ct = cp_model->add_constraints(); - LinearConstraintProto* arg = ct->mutable_linear(); - arg->add_domain(1); - arg->add_domain(std::numeric_limits::max()); - arg->add_vars(ranks[delivery]); - arg->add_coeffs(1); - arg->add_vars(ranks[pickup]); - arg->add_coeffs(-1); - } - { - // vehicles[pickup] == vehicles[delivery] - ConstraintProto* ct = cp_model->add_constraints(); - LinearConstraintProto* arg = ct->mutable_linear(); - arg->add_domain(0); - arg->add_domain(0); - arg->add_vars(vehicles[delivery]); - arg->add_coeffs(1); - arg->add_vars(vehicles[pickup]); - arg->add_coeffs(-1); - } + // ranks[pickup] + 1 <= ranks[delivery]. + AddLinearConstraint(cp_model, 1, std::numeric_limits::max(), + {{ranks[delivery], 1}, {ranks[pickup], -1}}); + // vehicles[pickup] == vehicles[delivery] + AddLinearConstraint(cp_model, 0, 0, + {{vehicles[delivery], 1}, {vehicles[pickup], -1}}); } } @@ -257,7 +248,7 @@ ArcVarMap PopulateMultiRouteModelFromRoutingModel(const RoutingModel& model, // as the CP solver will reject those. if (model.IsStart(head)) continue; const int head_index = model.IsEnd(head) ? depot : head; - if (head_index == tail_index) continue; + if (head_index == tail_index && head_index == depot) continue; const int64_t cost = tail != head ? model.GetHomogeneousCost(tail, head) : model.UnperformedPenalty(tail); if (cost == std::numeric_limits::max()) continue; @@ -270,80 +261,19 @@ ArcVarMap PopulateMultiRouteModelFromRoutingModel(const RoutingModel& model, } } - // The following flow constraints seem to be necessary with the Route - // constraint, greatly improving performance due to stronger LP relaxation - // (supposedly). - // TODO(user): Remove these constraints when the Route constraint handles - // LP relaxations properly. + // Limit the number of routes to the maximum number of vehicles. { - LinearConstraintProto* ct = cp_model->add_constraints()->mutable_linear(); - ct->add_domain(0); - ct->add_domain(0); - for (int node = 0; node < num_nodes; ++node) { - if (model.IsStart(node) || model.IsEnd(node)) continue; - int* const depot_node_var = gtl::FindOrNull(arc_vars, {depot, node}); - if (depot_node_var == nullptr) continue; - ct->add_vars(*depot_node_var); - ct->add_coeffs(1); - int* const node_depot_var = gtl::FindOrNull(arc_vars, {node, depot}); - if (node_depot_var == nullptr) continue; - ct->add_vars(*node_depot_var); - ct->add_coeffs(-1); - } - } - - { - LinearConstraintProto* ct = cp_model->add_constraints()->mutable_linear(); - ct->add_domain(0); - // Taking the min since as of 04/2020 fleet is homogeneous. - ct->add_domain( - std::min(model.vehicles(), model.GetMaximumNumberOfActiveVehicles())); + std::vector> variable_coeffs; for (int node = 0; node < num_nodes; ++node) { if (model.IsStart(node) || model.IsEnd(node)) continue; int* const var = gtl::FindOrNull(arc_vars, {depot, node}); if (var == nullptr) continue; - ct->add_vars(*var); - ct->add_coeffs(1); - } - } - - for (int tail = 0; tail < num_nodes; ++tail) { - if (model.IsStart(tail) || model.IsEnd(tail)) continue; - LinearConstraintProto* ct = cp_model->add_constraints()->mutable_linear(); - ct->add_domain(1); - ct->add_domain(1); - std::unique_ptr iter( - model.NextVar(tail)->MakeDomainIterator(false)); - bool depot_added = false; - for (int head : InitAndGetValues(iter.get())) { - if (model.IsStart(head)) continue; - if (tail == head) continue; - if (model.IsEnd(head)) { - if (depot_added) continue; - head = depot; - depot_added = true; - } - int* const var = gtl::FindOrNull(arc_vars, {tail, head}); - if (var == nullptr) continue; - ct->add_vars(*var); - ct->add_coeffs(1); - } - } - - for (int head = 0; head < num_nodes; ++head) { - if (model.IsStart(head) || model.IsEnd(head)) continue; - LinearConstraintProto* ct = cp_model->add_constraints()->mutable_linear(); - ct->add_domain(1); - ct->add_domain(1); - for (int tail = 0; tail < num_nodes; ++tail) { - if (model.IsEnd(head)) continue; - if (tail == head) continue; - if (model.IsStart(tail) && tail != depot) continue; - int* const var = gtl::FindOrNull(arc_vars, {tail, head}); - if (var == nullptr) continue; - ct->add_vars(*var); - ct->add_coeffs(1); + variable_coeffs.push_back({*var, 1}); } + AddLinearConstraint( + cp_model, 0, + std::min(model.vehicles(), model.GetMaximumNumberOfActiveVehicles()), + variable_coeffs); } AddPickupDeliveryConstraints(model, arc_vars, cp_model); @@ -471,6 +401,586 @@ bool ConvertToSolution(const CpSolverResponse& response, return true; } +// Adds dimensions to a CpModelProto for heterogeneous fleet. Adds path +// cumul constraints and cumul bounds. +void AddGeneralizedDimensions( + const RoutingModel& model, const ArcVarMap& arc_vars, + const std::vector>& vehicle_performs_node, + const std::vector>& + vehicle_class_performs_arc, + CpModelProto* cp_model) { + const int num_cp_nodes = model.Nexts().size() + model.vehicles() + 1; + for (const RoutingDimension* dimension : model.GetDimensions()) { + // Initialize cumuls. + std::vector cumuls(num_cp_nodes, -1); + for (int cp_node = 1; cp_node < num_cp_nodes; ++cp_node) { + const int node = cp_node - 1; + int64_t cumul_min = dimension->cumuls()[node]->Min(); + int64_t cumul_max = dimension->cumuls()[node]->Max(); + if (model.IsStart(node) || model.IsEnd(node)) { + const int vehicle = model.VehicleIndex(node); + cumul_max = + std::min(cumul_max, dimension->vehicle_capacities()[vehicle]); + } + cumuls[cp_node] = AddVariable(cp_model, cumul_min, cumul_max); + } + + // Constrain cumuls with vehicle capacities. + for (int vehicle = 0; vehicle < model.vehicles(); vehicle++) { + for (int cp_node = 1; cp_node < num_cp_nodes; cp_node++) { + if (!vehicle_performs_node[vehicle].contains(cp_node)) continue; + const int64_t vehicle_capacity = + dimension->vehicle_capacities()[vehicle]; + AddLinearConstraint(cp_model, std::numeric_limits::min(), + vehicle_capacity, {{cumuls[cp_node], 1}}, + {vehicle_performs_node[vehicle].at(cp_node)}); + } + } + + for (auto vehicle_class = RoutingVehicleClassIndex(0); + vehicle_class < model.GetVehicleClassesCount(); vehicle_class++) { + std::vector slack(num_cp_nodes, -1); + const int64_t span_cost = + dimension->GetSpanCostCoefficientForVehicleClass(vehicle_class); + for (const auto [arc, arc_var] : arc_vars) { + const auto [cp_tail, cp_head] = arc; + if (cp_tail == cp_head || cp_tail == 0 || cp_head == 0) continue; + if (!vehicle_class_performs_arc[vehicle_class.value()].contains( + arc_var)) { + continue; + } + // Create slack variable and add span cost to the objective. + if (slack[cp_tail] == -1) { + const int64_t slack_max = + cp_tail - 1 < dimension->slacks().size() + ? dimension->slacks()[cp_tail - 1]->Max() + : 0; + slack[cp_tail] = AddVariable(cp_model, 0, slack_max); + if (slack_max > 0 && span_cost > 0) { + cp_model->mutable_objective()->add_vars(slack[cp_tail]); + cp_model->mutable_objective()->add_coeffs(span_cost); + } + } + const int64_t transit = dimension->class_transit_evaluator( + vehicle_class)(cp_tail - 1, cp_head - 1); + // vehicle_class_performs_arc[vehicle][arc_var] = 1 -> + // cumuls[cp_head] - cumuls[cp_tail] - slack[cp_tail] = transit + AddLinearConstraint( + cp_model, transit, transit, + {{cumuls[cp_head], 1}, {cumuls[cp_tail], -1}, {slack[cp_tail], -1}}, + {vehicle_class_performs_arc[vehicle_class.value()].at(arc_var)}); + } + } + + // Constrain cumuls with span limits. + for (int vehicle = 0; vehicle < model.vehicles(); vehicle++) { + const int64_t span_limit = + dimension->vehicle_span_upper_bounds()[vehicle]; + if (span_limit == std::numeric_limits::max()) continue; + int cp_start = model.Start(vehicle) + 1; + int cp_end = model.End(vehicle) + 1; + AddLinearConstraint(cp_model, std::numeric_limits::min(), + span_limit, + {{cumuls[cp_end], 1}, {cumuls[cp_start], -1}}); + } + + // Set soft span upper bound costs. + if (dimension->HasSoftSpanUpperBounds()) { + for (int vehicle = 0; vehicle < model.vehicles(); vehicle++) { + const auto [bound, cost] = + dimension->GetSoftSpanUpperBoundForVehicle(vehicle); + const int cp_start = model.Start(vehicle) + 1; + const int cp_end = model.End(vehicle) + 1; + const int extra = + AddVariable(cp_model, 0, + std::min(dimension->cumuls()[model.End(vehicle)]->Max(), + dimension->vehicle_capacities()[vehicle])); + // -inf <= cumuls[cp_end] - cumuls[cp_start] - extra <= bound + AddLinearConstraint( + cp_model, std::numeric_limits::min(), bound, + {{cumuls[cp_end], 1}, {cumuls[cp_start], -1}, {extra, -1}}); + // Add extra * cost to objective. + cp_model->mutable_objective()->add_vars(extra); + cp_model->mutable_objective()->add_coeffs(cost); + } + } + } +} + +std::vector CreateGeneralizedRanks(const RoutingModel& model, + const ArcVarMap& arc_vars, + const std::vector& is_unperformed, + CpModelProto* cp_model) { + const int depot = 0; + const int num_cp_nodes = model.Nexts().size() + model.vehicles() + 1; + // Maximum length of a single route (excluding the depot & vehicle end nodes). + const int max_rank = num_cp_nodes - 2 * model.vehicles(); + std::vector ranks(num_cp_nodes, -1); + ranks[depot] = AddVariable(cp_model, 0, 0); + for (int cp_node = 1; cp_node < num_cp_nodes; cp_node++) { + if (model.IsEnd(cp_node - 1)) continue; + ranks[cp_node] = AddVariable(cp_model, 0, max_rank); + // For unperformed nodes rank is 0. + AddLinearConstraint(cp_model, 0, 0, {{ranks[cp_node], 1}}, + {is_unperformed[cp_node]}); + } + for (const auto [arc, arc_var] : arc_vars) { + const auto [cp_tail, cp_head] = arc; + if (model.IsEnd(cp_head - 1)) continue; + if (cp_tail == cp_head || cp_head == depot) continue; + // arc[tail][head] -> ranks[head] == ranks[tail] + 1. + AddLinearConstraint(cp_model, 1, 1, + {{ranks[cp_head], 1}, {ranks[cp_tail], -1}}, {arc_var}); + } + return ranks; +} + +void AddGeneralizedPickupDeliveryConstraints( + const RoutingModel& model, const ArcVarMap& arc_vars, + const std::vector>& vehicle_performs_node, + const std::vector& is_unperformed, CpModelProto* cp_model) { + if (model.GetPickupAndDeliveryPairs().empty()) return; + const std::vector ranks = + CreateGeneralizedRanks(model, arc_vars, is_unperformed, cp_model); + for (const auto& pairs : model.GetPickupAndDeliveryPairs()) { + for (const int delivery : pairs.second) { + const int cp_delivery = delivery + 1; + for (int vehicle = 0; vehicle < model.vehicles(); vehicle++) { + const Arc vehicle_start_delivery_arc = { + static_cast(model.Start(vehicle) + 1), cp_delivery}; + if (gtl::ContainsKey(arc_vars, vehicle_start_delivery_arc)) { + // Forbid vehicle_start -> delivery arc. + AddLinearConstraint(cp_model, 0, 0, + {{arc_vars.at(vehicle_start_delivery_arc), 1}}); + } + } + + for (const int pickup : pairs.first) { + const int cp_pickup = pickup + 1; + const Arc delivery_pickup_arc = {cp_delivery, cp_pickup}; + if (gtl::ContainsKey(arc_vars, delivery_pickup_arc)) { + // Forbid delivery -> pickup arc. + AddLinearConstraint(cp_model, 0, 0, + {{arc_vars.at(delivery_pickup_arc), 1}}); + } + + DCHECK_GE(is_unperformed[cp_delivery], 0); + DCHECK_GE(is_unperformed[cp_pickup], 0); + // A negative index i refers to NOT the literal at index -i - 1. + // -i - 1 ~ NOT i, if value of i in [0, 1] (boolean). + const int delivery_performed = -is_unperformed[cp_delivery] - 1; + const int pickup_performed = -is_unperformed[cp_pickup] - 1; + // The same vehicle performs pickup and delivery. + for (int vehicle = 0; vehicle < model.vehicles(); vehicle++) { + // delivery_performed & pickup_performed -> + // vehicle_performs_node[vehicle][cp_delivery] - + // vehicle_performs_node[vehicle][cp_pickup] = 0 + AddLinearConstraint( + cp_model, 0, 0, + {{vehicle_performs_node[vehicle].at(cp_delivery), 1}, + {vehicle_performs_node[vehicle].at(cp_pickup), -1}}, + {delivery_performed, pickup_performed}); + } + } + } + + std::vector> ranks_difference; + // -SUM(pickup)ranks[pickup]. + for (const int pickup : pairs.first) { + const int cp_pickup = pickup + 1; + ranks_difference.push_back({ranks[cp_pickup], -1}); + } + // SUM(delivery)ranks[delivery]. + for (const int delivery : pairs.second) { + const int cp_delivery = delivery + 1; + ranks_difference.push_back({ranks[cp_delivery], 1}); + } + // SUM(delivery)ranks[delivery] - SUM(pickup)ranks[pickup] >= 1 + AddLinearConstraint(cp_model, 1, std::numeric_limits::max(), + ranks_difference); + } +} + +// Converts a RoutingModel to CpModelProto for models with multiple +// vehicles. The node 0 is depot. All nodes in CpModel have index increased +// by 1 in comparison to the RoutingModel. Each start node has only 1 +// incoming arc (from depot), each end node has only 1 outgoing arc (to +// depot). The mapping from CPModelProto arcs to their corresponding arc +// variable is returned. +ArcVarMap PopulateGeneralizedRouteModelFromRoutingModel( + const RoutingModel& model, CpModelProto* cp_model) { + ArcVarMap arc_vars; + const int depot = 0; + const int num_nodes = model.Nexts().size(); + const int num_cp_nodes = num_nodes + model.vehicles() + 1; + // vehicle_performs_node[vehicle][node] equals to 1 if the vehicle performs + // the node, and 0 otherwise. + std::vector> vehicle_performs_node( + model.vehicles()); + // Connect vehicles start and end nodes to depot. + for (int vehicle = 0; vehicle < model.vehicles(); vehicle++) { + const int cp_start = model.Start(vehicle) + 1; + const Arc start_arc = {depot, cp_start}; + const int start_arc_var = AddVariable(cp_model, 1, 1); + DCHECK(!gtl::ContainsKey(arc_vars, start_arc)); + arc_vars.insert({start_arc, start_arc_var}); + + const int cp_end = model.End(vehicle) + 1; + const Arc end_arc = {cp_end, depot}; + const int end_arc_var = AddVariable(cp_model, 1, 1); + DCHECK(!gtl::ContainsKey(arc_vars, end_arc)); + arc_vars.insert({end_arc, end_arc_var}); + + vehicle_performs_node[vehicle][cp_start] = start_arc_var; + vehicle_performs_node[vehicle][cp_end] = end_arc_var; + } + + // is_unperformed[node] variable equals to 1 if visit is unperformed, and 0 + // otherwise. + std::vector is_unperformed(num_cp_nodes, -1); + // Initialize is_unperformed variables for nodes that must be performed. + for (int node = 0; node < num_nodes; node++) { + const int cp_node = node + 1; + // Forced active and nodes that are not involved in any disjunctions are + // always performed. + const std::vector& disjunction_indices = + model.GetDisjunctionIndices(node); + if (disjunction_indices.empty() || model.ActiveVar(node)->Min() == 1) { + is_unperformed[cp_node] = AddVariable(cp_model, 0, 0); + continue; + } + // Check if the node is in a forced active disjunction. + for (RoutingDisjunctionIndex disjunction_index : disjunction_indices) { + const int num_nodes = + model.GetDisjunctionNodeIndices(disjunction_index).size(); + const int64_t penalty = model.GetDisjunctionPenalty(disjunction_index); + const int64_t max_cardinality = + model.GetDisjunctionMaxCardinality(disjunction_index); + if (num_nodes == max_cardinality && + (penalty < 0 || penalty == std::numeric_limits::max())) { + // Nodes in this disjunction are forced active. + is_unperformed[cp_node] = AddVariable(cp_model, 0, 0); + break; + } + } + } + // Add alternative visits. Create self-looped arc variables. Set penalty for + // not performing disjunctions. + for (RoutingDisjunctionIndex disjunction_index(0); + disjunction_index < model.GetNumberOfDisjunctions(); + disjunction_index++) { + const std::vector& disjunction_indices = + model.GetDisjunctionNodeIndices(disjunction_index); + const int disjunction_size = disjunction_indices.size(); + const int64_t penalty = model.GetDisjunctionPenalty(disjunction_index); + const int64_t max_cardinality = + model.GetDisjunctionMaxCardinality(disjunction_index); + // Case when disjunction involves only 1 node, the node is only present in + // this disjunction, and the node can be unperformed. + if (disjunction_size == 1 && + model.GetDisjunctionIndices(disjunction_indices[0]).size() == 1 && + is_unperformed[disjunction_indices[0] + 1] == -1) { + const int cp_node = disjunction_indices[0] + 1; + const Arc arc = {cp_node, cp_node}; + DCHECK(!gtl::ContainsKey(arc_vars, arc)); + is_unperformed[cp_node] = AddVariable(cp_model, 0, 1); + arc_vars.insert({arc, is_unperformed[cp_node]}); + cp_model->mutable_objective()->add_vars(is_unperformed[cp_node]); + cp_model->mutable_objective()->add_coeffs(penalty); + continue; + } + // num_performed + SUM(node)is_unperformed[node] = disjunction_size + const int num_performed = AddVariable(cp_model, 0, max_cardinality); + std::vector> var_coeffs; + var_coeffs.push_back({num_performed, 1}); + for (const int node : disjunction_indices) { + const int cp_node = node + 1; + // Node can be unperformed. + if (is_unperformed[cp_node] == -1) { + const Arc arc = {cp_node, cp_node}; + DCHECK(!gtl::ContainsKey(arc_vars, arc)); + is_unperformed[cp_node] = AddVariable(cp_model, 0, 1); + arc_vars.insert({arc, is_unperformed[cp_node]}); + } + var_coeffs.push_back({is_unperformed[cp_node], 1}); + } + AddLinearConstraint(cp_model, disjunction_size, disjunction_size, + var_coeffs); + // When penalty is negative or max int64_t (forced active), num_violated is + // 0. + if (penalty < 0 || penalty == std::numeric_limits::max()) { + AddLinearConstraint(cp_model, max_cardinality, max_cardinality, + {{num_performed, 1}}); + continue; + } + // If number of active indices is less than max_cardinality, then for each + // violated index 'penalty' is paid. + const int num_violated = AddVariable(cp_model, 0, max_cardinality); + cp_model->mutable_objective()->add_vars(num_violated); + cp_model->mutable_objective()->add_coeffs(penalty); + // num_performed + num_violated = max_cardinality + AddLinearConstraint(cp_model, max_cardinality, max_cardinality, + {{num_performed, 1}, {num_violated, 1}}); + } + // Create "arc" variables. + for (int tail = 0; tail < num_nodes; ++tail) { + const int cp_tail = tail + 1; + std::unique_ptr iter( + model.NextVar(tail)->MakeDomainIterator(false)); + for (int head : InitAndGetValues(iter.get())) { + const int cp_head = head + 1; + if (model.IsStart(head)) continue; + // Arcs for unperformed visits have already been created. + if (tail == head) continue; + // Direct arcs from start to end nodes should exist only if they are + // for the same vehicle. + if (model.IsStart(tail) && model.IsEnd(head) && + model.VehicleIndex(tail) != model.VehicleIndex(head)) { + continue; + } + + bool feasible = false; + for (int vehicle = 0; vehicle < model.vehicles(); vehicle++) { + if (model.GetArcCostForVehicle(tail, head, vehicle) != + std::numeric_limits::max()) { + feasible = true; + break; + } + } + if (!feasible) continue; + + const Arc arc = {cp_tail, cp_head}; + DCHECK(!gtl::ContainsKey(arc_vars, arc)); + const int arc_var = AddVariable(cp_model, 0, 1); + arc_vars.insert({arc, arc_var}); + } + } + + // Set literals for vehicle performing node. + for (int cp_node = 1; cp_node < num_cp_nodes; cp_node++) { + // For starts and ends nodes vehicle_performs_node variables already set. + if (model.IsStart(cp_node - 1) || model.IsEnd(cp_node - 1)) continue; + // Each node should be performed by 1 vehicle, or be unperformed. + // SUM(vehicle)(vehicle_performs_node[vehicle][cp_node]) + loop(cp_node) = 1 + std::vector> var_coeffs; + for (int vehicle = 0; vehicle < model.vehicles(); vehicle++) { + vehicle_performs_node[vehicle][cp_node] = AddVariable(cp_model, 0, 1); + var_coeffs.push_back({vehicle_performs_node[vehicle][cp_node], 1}); + } + var_coeffs.push_back({is_unperformed[cp_node], 1}); + AddLinearConstraint(cp_model, 1, 1, var_coeffs); + } + const int num_vehicle_classes = model.GetVehicleClassesCount(); + // vehicle_class_performs_node[vehicle_class][node] equals to 1 if the + // vehicle of vehicle_class performs the node, and 0 otherwise. + std::vector> vehicle_class_performs_node( + num_vehicle_classes); + for (int cp_node = 1; cp_node < num_cp_nodes; cp_node++) { + const int node = cp_node - 1; + for (int vehicle_class = 0; vehicle_class < num_vehicle_classes; + vehicle_class++) { + if (model.IsStart(node) || model.IsEnd(node)) { + const int vehicle = model.VehicleIndex(node); + vehicle_class_performs_node[vehicle_class][cp_node] = + vehicle_class == + model.GetVehicleClassIndexOfVehicle(vehicle).value() + ? AddVariable(cp_model, 1, 1) + : AddVariable(cp_model, 0, 0); + continue; + } + vehicle_class_performs_node[vehicle_class][cp_node] = + AddVariable(cp_model, 0, 1); + std::vector> var_coeffs; + for (int vehicle = 0; vehicle < model.vehicles(); vehicle++) { + if (model.GetVehicleClassIndexOfVehicle(vehicle).value() == + vehicle_class) { + var_coeffs.push_back({vehicle_performs_node[vehicle][cp_node], 1}); + // vehicle_performs_node -> vehicle_class_performs_node + AddLinearConstraint( + cp_model, 1, 1, + {{vehicle_class_performs_node[vehicle_class][cp_node], 1}}, + {vehicle_performs_node[vehicle][cp_node]}); + } + } + // vehicle_class_performs_node -> exactly one vehicle from this class + // performs node. + AddLinearConstraint( + cp_model, 1, 1, var_coeffs, + {vehicle_class_performs_node[vehicle_class][cp_node]}); + } + } + // vehicle_class_performs_arc[vehicle_class][arc_var] equals to 1 if the + // vehicle of vehicle_class performs the arc, and 0 otherwise. + std::vector> vehicle_class_performs_arc( + num_vehicle_classes); + // Set "arc" costs. + for (const auto [arc, arc_var] : arc_vars) { + const auto [cp_tail, cp_head] = arc; + if (cp_tail == depot || cp_head == depot) continue; + const int tail = cp_tail - 1; + const int head = cp_head - 1; + // Costs for unperformed arcs have already been set. + if (tail == head) continue; + for (int vehicle = 0; vehicle < model.vehicles(); vehicle++) { + // The arc can't be performed by the vehicle when vehicle can't perform + // arc nodes. + if (!vehicle_performs_node[vehicle].contains(cp_tail) || + !vehicle_performs_node[vehicle].contains(cp_head)) { + continue; + } + int64_t cost = model.GetArcCostForVehicle(tail, head, vehicle); + // Arcs with int64_t's max cost are infeasible. + if (cost == std::numeric_limits::max()) continue; + const int vehicle_class = + model.GetVehicleClassIndexOfVehicle(vehicle).value(); + if (!vehicle_class_performs_arc[vehicle_class].contains(arc_var)) { + vehicle_class_performs_arc[vehicle_class][arc_var] = + AddVariable(cp_model, 0, 1); + // Create constraints to set vehicle_class_performs_arc. + // vehicle_class_performs_arc -> + // vehicle_class_performs_tail & vehicle_class_performs_head & + // arc_is_performed + ConstraintProto* ct = cp_model->add_constraints(); + ct->add_enforcement_literal( + vehicle_class_performs_arc[vehicle_class][arc_var]); + BoolArgumentProto* bool_and = ct->mutable_bool_and(); + bool_and->add_literals( + vehicle_class_performs_node[vehicle_class][cp_tail]); + bool_and->add_literals( + vehicle_class_performs_node[vehicle_class][cp_head]); + bool_and->add_literals(arc_var); + // Don't add arcs with zero cost to the objective. + if (cost != 0) { + cp_model->mutable_objective()->add_vars( + vehicle_class_performs_arc[vehicle_class][arc_var]); + cp_model->mutable_objective()->add_coeffs(cost); + } + } + // (arc_is_performed & vehicle_performs_tail) -> + // (vehicle_class_performs_arc & vehicle_performs_head) + ConstraintProto* ct_arc_tail = cp_model->add_constraints(); + ct_arc_tail->add_enforcement_literal(arc_var); + ct_arc_tail->add_enforcement_literal( + vehicle_performs_node[vehicle][cp_tail]); + ct_arc_tail->mutable_bool_and()->add_literals( + vehicle_class_performs_arc[vehicle_class][arc_var]); + ct_arc_tail->mutable_bool_and()->add_literals( + vehicle_performs_node[vehicle][cp_head]); + // (arc_is_performed & vehicle_performs_head) -> + // (vehicle_class_performs_arc & vehicle_performs_tail) + ConstraintProto* ct_arc_head = cp_model->add_constraints(); + ct_arc_head->add_enforcement_literal(arc_var); + ct_arc_head->add_enforcement_literal( + vehicle_performs_node[vehicle][cp_head]); + ct_arc_head->mutable_bool_and()->add_literals( + vehicle_class_performs_arc[vehicle_class][arc_var]); + ct_arc_head->mutable_bool_and()->add_literals( + vehicle_performs_node[vehicle][cp_tail]); + } + } + + AddGeneralizedPickupDeliveryConstraints( + model, arc_vars, vehicle_performs_node, is_unperformed, cp_model); + + AddGeneralizedDimensions(model, arc_vars, vehicle_performs_node, + vehicle_class_performs_arc, cp_model); + + // Create Routes constraint, ensuring circuits from and to the depot. + RoutesConstraintProto* routes_ct = + cp_model->add_constraints()->mutable_routes(); + for (const auto [arc, arc_var] : arc_vars) { + const int tail = arc.tail; + const int head = arc.head; + routes_ct->add_tails(tail); + routes_ct->add_heads(head); + routes_ct->add_literals(arc_var); + } + + // Add demands and capacities to improve the LP relaxation and cuts. These + // are based on the first "unary" dimension in the model if it exists. + // TODO(user): We might want to try to get demand lower bounds from + // non-unary dimensions if no unary exist. + const RoutingDimension* master_dimension = nullptr; + for (const RoutingDimension* dimension : model.GetDimensions()) { + bool is_unary = true; + for (int vehicle = 0; vehicle < model.vehicles(); vehicle++) { + if (dimension->GetUnaryTransitEvaluator(vehicle) == nullptr) { + is_unary = false; + break; + } + } + if (is_unary) { + master_dimension = dimension; + break; + } + } + if (master_dimension != nullptr) { + for (int cp_node = 0; cp_node < num_cp_nodes; ++cp_node) { + int64_t min_transit = std::numeric_limits::max(); + if (cp_node != 0 && !model.IsEnd(cp_node - 1)) { + for (int vehicle = 0; vehicle < model.vehicles(); vehicle++) { + const RoutingModel::TransitCallback1& transit = + master_dimension->GetUnaryTransitEvaluator(vehicle); + min_transit = std::min(min_transit, transit(cp_node - 1)); + } + } else { + min_transit = 0; + } + routes_ct->add_demands(min_transit); + } + DCHECK_EQ(routes_ct->demands_size(), num_cp_nodes); + int64_t max_capacity = std::numeric_limits::min(); + for (int vehicle = 0; vehicle < model.vehicles(); vehicle++) { + max_capacity = std::max(max_capacity, + master_dimension->vehicle_capacities()[vehicle]); + } + routes_ct->set_capacity(max_capacity); + } + return arc_vars; +} + +// Converts a CpSolverResponse to an Assignment containing next variables. +bool ConvertGeneralizedResponseToSolution(const CpSolverResponse& response, + const RoutingModel& model, + const ArcVarMap& arc_vars, + Assignment* solution) { + if (response.status() != CpSolverStatus::OPTIMAL && + response.status() != CpSolverStatus::FEASIBLE) { + return false; + } + const int depot = 0; + for (const auto [arc, arc_var] : arc_vars) { + if (response.solution(arc_var) == 0) continue; + const auto [tail, head] = arc; + if (head == depot || tail == depot) continue; + solution->Add(model.NextVar(tail - 1))->SetValue(head - 1); + } + return true; +} + +// Uses CP solution as hint for CP-SAT. +void AddSolutionAsHintToGeneralizedModel(const Assignment* solution, + const RoutingModel& model, + const ArcVarMap& arc_vars, + CpModelProto* cp_model) { + if (solution == nullptr) return; + PartialVariableAssignment* const hint = cp_model->mutable_solution_hint(); + hint->Clear(); + const int num_nodes = model.Nexts().size(); + for (int tail = 0; tail < num_nodes; ++tail) { + const int cp_tail = tail + 1; + const int cp_head = solution->Value(model.NextVar(tail)) + 1; + const int* const arc_var = gtl::FindOrNull(arc_vars, {cp_tail, cp_head}); + // Arcs with a cost of max int64_t are not added to the model (considered as + // infeasible). In some rare cases CP solutions might contain such arcs in + // which case they are skipped here and a partial solution is used as a + // hint. + if (arc_var == nullptr) continue; + hint->add_vars(*arc_var); + hint->add_values(1); + } +} + void AddSolutionAsHintToModel(const Assignment* solution, const RoutingModel& model, const ArcVarMap& arc_vars, @@ -501,14 +1011,20 @@ void AddSolutionAsHintToModel(const Assignment* solution, // Returns the response of the search. CpSolverResponse SolveRoutingModel( const CpModelProto& cp_model, absl::Duration remaining_time, + const RoutingSearchParameters& search_parameters, const std::function& observer) { - // TODO(user): Add CP-SAT parameters to routing parameters. - SatParameters parameters; - parameters.set_linearization_level(2); - parameters.set_max_time_in_seconds(absl::ToDoubleSeconds(remaining_time)); - parameters.set_num_search_workers(1); + // Copying to set remaining time. + SatParameters sat_parameters = search_parameters.sat_parameters(); + if (!sat_parameters.has_max_time_in_seconds()) { + sat_parameters.set_max_time_in_seconds( + absl::ToDoubleSeconds(remaining_time)); + } else { + sat_parameters.set_max_time_in_seconds( + std::min(absl::ToDoubleSeconds(remaining_time), + sat_parameters.max_time_in_seconds())); + } Model model; - model.Add(NewSatParameters(parameters)); + model.Add(NewSatParameters(sat_parameters)); if (observer != nullptr) { model.Add(NewFeasibleSolutionObserver(observer)); } @@ -517,6 +1033,20 @@ CpSolverResponse SolveRoutingModel( return SolveCpModel(cp_model, &model); } +// Check if all the nodes are present in arcs. Otherwise, CP-SAT solver may +// fail. +bool IsFeasibleArcVarMap(const ArcVarMap& arc_vars, int max_node_index) { + Bitset64<> present_in_arcs(max_node_index + 1); + for (const auto [arc, _] : arc_vars) { + present_in_arcs.Set(arc.head); + present_in_arcs.Set(arc.tail); + } + for (int i = 0; i <= max_node_index; i++) { + if (!present_in_arcs[i]) return false; + } + return true; +} + } // namespace } // namespace sat @@ -526,17 +1056,30 @@ bool SolveModelWithSat(const RoutingModel& model, const RoutingSearchParameters& search_parameters, const Assignment* initial_solution, Assignment* solution) { - if (!sat::RoutingModelCanBeSolvedBySat(model)) return false; sat::CpModelProto cp_model; cp_model.mutable_objective()->set_scaling_factor( search_parameters.log_cost_scaling_factor()); cp_model.mutable_objective()->set_offset(search_parameters.log_cost_offset()); + if (search_parameters.use_generalized_cp_sat() == BOOL_TRUE) { + const sat::ArcVarMap arc_vars = + sat::PopulateGeneralizedRouteModelFromRoutingModel(model, &cp_model); + const int max_node_index = model.Nexts().size() + model.vehicles(); + if (!sat::IsFeasibleArcVarMap(arc_vars, max_node_index)) return false; + sat::AddSolutionAsHintToGeneralizedModel(initial_solution, model, arc_vars, + &cp_model); + return sat::ConvertGeneralizedResponseToSolution( + sat::SolveRoutingModel(cp_model, model.RemainingTime(), + search_parameters, nullptr), + model, arc_vars, solution); + } + if (!sat::RoutingModelCanBeSolvedBySat(model)) return false; const sat::ArcVarMap arc_vars = sat::PopulateModelFromRoutingModel(model, &cp_model); sat::AddSolutionAsHintToModel(initial_solution, model, arc_vars, &cp_model); return sat::ConvertToSolution( - sat::SolveRoutingModel(cp_model, model.RemainingTime(), nullptr), model, - arc_vars, solution); + sat::SolveRoutingModel(cp_model, model.RemainingTime(), search_parameters, + nullptr), + model, arc_vars, solution); } } // namespace operations_research