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"
25
26namespace operations_research {
27
28namespace {
29
30// The following sets of parameters give the fastest response time without
31// impacting solutions found negatively.
32glop::GlopParameters GetGlopParametersForLocalLP() {
33 glop::GlopParameters parameters;
34 parameters.set_use_dual_simplex(true);
35 parameters.set_use_preprocessing(false);
36 return parameters;
37}
38
39glop::GlopParameters GetGlopParametersForGlobalLP() {
40 glop::GlopParameters parameters;
41 parameters.set_use_dual_simplex(true);
42 return parameters;
43}
44
45bool 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
73int64 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
83int64 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.
98void 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
222const int CumulBoundsPropagator::kNoParent = -2;
223const 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
235void 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
247bool 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
341bool 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
368bool 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
623DimensionSchedulingStatus 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) {
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.
671 index_to_cumul_variable_[model->Start(vehicle)], -1);
672 }
673 return solver->Solve(model->RemainingTime());
674}
675
676void 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
684bool 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
747bool 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
1259void 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
1291void 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
const RoutingDimension & dimension() const
bool PropagateCumulBounds(const std::function< int64(int64)> &next_accessor, int64 cumul_offset)
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
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
const std::string & name() const
Returns the name of the dimension.
Definition: routing.h:2619
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 GetSpanCostCoefficientForVehicle(int vehicle) const
Definition: routing.h:2673
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
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
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
IntVar * SlackVar(int64 index) const
Definition: routing.h:2389
int64 GetLocalOptimizerOffsetForVehicle(int vehicle) const
Definition: routing.h:2689
int GetPreTravelEvaluatorOfVehicle(int vehicle) const
!defined(SWIGPYTHON)
Definition: routing.cc:6911
RoutingModel * model() const
Returns the model on which the dimension was created.
Definition: routing.h:2372
const std::vector< IntervalVar * > & GetBreakIntervalsOfVehicle(int vehicle) const
Returns the break intervals set by SetBreakIntervalsOfVehicle().
Definition: routing.cc:6904
const std::vector< NodePrecedence > & GetNodePrecedences() const
Definition: routing.h:2656
const std::vector< SortedDisjointIntervalList > & forbidden_intervals() const
Returns forbidden intervals for each node.
Definition: routing.h:2400
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
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
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
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 CapAdd(int64 x, int64 y)
int64 CapProd(int64 x, int64 y)
int64 CapSub(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