1 // Copyright 2010-2021 Google LLC 2 // Licensed under the Apache License, Version 2.0 (the "License"); 3 // you may not use this file except in compliance with the License. 4 // You may obtain a copy of the License at 5 // 6 // http://www.apache.org/licenses/LICENSE-2.0 7 // 8 // Unless required by applicable law or agreed to in writing, software 9 // distributed under the License is distributed on an "AS IS" BASIS, 10 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 11 // See the License for the specific language governing permissions and 12 // limitations under the License. 13 14 #ifndef OR_TOOLS_CONSTRAINT_SOLVER_ROUTING_LP_SCHEDULING_H_ 15 #define OR_TOOLS_CONSTRAINT_SOLVER_ROUTING_LP_SCHEDULING_H_ 16 17 #include <algorithm> 18 #include <cstdint> 19 #include <deque> 20 #include <functional> 21 #include <limits> 22 #include <map> 23 #include <memory> 24 #include <string> 25 #include <utility> 26 #include <vector> 27 28 #include "absl/container/flat_hash_map.h" 29 #include "absl/memory/memory.h" 30 #include "absl/time/time.h" 31 #include "ortools/base/logging.h" 32 #include "ortools/base/mathutil.h" 33 #include "ortools/constraint_solver/routing.h" 34 #include "ortools/constraint_solver/routing_parameters.pb.h" 35 #include "ortools/glop/lp_solver.h" 36 #include "ortools/glop/parameters.pb.h" 37 #include "ortools/lp_data/lp_data.h" 38 #include "ortools/lp_data/lp_types.h" 39 #include "ortools/sat/cp_model.pb.h" 40 #include "ortools/sat/cp_model_solver.h" 41 #include "ortools/sat/lp_utils.h" 42 #include "ortools/sat/model.h" 43 #include "ortools/sat/sat_parameters.pb.h" 44 #include "ortools/util/saturated_arithmetic.h" 45 #include "ortools/util/sorted_interval_list.h" 46 47 namespace operations_research { 48 49 // Classes to solve dimension cumul placement (aka scheduling) problems using 50 // linear programming. 51 52 // Utility class used in the core optimizer to tighten the cumul bounds as much 53 // as possible based on the model precedences. 54 class CumulBoundsPropagator { 55 public: 56 explicit CumulBoundsPropagator(const RoutingDimension* dimension); 57 58 // Tightens the cumul bounds starting from the current cumul var min/max, 59 // and propagating the precedences resulting from the next_accessor, and the 60 // dimension's precedence rules. 61 // Returns false iff the precedences are infeasible with the given routes. 62 // Otherwise, the user can call CumulMin() and CumulMax() to retrieve the new 63 // bounds of an index. 64 bool PropagateCumulBounds( 65 const std::function<int64_t(int64_t)>& next_accessor, 66 int64_t cumul_offset); 67 CumulMin(int index)68 int64_t CumulMin(int index) const { 69 return propagated_bounds_[PositiveNode(index)]; 70 } 71 CumulMax(int index)72 int64_t CumulMax(int index) const { 73 const int64_t negated_upper_bound = propagated_bounds_[NegativeNode(index)]; 74 return negated_upper_bound == std::numeric_limits<int64_t>::min() 75 ? std::numeric_limits<int64_t>::max() 76 : -negated_upper_bound; 77 } 78 dimension()79 const RoutingDimension& dimension() const { return dimension_; } 80 81 private: 82 // An arc "tail --offset--> head" represents the relation 83 // tail + offset <= head. 84 // As arcs are stored by tail, we don't store it in the struct. 85 struct ArcInfo { 86 int head; 87 int64_t offset; 88 }; 89 static const int kNoParent; 90 static const int kParentToBePropagated; 91 92 // Return the node corresponding to the lower bound of the cumul of index and 93 // -index respectively. PositiveNode(int index)94 int PositiveNode(int index) const { return 2 * index; } NegativeNode(int index)95 int NegativeNode(int index) const { return 2 * index + 1; } 96 AddNodeToQueue(int node)97 void AddNodeToQueue(int node) { 98 if (!node_in_queue_[node]) { 99 bf_queue_.push_back(node); 100 node_in_queue_[node] = true; 101 } 102 } 103 104 // Adds the relation first_index + offset <= second_index, by adding arcs 105 // first_index --offset--> second_index and 106 // -second_index --offset--> -first_index. 107 void AddArcs(int first_index, int second_index, int64_t offset); 108 109 bool InitializeArcsAndBounds( 110 const std::function<int64_t(int64_t)>& next_accessor, 111 int64_t cumul_offset); 112 113 bool UpdateCurrentLowerBoundOfNode(int node, int64_t new_lb, int64_t offset); 114 115 bool DisassembleSubtree(int source, int target); 116 CleanupAndReturnFalse()117 bool CleanupAndReturnFalse() { 118 // We clean-up node_in_queue_ for future calls, and return false. 119 for (int node_to_cleanup : bf_queue_) { 120 node_in_queue_[node_to_cleanup] = false; 121 } 122 bf_queue_.clear(); 123 return false; 124 } 125 126 const RoutingDimension& dimension_; 127 const int64_t num_nodes_; 128 129 // TODO(user): Investigate if all arcs for a given tail can be created 130 // at the same time, in which case outgoing_arcs_ could point to an absl::Span 131 // for each tail index. 132 std::vector<std::vector<ArcInfo>> outgoing_arcs_; 133 134 std::deque<int> bf_queue_; 135 std::vector<bool> node_in_queue_; 136 std::vector<int> tree_parent_node_of_; 137 // After calling PropagateCumulBounds(), for each node index n, 138 // propagated_bounds_[2*n] and -propagated_bounds_[2*n+1] respectively contain 139 // the propagated lower and upper bounds of n's cumul variable. 140 std::vector<int64_t> propagated_bounds_; 141 142 // Vector used in DisassembleSubtree() to avoid memory reallocation. 143 std::vector<int> tmp_dfs_stack_; 144 145 // Used to store the pickup/delivery pairs encountered on the routes. 146 std::vector<std::pair<int64_t, int64_t>> 147 visited_pickup_delivery_indices_for_pair_; 148 }; 149 150 enum class DimensionSchedulingStatus { 151 // An optimal solution was found respecting all constraints. 152 OPTIMAL, 153 // An optimal solution was found, however constraints which were relaxed were 154 // violated. 155 RELAXED_OPTIMAL_ONLY, 156 // A solution could not be found. 157 INFEASIBLE 158 }; 159 160 class RoutingLinearSolverWrapper { 161 public: ~RoutingLinearSolverWrapper()162 virtual ~RoutingLinearSolverWrapper() {} 163 virtual void Clear() = 0; 164 virtual int CreateNewPositiveVariable() = 0; 165 virtual bool SetVariableBounds(int index, int64_t lower_bound, 166 int64_t upper_bound) = 0; 167 virtual void SetVariableDisjointBounds(int index, 168 const std::vector<int64_t>& starts, 169 const std::vector<int64_t>& ends) = 0; 170 virtual int64_t GetVariableLowerBound(int index) const = 0; 171 virtual int64_t GetVariableUpperBound(int index) const = 0; 172 virtual void SetObjectiveCoefficient(int index, double coefficient) = 0; 173 virtual double GetObjectiveCoefficient(int index) const = 0; 174 virtual void ClearObjective() = 0; 175 virtual int NumVariables() const = 0; 176 virtual int CreateNewConstraint(int64_t lower_bound, int64_t upper_bound) = 0; 177 virtual void SetCoefficient(int ct, int index, double coefficient) = 0; 178 virtual bool IsCPSATSolver() = 0; 179 virtual void AddObjectiveConstraint() = 0; 180 virtual void AddMaximumConstraint(int max_var, std::vector<int> vars) = 0; 181 virtual void AddProductConstraint(int product_var, std::vector<int> vars) = 0; 182 virtual void SetEnforcementLiteral(int ct, int condition) = 0; 183 virtual DimensionSchedulingStatus Solve(absl::Duration duration_limit) = 0; 184 virtual int64_t GetObjectiveValue() const = 0; 185 virtual double GetValue(int index) const = 0; 186 virtual bool SolutionIsInteger() const = 0; 187 188 // Adds a variable with bounds [lower_bound, upper_bound]. AddVariable(int64_t lower_bound,int64_t upper_bound)189 int AddVariable(int64_t lower_bound, int64_t upper_bound) { 190 CHECK_LE(lower_bound, upper_bound); 191 const int variable = CreateNewPositiveVariable(); 192 SetVariableBounds(variable, lower_bound, upper_bound); 193 return variable; 194 } 195 // Adds a linear constraint, enforcing 196 // lower_bound <= sum variable * coeff <= upper_bound, 197 // and returns the identifier of that constraint. AddLinearConstraint(int64_t lower_bound,int64_t upper_bound,const std::vector<std::pair<int,double>> & variable_coeffs)198 int AddLinearConstraint( 199 int64_t lower_bound, int64_t upper_bound, 200 const std::vector<std::pair<int, double>>& variable_coeffs) { 201 CHECK_LE(lower_bound, upper_bound); 202 const int ct = CreateNewConstraint(lower_bound, upper_bound); 203 for (const auto& variable_coeff : variable_coeffs) { 204 SetCoefficient(ct, variable_coeff.first, variable_coeff.second); 205 } 206 return ct; 207 } 208 // Adds a linear constraint and a 0/1 variable that is true iff 209 // lower_bound <= sum variable * coeff <= upper_bound, 210 // and returns the identifier of that variable. AddReifiedLinearConstraint(int64_t lower_bound,int64_t upper_bound,const std::vector<std::pair<int,double>> & weighted_variables)211 int AddReifiedLinearConstraint( 212 int64_t lower_bound, int64_t upper_bound, 213 const std::vector<std::pair<int, double>>& weighted_variables) { 214 const int reification_ct = AddLinearConstraint(1, 1, {}); 215 if (std::numeric_limits<int64_t>::min() < lower_bound) { 216 const int under_lower_bound = AddVariable(0, 1); 217 SetCoefficient(reification_ct, under_lower_bound, 1); 218 const int under_lower_bound_ct = 219 AddLinearConstraint(std::numeric_limits<int64_t>::min(), 220 lower_bound - 1, weighted_variables); 221 SetEnforcementLiteral(under_lower_bound_ct, under_lower_bound); 222 } 223 if (upper_bound < std::numeric_limits<int64_t>::max()) { 224 const int above_upper_bound = AddVariable(0, 1); 225 SetCoefficient(reification_ct, above_upper_bound, 1); 226 const int above_upper_bound_ct = AddLinearConstraint( 227 upper_bound + 1, std::numeric_limits<int64_t>::max(), 228 weighted_variables); 229 SetEnforcementLiteral(above_upper_bound_ct, above_upper_bound); 230 } 231 const int within_bounds = AddVariable(0, 1); 232 SetCoefficient(reification_ct, within_bounds, 1); 233 const int within_bounds_ct = 234 AddLinearConstraint(lower_bound, upper_bound, weighted_variables); 235 SetEnforcementLiteral(within_bounds_ct, within_bounds); 236 return within_bounds; 237 } 238 }; 239 240 class RoutingGlopWrapper : public RoutingLinearSolverWrapper { 241 public: RoutingGlopWrapper(bool is_relaxation,const glop::GlopParameters & parameters)242 RoutingGlopWrapper(bool is_relaxation, const glop::GlopParameters& parameters) 243 : is_relaxation_(is_relaxation) { 244 lp_solver_.SetParameters(parameters); 245 linear_program_.SetMaximizationProblem(false); 246 } Clear()247 void Clear() override { 248 linear_program_.Clear(); 249 linear_program_.SetMaximizationProblem(false); 250 allowed_intervals_.clear(); 251 } CreateNewPositiveVariable()252 int CreateNewPositiveVariable() override { 253 return linear_program_.CreateNewVariable().value(); 254 } SetVariableBounds(int index,int64_t lower_bound,int64_t upper_bound)255 bool SetVariableBounds(int index, int64_t lower_bound, 256 int64_t upper_bound) override { 257 DCHECK_GE(lower_bound, 0); 258 // When variable upper bounds are greater than this threshold, precision 259 // issues arise in GLOP. In this case we are just going to suppose that 260 // these high bound values are infinite and not set the upper bound. 261 const int64_t kMaxValue = 1e10; 262 const double lp_min = lower_bound; 263 const double lp_max = 264 (upper_bound > kMaxValue) ? glop::kInfinity : upper_bound; 265 if (lp_min <= lp_max) { 266 linear_program_.SetVariableBounds(glop::ColIndex(index), lp_min, lp_max); 267 return true; 268 } 269 // The linear_program would not be feasible, and it cannot handle the 270 // lp_min > lp_max case, so we must detect infeasibility here. 271 return false; 272 } SetVariableDisjointBounds(int index,const std::vector<int64_t> & starts,const std::vector<int64_t> & ends)273 void SetVariableDisjointBounds(int index, const std::vector<int64_t>& starts, 274 const std::vector<int64_t>& ends) override { 275 // TODO(user): Investigate if we can avoid rebuilding the interval list 276 // each time (we could keep a reference to the forbidden interval list in 277 // RoutingDimension but we would need to store cumul offsets and use them 278 // when checking intervals). 279 allowed_intervals_[index] = 280 absl::make_unique<SortedDisjointIntervalList>(starts, ends); 281 } GetVariableLowerBound(int index)282 int64_t GetVariableLowerBound(int index) const override { 283 return linear_program_.variable_lower_bounds()[glop::ColIndex(index)]; 284 } GetVariableUpperBound(int index)285 int64_t GetVariableUpperBound(int index) const override { 286 const double upper_bound = 287 linear_program_.variable_upper_bounds()[glop::ColIndex(index)]; 288 DCHECK_GE(upper_bound, 0); 289 return upper_bound == glop::kInfinity ? std::numeric_limits<int64_t>::max() 290 : static_cast<int64_t>(upper_bound); 291 } SetObjectiveCoefficient(int index,double coefficient)292 void SetObjectiveCoefficient(int index, double coefficient) override { 293 linear_program_.SetObjectiveCoefficient(glop::ColIndex(index), coefficient); 294 } GetObjectiveCoefficient(int index)295 double GetObjectiveCoefficient(int index) const override { 296 return linear_program_.objective_coefficients()[glop::ColIndex(index)]; 297 } ClearObjective()298 void ClearObjective() override { 299 for (glop::ColIndex i(0); i < linear_program_.num_variables(); ++i) { 300 linear_program_.SetObjectiveCoefficient(i, 0); 301 } 302 } NumVariables()303 int NumVariables() const override { 304 return linear_program_.num_variables().value(); 305 } CreateNewConstraint(int64_t lower_bound,int64_t upper_bound)306 int CreateNewConstraint(int64_t lower_bound, int64_t upper_bound) override { 307 const glop::RowIndex ct = linear_program_.CreateNewConstraint(); 308 linear_program_.SetConstraintBounds( 309 ct, 310 (lower_bound == std::numeric_limits<int64_t>::min()) ? -glop::kInfinity 311 : lower_bound, 312 (upper_bound == std::numeric_limits<int64_t>::max()) ? glop::kInfinity 313 : upper_bound); 314 return ct.value(); 315 } SetCoefficient(int ct,int index,double coefficient)316 void SetCoefficient(int ct, int index, double coefficient) override { 317 linear_program_.SetCoefficient(glop::RowIndex(ct), glop::ColIndex(index), 318 coefficient); 319 } IsCPSATSolver()320 bool IsCPSATSolver() override { return false; } AddObjectiveConstraint()321 void AddObjectiveConstraint() override { 322 const int ct = CreateNewConstraint(0, GetObjectiveValue()); 323 for (int variable = 0; variable < NumVariables(); variable++) { 324 const double coefficient = GetObjectiveCoefficient(variable); 325 if (coefficient != 0) { 326 SetCoefficient(ct, variable, coefficient); 327 } 328 } 329 } AddMaximumConstraint(int max_var,std::vector<int> vars)330 void AddMaximumConstraint(int max_var, std::vector<int> vars) override {} AddProductConstraint(int product_var,std::vector<int> vars)331 void AddProductConstraint(int product_var, std::vector<int> vars) override {} SetEnforcementLiteral(int ct,int condition)332 void SetEnforcementLiteral(int ct, int condition) override{}; Solve(absl::Duration duration_limit)333 DimensionSchedulingStatus Solve(absl::Duration duration_limit) override { 334 lp_solver_.GetMutableParameters()->set_max_time_in_seconds( 335 absl::ToDoubleSeconds(duration_limit)); 336 337 // Because we construct the lp one constraint at a time and we never call 338 // SetCoefficient() on the same variable twice for a constraint, we know 339 // that the columns do not contain duplicates and are already ordered by 340 // constraint so we do not need to call linear_program->CleanUp() which can 341 // be costly. Note that the assumptions are DCHECKed() in the call below. 342 linear_program_.NotifyThatColumnsAreClean(); 343 VLOG(2) << linear_program_.Dump(); 344 const glop::ProblemStatus status = lp_solver_.Solve(linear_program_); 345 if (status != glop::ProblemStatus::OPTIMAL && 346 status != glop::ProblemStatus::IMPRECISE) { 347 return DimensionSchedulingStatus::INFEASIBLE; 348 } 349 if (is_relaxation_) { 350 return DimensionSchedulingStatus::RELAXED_OPTIMAL_ONLY; 351 } 352 for (const auto& allowed_interval : allowed_intervals_) { 353 const double value_double = GetValue(allowed_interval.first); 354 const int64_t value = 355 (value_double >= std::numeric_limits<int64_t>::max()) 356 ? std::numeric_limits<int64_t>::max() 357 : MathUtil::FastInt64Round(value_double); 358 const SortedDisjointIntervalList* const interval_list = 359 allowed_interval.second.get(); 360 const auto it = interval_list->FirstIntervalGreaterOrEqual(value); 361 if (it == interval_list->end() || value < it->start) { 362 return DimensionSchedulingStatus::RELAXED_OPTIMAL_ONLY; 363 } 364 } 365 return DimensionSchedulingStatus::OPTIMAL; 366 } GetObjectiveValue()367 int64_t GetObjectiveValue() const override { 368 return MathUtil::FastInt64Round(lp_solver_.GetObjectiveValue()); 369 } GetValue(int index)370 double GetValue(int index) const override { 371 return lp_solver_.variable_values()[glop::ColIndex(index)]; 372 } SolutionIsInteger()373 bool SolutionIsInteger() const override { 374 return linear_program_.SolutionIsInteger(lp_solver_.variable_values(), 375 /*absolute_tolerance*/ 1e-3); 376 } 377 378 private: 379 const bool is_relaxation_; 380 glop::LinearProgram linear_program_; 381 glop::LPSolver lp_solver_; 382 absl::flat_hash_map<int, std::unique_ptr<SortedDisjointIntervalList>> 383 allowed_intervals_; 384 }; 385 386 class RoutingCPSatWrapper : public RoutingLinearSolverWrapper { 387 public: RoutingCPSatWrapper()388 RoutingCPSatWrapper() { 389 parameters_.set_num_search_workers(1); 390 // Keeping presolve but with 0 iterations; as of 11/2019 it is 391 // significantly faster than both full presolve and no presolve. 392 parameters_.set_cp_model_presolve(true); 393 parameters_.set_max_presolve_iterations(0); 394 parameters_.set_catch_sigint_signal(false); 395 parameters_.set_mip_max_bound(1e8); 396 parameters_.set_search_branching(sat::SatParameters::LP_SEARCH); 397 } ~RoutingCPSatWrapper()398 ~RoutingCPSatWrapper() override {} Clear()399 void Clear() override { 400 model_.Clear(); 401 response_.Clear(); 402 objective_coefficients_.clear(); 403 } CreateNewPositiveVariable()404 int CreateNewPositiveVariable() override { 405 const int index = model_.variables_size(); 406 sat::IntegerVariableProto* const variable = model_.add_variables(); 407 variable->add_domain(0); 408 variable->add_domain(static_cast<int64_t>(parameters_.mip_max_bound())); 409 return index; 410 } SetVariableBounds(int index,int64_t lower_bound,int64_t upper_bound)411 bool SetVariableBounds(int index, int64_t lower_bound, 412 int64_t upper_bound) override { 413 DCHECK_GE(lower_bound, 0); 414 const int64_t capped_upper_bound = 415 std::min<int64_t>(upper_bound, parameters_.mip_max_bound()); 416 if (lower_bound > capped_upper_bound) return false; 417 sat::IntegerVariableProto* const variable = model_.mutable_variables(index); 418 variable->set_domain(0, lower_bound); 419 variable->set_domain(1, capped_upper_bound); 420 return true; 421 } SetVariableDisjointBounds(int index,const std::vector<int64_t> & starts,const std::vector<int64_t> & ends)422 void SetVariableDisjointBounds(int index, const std::vector<int64_t>& starts, 423 const std::vector<int64_t>& ends) override { 424 DCHECK_EQ(starts.size(), ends.size()); 425 const int ct = CreateNewConstraint(1, 1); 426 for (int i = 0; i < starts.size(); ++i) { 427 const int variable = CreateNewPositiveVariable(); 428 SetVariableBounds(variable, 0, 1); 429 SetCoefficient(ct, variable, 1); 430 const int window_ct = CreateNewConstraint(starts[i], ends[i]); 431 SetCoefficient(window_ct, index, 1); 432 model_.mutable_constraints(window_ct)->add_enforcement_literal(variable); 433 } 434 } GetVariableLowerBound(int index)435 int64_t GetVariableLowerBound(int index) const override { 436 return model_.variables(index).domain(0); 437 } GetVariableUpperBound(int index)438 int64_t GetVariableUpperBound(int index) const override { 439 const auto& domain = model_.variables(index).domain(); 440 return domain[domain.size() - 1]; 441 } SetObjectiveCoefficient(int index,double coefficient)442 void SetObjectiveCoefficient(int index, double coefficient) override { 443 if (index >= objective_coefficients_.size()) { 444 objective_coefficients_.resize(index + 1, 0); 445 } 446 objective_coefficients_[index] = coefficient; 447 sat::FloatObjectiveProto* const objective = 448 model_.mutable_floating_point_objective(); 449 objective->add_vars(index); 450 objective->add_coeffs(coefficient); 451 } GetObjectiveCoefficient(int index)452 double GetObjectiveCoefficient(int index) const override { 453 return (index < objective_coefficients_.size()) 454 ? objective_coefficients_[index] 455 : 0; 456 } ClearObjective()457 void ClearObjective() override { 458 model_.mutable_floating_point_objective()->Clear(); 459 } NumVariables()460 int NumVariables() const override { return model_.variables_size(); } CreateNewConstraint(int64_t lower_bound,int64_t upper_bound)461 int CreateNewConstraint(int64_t lower_bound, int64_t upper_bound) override { 462 sat::LinearConstraintProto* const ct = 463 model_.add_constraints()->mutable_linear(); 464 ct->add_domain(lower_bound); 465 ct->add_domain(upper_bound); 466 return model_.constraints_size() - 1; 467 } SetCoefficient(int ct_index,int index,double coefficient)468 void SetCoefficient(int ct_index, int index, double coefficient) override { 469 sat::LinearConstraintProto* const ct = 470 model_.mutable_constraints(ct_index)->mutable_linear(); 471 ct->add_vars(index); 472 const int64_t integer_coefficient = coefficient; 473 ct->add_coeffs(integer_coefficient); 474 } IsCPSATSolver()475 bool IsCPSATSolver() override { return true; } AddObjectiveConstraint()476 void AddObjectiveConstraint() override { 477 // Scale the objective to get its integer representation. 478 const sat::FloatObjectiveProto& fp_objective = 479 model_.floating_point_objective(); 480 std::vector<std::pair<int, double>> terms; 481 for (int i = 0; i < fp_objective.vars_size(); ++i) { 482 terms.push_back({fp_objective.vars(i), fp_objective.coeffs(i)}); 483 } 484 SolverLogger logger; 485 const bool status = 486 sat::ScaleAndSetObjective(parameters_, terms, fp_objective.offset(), 487 fp_objective.maximize(), &model_, &logger); 488 DCHECK(status); 489 const sat::CpObjectiveProto& objective = model_.objective(); 490 int64_t activity = 0; 491 for (int i = 0; i < objective.vars_size(); ++i) { 492 activity += response_.solution(objective.vars(i)) * objective.coeffs(i); 493 } 494 const int ct = CreateNewConstraint(0, activity); 495 for (int i = 0; i < objective.vars_size(); ++i) { 496 SetCoefficient(ct, objective.vars(i), objective.coeffs(i)); 497 } 498 model_.clear_objective(); 499 } AddMaximumConstraint(int max_var,std::vector<int> vars)500 void AddMaximumConstraint(int max_var, std::vector<int> vars) override { 501 sat::LinearArgumentProto* const ct = 502 model_.add_constraints()->mutable_lin_max(); 503 ct->mutable_target()->add_vars(max_var); 504 ct->mutable_target()->add_coeffs(1); 505 for (const int var : vars) { 506 sat::LinearExpressionProto* const expr = ct->add_exprs(); 507 expr->add_vars(var); 508 expr->add_coeffs(1); 509 } 510 } AddProductConstraint(int product_var,std::vector<int> vars)511 void AddProductConstraint(int product_var, std::vector<int> vars) override { 512 sat::LinearArgumentProto* const ct = 513 model_.add_constraints()->mutable_int_prod(); 514 ct->mutable_target()->add_vars(product_var); 515 ct->mutable_target()->add_coeffs(1); 516 for (const int var : vars) { 517 sat::LinearExpressionProto* expr = ct->add_exprs(); 518 expr->add_vars(var); 519 expr->add_coeffs(1); 520 } 521 } SetEnforcementLiteral(int ct,int condition)522 void SetEnforcementLiteral(int ct, int condition) override { 523 DCHECK_LT(ct, model_.constraints_size()); 524 model_.mutable_constraints(ct)->add_enforcement_literal(condition); 525 } Solve(absl::Duration duration_limit)526 DimensionSchedulingStatus Solve(absl::Duration duration_limit) override { 527 parameters_.set_max_time_in_seconds(absl::ToDoubleSeconds(duration_limit)); 528 VLOG(2) << model_.DebugString(); 529 if (hint_.vars_size() == model_.variables_size()) { 530 *model_.mutable_solution_hint() = hint_; 531 } 532 sat::Model model; 533 model.Add(sat::NewSatParameters(parameters_)); 534 response_ = sat::SolveCpModel(model_, &model); 535 VLOG(2) << response_.DebugString(); 536 if (response_.status() == sat::CpSolverStatus::OPTIMAL || 537 (response_.status() == sat::CpSolverStatus::FEASIBLE && 538 !model_.has_floating_point_objective())) { 539 hint_.Clear(); 540 for (int i = 0; i < response_.solution_size(); ++i) { 541 hint_.add_vars(i); 542 hint_.add_values(response_.solution(i)); 543 } 544 return DimensionSchedulingStatus::OPTIMAL; 545 } 546 return DimensionSchedulingStatus::INFEASIBLE; 547 } GetObjectiveValue()548 int64_t GetObjectiveValue() const override { 549 return MathUtil::FastInt64Round(response_.objective_value()); 550 } GetValue(int index)551 double GetValue(int index) const override { 552 return response_.solution(index); 553 } SolutionIsInteger()554 bool SolutionIsInteger() const override { return true; } 555 556 private: 557 sat::CpModelProto model_; 558 sat::CpSolverResponse response_; 559 sat::SatParameters parameters_; 560 std::vector<double> objective_coefficients_; 561 sat::PartialVariableAssignment hint_; 562 }; 563 564 // Utility class used in Local/GlobalDimensionCumulOptimizer to set the linear 565 // solver constraints and solve the problem. 566 class DimensionCumulOptimizerCore { 567 public: 568 DimensionCumulOptimizerCore(const RoutingDimension* dimension, 569 bool use_precedence_propagator); 570 571 // In the OptimizeSingleRoute() and Optimize() methods, if both "cumul_values" 572 // and "cost" parameters are null, we don't optimize the cost and stop at the 573 // first feasible solution in the linear solver (since in this case only 574 // feasibility is of interest). 575 DimensionSchedulingStatus OptimizeSingleRoute( 576 int vehicle, const std::function<int64_t(int64_t)>& next_accessor, 577 RoutingLinearSolverWrapper* solver, std::vector<int64_t>* cumul_values, 578 std::vector<int64_t>* break_values, int64_t* cost, int64_t* transit_cost, 579 bool clear_lp = true); 580 581 std::vector<DimensionSchedulingStatus> OptimizeSingleRouteWithResources( 582 int vehicle, const std::function<int64_t(int64_t)>& next_accessor, 583 const std::vector<RoutingModel::ResourceGroup::Resource>& resources, 584 const std::vector<int>& resource_indices, bool optimize_vehicle_costs, 585 RoutingLinearSolverWrapper* solver, 586 std::vector<int64_t>* costs_without_transits, 587 std::vector<std::vector<int64_t>>* cumul_values, 588 std::vector<std::vector<int64_t>>* break_values, bool clear_lp = true); 589 590 DimensionSchedulingStatus Optimize( 591 const std::function<int64_t(int64_t)>& next_accessor, 592 RoutingLinearSolverWrapper* solver, std::vector<int64_t>* cumul_values, 593 std::vector<int64_t>* break_values, 594 std::vector<std::vector<int>>* resource_indices_per_group, int64_t* cost, 595 int64_t* transit_cost, bool clear_lp = true); 596 597 DimensionSchedulingStatus OptimizeAndPack( 598 const std::function<int64_t(int64_t)>& next_accessor, 599 RoutingLinearSolverWrapper* solver, std::vector<int64_t>* cumul_values, 600 std::vector<int64_t>* break_values, 601 std::vector<std::vector<int>>* resource_indices_per_group); 602 603 DimensionSchedulingStatus OptimizeAndPackSingleRoute( 604 int vehicle, const std::function<int64_t(int64_t)>& next_accessor, 605 const RoutingModel::ResourceGroup::Resource* resource, 606 RoutingLinearSolverWrapper* solver, std::vector<int64_t>* cumul_values, 607 std::vector<int64_t>* break_values); 608 dimension()609 const RoutingDimension* dimension() const { return dimension_; } 610 611 private: 612 // Initializes the containers and given solver. Must be called prior to 613 // setting any contraints and solving. 614 void InitOptimizer(RoutingLinearSolverWrapper* solver); 615 616 // Computes the minimum/maximum of cumuls for nodes on "route", and sets them 617 // in current_route_[min|max]_cumuls_ respectively. 618 // If the propagator_ is not null, uses the bounds tightened by the 619 // propagator. 620 // Otherwise, the bounds are computed by going over the nodes on the route 621 // using the CP bounds, and the fixed transits are used to tighten them. 622 bool ComputeRouteCumulBounds(const std::vector<int64_t>& route, 623 const std::vector<int64_t>& fixed_transits, 624 int64_t cumul_offset); 625 626 // Sets the constraints for all nodes on "vehicle"'s route according to 627 // "next_accessor". If optimize_costs is true, also sets the objective 628 // coefficients for the LP. 629 // Returns false if some infeasibility was detected, true otherwise. 630 bool SetRouteCumulConstraints( 631 int vehicle, const std::function<int64_t(int64_t)>& next_accessor, 632 int64_t cumul_offset, bool optimize_costs, 633 RoutingLinearSolverWrapper* solver, int64_t* route_transit_cost, 634 int64_t* route_cost_offset); 635 636 // Sets the global constraints on the dimension, and adds global objective 637 // cost coefficients if optimize_costs is true. 638 // NOTE: When called, the call to this function MUST come after 639 // SetRouteCumulConstraints() has been called on all routes, so that 640 // index_to_cumul_variable_ and min_start/max_end_cumul_ are correctly 641 // initialized. 642 // Returns false if some infeasibility was detected, true otherwise. 643 bool SetGlobalConstraints( 644 const std::function<int64_t(int64_t)>& next_accessor, 645 int64_t cumul_offset, bool optimize_costs, 646 RoutingLinearSolverWrapper* solver); 647 648 void SetValuesFromLP(const std::vector<int>& lp_variables, int64_t offset, 649 RoutingLinearSolverWrapper* solver, 650 std::vector<int64_t>* lp_values) const; 651 652 void SetResourceIndices( 653 RoutingLinearSolverWrapper* solver, 654 std::vector<std::vector<int>>* resource_indices_per_group) const; 655 656 // This function packs the routes of the given vehicles while keeping the cost 657 // of the LP lower than its current (supposed optimal) objective value. 658 // It does so by setting the current objective variables' coefficient to 0 and 659 // setting the coefficient of the route ends to 1, to first minimize the route 660 // ends' cumuls, and then maximizes the starts' cumuls without increasing the 661 // ends. 662 DimensionSchedulingStatus PackRoutes(std::vector<int> vehicles, 663 RoutingLinearSolverWrapper* solver); 664 665 std::unique_ptr<CumulBoundsPropagator> propagator_; 666 std::vector<int64_t> current_route_min_cumuls_; 667 std::vector<int64_t> current_route_max_cumuls_; 668 const RoutingDimension* const dimension_; 669 // Scheduler variables for current route cumuls and for all nodes cumuls. 670 std::vector<int> current_route_cumul_variables_; 671 std::vector<int> index_to_cumul_variable_; 672 // Scheduler variables for current route breaks and all vehicle breaks. 673 // There are two variables for each break: start and end. 674 // current_route_break_variables_ has variables corresponding to 675 // break[0] start, break[0] end, break[1] start, break[1] end, etc. 676 std::vector<int> current_route_break_variables_; 677 // Vector all_break_variables contains the break variables of all vehicles, 678 // in the same format as current_route_break_variables. 679 // It is the concatenation of break variables of vehicles in [0, #vehicles). 680 std::vector<int> all_break_variables_; 681 // Allows to retrieve break variables of a given vehicle: those go from 682 // all_break_variables_[vehicle_to_all_break_variables_offset_[vehicle]] to 683 // all_break_variables[vehicle_to_all_break_variables_offset_[vehicle+1]-1]. 684 std::vector<int> vehicle_to_all_break_variables_offset_; 685 // The following vector contains indices of resource-to-vehicle assignment 686 // variables. For every resource group, stores indices of 687 // num_resources*num_vehicles boolean variables indicating whether resource #r 688 // is assigned to vehicle #v. 689 std::vector<std::vector<int>> 690 resource_group_to_resource_to_vehicle_assignment_variables_; 691 692 int max_end_cumul_; 693 int min_start_cumul_; 694 std::vector<std::pair<int64_t, int64_t>> 695 visited_pickup_delivery_indices_for_pair_; 696 }; 697 698 // Class used to compute optimal values for dimension cumuls of routes, 699 // minimizing cumul soft lower and upper bound costs, and vehicle span costs of 700 // a route. 701 // In its methods, next_accessor is a callback returning the next node of a 702 // given node on a route. 703 class LocalDimensionCumulOptimizer { 704 public: 705 LocalDimensionCumulOptimizer( 706 const RoutingDimension* dimension, 707 RoutingSearchParameters::SchedulingSolver solver_type); 708 709 // If feasible, computes the optimal cost of the route performed by a vehicle, 710 // minimizing cumul soft lower and upper bound costs and vehicle span costs, 711 // and stores it in "optimal_cost" (if not null). 712 // Returns true iff the route respects all constraints. 713 DimensionSchedulingStatus ComputeRouteCumulCost( 714 int vehicle, const std::function<int64_t(int64_t)>& next_accessor, 715 int64_t* optimal_cost); 716 717 // Same as ComputeRouteCumulCost, but the cost computed does not contain 718 // the part of the vehicle span cost due to fixed transits. 719 DimensionSchedulingStatus ComputeRouteCumulCostWithoutFixedTransits( 720 int vehicle, const std::function<int64_t(int64_t)>& next_accessor, 721 int64_t* optimal_cost_without_transits); 722 723 std::vector<DimensionSchedulingStatus> 724 ComputeRouteCumulCostsForResourcesWithoutFixedTransits( 725 int vehicle, const std::function<int64_t(int64_t)>& next_accessor, 726 const std::vector<RoutingModel::ResourceGroup::Resource>& resources, 727 const std::vector<int>& resource_indices, bool optimize_vehicle_costs, 728 std::vector<int64_t>* optimal_costs_without_transits, 729 std::vector<std::vector<int64_t>>* optimal_cumuls, 730 std::vector<std::vector<int64_t>>* optimal_breaks); 731 732 // If feasible, computes the optimal values for cumul and break variables 733 // of the route performed by a vehicle, minimizing cumul soft lower, upper 734 // bound costs and vehicle span costs, stores them in "optimal_cumuls" 735 // (if not null), and optimal_breaks, and returns true. 736 // Returns false if the route is not feasible. 737 DimensionSchedulingStatus ComputeRouteCumuls( 738 int vehicle, const std::function<int64_t(int64_t)>& next_accessor, 739 std::vector<int64_t>* optimal_cumuls, 740 std::vector<int64_t>* optimal_breaks); 741 742 // Similar to ComputeRouteCumuls, but also tries to pack the cumul values on 743 // the route, such that the cost remains the same, the cumul of route end is 744 // minimized, and then the cumul of the start of the route is maximized. 745 // If 'resource' is non-null, the packed route must also respect its start/end 746 // time window. 747 DimensionSchedulingStatus ComputePackedRouteCumuls( 748 int vehicle, const std::function<int64_t(int64_t)>& next_accessor, 749 const RoutingModel::ResourceGroup::Resource* resource, 750 std::vector<int64_t>* packed_cumuls, std::vector<int64_t>* packed_breaks); 751 dimension()752 const RoutingDimension* dimension() const { 753 return optimizer_core_.dimension(); 754 } 755 756 private: 757 std::vector<std::unique_ptr<RoutingLinearSolverWrapper>> solver_; 758 DimensionCumulOptimizerCore optimizer_core_; 759 }; 760 761 class GlobalDimensionCumulOptimizer { 762 public: 763 GlobalDimensionCumulOptimizer( 764 const RoutingDimension* dimension, 765 RoutingSearchParameters::SchedulingSolver solver_type); 766 // If feasible, computes the optimal cost of the entire model with regards to 767 // the optimizer_core_'s dimension costs, minimizing cumul soft lower/upper 768 // bound costs and vehicle/global span costs, and stores it in "optimal_cost" 769 // (if not null). 770 // Returns true iff all the constraints can be respected. 771 DimensionSchedulingStatus ComputeCumulCostWithoutFixedTransits( 772 const std::function<int64_t(int64_t)>& next_accessor, 773 int64_t* optimal_cost_without_transits); 774 // If feasible, computes the optimal values for cumul, break and resource 775 // variables, minimizing cumul soft lower/upper bound costs and vehicle/global 776 // span costs, stores them in "optimal_cumuls" (if not null), "optimal_breaks" 777 // and "optimal_resource_indices_per_group", and returns true. 778 // Returns false if the routes are not feasible. 779 DimensionSchedulingStatus ComputeCumuls( 780 const std::function<int64_t(int64_t)>& next_accessor, 781 std::vector<int64_t>* optimal_cumuls, 782 std::vector<int64_t>* optimal_breaks, 783 std::vector<std::vector<int>>* optimal_resource_indices_per_group); 784 785 // Similar to ComputeCumuls, but also tries to pack the cumul values on all 786 // routes, such that the cost remains the same, the cumuls of route ends are 787 // minimized, and then the cumuls of the starts of the routes are maximized. 788 DimensionSchedulingStatus ComputePackedCumuls( 789 const std::function<int64_t(int64_t)>& next_accessor, 790 std::vector<int64_t>* packed_cumuls, std::vector<int64_t>* packed_breaks, 791 std::vector<std::vector<int>>* resource_indices_per_group); 792 dimension()793 const RoutingDimension* dimension() const { 794 return optimizer_core_.dimension(); 795 } 796 797 private: 798 std::unique_ptr<RoutingLinearSolverWrapper> solver_; 799 DimensionCumulOptimizerCore optimizer_core_; 800 }; 801 802 class ResourceAssignmentOptimizer { 803 public: 804 ResourceAssignmentOptimizer(const RoutingModel::ResourceGroup* resource_group, 805 LocalDimensionCumulOptimizer* optimizer, 806 LocalDimensionCumulOptimizer* mp_optimizer); 807 808 // Returns the cost resulting from the min-cost assignment of resources to 809 // (used) vehicles, or -1 if the assignment is impossible. 810 // For each vehicle v and resource r, the cost of assigning r to v is equal to 811 // - primary_vehicle_to_resource_assignment_costs[v][r] if 812 // primary_vehicle_to_resource_assignment_costs[v] is not empty, 813 // - secondary_vehicle_to_resource_assignment_costs[v][r] otherwise 814 // (secondary_vehicle_to_resource_assignment_costs[v] can never be empty). 815 // If non-null, 'resource_indices' contains the index of the resource assigned 816 // to each vehicle. 817 int64_t ComputeBestAssignmentCost( 818 const std::vector<std::vector<int64_t>>& 819 primary_vehicle_to_resource_assignment_costs, 820 const std::vector<std::vector<int64_t>>& 821 secondary_vehicle_to_resource_assignment_costs, 822 const std::function<bool(int)>& use_primary_for_vehicle, 823 std::vector<int>* resource_indices) const; 824 825 // Computes the vehicle-to-resource assignment costs for the given vehicle to 826 // all resources in the group, and sets these costs in assignment_costs (if 827 // non-null). 828 // optimize_vehicle_costs indicates if the costs should be optimized or if 829 // we merely care about feasibility (cost of 0) and infeasibility (cost of -1) 830 // of the assignments. 831 // The cumul and break values corresponding to the assignment of each resource 832 // are also set in cumul_values and break_values, if non-null. 833 bool ComputeAssignmentCostsForVehicle( 834 int v, const std::function<int64_t(int64_t)>& next_accessor, 835 bool optimize_vehicle_costs, std::vector<int64_t>* assignment_costs, 836 std::vector<std::vector<int64_t>>* cumul_values, 837 std::vector<std::vector<int64_t>>* break_values); 838 dimension()839 const RoutingDimension* const dimension() const { 840 return optimizer_.dimension(); 841 } 842 843 private: 844 LocalDimensionCumulOptimizer& optimizer_; 845 LocalDimensionCumulOptimizer& mp_optimizer_; 846 const RoutingModel& model_; 847 const RoutingModel::ResourceGroup& resource_group_; 848 }; 849 850 } // namespace operations_research 851 852 #endif // OR_TOOLS_CONSTRAINT_SOLVER_ROUTING_LP_SCHEDULING_H_ 853