OR-Tools  8.2
routing_lp_scheduling.cc
Go to the documentation of this file.
1 // Copyright 2010-2018 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 
15 
16 #include <numeric>
17 
18 #include "absl/container/flat_hash_set.h"
19 #include "absl/time/time.h"
21 #include "ortools/base/map_util.h"
23 #include "ortools/glop/lp_solver.h"
25 
26 namespace operations_research {
27 
28 namespace {
29 
30 // The following sets of parameters give the fastest response time without
31 // impacting solutions found negatively.
32 glop::GlopParameters GetGlopParametersForLocalLP() {
33  glop::GlopParameters parameters;
34  parameters.set_use_dual_simplex(true);
35  parameters.set_use_preprocessing(false);
36  return parameters;
37 }
38 
39 glop::GlopParameters GetGlopParametersForGlobalLP() {
40  glop::GlopParameters parameters;
41  parameters.set_use_dual_simplex(true);
42  return parameters;
43 }
44 
45 bool GetCumulBoundsWithOffset(const RoutingDimension& dimension,
46  int64 node_index, int64 cumul_offset,
47  int64* lower_bound, int64* upper_bound) {
48  DCHECK(lower_bound != nullptr);
49  DCHECK(upper_bound != nullptr);
50 
51  const IntVar& cumul_var = *dimension.CumulVar(node_index);
52  *upper_bound = cumul_var.Max();
53  if (*upper_bound < cumul_offset) {
54  return false;
55  }
56 
57  const int64 first_after_offset =
58  std::max(dimension.GetFirstPossibleGreaterOrEqualValueForNode(
59  node_index, cumul_offset),
60  cumul_var.Min());
61  DCHECK_LT(first_after_offset, kint64max);
62  *lower_bound = CapSub(first_after_offset, cumul_offset);
63  DCHECK_GE(*lower_bound, 0);
64 
65  if (*upper_bound == kint64max) {
66  return true;
67  }
68  *upper_bound = CapSub(*upper_bound, cumul_offset);
69  DCHECK_GE(*upper_bound, *lower_bound);
70  return true;
71 }
72 
73 int64 GetFirstPossibleValueForCumulWithOffset(const RoutingDimension& dimension,
74  int64 node_index,
75  int64 lower_bound_without_offset,
76  int64 cumul_offset) {
77  return CapSub(
78  dimension.GetFirstPossibleGreaterOrEqualValueForNode(
79  node_index, CapAdd(lower_bound_without_offset, cumul_offset)),
80  cumul_offset);
81 }
82 
83 int64 GetLastPossibleValueForCumulWithOffset(const RoutingDimension& dimension,
84  int64 node_index,
85  int64 upper_bound_without_offset,
86  int64 cumul_offset) {
87  return CapSub(
88  dimension.GetLastPossibleLessOrEqualValueForNode(
89  node_index, CapAdd(upper_bound_without_offset, cumul_offset)),
90  cumul_offset);
91 }
92 
93 // Finds the pickup/delivery pairs of nodes on a given vehicle's route.
94 // Returns the vector of visited pair indices, and stores the corresponding
95 // pickup/delivery indices in visited_pickup_delivery_indices_for_pair_.
96 // NOTE: Supposes that visited_pickup_delivery_indices_for_pair is correctly
97 // sized and initialized to {-1, -1} for all pairs.
98 void StoreVisitedPickupDeliveryPairsOnRoute(
99  const RoutingDimension& dimension, int vehicle,
100  const std::function<int64(int64)>& next_accessor,
101  std::vector<int>* visited_pairs,
102  std::vector<std::pair<int64, int64>>*
103  visited_pickup_delivery_indices_for_pair) {
104  // visited_pickup_delivery_indices_for_pair must be all {-1, -1}.
105  DCHECK_EQ(visited_pickup_delivery_indices_for_pair->size(),
106  dimension.model()->GetPickupAndDeliveryPairs().size());
107  DCHECK(std::all_of(visited_pickup_delivery_indices_for_pair->begin(),
108  visited_pickup_delivery_indices_for_pair->end(),
109  [](std::pair<int64, int64> p) {
110  return p.first == -1 && p.second == -1;
111  }));
112  visited_pairs->clear();
113  if (!dimension.HasPickupToDeliveryLimits()) {
114  return;
115  }
116  const RoutingModel& model = *dimension.model();
117 
118  int64 node_index = model.Start(vehicle);
119  while (!model.IsEnd(node_index)) {
120  const std::vector<std::pair<int, int>>& pickup_index_pairs =
121  model.GetPickupIndexPairs(node_index);
122  const std::vector<std::pair<int, int>>& delivery_index_pairs =
123  model.GetDeliveryIndexPairs(node_index);
124  if (!pickup_index_pairs.empty()) {
125  // The current node is a pickup. We verify that it belongs to a single
126  // pickup index pair and that it's not a delivery, and store the index.
127  DCHECK(delivery_index_pairs.empty());
128  DCHECK_EQ(pickup_index_pairs.size(), 1);
129  (*visited_pickup_delivery_indices_for_pair)[pickup_index_pairs[0].first]
130  .first = node_index;
131  visited_pairs->push_back(pickup_index_pairs[0].first);
132  } else if (!delivery_index_pairs.empty()) {
133  // The node is a delivery. We verify that it belongs to a single
134  // delivery pair, and set the limit with its pickup if one has been
135  // visited for this pair.
136  DCHECK_EQ(delivery_index_pairs.size(), 1);
137  const int pair_index = delivery_index_pairs[0].first;
138  std::pair<int64, int64>& pickup_delivery_index =
139  (*visited_pickup_delivery_indices_for_pair)[pair_index];
140  if (pickup_delivery_index.first < 0) {
141  // This case should not happen, as a delivery must have its pickup
142  // on the route, but we ignore it here.
143  node_index = next_accessor(node_index);
144  continue;
145  }
146  pickup_delivery_index.second = node_index;
147  }
148  node_index = next_accessor(node_index);
149  }
150 }
151 
152 } // namespace
153 
155  const RoutingDimension* dimension,
156  RoutingSearchParameters::SchedulingSolver solver_type)
157  : optimizer_core_(dimension, /*use_precedence_propagator=*/false) {
158  // Using one solver per vehicle in the hope that if routes don't change this
159  // will be faster.
160  const int vehicles = dimension->model()->vehicles();
161  solver_.resize(vehicles);
162  switch (solver_type) {
163  case RoutingSearchParameters::GLOP: {
164  const glop::GlopParameters parameters = GetGlopParametersForLocalLP();
165  for (int vehicle = 0; vehicle < vehicles; ++vehicle) {
166  solver_[vehicle] = absl::make_unique<RoutingGlopWrapper>(parameters);
167  }
168  break;
169  }
170  case RoutingSearchParameters::CP_SAT: {
171  for (int vehicle = 0; vehicle < vehicles; ++vehicle) {
172  solver_[vehicle] = absl::make_unique<RoutingCPSatWrapper>();
173  }
174  break;
175  }
176  default:
177  LOG(DFATAL) << "Unrecognized solver type: " << solver_type;
178  }
179 }
180 
182  int vehicle, const std::function<int64(int64)>& next_accessor,
183  int64* optimal_cost) {
184  return optimizer_core_.OptimizeSingleRoute(vehicle, next_accessor,
185  solver_[vehicle].get(), nullptr,
186  nullptr, optimal_cost, nullptr);
187 }
188 
191  int vehicle, const std::function<int64(int64)>& next_accessor,
192  int64* optimal_cost_without_transits) {
193  int64 cost = 0;
194  int64 transit_cost = 0;
195  const DimensionSchedulingStatus status = optimizer_core_.OptimizeSingleRoute(
196  vehicle, next_accessor, solver_[vehicle].get(), nullptr, nullptr, &cost,
197  &transit_cost);
199  optimal_cost_without_transits != nullptr) {
200  *optimal_cost_without_transits = CapSub(cost, transit_cost);
201  }
202  return status;
203 }
204 
206  int vehicle, const std::function<int64(int64)>& next_accessor,
207  std::vector<int64>* optimal_cumuls, std::vector<int64>* optimal_breaks) {
208  return optimizer_core_.OptimizeSingleRoute(
209  vehicle, next_accessor, solver_[vehicle].get(), optimal_cumuls,
210  optimal_breaks, nullptr, nullptr);
211 }
212 
215  int vehicle, const std::function<int64(int64)>& next_accessor,
216  std::vector<int64>* packed_cumuls, std::vector<int64>* packed_breaks) {
217  return optimizer_core_.OptimizeAndPackSingleRoute(
218  vehicle, next_accessor, solver_[vehicle].get(), packed_cumuls,
219  packed_breaks);
220 }
221 
222 const int CumulBoundsPropagator::kNoParent = -2;
223 const int CumulBoundsPropagator::kParentToBePropagated = -1;
224 
226  : dimension_(*dimension), num_nodes_(2 * dimension->cumuls().size()) {
227  outgoing_arcs_.resize(num_nodes_);
228  node_in_queue_.resize(num_nodes_, false);
229  tree_parent_node_of_.resize(num_nodes_, kNoParent);
230  propagated_bounds_.resize(num_nodes_);
231  visited_pickup_delivery_indices_for_pair_.resize(
232  dimension->model()->GetPickupAndDeliveryPairs().size(), {-1, -1});
233 }
234 
235 void CumulBoundsPropagator::AddArcs(int first_index, int second_index,
236  int64 offset) {
237  // Add arc first_index + offset <= second_index
238  outgoing_arcs_[PositiveNode(first_index)].push_back(
239  {PositiveNode(second_index), offset});
240  AddNodeToQueue(PositiveNode(first_index));
241  // Add arc -second_index + transit <= -first_index
242  outgoing_arcs_[NegativeNode(second_index)].push_back(
243  {NegativeNode(first_index), offset});
244  AddNodeToQueue(NegativeNode(second_index));
245 }
246 
247 bool CumulBoundsPropagator::InitializeArcsAndBounds(
248  const std::function<int64(int64)>& next_accessor, int64 cumul_offset) {
249  propagated_bounds_.assign(num_nodes_, kint64min);
250 
251  for (std::vector<ArcInfo>& arcs : outgoing_arcs_) {
252  arcs.clear();
253  }
254 
255  RoutingModel* const model = dimension_.model();
256  std::vector<int64>& lower_bounds = propagated_bounds_;
257 
258  for (int vehicle = 0; vehicle < model->vehicles(); vehicle++) {
259  const std::function<int64(int64, int64)>& transit_accessor =
260  dimension_.transit_evaluator(vehicle);
261 
262  int node = model->Start(vehicle);
263  while (true) {
264  int64 cumul_lb, cumul_ub;
265  if (!GetCumulBoundsWithOffset(dimension_, node, cumul_offset, &cumul_lb,
266  &cumul_ub)) {
267  return false;
268  }
269  lower_bounds[PositiveNode(node)] = cumul_lb;
270  if (cumul_ub < kint64max) {
271  lower_bounds[NegativeNode(node)] = -cumul_ub;
272  }
273 
274  if (model->IsEnd(node)) {
275  break;
276  }
277 
278  const int next = next_accessor(node);
279  const int64 transit = transit_accessor(node, next);
280  const IntVar& slack_var = *dimension_.SlackVar(node);
281  // node + transit + slack_var == next
282  // Add arcs for node + transit + slack_min <= next
283  AddArcs(node, next, CapAdd(transit, slack_var.Min()));
284  if (slack_var.Max() < kint64max) {
285  // Add arcs for node + transit + slack_max >= next.
286  AddArcs(next, node, CapSub(-slack_var.Max(), transit));
287  }
288 
289  node = next;
290  }
291 
292  // Add vehicle span upper bound: end - span_ub <= start.
293  const int64 span_ub = dimension_.GetSpanUpperBoundForVehicle(vehicle);
294  if (span_ub < kint64max) {
295  AddArcs(model->End(vehicle), model->Start(vehicle), -span_ub);
296  }
297 
298  // Set pickup/delivery limits on route.
299  std::vector<int> visited_pairs;
300  StoreVisitedPickupDeliveryPairsOnRoute(
301  dimension_, vehicle, next_accessor, &visited_pairs,
302  &visited_pickup_delivery_indices_for_pair_);
303  for (int pair_index : visited_pairs) {
304  const int64 pickup_index =
305  visited_pickup_delivery_indices_for_pair_[pair_index].first;
306  const int64 delivery_index =
307  visited_pickup_delivery_indices_for_pair_[pair_index].second;
308  visited_pickup_delivery_indices_for_pair_[pair_index] = {-1, -1};
309 
310  DCHECK_GE(pickup_index, 0);
311  if (delivery_index < 0) {
312  // We didn't encounter a delivery for this pickup.
313  continue;
314  }
315 
316  const int64 limit = dimension_.GetPickupToDeliveryLimitForPair(
317  pair_index, model->GetPickupIndexPairs(pickup_index)[0].second,
318  model->GetDeliveryIndexPairs(delivery_index)[0].second);
319  if (limit < kint64max) {
320  // delivery_cumul - limit <= pickup_cumul.
321  AddArcs(delivery_index, pickup_index, -limit);
322  }
323  }
324  }
325 
326  for (const RoutingDimension::NodePrecedence& precedence :
327  dimension_.GetNodePrecedences()) {
328  const int first_index = precedence.first_node;
329  const int second_index = precedence.second_node;
330  if (lower_bounds[PositiveNode(first_index)] == kint64min ||
331  lower_bounds[PositiveNode(second_index)] == kint64min) {
332  // One of the nodes is unperformed, so the precedence rule doesn't apply.
333  continue;
334  }
335  AddArcs(first_index, second_index, precedence.offset);
336  }
337 
338  return true;
339 }
340 
341 bool CumulBoundsPropagator::UpdateCurrentLowerBoundOfNode(int node,
342  int64 new_lb,
343  int64 offset) {
344  const int cumul_var_index = node / 2;
345 
346  if (node == PositiveNode(cumul_var_index)) {
347  // new_lb is a lower bound of the cumul of variable 'cumul_var_index'.
348  propagated_bounds_[node] = GetFirstPossibleValueForCumulWithOffset(
349  dimension_, cumul_var_index, new_lb, offset);
350  } else {
351  // -new_lb is an upper bound of the cumul of variable 'cumul_var_index'.
352  const int64 new_ub = CapSub(0, new_lb);
353  propagated_bounds_[node] =
354  CapSub(0, GetLastPossibleValueForCumulWithOffset(
355  dimension_, cumul_var_index, new_ub, offset));
356  }
357 
358  // Test that the lower/upper bounds do not cross each other.
359  const int64 cumul_lower_bound =
360  propagated_bounds_[PositiveNode(cumul_var_index)];
361 
362  const int64 negated_cumul_upper_bound =
363  propagated_bounds_[NegativeNode(cumul_var_index)];
364 
365  return CapAdd(negated_cumul_upper_bound, cumul_lower_bound) <= 0;
366 }
367 
368 bool CumulBoundsPropagator::DisassembleSubtree(int source, int target) {
369  tmp_dfs_stack_.clear();
370  tmp_dfs_stack_.push_back(source);
371  while (!tmp_dfs_stack_.empty()) {
372  const int tail = tmp_dfs_stack_.back();
373  tmp_dfs_stack_.pop_back();
374  for (const ArcInfo& arc : outgoing_arcs_[tail]) {
375  const int child_node = arc.head;
376  if (tree_parent_node_of_[child_node] != tail) continue;
377  if (child_node == target) return false;
378  tree_parent_node_of_[child_node] = kParentToBePropagated;
379  tmp_dfs_stack_.push_back(child_node);
380  }
381  }
382  return true;
383 }
384 
386  const std::function<int64(int64)>& next_accessor, int64 cumul_offset) {
387  tree_parent_node_of_.assign(num_nodes_, kNoParent);
388  DCHECK(std::none_of(node_in_queue_.begin(), node_in_queue_.end(),
389  [](bool b) { return b; }));
390  DCHECK(bf_queue_.empty());
391 
392  if (!InitializeArcsAndBounds(next_accessor, cumul_offset)) {
393  return CleanupAndReturnFalse();
394  }
395 
396  std::vector<int64>& current_lb = propagated_bounds_;
397 
398  // Bellman-Ford-Tarjan algorithm.
399  while (!bf_queue_.empty()) {
400  const int node = bf_queue_.front();
401  bf_queue_.pop_front();
402  node_in_queue_[node] = false;
403 
404  if (tree_parent_node_of_[node] == kParentToBePropagated) {
405  // The parent of this node is still in the queue, so no need to process
406  // node now, since it will be re-enqued when its parent is processed.
407  continue;
408  }
409 
410  const int64 lower_bound = current_lb[node];
411  for (const ArcInfo& arc : outgoing_arcs_[node]) {
412  // NOTE: kint64min as a lower bound means no lower bound at all, so we
413  // don't use this value to propagate.
414  const int64 induced_lb = (lower_bound == kint64min)
415  ? kint64min
416  : CapAdd(lower_bound, arc.offset);
417 
418  const int head_node = arc.head;
419  if (induced_lb <= current_lb[head_node]) {
420  // No update necessary for the head_node, continue to next children of
421  // node.
422  continue;
423  }
424  if (!UpdateCurrentLowerBoundOfNode(head_node, induced_lb, cumul_offset) ||
425  !DisassembleSubtree(head_node, node)) {
426  // The new lower bound is infeasible, or a positive cycle was detected
427  // in the precedence graph by DisassembleSubtree().
428  return CleanupAndReturnFalse();
429  }
430 
431  tree_parent_node_of_[head_node] = node;
432  AddNodeToQueue(head_node);
433  }
434  }
435  return true;
436 }
437 
439  const RoutingDimension* dimension, bool use_precedence_propagator)
440  : dimension_(dimension),
441  visited_pickup_delivery_indices_for_pair_(
442  dimension->model()->GetPickupAndDeliveryPairs().size(), {-1, -1}) {
443  if (use_precedence_propagator) {
444  propagator_ = absl::make_unique<CumulBoundsPropagator>(dimension);
445  }
446  if (dimension_->HasBreakConstraints()) {
447  // Initialize vehicle_to_first_index_ so the variables of the breaks of
448  // vehicle v are stored from vehicle_to_first_index_[v] to
449  // vehicle_to_first_index_[v+1] - 1.
450  const int num_vehicles = dimension_->model()->vehicles();
451  vehicle_to_all_break_variables_offset_.reserve(num_vehicles);
452  int num_break_vars = 0;
453  for (int vehicle = 0; vehicle < num_vehicles; ++vehicle) {
454  vehicle_to_all_break_variables_offset_.push_back(num_break_vars);
455  const auto& intervals = dimension_->GetBreakIntervalsOfVehicle(vehicle);
456  num_break_vars += 2 * intervals.size(); // 2 variables per break.
457  }
458  all_break_variables_.resize(num_break_vars, -1);
459  }
460 }
461 
463  int vehicle, const std::function<int64(int64)>& next_accessor,
464  RoutingLinearSolverWrapper* solver, std::vector<int64>* cumul_values,
465  std::vector<int64>* break_values, int64* cost, int64* transit_cost,
466  bool clear_lp) {
467  InitOptimizer(solver);
468  // Make sure SetRouteCumulConstraints will properly set the cumul bounds by
469  // looking at this route only.
470  DCHECK(propagator_ == nullptr);
471 
472  const RoutingModel* const model = dimension()->model();
473  const bool optimize_vehicle_costs =
474  (cumul_values != nullptr || cost != nullptr) &&
475  (!model->IsEnd(next_accessor(model->Start(vehicle))) ||
476  model->AreEmptyRouteCostsConsideredForVehicle(vehicle));
477  const int64 cumul_offset =
478  dimension_->GetLocalOptimizerOffsetForVehicle(vehicle);
479  int64 cost_offset = 0;
480  if (!SetRouteCumulConstraints(vehicle, next_accessor, cumul_offset,
481  optimize_vehicle_costs, solver, transit_cost,
482  &cost_offset)) {
484  }
485  const DimensionSchedulingStatus status =
486  solver->Solve(model->RemainingTime());
488  return status;
489  }
490 
491  SetValuesFromLP(current_route_cumul_variables_, cumul_offset, solver,
492  cumul_values);
493  SetValuesFromLP(current_route_break_variables_, cumul_offset, solver,
494  break_values);
495  if (cost != nullptr) {
496  *cost = CapAdd(cost_offset, solver->GetObjectiveValue());
497  }
498 
499  if (clear_lp) {
500  solver->Clear();
501  }
502  return status;
503 }
504 
506  const std::function<int64(int64)>& next_accessor,
507  RoutingLinearSolverWrapper* solver, std::vector<int64>* cumul_values,
508  std::vector<int64>* break_values, int64* cost, int64* transit_cost,
509  bool clear_lp) {
510  InitOptimizer(solver);
511 
512  // If both "cumul_values" and "cost" parameters are null, we don't try to
513  // optimize the cost and stop at the first feasible solution.
514  const bool optimize_costs = (cumul_values != nullptr) || (cost != nullptr);
515  bool has_vehicles_being_optimized = false;
516 
517  const int64 cumul_offset = dimension_->GetGlobalOptimizerOffset();
518 
519  if (propagator_ != nullptr &&
520  !propagator_->PropagateCumulBounds(next_accessor, cumul_offset)) {
521  return false;
522  }
523 
524  int64 total_transit_cost = 0;
525  int64 total_cost_offset = 0;
526  const RoutingModel* model = dimension()->model();
527  for (int vehicle = 0; vehicle < model->vehicles(); vehicle++) {
528  int64 route_transit_cost = 0;
529  int64 route_cost_offset = 0;
530  const bool optimize_vehicle_costs =
531  optimize_costs &&
532  (!model->IsEnd(next_accessor(model->Start(vehicle))) ||
533  model->AreEmptyRouteCostsConsideredForVehicle(vehicle));
534  if (!SetRouteCumulConstraints(vehicle, next_accessor, cumul_offset,
535  optimize_vehicle_costs, solver,
536  &route_transit_cost, &route_cost_offset)) {
537  return false;
538  }
539  total_transit_cost = CapAdd(total_transit_cost, route_transit_cost);
540  total_cost_offset = CapAdd(total_cost_offset, route_cost_offset);
541  has_vehicles_being_optimized |= optimize_vehicle_costs;
542  }
543  if (transit_cost != nullptr) {
544  *transit_cost = total_transit_cost;
545  }
546 
547  SetGlobalConstraints(has_vehicles_being_optimized, solver);
548 
549  if (solver->Solve(model->RemainingTime()) ==
551  return false;
552  }
553 
554  SetValuesFromLP(index_to_cumul_variable_, cumul_offset, solver, cumul_values);
555  SetValuesFromLP(all_break_variables_, cumul_offset, solver, break_values);
556 
557  if (cost != nullptr) {
558  *cost = CapAdd(solver->GetObjectiveValue(), total_cost_offset);
559  }
560 
561  if (clear_lp) {
562  solver->Clear();
563  }
564  return true;
565 }
566 
568  const std::function<int64(int64)>& next_accessor,
569  RoutingLinearSolverWrapper* solver, std::vector<int64>* cumul_values,
570  std::vector<int64>* break_values) {
571  // Note: We pass a non-nullptr cost to the Optimize() method so the costs are
572  // optimized by the LP.
573  int64 cost = 0;
574  if (!Optimize(next_accessor, solver,
575  /*cumul_values=*/nullptr, /*break_values=*/nullptr, &cost,
576  /*transit_cost=*/nullptr, /*clear_lp=*/false)) {
577  return false;
578  }
579 
580  std::vector<int> vehicles(dimension()->model()->vehicles());
581  std::iota(vehicles.begin(), vehicles.end(), 0);
582  if (PackRoutes(std::move(vehicles), solver) ==
584  return false;
585  }
586  const int64 global_offset = dimension_->GetGlobalOptimizerOffset();
587  SetValuesFromLP(index_to_cumul_variable_, global_offset, solver,
588  cumul_values);
589  SetValuesFromLP(all_break_variables_, global_offset, solver, break_values);
590  solver->Clear();
591  return true;
592 }
593 
596  int vehicle, const std::function<int64(int64)>& next_accessor,
597  RoutingLinearSolverWrapper* solver, std::vector<int64>* cumul_values,
598  std::vector<int64>* break_values) {
599  // Note: We pass a non-nullptr cost to the OptimizeSingleRoute() method so the
600  // costs are optimized by the LP.
601  int64 cost = 0;
602  if (OptimizeSingleRoute(vehicle, next_accessor, solver,
603  /*cumul_values=*/nullptr, /*break_values=*/nullptr,
604  &cost, /*transit_cost=*/nullptr,
605  /*clear_lp=*/false) ==
608  }
609  const DimensionSchedulingStatus status = PackRoutes({vehicle}, solver);
612  }
613  const int64 local_offset =
614  dimension_->GetLocalOptimizerOffsetForVehicle(vehicle);
615  SetValuesFromLP(current_route_cumul_variables_, local_offset, solver,
616  cumul_values);
617  SetValuesFromLP(current_route_break_variables_, local_offset, solver,
618  break_values);
619  solver->Clear();
620  return status;
621 }
622 
623 DimensionSchedulingStatus DimensionCumulOptimizerCore::PackRoutes(
624  std::vector<int> vehicles, RoutingLinearSolverWrapper* solver) {
625  const RoutingModel* model = dimension_->model();
626 
627  // NOTE(user): Given our constraint matrix, our problem *should* always
628  // have an integer optimal solution, in which case we can round to the nearest
629  // integer both for the objective constraint bound (returned by
630  // GetObjectiveValue()) and the end cumul variable bound after minimizing
631  // (see b/154381899 showcasing an example where std::ceil leads to an
632  // "imperfect" packing due to rounding precision errors).
633  // If this DCHECK ever fails, it can be removed but the code below should be
634  // adapted to have a 2-phase approach, solving once with the rounded value as
635  // bound and if this fails, solve again using std::ceil.
636  DCHECK(solver->SolutionIsInteger());
637 
638  // Minimize the route end times without increasing the cost.
639  const int objective_ct =
640  solver->CreateNewConstraint(0, solver->GetObjectiveValue());
641 
642  for (int variable = 0; variable < solver->NumVariables(); variable++) {
643  const double coefficient = solver->GetObjectiveCoefficient(variable);
644  if (coefficient != 0) {
645  solver->SetCoefficient(objective_ct, variable, coefficient);
646  }
647  }
648  solver->ClearObjective();
649  for (int vehicle : vehicles) {
650  solver->SetObjectiveCoefficient(
651  index_to_cumul_variable_[model->End(vehicle)], 1);
652  }
653 
654  if (solver->Solve(model->RemainingTime()) ==
657  }
658 
659  // Maximize the route start times without increasing the cost or the route end
660  // times.
661  solver->ClearObjective();
662  for (int vehicle : vehicles) {
663  const int end_cumul_var = index_to_cumul_variable_[model->End(vehicle)];
664  // end_cumul_var <= solver.GetValue(end_cumul_var)
665  solver->SetVariableBounds(
666  end_cumul_var, solver->GetVariableLowerBound(end_cumul_var),
667  MathUtil::FastInt64Round(solver->GetValue(end_cumul_var)));
668 
669  // Maximize the starts of the routes.
670  solver->SetObjectiveCoefficient(
671  index_to_cumul_variable_[model->Start(vehicle)], -1);
672  }
673  return solver->Solve(model->RemainingTime());
674 }
675 
676 void DimensionCumulOptimizerCore::InitOptimizer(
677  RoutingLinearSolverWrapper* solver) {
678  solver->Clear();
679  index_to_cumul_variable_.assign(dimension_->cumuls().size(), -1);
680  max_end_cumul_ = solver->CreateNewPositiveVariable();
681  min_start_cumul_ = solver->CreateNewPositiveVariable();
682 }
683 
684 bool DimensionCumulOptimizerCore::ComputeRouteCumulBounds(
685  const std::vector<int64>& route, const std::vector<int64>& fixed_transits,
686  int64 cumul_offset) {
687  const int route_size = route.size();
688  current_route_min_cumuls_.resize(route_size);
689  current_route_max_cumuls_.resize(route_size);
690  if (propagator_ != nullptr) {
691  for (int pos = 0; pos < route_size; pos++) {
692  const int64 node = route[pos];
693  current_route_min_cumuls_[pos] = propagator_->CumulMin(node);
694  DCHECK_GE(current_route_min_cumuls_[pos], 0);
695  current_route_max_cumuls_[pos] = propagator_->CumulMax(node);
696  DCHECK_GE(current_route_max_cumuls_[pos], current_route_min_cumuls_[pos]);
697  }
698  return true;
699  }
700 
701  // Extract cumul min/max and fixed transits from CP.
702  for (int pos = 0; pos < route_size; ++pos) {
703  if (!GetCumulBoundsWithOffset(*dimension_, route[pos], cumul_offset,
704  &current_route_min_cumuls_[pos],
705  &current_route_max_cumuls_[pos])) {
706  return false;
707  }
708  }
709 
710  // Refine cumul bounds using
711  // cumul[i+1] >= cumul[i] + fixed_transit[i] + slack[i].
712  for (int pos = 1; pos < route_size; ++pos) {
713  const int64 slack_min = dimension_->SlackVar(route[pos - 1])->Min();
714  current_route_min_cumuls_[pos] = std::max(
715  current_route_min_cumuls_[pos],
716  CapAdd(
717  CapAdd(current_route_min_cumuls_[pos - 1], fixed_transits[pos - 1]),
718  slack_min));
719  current_route_min_cumuls_[pos] = GetFirstPossibleValueForCumulWithOffset(
720  *dimension_, route[pos], current_route_min_cumuls_[pos], cumul_offset);
721  if (current_route_min_cumuls_[pos] > current_route_max_cumuls_[pos]) {
722  return false;
723  }
724  }
725 
726  for (int pos = route_size - 2; pos >= 0; --pos) {
727  // If cumul_max[pos+1] is kint64max, it will be translated to
728  // double +infinity, so it must not constrain cumul_max[pos].
729  if (current_route_max_cumuls_[pos + 1] < kint64max) {
730  const int64 slack_min = dimension_->SlackVar(route[pos])->Min();
731  current_route_max_cumuls_[pos] = std::min(
732  current_route_max_cumuls_[pos],
733  CapSub(
734  CapSub(current_route_max_cumuls_[pos + 1], fixed_transits[pos]),
735  slack_min));
736  current_route_max_cumuls_[pos] = GetLastPossibleValueForCumulWithOffset(
737  *dimension_, route[pos], current_route_max_cumuls_[pos],
738  cumul_offset);
739  if (current_route_max_cumuls_[pos] < current_route_min_cumuls_[pos]) {
740  return false;
741  }
742  }
743  }
744  return true;
745 }
746 
747 bool DimensionCumulOptimizerCore::SetRouteCumulConstraints(
748  int vehicle, const std::function<int64(int64)>& next_accessor,
749  int64 cumul_offset, bool optimize_costs, RoutingLinearSolverWrapper* solver,
750  int64* route_transit_cost, int64* route_cost_offset) {
751  RoutingModel* const model = dimension_->model();
752  // Extract the vehicle's path from next_accessor.
753  std::vector<int64> path;
754  {
755  int node = model->Start(vehicle);
756  path.push_back(node);
757  while (!model->IsEnd(node)) {
758  node = next_accessor(node);
759  path.push_back(node);
760  }
761  DCHECK_GE(path.size(), 2);
762  }
763  const int path_size = path.size();
764 
765  std::vector<int64> fixed_transit(path_size - 1);
766  {
767  const std::function<int64(int64, int64)>& transit_accessor =
768  dimension_->transit_evaluator(vehicle);
769  for (int pos = 1; pos < path_size; ++pos) {
770  fixed_transit[pos - 1] = transit_accessor(path[pos - 1], path[pos]);
771  }
772  }
773 
774  if (!ComputeRouteCumulBounds(path, fixed_transit, cumul_offset)) {
775  return false;
776  }
777 
778  // LP Model variables, current_route_cumul_variables_ and lp_slacks.
779  // Create LP variables for cumuls.
780  std::vector<int>& lp_cumuls = current_route_cumul_variables_;
781  lp_cumuls.assign(path_size, -1);
782  for (int pos = 0; pos < path_size; ++pos) {
783  const int lp_cumul = solver->CreateNewPositiveVariable();
784  index_to_cumul_variable_[path[pos]] = lp_cumul;
785  lp_cumuls[pos] = lp_cumul;
786  if (!solver->SetVariableBounds(lp_cumul, current_route_min_cumuls_[pos],
787  current_route_max_cumuls_[pos])) {
788  return false;
789  }
790  const SortedDisjointIntervalList& forbidden =
791  dimension_->forbidden_intervals()[path[pos]];
792  if (forbidden.NumIntervals() > 0) {
793  std::vector<int64> starts;
794  std::vector<int64> ends;
795  for (const ClosedInterval interval :
796  dimension_->GetAllowedIntervalsInRange(
797  path[pos], CapAdd(current_route_min_cumuls_[pos], cumul_offset),
798  CapAdd(current_route_max_cumuls_[pos], cumul_offset))) {
799  starts.push_back(CapSub(interval.start, cumul_offset));
800  ends.push_back(CapSub(interval.end, cumul_offset));
801  }
802  solver->SetVariableDisjointBounds(lp_cumul, starts, ends);
803  }
804  }
805  // Create LP variables for slacks.
806  std::vector<int> lp_slacks(path_size - 1, -1);
807  for (int pos = 0; pos < path_size - 1; ++pos) {
808  const IntVar* cp_slack = dimension_->SlackVar(path[pos]);
809  lp_slacks[pos] = solver->CreateNewPositiveVariable();
810  if (!solver->SetVariableBounds(lp_slacks[pos], cp_slack->Min(),
811  cp_slack->Max())) {
812  return false;
813  }
814  }
815 
816  // LP Model constraints and costs.
817  // Add all path constraints to LP:
818  // cumul[i] + fixed_transit[i] + slack[i] == cumul[i+1]
819  // <=> fixed_transit[i] == cumul[i+1] - cumul[i] - slack[i].
820  for (int pos = 0; pos < path_size - 1; ++pos) {
821  const int ct =
822  solver->CreateNewConstraint(fixed_transit[pos], fixed_transit[pos]);
823  solver->SetCoefficient(ct, lp_cumuls[pos + 1], 1);
824  solver->SetCoefficient(ct, lp_cumuls[pos], -1);
825  solver->SetCoefficient(ct, lp_slacks[pos], -1);
826  }
827  if (route_cost_offset != nullptr) *route_cost_offset = 0;
828  if (optimize_costs) {
829  // Add soft upper bounds.
830  for (int pos = 0; pos < path_size; ++pos) {
831  if (!dimension_->HasCumulVarSoftUpperBound(path[pos])) continue;
832  const int64 coef =
833  dimension_->GetCumulVarSoftUpperBoundCoefficient(path[pos]);
834  if (coef == 0) continue;
835  int64 bound = dimension_->GetCumulVarSoftUpperBound(path[pos]);
836  if (bound < cumul_offset && route_cost_offset != nullptr) {
837  // Add coef * (cumul_offset - bound) to the cost offset.
838  *route_cost_offset = CapAdd(*route_cost_offset,
839  CapProd(CapSub(cumul_offset, bound), coef));
840  }
841  bound = std::max<int64>(0, CapSub(bound, cumul_offset));
842  if (current_route_max_cumuls_[pos] <= bound) {
843  // constraint is never violated.
844  continue;
845  }
846  const int soft_ub_diff = solver->CreateNewPositiveVariable();
847  solver->SetObjectiveCoefficient(soft_ub_diff, coef);
848  // cumul - soft_ub_diff <= bound.
849  const int ct = solver->CreateNewConstraint(kint64min, bound);
850  solver->SetCoefficient(ct, lp_cumuls[pos], 1);
851  solver->SetCoefficient(ct, soft_ub_diff, -1);
852  }
853  // Add soft lower bounds.
854  for (int pos = 0; pos < path_size; ++pos) {
855  if (!dimension_->HasCumulVarSoftLowerBound(path[pos])) continue;
856  const int64 coef =
857  dimension_->GetCumulVarSoftLowerBoundCoefficient(path[pos]);
858  if (coef == 0) continue;
859  const int64 bound = std::max<int64>(
860  0, CapSub(dimension_->GetCumulVarSoftLowerBound(path[pos]),
861  cumul_offset));
862  if (current_route_min_cumuls_[pos] >= bound) {
863  // constraint is never violated.
864  continue;
865  }
866  const int soft_lb_diff = solver->CreateNewPositiveVariable();
867  solver->SetObjectiveCoefficient(soft_lb_diff, coef);
868  // bound - cumul <= soft_lb_diff
869  const int ct = solver->CreateNewConstraint(bound, kint64max);
870  solver->SetCoefficient(ct, lp_cumuls[pos], 1);
871  solver->SetCoefficient(ct, soft_lb_diff, 1);
872  }
873  }
874  // Add pickup and delivery limits.
875  std::vector<int> visited_pairs;
876  StoreVisitedPickupDeliveryPairsOnRoute(
877  *dimension_, vehicle, next_accessor, &visited_pairs,
878  &visited_pickup_delivery_indices_for_pair_);
879  for (int pair_index : visited_pairs) {
880  const int64 pickup_index =
881  visited_pickup_delivery_indices_for_pair_[pair_index].first;
882  const int64 delivery_index =
883  visited_pickup_delivery_indices_for_pair_[pair_index].second;
884  visited_pickup_delivery_indices_for_pair_[pair_index] = {-1, -1};
885 
886  DCHECK_GE(pickup_index, 0);
887  if (delivery_index < 0) {
888  // We didn't encounter a delivery for this pickup.
889  continue;
890  }
891 
892  const int64 limit = dimension_->GetPickupToDeliveryLimitForPair(
893  pair_index, model->GetPickupIndexPairs(pickup_index)[0].second,
894  model->GetDeliveryIndexPairs(delivery_index)[0].second);
895  if (limit < kint64max) {
896  // delivery_cumul - pickup_cumul <= limit.
897  const int ct = solver->CreateNewConstraint(kint64min, limit);
898  solver->SetCoefficient(ct, index_to_cumul_variable_[delivery_index], 1);
899  solver->SetCoefficient(ct, index_to_cumul_variable_[pickup_index], -1);
900  }
901  }
902 
903  // Add span bound constraint.
904  const int64 span_bound = dimension_->GetSpanUpperBoundForVehicle(vehicle);
905  if (span_bound < kint64max) {
906  // end_cumul - start_cumul <= bound
907  const int ct = solver->CreateNewConstraint(kint64min, span_bound);
908  solver->SetCoefficient(ct, lp_cumuls.back(), 1);
909  solver->SetCoefficient(ct, lp_cumuls.front(), -1);
910  }
911  // Add span cost.
912  const int64 span_cost_coef =
913  dimension_->GetSpanCostCoefficientForVehicle(vehicle);
914  if (optimize_costs && span_cost_coef > 0) {
915  solver->SetObjectiveCoefficient(lp_cumuls.back(), span_cost_coef);
916  solver->SetObjectiveCoefficient(lp_cumuls.front(), -span_cost_coef);
917  }
918  // Add soft span cost.
919  if (optimize_costs && dimension_->HasSoftSpanUpperBounds()) {
920  SimpleBoundCosts::BoundCost bound_cost =
921  dimension_->GetSoftSpanUpperBoundForVehicle(vehicle);
922  if (bound_cost.bound < kint64max && bound_cost.cost > 0) {
923  const int span_violation = solver->CreateNewPositiveVariable();
924  // end - start <= bound + span_violation
925  const int violation =
926  solver->CreateNewConstraint(kint64min, bound_cost.bound);
927  solver->SetCoefficient(violation, lp_cumuls.back(), 1.0);
928  solver->SetCoefficient(violation, lp_cumuls.front(), -1.0);
929  solver->SetCoefficient(violation, span_violation, -1.0);
930  // Add span_violation * cost to objective.
931  solver->SetObjectiveCoefficient(span_violation, bound_cost.cost);
932  }
933  }
934  // Add global span constraint.
935  if (optimize_costs && dimension_->global_span_cost_coefficient() > 0) {
936  // min_start_cumul_ <= cumuls[start]
937  int ct = solver->CreateNewConstraint(kint64min, 0);
938  solver->SetCoefficient(ct, min_start_cumul_, 1);
939  solver->SetCoefficient(ct, lp_cumuls.front(), -1);
940  // max_end_cumul_ >= cumuls[end]
941  ct = solver->CreateNewConstraint(0, kint64max);
942  solver->SetCoefficient(ct, max_end_cumul_, 1);
943  solver->SetCoefficient(ct, lp_cumuls.back(), -1);
944  }
945  // Fill transit cost if specified.
946  if (route_transit_cost != nullptr) {
947  if (optimize_costs && span_cost_coef > 0) {
948  const int64 total_fixed_transit = std::accumulate(
949  fixed_transit.begin(), fixed_transit.end(), 0, CapAdd);
950  *route_transit_cost = CapProd(total_fixed_transit, span_cost_coef);
951  } else {
952  *route_transit_cost = 0;
953  }
954  }
955  // For every break that must be inside the route, the duration of that break
956  // must be flowed in the slacks of arcs that can intersect the break.
957  // This LP modelization is correct but not complete:
958  // can miss some cases where the breaks cannot fit.
959  // TODO(user): remove the need for returns in the code below.
960  current_route_break_variables_.clear();
961  if (!dimension_->HasBreakConstraints()) return true;
962  const std::vector<IntervalVar*>& breaks =
963  dimension_->GetBreakIntervalsOfVehicle(vehicle);
964  const int num_breaks = breaks.size();
965  // When there are no breaks, only break distance needs to be modeled,
966  // and it reduces to a span maximum.
967  // TODO(user): Also add the case where no breaks can intersect the route.
968  if (num_breaks == 0) {
969  int64 maximum_route_span = kint64max;
970  for (const auto& distance_duration :
971  dimension_->GetBreakDistanceDurationOfVehicle(vehicle)) {
972  maximum_route_span =
973  std::min(maximum_route_span, distance_duration.first);
974  }
975  if (maximum_route_span < kint64max) {
976  const int ct = solver->CreateNewConstraint(kint64min, maximum_route_span);
977  solver->SetCoefficient(ct, lp_cumuls.back(), 1);
978  solver->SetCoefficient(ct, lp_cumuls.front(), -1);
979  }
980  return true;
981  }
982  // Gather visit information: the visit of node i has [start, end) =
983  // [cumul[i] - post_travel[i-1], cumul[i] + pre_travel[i]).
984  // Breaks cannot overlap those visit intervals.
985  std::vector<int64> pre_travel(path_size - 1, 0);
986  std::vector<int64> post_travel(path_size - 1, 0);
987  {
988  const int pre_travel_index =
989  dimension_->GetPreTravelEvaluatorOfVehicle(vehicle);
990  if (pre_travel_index != -1) {
991  FillPathEvaluation(path, model->TransitCallback(pre_travel_index),
992  &pre_travel);
993  }
994  const int post_travel_index =
995  dimension_->GetPostTravelEvaluatorOfVehicle(vehicle);
996  if (post_travel_index != -1) {
997  FillPathEvaluation(path, model->TransitCallback(post_travel_index),
998  &post_travel);
999  }
1000  }
1001  // If the solver is CPSAT, it will need to represent the times at which
1002  // breaks are scheduled, those variables are used both in the pure breaks
1003  // part and in the break distance part of the model.
1004  // Otherwise, it doesn't need the variables and they are not created.
1005  std::vector<int> lp_break_start;
1006  std::vector<int> lp_break_duration;
1007  std::vector<int> lp_break_end;
1008  if (solver->IsCPSATSolver()) {
1009  lp_break_start.resize(num_breaks, -1);
1010  lp_break_duration.resize(num_breaks, -1);
1011  lp_break_end.resize(num_breaks, -1);
1012  }
1013 
1014  std::vector<int> slack_exact_lower_bound_ct(path_size - 1, -1);
1015  std::vector<int> slack_linear_lower_bound_ct(path_size - 1, -1);
1016 
1017  const int64 vehicle_start_min = current_route_min_cumuls_.front();
1018  const int64 vehicle_start_max = current_route_max_cumuls_.front();
1019  const int64 vehicle_end_min = current_route_min_cumuls_.back();
1020  const int64 vehicle_end_max = current_route_max_cumuls_.back();
1021  const int all_break_variables_offset =
1022  vehicle_to_all_break_variables_offset_[vehicle];
1023  for (int br = 0; br < num_breaks; ++br) {
1024  const IntervalVar& break_var = *breaks[br];
1025  if (!break_var.MustBePerformed()) continue;
1026  const int64 break_start_min = CapSub(break_var.StartMin(), cumul_offset);
1027  const int64 break_start_max = CapSub(break_var.StartMax(), cumul_offset);
1028  const int64 break_end_min = CapSub(break_var.EndMin(), cumul_offset);
1029  const int64 break_end_max = CapSub(break_var.EndMax(), cumul_offset);
1030  const int64 break_duration_min = break_var.DurationMin();
1031  const int64 break_duration_max = break_var.DurationMax();
1032  // The CPSAT solver encodes all breaks that can intersect the route,
1033  // the LP solver only encodes the breaks that must intersect the route.
1034  if (solver->IsCPSATSolver()) {
1035  if (break_end_max <= vehicle_start_min ||
1036  vehicle_end_max <= break_start_min) {
1037  all_break_variables_[all_break_variables_offset + 2 * br] = -1;
1038  all_break_variables_[all_break_variables_offset + 2 * br + 1] = -1;
1039  current_route_break_variables_.push_back(-1);
1040  current_route_break_variables_.push_back(-1);
1041  continue;
1042  }
1043  lp_break_start[br] =
1044  solver->AddVariable(break_start_min, break_start_max);
1045  lp_break_end[br] = solver->AddVariable(break_end_min, break_end_max);
1046  lp_break_duration[br] =
1047  solver->AddVariable(break_duration_min, break_duration_max);
1048  // start + duration = end.
1049  solver->AddLinearConstraint(0, 0,
1050  {{lp_break_end[br], 1},
1051  {lp_break_start[br], -1},
1052  {lp_break_duration[br], -1}});
1053  // Record index of variables
1054  all_break_variables_[all_break_variables_offset + 2 * br] =
1055  lp_break_start[br];
1056  all_break_variables_[all_break_variables_offset + 2 * br + 1] =
1057  lp_break_end[br];
1058  current_route_break_variables_.push_back(lp_break_start[br]);
1059  current_route_break_variables_.push_back(lp_break_end[br]);
1060  } else {
1061  if (break_end_min <= vehicle_start_max ||
1062  vehicle_end_min <= break_start_max) {
1063  continue;
1064  }
1065  }
1066 
1067  // Create a constraint for every break, that forces it to be scheduled
1068  // in exactly one place, i.e. one slack or before/after the route.
1069  // sum_i break_in_slack_i == 1.
1070  const int break_in_one_slack_ct = solver->CreateNewConstraint(1, 1);
1071 
1072  if (solver->IsCPSATSolver()) {
1073  // Break can be before route.
1074  if (break_end_min <= vehicle_start_max) {
1075  const int ct = solver->AddLinearConstraint(
1076  0, kint64max, {{lp_cumuls.front(), 1}, {lp_break_end[br], -1}});
1077  const int break_is_before_route = solver->AddVariable(0, 1);
1078  solver->SetEnforcementLiteral(ct, break_is_before_route);
1079  solver->SetCoefficient(break_in_one_slack_ct, break_is_before_route, 1);
1080  }
1081  // Break can be after route.
1082  if (vehicle_end_min <= break_start_max) {
1083  const int ct = solver->AddLinearConstraint(
1084  0, kint64max, {{lp_break_start[br], 1}, {lp_cumuls.back(), -1}});
1085  const int break_is_after_route = solver->AddVariable(0, 1);
1086  solver->SetEnforcementLiteral(ct, break_is_after_route);
1087  solver->SetCoefficient(break_in_one_slack_ct, break_is_after_route, 1);
1088  }
1089  }
1090 
1091  // Add the possibility of fitting the break during each slack where it can.
1092  for (int pos = 0; pos < path_size - 1; ++pos) {
1093  // Pass on slacks that cannot start before, cannot end after,
1094  // or are not long enough to contain the break.
1095  const int64 slack_start_min =
1096  CapAdd(current_route_min_cumuls_[pos], pre_travel[pos]);
1097  if (slack_start_min > break_start_max) break;
1098  const int64 slack_end_max =
1099  CapSub(current_route_max_cumuls_[pos + 1], post_travel[pos]);
1100  if (break_end_min > slack_end_max) continue;
1101  const int64 slack_duration_max =
1102  std::min(CapSub(CapSub(current_route_max_cumuls_[pos + 1],
1103  current_route_min_cumuls_[pos]),
1104  fixed_transit[pos]),
1105  dimension_->SlackVar(path[pos])->Max());
1106  if (slack_duration_max < break_duration_min) continue;
1107 
1108  // Break can fit into slack: make LP variable, add to break and slack
1109  // constraints.
1110  // Make a linearized slack lower bound (lazily), that represents
1111  // sum_br break_duration_min(br) * break_in_slack(br, pos) <=
1112  // lp_slacks(pos).
1113  const int break_in_slack = solver->AddVariable(0, 1);
1114  solver->SetCoefficient(break_in_one_slack_ct, break_in_slack, 1);
1115  if (slack_linear_lower_bound_ct[pos] == -1) {
1116  slack_linear_lower_bound_ct[pos] =
1117  solver->AddLinearConstraint(kint64min, 0, {{lp_slacks[pos], -1}});
1118  }
1119  solver->SetCoefficient(slack_linear_lower_bound_ct[pos], break_in_slack,
1120  break_duration_min);
1121  if (solver->IsCPSATSolver()) {
1122  // Exact relation between breaks, slacks and cumul variables.
1123  // Make an exact slack lower bound (lazily), that represents
1124  // sum_br break_duration(br) * break_in_slack(br, pos) <=
1125  // lp_slacks(pos).
1126  const int break_duration_in_slack =
1127  solver->AddVariable(0, slack_duration_max);
1128  solver->AddProductConstraint(break_duration_in_slack,
1129  {break_in_slack, lp_break_duration[br]});
1130  if (slack_exact_lower_bound_ct[pos] == -1) {
1131  slack_exact_lower_bound_ct[pos] =
1132  solver->AddLinearConstraint(kint64min, 0, {{lp_slacks[pos], -1}});
1133  }
1134  solver->SetCoefficient(slack_exact_lower_bound_ct[pos],
1135  break_duration_in_slack, 1);
1136  // If break_in_slack_i == 1, then
1137  // 1) break_start >= cumul[pos] + pre_travel[pos]
1138  const int break_start_after_current_ct = solver->AddLinearConstraint(
1139  pre_travel[pos], kint64max,
1140  {{lp_break_start[br], 1}, {lp_cumuls[pos], -1}});
1141  solver->SetEnforcementLiteral(break_start_after_current_ct,
1142  break_in_slack);
1143  // 2) break_end <= cumul[pos+1] - post_travel[pos]
1144  const int break_ends_before_next_ct = solver->AddLinearConstraint(
1145  post_travel[pos], kint64max,
1146  {{lp_cumuls[pos + 1], 1}, {lp_break_end[br], -1}});
1147  solver->SetEnforcementLiteral(break_ends_before_next_ct,
1148  break_in_slack);
1149  }
1150  }
1151  }
1152 
1153  if (!solver->IsCPSATSolver()) return true;
1154  if (!dimension_->GetBreakDistanceDurationOfVehicle(vehicle).empty()) {
1155  // If there is an optional interval, the following model would be wrong.
1156  // TODO(user): support optional intervals.
1157  for (const IntervalVar* interval :
1158  dimension_->GetBreakIntervalsOfVehicle(vehicle)) {
1159  if (!interval->MustBePerformed()) return true;
1160  }
1161  // When this feature is used, breaks are in sorted order.
1162  for (int br = 1; br < num_breaks; ++br) {
1163  solver->AddLinearConstraint(
1164  0, kint64max, {{lp_break_end[br - 1], -1}, {lp_break_start[br], 1}});
1165  }
1166  }
1167  for (const auto& distance_duration :
1168  dimension_->GetBreakDistanceDurationOfVehicle(vehicle)) {
1169  const int64 limit = distance_duration.first;
1170  const int64 min_break_duration = distance_duration.second;
1171  // Interbreak limit constraint: breaks are interpreted as being in sorted
1172  // order, and the maximum duration between two consecutive
1173  // breaks of duration more than 'min_break_duration' is 'limit'. This
1174  // considers the time until start of route and after end of route to be
1175  // infinite breaks.
1176  // The model for this constraint adds some 'cover_i' variables, such that
1177  // the breaks up to i and the start of route allows to go without a break.
1178  // With s_i the start of break i and e_i its end:
1179  // - the route start covers time from start to start + limit:
1180  // cover_0 = route_start + limit
1181  // - the coverage up to a given break is the largest of the coverage of the
1182  // previous break and if the break is long enough, break end + limit:
1183  // cover_{i+1} = max(cover_i,
1184  // e_i - s_i >= min_break_duration ? e_i + limit : -inf)
1185  // - the coverage of the last break must be at least the route end,
1186  // to ensure the time point route_end-1 is covered:
1187  // cover_{num_breaks} >= route_end
1188  // - similarly, time point s_i-1 must be covered by breaks up to i-1,
1189  // but only if the cover has not reached the route end.
1190  // For instance, a vehicle could have a choice between two days,
1191  // with a potential break on day 1 and a potential break on day 2,
1192  // but the break of day 1 does not have to cover that of day 2!
1193  // cover_{i-1} < route_end => s_i <= cover_{i-1}
1194  // This is sufficient to ensure that the union of the intervals
1195  // (-infinity, route_start], [route_end, +infinity) and all
1196  // [s_i, e_i+limit) where e_i - s_i >= min_break_duration is
1197  // the whole timeline (-infinity, +infinity).
1198  int previous_cover = solver->AddVariable(CapAdd(vehicle_start_min, limit),
1199  CapAdd(vehicle_start_max, limit));
1200  solver->AddLinearConstraint(limit, limit,
1201  {{previous_cover, 1}, {lp_cumuls.front(), -1}});
1202  for (int br = 0; br < num_breaks; ++br) {
1203  if (lp_break_start[br] == -1) continue;
1204  const int64 break_end_min = CapSub(breaks[br]->EndMin(), cumul_offset);
1205  const int64 break_end_max = CapSub(breaks[br]->EndMax(), cumul_offset);
1206  // break_is_eligible <=>
1207  // break_end - break_start >= break_minimum_duration.
1208  const int break_is_eligible = solver->AddVariable(0, 1);
1209  const int break_is_not_eligible = solver->AddVariable(0, 1);
1210  {
1211  solver->AddLinearConstraint(
1212  1, 1, {{break_is_eligible, 1}, {break_is_not_eligible, 1}});
1213  const int positive_ct = solver->AddLinearConstraint(
1214  min_break_duration, kint64max,
1215  {{lp_break_end[br], 1}, {lp_break_start[br], -1}});
1216  solver->SetEnforcementLiteral(positive_ct, break_is_eligible);
1217  const int negative_ct = solver->AddLinearConstraint(
1218  kint64min, min_break_duration - 1,
1219  {{lp_break_end[br], 1}, {lp_break_start[br], -1}});
1220  solver->SetEnforcementLiteral(negative_ct, break_is_not_eligible);
1221  }
1222  // break_is_eligible => break_cover == break_end + limit.
1223  // break_is_not_eligible => break_cover == vehicle_start_min + limit.
1224  // break_cover's initial domain is the smallest interval that contains the
1225  // union of sets {vehicle_start_min+limit} and
1226  // [break_end_min+limit, break_end_max+limit).
1227  const int break_cover = solver->AddVariable(
1228  CapAdd(std::min(vehicle_start_min, break_end_min), limit),
1229  CapAdd(std::max(vehicle_start_min, break_end_max), limit));
1230  const int limit_cover_ct = solver->AddLinearConstraint(
1231  limit, limit, {{break_cover, 1}, {lp_break_end[br], -1}});
1232  solver->SetEnforcementLiteral(limit_cover_ct, break_is_eligible);
1233  const int empty_cover_ct = solver->AddLinearConstraint(
1234  CapAdd(vehicle_start_min, limit), CapAdd(vehicle_start_min, limit),
1235  {{break_cover, 1}});
1236  solver->SetEnforcementLiteral(empty_cover_ct, break_is_not_eligible);
1237 
1238  const int cover =
1239  solver->AddVariable(CapAdd(vehicle_start_min, limit), kint64max);
1240  solver->AddMaximumConstraint(cover, {previous_cover, break_cover});
1241  // Cover chaining. If route end is not covered, break start must be:
1242  // cover_{i-1} < route_end => s_i <= cover_{i-1}
1243  const int route_end_is_not_covered = solver->AddReifiedLinearConstraint(
1244  1, kint64max, {{lp_cumuls.back(), 1}, {previous_cover, -1}});
1245  const int break_start_cover_ct = solver->AddLinearConstraint(
1246  0, kint64max, {{previous_cover, 1}, {lp_break_start[br], -1}});
1247  solver->SetEnforcementLiteral(break_start_cover_ct,
1248  route_end_is_not_covered);
1249 
1250  previous_cover = cover;
1251  }
1252  solver->AddLinearConstraint(0, kint64max,
1253  {{previous_cover, 1}, {lp_cumuls.back(), -1}});
1254  }
1255 
1256  return true;
1257 }
1258 
1259 void DimensionCumulOptimizerCore::SetGlobalConstraints(
1260  bool optimize_costs, RoutingLinearSolverWrapper* solver) {
1261  // Global span cost =
1262  // global_span_cost_coefficient * (max_end_cumul - min_start_cumul).
1263  const int64 global_span_coeff = dimension_->global_span_cost_coefficient();
1264  if (optimize_costs && global_span_coeff > 0) {
1265  solver->SetObjectiveCoefficient(max_end_cumul_, global_span_coeff);
1266  solver->SetObjectiveCoefficient(min_start_cumul_, -global_span_coeff);
1267  }
1268 
1269  // Node precedence constraints, set when both nodes are visited.
1270  for (const RoutingDimension::NodePrecedence& precedence :
1271  dimension_->GetNodePrecedences()) {
1272  const int first_cumul_var = index_to_cumul_variable_[precedence.first_node];
1273  const int second_cumul_var =
1274  index_to_cumul_variable_[precedence.second_node];
1275  if (first_cumul_var < 0 || second_cumul_var < 0) {
1276  // At least one of the nodes is not on any route, skip this precedence
1277  // constraint.
1278  continue;
1279  }
1280  DCHECK_NE(first_cumul_var, second_cumul_var)
1281  << "Dimension " << dimension_->name()
1282  << " has a self-precedence on node " << precedence.first_node << ".";
1283 
1284  // cumul[second_node] - cumul[first_node] >= offset.
1285  const int ct = solver->CreateNewConstraint(precedence.offset, kint64max);
1286  solver->SetCoefficient(ct, second_cumul_var, 1);
1287  solver->SetCoefficient(ct, first_cumul_var, -1);
1288  }
1289 }
1290 
1291 void DimensionCumulOptimizerCore::SetValuesFromLP(
1292  const std::vector<int>& lp_variables, int64 offset,
1293  RoutingLinearSolverWrapper* solver, std::vector<int64>* lp_values) {
1294  if (lp_values == nullptr) return;
1295  lp_values->assign(lp_variables.size(), kint64min);
1296  for (int i = 0; i < lp_variables.size(); i++) {
1297  const int cumul_var = lp_variables[i];
1298  if (cumul_var < 0) continue; // Keep default value, kint64min.
1299  const double lp_value_double = solver->GetValue(cumul_var);
1300  const int64 lp_value_int64 =
1301  (lp_value_double >= kint64max)
1302  ? kint64max
1303  : MathUtil::FastInt64Round(lp_value_double);
1304  (*lp_values)[i] = CapAdd(lp_value_int64, offset);
1305  }
1306 }
1307 
1309  const RoutingDimension* dimension)
1310  : optimizer_core_(dimension,
1311  /*use_precedence_propagator=*/
1312  !dimension->GetNodePrecedences().empty()) {
1313  solver_ =
1314  absl::make_unique<RoutingGlopWrapper>(GetGlopParametersForGlobalLP());
1315 }
1316 
1318  const std::function<int64(int64)>& next_accessor,
1319  int64* optimal_cost_without_transits) {
1320  int64 cost = 0;
1321  int64 transit_cost = 0;
1322  if (optimizer_core_.Optimize(next_accessor, solver_.get(), nullptr, nullptr,
1323  &cost, &transit_cost)) {
1324  if (optimal_cost_without_transits != nullptr) {
1325  *optimal_cost_without_transits = CapSub(cost, transit_cost);
1326  }
1327  return true;
1328  }
1329  return false;
1330 }
1331 
1333  const std::function<int64(int64)>& next_accessor,
1334  std::vector<int64>* optimal_cumuls, std::vector<int64>* optimal_breaks) {
1335  return optimizer_core_.Optimize(next_accessor, solver_.get(), optimal_cumuls,
1336  optimal_breaks, nullptr, nullptr);
1337 }
1338 
1340  const std::function<int64(int64)>& next_accessor) {
1341  return optimizer_core_.Optimize(next_accessor, solver_.get(), nullptr,
1342  nullptr, nullptr, nullptr);
1343 }
1344 
1346  const std::function<int64(int64)>& next_accessor,
1347  std::vector<int64>* packed_cumuls, std::vector<int64>* packed_breaks) {
1348  return optimizer_core_.OptimizeAndPack(next_accessor, solver_.get(),
1349  packed_cumuls, packed_breaks);
1350 }
1351 } // namespace operations_research
int64 min
Definition: alldiff_cst.cc:138
int64 max
Definition: alldiff_cst.cc:139
#define DCHECK_NE(val1, val2)
Definition: base/logging.h:886
#define DCHECK_GE(val1, val2)
Definition: base/logging.h:889
#define DCHECK_LT(val1, val2)
Definition: base/logging.h:888
#define LOG(severity)
Definition: base/logging.h:420
#define DCHECK(condition)
Definition: base/logging.h:884
#define DCHECK_EQ(val1, val2)
Definition: base/logging.h:885
bool PropagateCumulBounds(const std::function< int64(int64)> &next_accessor, int64 cumul_offset)
const RoutingDimension & dimension() const
CumulBoundsPropagator(const RoutingDimension *dimension)
bool OptimizeAndPack(const std::function< int64(int64)> &next_accessor, RoutingLinearSolverWrapper *solver, std::vector< int64 > *cumul_values, std::vector< int64 > *break_values)
DimensionCumulOptimizerCore(const RoutingDimension *dimension, bool use_precedence_propagator)
bool Optimize(const std::function< int64(int64)> &next_accessor, RoutingLinearSolverWrapper *solver, std::vector< int64 > *cumul_values, std::vector< int64 > *break_values, int64 *cost, int64 *transit_cost, bool clear_lp=true)
DimensionSchedulingStatus OptimizeAndPackSingleRoute(int vehicle, const std::function< int64(int64)> &next_accessor, RoutingLinearSolverWrapper *solver, std::vector< int64 > *cumul_values, std::vector< int64 > *break_values)
DimensionSchedulingStatus OptimizeSingleRoute(int vehicle, const std::function< int64(int64)> &next_accessor, RoutingLinearSolverWrapper *solver, std::vector< int64 > *cumul_values, std::vector< int64 > *break_values, int64 *cost, int64 *transit_cost, bool clear_lp=true)
bool ComputeCumuls(const std::function< int64(int64)> &next_accessor, std::vector< int64 > *optimal_cumuls, std::vector< int64 > *optimal_breaks)
bool IsFeasible(const std::function< int64(int64)> &next_accessor)
GlobalDimensionCumulOptimizer(const RoutingDimension *dimension)
bool ComputeCumulCostWithoutFixedTransits(const std::function< int64(int64)> &next_accessor, int64 *optimal_cost_without_transits)
bool ComputePackedCumuls(const std::function< int64(int64)> &next_accessor, std::vector< int64 > *packed_cumuls, std::vector< int64 > *packed_breaks)
virtual int64 Max() const =0
virtual int64 Min() const =0
DimensionSchedulingStatus ComputePackedRouteCumuls(int vehicle, const std::function< int64(int64)> &next_accessor, std::vector< int64 > *packed_cumuls, std::vector< int64 > *packed_breaks)
LocalDimensionCumulOptimizer(const RoutingDimension *dimension, RoutingSearchParameters::SchedulingSolver solver_type)
DimensionSchedulingStatus ComputeRouteCumuls(int vehicle, const std::function< int64(int64)> &next_accessor, std::vector< int64 > *optimal_cumuls, std::vector< int64 > *optimal_breaks)
DimensionSchedulingStatus ComputeRouteCumulCost(int vehicle, const std::function< int64(int64)> &next_accessor, int64 *optimal_cost)
DimensionSchedulingStatus ComputeRouteCumulCostWithoutFixedTransits(int vehicle, const std::function< int64(int64)> &next_accessor, int64 *optimal_cost_without_transits)
static int64 FastInt64Round(double x)
Definition: mathutil.h:138
Dimensions represent quantities accumulated at nodes along the routes.
Definition: routing.h:2368
SimpleBoundCosts::BoundCost GetSoftSpanUpperBoundForVehicle(int vehicle) const
Definition: routing.h:2710
const std::vector< IntVar * > & cumuls() const
Like CumulVar(), TransitVar(), SlackVar() but return the whole variable vectors instead (indexed by i...
Definition: routing.h:2394
int64 global_span_cost_coefficient() const
Definition: routing.h:2681
const std::vector< std::pair< int64, int64 > > & GetBreakDistanceDurationOfVehicle(int vehicle) const
Returns the pairs (distance, duration) specified by break distance constraints.
Definition: routing.cc:6939
int64 GetSpanCostCoefficientForVehicle(int vehicle) const
Definition: routing.h:2673
RoutingModel * model() const
Returns the model on which the dimension was created.
Definition: routing.h:2372
int64 GetCumulVarSoftLowerBoundCoefficient(int64 index) const
Returns the cost coefficient of the soft lower bound of a cumul variable for a given variable index.
Definition: routing.cc:6751
bool HasCumulVarSoftLowerBound(int64 index) const
Returns true if a soft lower bound has been set for a given variable index.
Definition: routing.cc:6738
int64 GetCumulVarSoftLowerBound(int64 index) const
Returns the soft lower bound of a cumul variable for a given variable index.
Definition: routing.cc:6743
bool HasBreakConstraints() const
Returns true if any break interval or break distance was defined.
Definition: routing.cc:6900
int64 GetLocalOptimizerOffsetForVehicle(int vehicle) const
Definition: routing.h:2689
int GetPreTravelEvaluatorOfVehicle(int vehicle) const
!defined(SWIGPYTHON)
Definition: routing.cc:6911
const std::vector< IntervalVar * > & GetBreakIntervalsOfVehicle(int vehicle) const
Returns the break intervals set by SetBreakIntervalsOfVehicle().
Definition: routing.cc:6904
int64 GetPickupToDeliveryLimitForPair(int pair_index, int pickup, int delivery) const
Definition: routing.cc:6959
int64 GetCumulVarSoftUpperBound(int64 index) const
Returns the soft upper bound of a cumul variable for a given variable index.
Definition: routing.cc:6692
const RoutingModel::TransitCallback2 & transit_evaluator(int vehicle) const
Returns the callback evaluating the transit value between two node indices for a given vehicle.
Definition: routing.h:2449
IntVar * SlackVar(int64 index) const
Definition: routing.h:2389
SortedDisjointIntervalList GetAllowedIntervalsInRange(int64 index, int64 min_value, int64 max_value) const
Returns allowed intervals for a given node in a given interval.
Definition: routing.cc:6563
const std::string & name() const
Returns the name of the dimension.
Definition: routing.h:2619
const std::vector< NodePrecedence > & GetNodePrecedences() const
Definition: routing.h:2656
bool HasCumulVarSoftUpperBound(int64 index) const
Returns true if a soft upper bound has been set for a given variable index.
Definition: routing.cc:6687
int64 GetCumulVarSoftUpperBoundCoefficient(int64 index) const
Returns the cost coefficient of the soft upper bound of a cumul variable for a given variable index.
Definition: routing.cc:6700
int GetPostTravelEvaluatorOfVehicle(int vehicle) const
Definition: routing.cc:6917
const std::vector< SortedDisjointIntervalList > & forbidden_intervals() const
Returns forbidden intervals for each node.
Definition: routing.h:2400
int64 GetSpanUpperBoundForVehicle(int vehicle) const
Definition: routing.h:2665
virtual void SetCoefficient(int ct, int index, double coefficient)=0
virtual int CreateNewConstraint(int64 lower_bound, int64 upper_bound)=0
virtual double GetObjectiveCoefficient(int index) const =0
virtual int64 GetVariableLowerBound(int index) const =0
virtual bool SetVariableBounds(int index, int64 lower_bound, int64 upper_bound)=0
virtual double GetValue(int index) const =0
virtual DimensionSchedulingStatus Solve(absl::Duration duration_limit)=0
virtual void SetObjectiveCoefficient(int index, double coefficient)=0
const IndexPairs & GetPickupAndDeliveryPairs() const
Returns pickup and delivery pairs currently in the model.
Definition: routing.h:746
int vehicles() const
Returns the number of vehicle routes in the model.
Definition: routing.h:1345
Block * next
SatParameters parameters
const Constraint * ct
int64 coef
Definition: expr_array.cc:1859
GRBmodel * model
static const int64 kint64max
int64_t int64
static const int64 kint64min
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
int64 CapSub(int64 x, int64 y)
int64 CapAdd(int64 x, int64 y)
int64 CapProd(int64 x, int64 y)
void FillPathEvaluation(const std::vector< int64 > &path, const RoutingModel::TransitCallback2 &evaluator, std::vector< int64 > *values)
Definition: routing.cc:6215
IntervalVar * interval
Definition: resource.cc:98
int64 tail
int64 cost
int64 bound
int64 coefficient
std::vector< double > lower_bounds