1 // Copyright 2010-2021 Google LLC
2 // Licensed under the Apache License, Version 2.0 (the "License");
3 // you may not use this file except in compliance with the License.
4 // You may obtain a copy of the License at
5 //
6 // http://www.apache.org/licenses/LICENSE-2.0
7 //
8 // Unless required by applicable law or agreed to in writing, software
9 // distributed under the License is distributed on an "AS IS" BASIS,
10 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11 // See the License for the specific language governing permissions and
12 // limitations under the License.
13
14 #include "ortools/constraint_solver/routing_lp_scheduling.h"
15
16 #include <algorithm>
17 #include <cstdint>
18 #include <deque>
19 #include <functional>
20 #include <limits>
21 #include <memory>
22 #include <numeric>
23 #include <utility>
24 #include <vector>
25
26 #include "absl/memory/memory.h"
27 #include "absl/strings/str_join.h"
28 #include "absl/time/time.h"
29 #include "ortools/base/logging.h"
30 #include "ortools/base/mathutil.h"
31 #include "ortools/constraint_solver/constraint_solver.h"
32 #include "ortools/constraint_solver/routing.h"
33 #include "ortools/constraint_solver/routing_parameters.pb.h"
34 #include "ortools/glop/parameters.pb.h"
35 #include "ortools/graph/min_cost_flow.h"
36 #include "ortools/util/saturated_arithmetic.h"
37 #include "ortools/util/sorted_interval_list.h"
38
39 namespace operations_research {
40
41 namespace {
42
43 // The following sets of parameters give the fastest response time without
44 // impacting solutions found negatively.
GetGlopParametersForLocalLP()45 glop::GlopParameters GetGlopParametersForLocalLP() {
46 glop::GlopParameters parameters;
47 parameters.set_use_dual_simplex(true);
48 parameters.set_use_preprocessing(false);
49 return parameters;
50 }
51
GetGlopParametersForGlobalLP()52 glop::GlopParameters GetGlopParametersForGlobalLP() {
53 glop::GlopParameters parameters;
54 parameters.set_use_dual_simplex(true);
55 return parameters;
56 }
57
GetCumulBoundsWithOffset(const RoutingDimension & dimension,int64_t node_index,int64_t cumul_offset,int64_t * lower_bound,int64_t * upper_bound)58 bool GetCumulBoundsWithOffset(const RoutingDimension& dimension,
59 int64_t node_index, int64_t cumul_offset,
60 int64_t* lower_bound, int64_t* upper_bound) {
61 DCHECK(lower_bound != nullptr);
62 DCHECK(upper_bound != nullptr);
63
64 const IntVar& cumul_var = *dimension.CumulVar(node_index);
65 *upper_bound = cumul_var.Max();
66 if (*upper_bound < cumul_offset) {
67 return false;
68 }
69
70 const int64_t first_after_offset =
71 std::max(dimension.GetFirstPossibleGreaterOrEqualValueForNode(
72 node_index, cumul_offset),
73 cumul_var.Min());
74 DCHECK_LT(first_after_offset, std::numeric_limits<int64_t>::max());
75 *lower_bound = CapSub(first_after_offset, cumul_offset);
76 DCHECK_GE(*lower_bound, 0);
77
78 if (*upper_bound == std::numeric_limits<int64_t>::max()) {
79 return true;
80 }
81 *upper_bound = CapSub(*upper_bound, cumul_offset);
82 DCHECK_GE(*upper_bound, *lower_bound);
83 return true;
84 }
85
GetFirstPossibleValueForCumulWithOffset(const RoutingDimension & dimension,int64_t node_index,int64_t lower_bound_without_offset,int64_t cumul_offset)86 int64_t GetFirstPossibleValueForCumulWithOffset(
87 const RoutingDimension& dimension, int64_t node_index,
88 int64_t lower_bound_without_offset, int64_t cumul_offset) {
89 return CapSub(
90 dimension.GetFirstPossibleGreaterOrEqualValueForNode(
91 node_index, CapAdd(lower_bound_without_offset, cumul_offset)),
92 cumul_offset);
93 }
94
GetLastPossibleValueForCumulWithOffset(const RoutingDimension & dimension,int64_t node_index,int64_t upper_bound_without_offset,int64_t cumul_offset)95 int64_t GetLastPossibleValueForCumulWithOffset(
96 const RoutingDimension& dimension, int64_t node_index,
97 int64_t upper_bound_without_offset, int64_t cumul_offset) {
98 return CapSub(
99 dimension.GetLastPossibleLessOrEqualValueForNode(
100 node_index, CapAdd(upper_bound_without_offset, cumul_offset)),
101 cumul_offset);
102 }
103
104 // Finds the pickup/delivery pairs of nodes on a given vehicle's route.
105 // Returns the vector of visited pair indices, and stores the corresponding
106 // pickup/delivery indices in visited_pickup_delivery_indices_for_pair_.
107 // NOTE: Supposes that visited_pickup_delivery_indices_for_pair is correctly
108 // sized and initialized to {-1, -1} for all pairs.
StoreVisitedPickupDeliveryPairsOnRoute(const RoutingDimension & dimension,int vehicle,const std::function<int64_t (int64_t)> & next_accessor,std::vector<int> * visited_pairs,std::vector<std::pair<int64_t,int64_t>> * visited_pickup_delivery_indices_for_pair)109 void StoreVisitedPickupDeliveryPairsOnRoute(
110 const RoutingDimension& dimension, int vehicle,
111 const std::function<int64_t(int64_t)>& next_accessor,
112 std::vector<int>* visited_pairs,
113 std::vector<std::pair<int64_t, int64_t>>*
114 visited_pickup_delivery_indices_for_pair) {
115 // visited_pickup_delivery_indices_for_pair must be all {-1, -1}.
116 DCHECK_EQ(visited_pickup_delivery_indices_for_pair->size(),
117 dimension.model()->GetPickupAndDeliveryPairs().size());
118 DCHECK(std::all_of(visited_pickup_delivery_indices_for_pair->begin(),
119 visited_pickup_delivery_indices_for_pair->end(),
120 [](std::pair<int64_t, int64_t> p) {
121 return p.first == -1 && p.second == -1;
122 }));
123 visited_pairs->clear();
124 if (!dimension.HasPickupToDeliveryLimits()) {
125 return;
126 }
127 const RoutingModel& model = *dimension.model();
128
129 int64_t node_index = model.Start(vehicle);
130 while (!model.IsEnd(node_index)) {
131 const std::vector<std::pair<int, int>>& pickup_index_pairs =
132 model.GetPickupIndexPairs(node_index);
133 const std::vector<std::pair<int, int>>& delivery_index_pairs =
134 model.GetDeliveryIndexPairs(node_index);
135 if (!pickup_index_pairs.empty()) {
136 // The current node is a pickup. We verify that it belongs to a single
137 // pickup index pair and that it's not a delivery, and store the index.
138 DCHECK(delivery_index_pairs.empty());
139 DCHECK_EQ(pickup_index_pairs.size(), 1);
140 (*visited_pickup_delivery_indices_for_pair)[pickup_index_pairs[0].first]
141 .first = node_index;
142 visited_pairs->push_back(pickup_index_pairs[0].first);
143 } else if (!delivery_index_pairs.empty()) {
144 // The node is a delivery. We verify that it belongs to a single
145 // delivery pair, and set the limit with its pickup if one has been
146 // visited for this pair.
147 DCHECK_EQ(delivery_index_pairs.size(), 1);
148 const int pair_index = delivery_index_pairs[0].first;
149 std::pair<int64_t, int64_t>& pickup_delivery_index =
150 (*visited_pickup_delivery_indices_for_pair)[pair_index];
151 if (pickup_delivery_index.first < 0) {
152 // This case should not happen, as a delivery must have its pickup
153 // on the route, but we ignore it here.
154 node_index = next_accessor(node_index);
155 continue;
156 }
157 pickup_delivery_index.second = node_index;
158 }
159 node_index = next_accessor(node_index);
160 }
161 }
162
163 } // namespace
164
165 // LocalDimensionCumulOptimizer
166
LocalDimensionCumulOptimizer(const RoutingDimension * dimension,RoutingSearchParameters::SchedulingSolver solver_type)167 LocalDimensionCumulOptimizer::LocalDimensionCumulOptimizer(
168 const RoutingDimension* dimension,
169 RoutingSearchParameters::SchedulingSolver solver_type)
170 : optimizer_core_(dimension, /*use_precedence_propagator=*/false) {
171 // Using one solver per vehicle in the hope that if routes don't change this
172 // will be faster.
173 const int vehicles = dimension->model()->vehicles();
174 solver_.resize(vehicles);
175 switch (solver_type) {
176 case RoutingSearchParameters::GLOP: {
177 const glop::GlopParameters parameters = GetGlopParametersForLocalLP();
178 for (int vehicle = 0; vehicle < vehicles; ++vehicle) {
179 // TODO(user): Instead of passing false, detect if the relaxation
180 // will always violate the MIPL constraints.
181 solver_[vehicle] =
182 absl::make_unique<RoutingGlopWrapper>(false, parameters);
183 }
184 break;
185 }
186 case RoutingSearchParameters::CP_SAT: {
187 for (int vehicle = 0; vehicle < vehicles; ++vehicle) {
188 solver_[vehicle] = absl::make_unique<RoutingCPSatWrapper>();
189 }
190 break;
191 }
192 default:
193 LOG(DFATAL) << "Unrecognized solver type: " << solver_type;
194 }
195 }
196
ComputeRouteCumulCost(int vehicle,const std::function<int64_t (int64_t)> & next_accessor,int64_t * optimal_cost)197 DimensionSchedulingStatus LocalDimensionCumulOptimizer::ComputeRouteCumulCost(
198 int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
199 int64_t* optimal_cost) {
200 return optimizer_core_.OptimizeSingleRoute(vehicle, next_accessor,
201 solver_[vehicle].get(), nullptr,
202 nullptr, optimal_cost, nullptr);
203 }
204
205 DimensionSchedulingStatus
ComputeRouteCumulCostWithoutFixedTransits(int vehicle,const std::function<int64_t (int64_t)> & next_accessor,int64_t * optimal_cost_without_transits)206 LocalDimensionCumulOptimizer::ComputeRouteCumulCostWithoutFixedTransits(
207 int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
208 int64_t* optimal_cost_without_transits) {
209 int64_t cost = 0;
210 int64_t transit_cost = 0;
211 const DimensionSchedulingStatus status = optimizer_core_.OptimizeSingleRoute(
212 vehicle, next_accessor, solver_[vehicle].get(), nullptr, nullptr, &cost,
213 &transit_cost);
214 if (status != DimensionSchedulingStatus::INFEASIBLE &&
215 optimal_cost_without_transits != nullptr) {
216 *optimal_cost_without_transits = CapSub(cost, transit_cost);
217 }
218 return status;
219 }
220
221 std::vector<DimensionSchedulingStatus> LocalDimensionCumulOptimizer::
ComputeRouteCumulCostsForResourcesWithoutFixedTransits(int vehicle,const std::function<int64_t (int64_t)> & next_accessor,const std::vector<RoutingModel::ResourceGroup::Resource> & resources,const std::vector<int> & resource_indices,bool optimize_vehicle_costs,std::vector<int64_t> * optimal_costs_without_transits,std::vector<std::vector<int64_t>> * optimal_cumuls,std::vector<std::vector<int64_t>> * optimal_breaks)222 ComputeRouteCumulCostsForResourcesWithoutFixedTransits(
223 int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
224 const std::vector<RoutingModel::ResourceGroup::Resource>& resources,
225 const std::vector<int>& resource_indices, bool optimize_vehicle_costs,
226 std::vector<int64_t>* optimal_costs_without_transits,
227 std::vector<std::vector<int64_t>>* optimal_cumuls,
228 std::vector<std::vector<int64_t>>* optimal_breaks) {
229 return optimizer_core_.OptimizeSingleRouteWithResources(
230 vehicle, next_accessor, resources, resource_indices,
231 optimize_vehicle_costs, solver_[vehicle].get(),
232 optimal_costs_without_transits, optimal_cumuls, optimal_breaks);
233 }
234
ComputeRouteCumuls(int vehicle,const std::function<int64_t (int64_t)> & next_accessor,std::vector<int64_t> * optimal_cumuls,std::vector<int64_t> * optimal_breaks)235 DimensionSchedulingStatus LocalDimensionCumulOptimizer::ComputeRouteCumuls(
236 int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
237 std::vector<int64_t>* optimal_cumuls,
238 std::vector<int64_t>* optimal_breaks) {
239 return optimizer_core_.OptimizeSingleRoute(
240 vehicle, next_accessor, solver_[vehicle].get(), optimal_cumuls,
241 optimal_breaks, nullptr, nullptr);
242 }
243
244 DimensionSchedulingStatus
ComputePackedRouteCumuls(int vehicle,const std::function<int64_t (int64_t)> & next_accessor,const RoutingModel::ResourceGroup::Resource * resource,std::vector<int64_t> * packed_cumuls,std::vector<int64_t> * packed_breaks)245 LocalDimensionCumulOptimizer::ComputePackedRouteCumuls(
246 int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
247 const RoutingModel::ResourceGroup::Resource* resource,
248 std::vector<int64_t>* packed_cumuls, std::vector<int64_t>* packed_breaks) {
249 return optimizer_core_.OptimizeAndPackSingleRoute(
250 vehicle, next_accessor, resource, solver_[vehicle].get(), packed_cumuls,
251 packed_breaks);
252 }
253
254 const int CumulBoundsPropagator::kNoParent = -2;
255 const int CumulBoundsPropagator::kParentToBePropagated = -1;
256
CumulBoundsPropagator(const RoutingDimension * dimension)257 CumulBoundsPropagator::CumulBoundsPropagator(const RoutingDimension* dimension)
258 : dimension_(*dimension), num_nodes_(2 * dimension->cumuls().size()) {
259 outgoing_arcs_.resize(num_nodes_);
260 node_in_queue_.resize(num_nodes_, false);
261 tree_parent_node_of_.resize(num_nodes_, kNoParent);
262 propagated_bounds_.resize(num_nodes_);
263 visited_pickup_delivery_indices_for_pair_.resize(
264 dimension->model()->GetPickupAndDeliveryPairs().size(), {-1, -1});
265 }
266
AddArcs(int first_index,int second_index,int64_t offset)267 void CumulBoundsPropagator::AddArcs(int first_index, int second_index,
268 int64_t offset) {
269 // Add arc first_index + offset <= second_index
270 outgoing_arcs_[PositiveNode(first_index)].push_back(
271 {PositiveNode(second_index), offset});
272 AddNodeToQueue(PositiveNode(first_index));
273 // Add arc -second_index + transit <= -first_index
274 outgoing_arcs_[NegativeNode(second_index)].push_back(
275 {NegativeNode(first_index), offset});
276 AddNodeToQueue(NegativeNode(second_index));
277 }
278
InitializeArcsAndBounds(const std::function<int64_t (int64_t)> & next_accessor,int64_t cumul_offset)279 bool CumulBoundsPropagator::InitializeArcsAndBounds(
280 const std::function<int64_t(int64_t)>& next_accessor,
281 int64_t cumul_offset) {
282 propagated_bounds_.assign(num_nodes_, std::numeric_limits<int64_t>::min());
283
284 for (std::vector<ArcInfo>& arcs : outgoing_arcs_) {
285 arcs.clear();
286 }
287
288 RoutingModel* const model = dimension_.model();
289 std::vector<int64_t>& lower_bounds = propagated_bounds_;
290
291 for (int vehicle = 0; vehicle < model->vehicles(); vehicle++) {
292 const std::function<int64_t(int64_t, int64_t)>& transit_accessor =
293 dimension_.transit_evaluator(vehicle);
294
295 int node = model->Start(vehicle);
296 while (true) {
297 int64_t cumul_lb, cumul_ub;
298 if (!GetCumulBoundsWithOffset(dimension_, node, cumul_offset, &cumul_lb,
299 &cumul_ub)) {
300 return false;
301 }
302 lower_bounds[PositiveNode(node)] = cumul_lb;
303 if (cumul_ub < std::numeric_limits<int64_t>::max()) {
304 lower_bounds[NegativeNode(node)] = -cumul_ub;
305 }
306
307 if (model->IsEnd(node)) {
308 break;
309 }
310
311 const int next = next_accessor(node);
312 const int64_t transit = transit_accessor(node, next);
313 const IntVar& slack_var = *dimension_.SlackVar(node);
314 // node + transit + slack_var == next
315 // Add arcs for node + transit + slack_min <= next
316 AddArcs(node, next, CapAdd(transit, slack_var.Min()));
317 if (slack_var.Max() < std::numeric_limits<int64_t>::max()) {
318 // Add arcs for node + transit + slack_max >= next.
319 AddArcs(next, node, CapSub(-slack_var.Max(), transit));
320 }
321
322 node = next;
323 }
324
325 // Add vehicle span upper bound: end - span_ub <= start.
326 const int64_t span_ub = dimension_.GetSpanUpperBoundForVehicle(vehicle);
327 if (span_ub < std::numeric_limits<int64_t>::max()) {
328 AddArcs(model->End(vehicle), model->Start(vehicle), -span_ub);
329 }
330
331 // Set pickup/delivery limits on route.
332 std::vector<int> visited_pairs;
333 StoreVisitedPickupDeliveryPairsOnRoute(
334 dimension_, vehicle, next_accessor, &visited_pairs,
335 &visited_pickup_delivery_indices_for_pair_);
336 for (int pair_index : visited_pairs) {
337 const int64_t pickup_index =
338 visited_pickup_delivery_indices_for_pair_[pair_index].first;
339 const int64_t delivery_index =
340 visited_pickup_delivery_indices_for_pair_[pair_index].second;
341 visited_pickup_delivery_indices_for_pair_[pair_index] = {-1, -1};
342
343 DCHECK_GE(pickup_index, 0);
344 if (delivery_index < 0) {
345 // We didn't encounter a delivery for this pickup.
346 continue;
347 }
348
349 const int64_t limit = dimension_.GetPickupToDeliveryLimitForPair(
350 pair_index, model->GetPickupIndexPairs(pickup_index)[0].second,
351 model->GetDeliveryIndexPairs(delivery_index)[0].second);
352 if (limit < std::numeric_limits<int64_t>::max()) {
353 // delivery_cumul - limit <= pickup_cumul.
354 AddArcs(delivery_index, pickup_index, -limit);
355 }
356 }
357 }
358
359 for (const RoutingDimension::NodePrecedence& precedence :
360 dimension_.GetNodePrecedences()) {
361 const int first_index = precedence.first_node;
362 const int second_index = precedence.second_node;
363 if (lower_bounds[PositiveNode(first_index)] ==
364 std::numeric_limits<int64_t>::min() ||
365 lower_bounds[PositiveNode(second_index)] ==
366 std::numeric_limits<int64_t>::min()) {
367 // One of the nodes is unperformed, so the precedence rule doesn't apply.
368 continue;
369 }
370 AddArcs(first_index, second_index, precedence.offset);
371 }
372
373 return true;
374 }
375
UpdateCurrentLowerBoundOfNode(int node,int64_t new_lb,int64_t offset)376 bool CumulBoundsPropagator::UpdateCurrentLowerBoundOfNode(int node,
377 int64_t new_lb,
378 int64_t offset) {
379 const int cumul_var_index = node / 2;
380
381 if (node == PositiveNode(cumul_var_index)) {
382 // new_lb is a lower bound of the cumul of variable 'cumul_var_index'.
383 propagated_bounds_[node] = GetFirstPossibleValueForCumulWithOffset(
384 dimension_, cumul_var_index, new_lb, offset);
385 } else {
386 // -new_lb is an upper bound of the cumul of variable 'cumul_var_index'.
387 const int64_t new_ub = CapSub(0, new_lb);
388 propagated_bounds_[node] =
389 CapSub(0, GetLastPossibleValueForCumulWithOffset(
390 dimension_, cumul_var_index, new_ub, offset));
391 }
392
393 // Test that the lower/upper bounds do not cross each other.
394 const int64_t cumul_lower_bound =
395 propagated_bounds_[PositiveNode(cumul_var_index)];
396
397 const int64_t negated_cumul_upper_bound =
398 propagated_bounds_[NegativeNode(cumul_var_index)];
399
400 return CapAdd(negated_cumul_upper_bound, cumul_lower_bound) <= 0;
401 }
402
DisassembleSubtree(int source,int target)403 bool CumulBoundsPropagator::DisassembleSubtree(int source, int target) {
404 tmp_dfs_stack_.clear();
405 tmp_dfs_stack_.push_back(source);
406 while (!tmp_dfs_stack_.empty()) {
407 const int tail = tmp_dfs_stack_.back();
408 tmp_dfs_stack_.pop_back();
409 for (const ArcInfo& arc : outgoing_arcs_[tail]) {
410 const int child_node = arc.head;
411 if (tree_parent_node_of_[child_node] != tail) continue;
412 if (child_node == target) return false;
413 tree_parent_node_of_[child_node] = kParentToBePropagated;
414 tmp_dfs_stack_.push_back(child_node);
415 }
416 }
417 return true;
418 }
419
PropagateCumulBounds(const std::function<int64_t (int64_t)> & next_accessor,int64_t cumul_offset)420 bool CumulBoundsPropagator::PropagateCumulBounds(
421 const std::function<int64_t(int64_t)>& next_accessor,
422 int64_t cumul_offset) {
423 tree_parent_node_of_.assign(num_nodes_, kNoParent);
424 DCHECK(std::none_of(node_in_queue_.begin(), node_in_queue_.end(),
425 [](bool b) { return b; }));
426 DCHECK(bf_queue_.empty());
427
428 if (!InitializeArcsAndBounds(next_accessor, cumul_offset)) {
429 return CleanupAndReturnFalse();
430 }
431
432 std::vector<int64_t>& current_lb = propagated_bounds_;
433
434 // Bellman-Ford-Tarjan algorithm.
435 while (!bf_queue_.empty()) {
436 const int node = bf_queue_.front();
437 bf_queue_.pop_front();
438 node_in_queue_[node] = false;
439
440 if (tree_parent_node_of_[node] == kParentToBePropagated) {
441 // The parent of this node is still in the queue, so no need to process
442 // node now, since it will be re-enqued when its parent is processed.
443 continue;
444 }
445
446 const int64_t lower_bound = current_lb[node];
447 for (const ArcInfo& arc : outgoing_arcs_[node]) {
448 // NOTE: kint64min as a lower bound means no lower bound at all, so we
449 // don't use this value to propagate.
450 const int64_t induced_lb =
451 (lower_bound == std::numeric_limits<int64_t>::min())
452 ? std::numeric_limits<int64_t>::min()
453 : CapAdd(lower_bound, arc.offset);
454
455 const int head_node = arc.head;
456 if (induced_lb <= current_lb[head_node]) {
457 // No update necessary for the head_node, continue to next children of
458 // node.
459 continue;
460 }
461 if (!UpdateCurrentLowerBoundOfNode(head_node, induced_lb, cumul_offset) ||
462 !DisassembleSubtree(head_node, node)) {
463 // The new lower bound is infeasible, or a positive cycle was detected
464 // in the precedence graph by DisassembleSubtree().
465 return CleanupAndReturnFalse();
466 }
467
468 tree_parent_node_of_[head_node] = node;
469 AddNodeToQueue(head_node);
470 }
471 }
472 return true;
473 }
474
DimensionCumulOptimizerCore(const RoutingDimension * dimension,bool use_precedence_propagator)475 DimensionCumulOptimizerCore::DimensionCumulOptimizerCore(
476 const RoutingDimension* dimension, bool use_precedence_propagator)
477 : dimension_(dimension),
478 visited_pickup_delivery_indices_for_pair_(
479 dimension->model()->GetPickupAndDeliveryPairs().size(), {-1, -1}) {
480 if (use_precedence_propagator) {
481 propagator_ = absl::make_unique<CumulBoundsPropagator>(dimension);
482 }
483 const RoutingModel& model = *dimension_->model();
484 if (dimension_->HasBreakConstraints()) {
485 // Initialize vehicle_to_first_index_ so the variables of the breaks of
486 // vehicle v are stored from vehicle_to_first_index_[v] to
487 // vehicle_to_first_index_[v+1] - 1.
488 const int num_vehicles = model.vehicles();
489 vehicle_to_all_break_variables_offset_.reserve(num_vehicles);
490 int num_break_vars = 0;
491 for (int vehicle = 0; vehicle < num_vehicles; ++vehicle) {
492 vehicle_to_all_break_variables_offset_.push_back(num_break_vars);
493 const auto& intervals = dimension_->GetBreakIntervalsOfVehicle(vehicle);
494 num_break_vars += 2 * intervals.size(); // 2 variables per break.
495 }
496 all_break_variables_.resize(num_break_vars, -1);
497 }
498 if (!model.GetDimensionResourceGroupIndices(dimension_).empty()) {
499 resource_group_to_resource_to_vehicle_assignment_variables_.resize(
500 model.GetResourceGroups().size());
501 }
502 }
503
OptimizeSingleRoute(int vehicle,const std::function<int64_t (int64_t)> & next_accessor,RoutingLinearSolverWrapper * solver,std::vector<int64_t> * cumul_values,std::vector<int64_t> * break_values,int64_t * cost,int64_t * transit_cost,bool clear_lp)504 DimensionSchedulingStatus DimensionCumulOptimizerCore::OptimizeSingleRoute(
505 int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
506 RoutingLinearSolverWrapper* solver, std::vector<int64_t>* cumul_values,
507 std::vector<int64_t>* break_values, int64_t* cost, int64_t* transit_cost,
508 bool clear_lp) {
509 InitOptimizer(solver);
510 // Make sure SetRouteCumulConstraints will properly set the cumul bounds by
511 // looking at this route only.
512 DCHECK_EQ(propagator_.get(), nullptr);
513
514 RoutingModel* const model = dimension()->model();
515 const bool optimize_vehicle_costs =
516 (cumul_values != nullptr || cost != nullptr) &&
517 (!model->IsEnd(next_accessor(model->Start(vehicle))) ||
518 model->IsVehicleUsedWhenEmpty(vehicle));
519 const int64_t cumul_offset =
520 dimension_->GetLocalOptimizerOffsetForVehicle(vehicle);
521 int64_t cost_offset = 0;
522 if (!SetRouteCumulConstraints(vehicle, next_accessor, cumul_offset,
523 optimize_vehicle_costs, solver, transit_cost,
524 &cost_offset)) {
525 return DimensionSchedulingStatus::INFEASIBLE;
526 }
527 if (model->CheckLimit()) {
528 return DimensionSchedulingStatus::INFEASIBLE;
529 }
530 const DimensionSchedulingStatus status =
531 solver->Solve(model->RemainingTime());
532 if (status == DimensionSchedulingStatus::INFEASIBLE) {
533 solver->Clear();
534 return status;
535 }
536
537 SetValuesFromLP(current_route_cumul_variables_, cumul_offset, solver,
538 cumul_values);
539 SetValuesFromLP(current_route_break_variables_, cumul_offset, solver,
540 break_values);
541 if (cost != nullptr) {
542 *cost = CapAdd(cost_offset, solver->GetObjectiveValue());
543 }
544
545 if (clear_lp) {
546 solver->Clear();
547 }
548 return status;
549 }
550
551 namespace {
552
553 using ResourceGroup = RoutingModel::ResourceGroup;
554
GetDomainOffsetBounds(const Domain & domain,int64_t offset,ClosedInterval * interval)555 bool GetDomainOffsetBounds(const Domain& domain, int64_t offset,
556 ClosedInterval* interval) {
557 const int64_t lower_bound =
558 std::max<int64_t>(CapSub(domain.Min(), offset), 0);
559 const int64_t upper_bound =
560 domain.Max() == std::numeric_limits<int64_t>::max()
561 ? std::numeric_limits<int64_t>::max()
562 : CapSub(domain.Max(), offset);
563 if (lower_bound > upper_bound) return false;
564
565 *interval = ClosedInterval(lower_bound, upper_bound);
566 return true;
567 }
568
GetIntervalIntersectionWithOffsetDomain(const ClosedInterval & interval,const Domain & domain,int64_t offset,ClosedInterval * intersection)569 bool GetIntervalIntersectionWithOffsetDomain(const ClosedInterval& interval,
570 const Domain& domain,
571 int64_t offset,
572 ClosedInterval* intersection) {
573 ClosedInterval domain_bounds;
574 if (!GetDomainOffsetBounds(domain, offset, &domain_bounds)) {
575 return false;
576 }
577 const int64_t intersection_lb = std::max(interval.start, domain_bounds.start);
578 const int64_t intersection_ub = std::min(interval.end, domain_bounds.end);
579 if (intersection_lb > intersection_ub) return false;
580
581 *intersection = ClosedInterval(intersection_lb, intersection_ub);
582 return true;
583 }
584
GetVariableBounds(int index,const RoutingLinearSolverWrapper & solver)585 ClosedInterval GetVariableBounds(int index,
586 const RoutingLinearSolverWrapper& solver) {
587 return ClosedInterval(solver.GetVariableLowerBound(index),
588 solver.GetVariableUpperBound(index));
589 }
590
TightenStartEndVariableBoundsWithResource(const RoutingDimension & dimension,const ResourceGroup::Resource & resource,const ClosedInterval & start_bounds,int start_index,const ClosedInterval & end_bounds,int end_index,int64_t offset,RoutingLinearSolverWrapper * solver)591 bool TightenStartEndVariableBoundsWithResource(
592 const RoutingDimension& dimension, const ResourceGroup::Resource& resource,
593 const ClosedInterval& start_bounds, int start_index,
594 const ClosedInterval& end_bounds, int end_index, int64_t offset,
595 RoutingLinearSolverWrapper* solver) {
596 const ResourceGroup::Attributes& attributes =
597 resource.GetDimensionAttributes(&dimension);
598 ClosedInterval new_start_bounds;
599 ClosedInterval new_end_bounds;
600 return GetIntervalIntersectionWithOffsetDomain(start_bounds,
601 attributes.start_domain(),
602 offset, &new_start_bounds) &&
603 solver->SetVariableBounds(start_index, new_start_bounds.start,
604 new_start_bounds.end) &&
605 GetIntervalIntersectionWithOffsetDomain(
606 end_bounds, attributes.end_domain(), offset, &new_end_bounds) &&
607 solver->SetVariableBounds(end_index, new_end_bounds.start,
608 new_end_bounds.end);
609 }
610
611 } // namespace
612
613 std::vector<DimensionSchedulingStatus>
OptimizeSingleRouteWithResources(int vehicle,const std::function<int64_t (int64_t)> & next_accessor,const std::vector<RoutingModel::ResourceGroup::Resource> & resources,const std::vector<int> & resource_indices,bool optimize_vehicle_costs,RoutingLinearSolverWrapper * solver,std::vector<int64_t> * costs_without_transits,std::vector<std::vector<int64_t>> * cumul_values,std::vector<std::vector<int64_t>> * break_values,bool clear_lp)614 DimensionCumulOptimizerCore::OptimizeSingleRouteWithResources(
615 int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
616 const std::vector<RoutingModel::ResourceGroup::Resource>& resources,
617 const std::vector<int>& resource_indices, bool optimize_vehicle_costs,
618 RoutingLinearSolverWrapper* solver,
619 std::vector<int64_t>* costs_without_transits,
620 std::vector<std::vector<int64_t>>* cumul_values,
621 std::vector<std::vector<int64_t>>* break_values, bool clear_lp) {
622 if (resource_indices.empty()) return {};
623
624 InitOptimizer(solver);
625 // Make sure SetRouteCumulConstraints will properly set the cumul bounds by
626 // looking at this route only.
627 DCHECK_EQ(propagator_.get(), nullptr);
628 DCHECK_NE(costs_without_transits, nullptr);
629 costs_without_transits->clear();
630
631 RoutingModel* const model = dimension()->model();
632 if (model->IsEnd(next_accessor(model->Start(vehicle))) &&
633 !model->IsVehicleUsedWhenEmpty(vehicle)) {
634 // An unused empty vehicle doesn't require resources.
635 return {};
636 }
637
638 const int64_t cumul_offset =
639 dimension_->GetLocalOptimizerOffsetForVehicle(vehicle);
640 int64_t cost_offset = 0;
641 int64_t transit_cost = 0;
642 if (!SetRouteCumulConstraints(vehicle, next_accessor, cumul_offset,
643 optimize_vehicle_costs, solver, &transit_cost,
644 &cost_offset)) {
645 return {DimensionSchedulingStatus::INFEASIBLE};
646 }
647
648 costs_without_transits->assign(resource_indices.size(), -1);
649 if (cumul_values != nullptr) {
650 cumul_values->assign(resource_indices.size(), {});
651 }
652 if (break_values != nullptr) {
653 break_values->assign(resource_indices.size(), {});
654 }
655
656 DCHECK_GE(current_route_cumul_variables_.size(), 2);
657
658 const int start_cumul = current_route_cumul_variables_[0];
659 const ClosedInterval start_bounds = GetVariableBounds(start_cumul, *solver);
660 const int end_cumul = current_route_cumul_variables_.back();
661 const ClosedInterval end_bounds = GetVariableBounds(end_cumul, *solver);
662 std::vector<DimensionSchedulingStatus> statuses;
663 for (int i = 0; i < resource_indices.size(); i++) {
664 if (model->CheckLimit()) {
665 // The model's deadline has been reached, stop.
666 costs_without_transits->clear();
667 if (cumul_values != nullptr) {
668 cumul_values->clear();
669 }
670 if (break_values != nullptr) {
671 break_values->clear();
672 }
673 return {};
674 }
675 if (!TightenStartEndVariableBoundsWithResource(
676 *dimension_, resources[resource_indices[i]], start_bounds,
677 start_cumul, end_bounds, end_cumul, cumul_offset, solver)) {
678 // The resource attributes don't match this vehicle.
679 statuses.push_back(DimensionSchedulingStatus::INFEASIBLE);
680 continue;
681 }
682
683 statuses.push_back(solver->Solve(model->RemainingTime()));
684 if (statuses.back() == DimensionSchedulingStatus::INFEASIBLE) {
685 continue;
686 }
687 costs_without_transits->at(i) =
688 optimize_vehicle_costs
689 ? CapSub(CapAdd(cost_offset, solver->GetObjectiveValue()),
690 transit_cost)
691 : 0;
692
693 if (cumul_values != nullptr) {
694 SetValuesFromLP(current_route_cumul_variables_, cumul_offset, solver,
695 &cumul_values->at(i));
696 }
697 if (break_values != nullptr) {
698 SetValuesFromLP(current_route_break_variables_, cumul_offset, solver,
699 &break_values->at(i));
700 }
701 }
702
703 if (clear_lp) {
704 solver->Clear();
705 }
706 return statuses;
707 }
708
Optimize(const std::function<int64_t (int64_t)> & next_accessor,RoutingLinearSolverWrapper * solver,std::vector<int64_t> * cumul_values,std::vector<int64_t> * break_values,std::vector<std::vector<int>> * resource_indices_per_group,int64_t * cost,int64_t * transit_cost,bool clear_lp)709 DimensionSchedulingStatus DimensionCumulOptimizerCore::Optimize(
710 const std::function<int64_t(int64_t)>& next_accessor,
711 RoutingLinearSolverWrapper* solver, std::vector<int64_t>* cumul_values,
712 std::vector<int64_t>* break_values,
713 std::vector<std::vector<int>>* resource_indices_per_group, int64_t* cost,
714 int64_t* transit_cost, bool clear_lp) {
715 InitOptimizer(solver);
716
717 // If both "cumul_values" and "cost" parameters are null, we don't try to
718 // optimize the cost and stop at the first feasible solution.
719 const bool optimize_costs = (cumul_values != nullptr) || (cost != nullptr);
720 bool has_vehicles_being_optimized = false;
721
722 const int64_t cumul_offset = dimension_->GetGlobalOptimizerOffset();
723
724 if (propagator_ != nullptr &&
725 !propagator_->PropagateCumulBounds(next_accessor, cumul_offset)) {
726 return DimensionSchedulingStatus::INFEASIBLE;
727 }
728
729 int64_t total_transit_cost = 0;
730 int64_t total_cost_offset = 0;
731 const RoutingModel* model = dimension()->model();
732 for (int vehicle = 0; vehicle < model->vehicles(); vehicle++) {
733 int64_t route_transit_cost = 0;
734 int64_t route_cost_offset = 0;
735 const bool vehicle_is_used =
736 !model->IsEnd(next_accessor(model->Start(vehicle))) ||
737 model->IsVehicleUsedWhenEmpty(vehicle);
738 const bool optimize_vehicle_costs = optimize_costs && vehicle_is_used;
739 if (!SetRouteCumulConstraints(vehicle, next_accessor, cumul_offset,
740 optimize_vehicle_costs, solver,
741 &route_transit_cost, &route_cost_offset)) {
742 return DimensionSchedulingStatus::INFEASIBLE;
743 }
744 total_transit_cost = CapAdd(total_transit_cost, route_transit_cost);
745 total_cost_offset = CapAdd(total_cost_offset, route_cost_offset);
746 has_vehicles_being_optimized |= optimize_vehicle_costs;
747 }
748 if (transit_cost != nullptr) {
749 *transit_cost = total_transit_cost;
750 }
751
752 if (!SetGlobalConstraints(next_accessor, cumul_offset,
753 has_vehicles_being_optimized, solver)) {
754 return DimensionSchedulingStatus::INFEASIBLE;
755 }
756
757 const DimensionSchedulingStatus status =
758 solver->Solve(model->RemainingTime());
759 if (status == DimensionSchedulingStatus::INFEASIBLE) {
760 solver->Clear();
761 return status;
762 }
763
764 // TODO(user): In case the status is RELAXED_OPTIMAL_ONLY, check we can
765 // safely avoid filling variable and cost values.
766 SetValuesFromLP(index_to_cumul_variable_, cumul_offset, solver, cumul_values);
767 SetValuesFromLP(all_break_variables_, cumul_offset, solver, break_values);
768 SetResourceIndices(solver, resource_indices_per_group);
769
770 if (cost != nullptr) {
771 *cost = CapAdd(solver->GetObjectiveValue(), total_cost_offset);
772 }
773
774 if (clear_lp) {
775 solver->Clear();
776 }
777 return status;
778 }
779
OptimizeAndPack(const std::function<int64_t (int64_t)> & next_accessor,RoutingLinearSolverWrapper * solver,std::vector<int64_t> * cumul_values,std::vector<int64_t> * break_values,std::vector<std::vector<int>> * resource_indices_per_group)780 DimensionSchedulingStatus DimensionCumulOptimizerCore::OptimizeAndPack(
781 const std::function<int64_t(int64_t)>& next_accessor,
782 RoutingLinearSolverWrapper* solver, std::vector<int64_t>* cumul_values,
783 std::vector<int64_t>* break_values,
784 std::vector<std::vector<int>>* resource_indices_per_group) {
785 // Note: We pass a non-nullptr cost to the Optimize() method so the costs
786 // are optimized by the solver.
787 int64_t cost = 0;
788 if (Optimize(next_accessor, solver,
789 /*cumul_values=*/nullptr, /*break_values=*/nullptr,
790 /*resource_indices_per_group=*/nullptr, &cost,
791 /*transit_cost=*/nullptr,
792 /*clear_lp=*/false) == DimensionSchedulingStatus::INFEASIBLE) {
793 return DimensionSchedulingStatus::INFEASIBLE;
794 }
795 std::vector<int> vehicles(dimension()->model()->vehicles());
796 std::iota(vehicles.begin(), vehicles.end(), 0);
797 // Subtle: Even if the status was RELAXED_OPTIMAL_ONLY we try to pack just in
798 // case packing manages to make the solution completely feasible.
799 DimensionSchedulingStatus status = PackRoutes(vehicles, solver);
800 if (status == DimensionSchedulingStatus::INFEASIBLE) {
801 return status;
802 }
803 // TODO(user): In case the status is RELAXED_OPTIMAL_ONLY, check we can
804 // safely avoid filling variable values.
805 const int64_t global_offset = dimension_->GetGlobalOptimizerOffset();
806 SetValuesFromLP(index_to_cumul_variable_, global_offset, solver,
807 cumul_values);
808 SetValuesFromLP(all_break_variables_, global_offset, solver, break_values);
809 SetResourceIndices(solver, resource_indices_per_group);
810 solver->Clear();
811 return status;
812 }
813
814 DimensionSchedulingStatus
OptimizeAndPackSingleRoute(int vehicle,const std::function<int64_t (int64_t)> & next_accessor,const RoutingModel::ResourceGroup::Resource * resource,RoutingLinearSolverWrapper * solver,std::vector<int64_t> * cumul_values,std::vector<int64_t> * break_values)815 DimensionCumulOptimizerCore::OptimizeAndPackSingleRoute(
816 int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
817 const RoutingModel::ResourceGroup::Resource* resource,
818 RoutingLinearSolverWrapper* solver, std::vector<int64_t>* cumul_values,
819 std::vector<int64_t>* break_values) {
820 if (resource == nullptr) {
821 // Note: We pass a non-nullptr cost to the OptimizeSingleRoute() method so
822 // the costs are optimized by the LP.
823 int64_t cost = 0;
824 if (OptimizeSingleRoute(vehicle, next_accessor, solver,
825 /*cumul_values=*/nullptr, /*break_values=*/nullptr,
826 &cost, /*transit_cost=*/nullptr,
827 /*clear_lp=*/false) ==
828 DimensionSchedulingStatus::INFEASIBLE) {
829 return DimensionSchedulingStatus::INFEASIBLE;
830 }
831 } else {
832 std::vector<int64_t> costs_without_transits;
833 const std::vector<DimensionSchedulingStatus> statuses =
834 OptimizeSingleRouteWithResources(
835 vehicle, next_accessor, {*resource}, {0},
836 /*optimize_vehicle_costs=*/true, solver, &costs_without_transits,
837 /*cumul_values=*/nullptr,
838 /*break_values=*/nullptr,
839 /*clear_lp=*/false);
840 DCHECK_EQ(statuses.size(), 1);
841 if (statuses[0] == DimensionSchedulingStatus::INFEASIBLE) {
842 return DimensionSchedulingStatus::INFEASIBLE;
843 }
844 }
845
846 const DimensionSchedulingStatus status = PackRoutes({vehicle}, solver);
847 if (status == DimensionSchedulingStatus::INFEASIBLE) {
848 return DimensionSchedulingStatus::INFEASIBLE;
849 }
850 const int64_t local_offset =
851 dimension_->GetLocalOptimizerOffsetForVehicle(vehicle);
852 SetValuesFromLP(current_route_cumul_variables_, local_offset, solver,
853 cumul_values);
854 SetValuesFromLP(current_route_break_variables_, local_offset, solver,
855 break_values);
856 solver->Clear();
857 return status;
858 }
859
PackRoutes(std::vector<int> vehicles,RoutingLinearSolverWrapper * solver)860 DimensionSchedulingStatus DimensionCumulOptimizerCore::PackRoutes(
861 std::vector<int> vehicles, RoutingLinearSolverWrapper* solver) {
862 const RoutingModel* model = dimension_->model();
863
864 // NOTE(user): Given our constraint matrix, our problem *should* always
865 // have an integer optimal solution, in which case we can round to the nearest
866 // integer both for the objective constraint bound (returned by
867 // GetObjectiveValue()) and the end cumul variable bound after minimizing
868 // (see b/154381899 showcasing an example where std::ceil leads to an
869 // "imperfect" packing due to rounding precision errors).
870 // If this DCHECK ever fails, it can be removed but the code below should be
871 // adapted to have a 2-phase approach, solving once with the rounded value as
872 // bound and if this fails, solve again using std::ceil.
873 DCHECK(solver->SolutionIsInteger());
874
875 // Minimize the route end times without increasing the cost.
876 solver->AddObjectiveConstraint();
877 solver->ClearObjective();
878 for (int vehicle : vehicles) {
879 solver->SetObjectiveCoefficient(
880 index_to_cumul_variable_[model->End(vehicle)], 1);
881 }
882
883 if (solver->Solve(model->RemainingTime()) ==
884 DimensionSchedulingStatus::INFEASIBLE) {
885 return DimensionSchedulingStatus::INFEASIBLE;
886 }
887
888 // Maximize the route start times without increasing the cost or the route end
889 // times.
890 solver->ClearObjective();
891 for (int vehicle : vehicles) {
892 const int end_cumul_var = index_to_cumul_variable_[model->End(vehicle)];
893 // end_cumul_var <= solver.GetValue(end_cumul_var)
894 solver->SetVariableBounds(
895 end_cumul_var, solver->GetVariableLowerBound(end_cumul_var),
896 MathUtil::FastInt64Round(solver->GetValue(end_cumul_var)));
897
898 // Maximize the starts of the routes.
899 solver->SetObjectiveCoefficient(
900 index_to_cumul_variable_[model->Start(vehicle)], -1);
901 }
902 return solver->Solve(model->RemainingTime());
903 }
904
InitOptimizer(RoutingLinearSolverWrapper * solver)905 void DimensionCumulOptimizerCore::InitOptimizer(
906 RoutingLinearSolverWrapper* solver) {
907 solver->Clear();
908 index_to_cumul_variable_.assign(dimension_->cumuls().size(), -1);
909 max_end_cumul_ = solver->CreateNewPositiveVariable();
910 min_start_cumul_ = solver->CreateNewPositiveVariable();
911 }
912
ComputeRouteCumulBounds(const std::vector<int64_t> & route,const std::vector<int64_t> & fixed_transits,int64_t cumul_offset)913 bool DimensionCumulOptimizerCore::ComputeRouteCumulBounds(
914 const std::vector<int64_t>& route,
915 const std::vector<int64_t>& fixed_transits, int64_t cumul_offset) {
916 const int route_size = route.size();
917 current_route_min_cumuls_.resize(route_size);
918 current_route_max_cumuls_.resize(route_size);
919 if (propagator_ != nullptr) {
920 for (int pos = 0; pos < route_size; pos++) {
921 const int64_t node = route[pos];
922 current_route_min_cumuls_[pos] = propagator_->CumulMin(node);
923 DCHECK_GE(current_route_min_cumuls_[pos], 0);
924 current_route_max_cumuls_[pos] = propagator_->CumulMax(node);
925 DCHECK_GE(current_route_max_cumuls_[pos], current_route_min_cumuls_[pos]);
926 }
927 return true;
928 }
929
930 // Extract cumul min/max and fixed transits from CP.
931 for (int pos = 0; pos < route_size; ++pos) {
932 if (!GetCumulBoundsWithOffset(*dimension_, route[pos], cumul_offset,
933 ¤t_route_min_cumuls_[pos],
934 ¤t_route_max_cumuls_[pos])) {
935 return false;
936 }
937 }
938
939 // Refine cumul bounds using
940 // cumul[i+1] >= cumul[i] + fixed_transit[i] + slack[i].
941 for (int pos = 1; pos < route_size; ++pos) {
942 const int64_t slack_min = dimension_->SlackVar(route[pos - 1])->Min();
943 current_route_min_cumuls_[pos] = std::max(
944 current_route_min_cumuls_[pos],
945 CapAdd(
946 CapAdd(current_route_min_cumuls_[pos - 1], fixed_transits[pos - 1]),
947 slack_min));
948 current_route_min_cumuls_[pos] = GetFirstPossibleValueForCumulWithOffset(
949 *dimension_, route[pos], current_route_min_cumuls_[pos], cumul_offset);
950 if (current_route_min_cumuls_[pos] > current_route_max_cumuls_[pos]) {
951 return false;
952 }
953 }
954
955 for (int pos = route_size - 2; pos >= 0; --pos) {
956 // If cumul_max[pos+1] is kint64max, it will be translated to
957 // double +infinity, so it must not constrain cumul_max[pos].
958 if (current_route_max_cumuls_[pos + 1] <
959 std::numeric_limits<int64_t>::max()) {
960 const int64_t slack_min = dimension_->SlackVar(route[pos])->Min();
961 current_route_max_cumuls_[pos] = std::min(
962 current_route_max_cumuls_[pos],
963 CapSub(
964 CapSub(current_route_max_cumuls_[pos + 1], fixed_transits[pos]),
965 slack_min));
966 current_route_max_cumuls_[pos] = GetLastPossibleValueForCumulWithOffset(
967 *dimension_, route[pos], current_route_max_cumuls_[pos],
968 cumul_offset);
969 if (current_route_max_cumuls_[pos] < current_route_min_cumuls_[pos]) {
970 return false;
971 }
972 }
973 }
974 return true;
975 }
976
SetRouteCumulConstraints(int vehicle,const std::function<int64_t (int64_t)> & next_accessor,int64_t cumul_offset,bool optimize_costs,RoutingLinearSolverWrapper * solver,int64_t * route_transit_cost,int64_t * route_cost_offset)977 bool DimensionCumulOptimizerCore::SetRouteCumulConstraints(
978 int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
979 int64_t cumul_offset, bool optimize_costs,
980 RoutingLinearSolverWrapper* solver, int64_t* route_transit_cost,
981 int64_t* route_cost_offset) {
982 RoutingModel* const model = dimension_->model();
983 // Extract the vehicle's path from next_accessor.
984 std::vector<int64_t> path;
985 {
986 int node = model->Start(vehicle);
987 path.push_back(node);
988 while (!model->IsEnd(node)) {
989 node = next_accessor(node);
990 path.push_back(node);
991 }
992 DCHECK_GE(path.size(), 2);
993 }
994 const int path_size = path.size();
995
996 std::vector<int64_t> fixed_transit(path_size - 1);
997 {
998 const std::function<int64_t(int64_t, int64_t)>& transit_accessor =
999 dimension_->transit_evaluator(vehicle);
1000 for (int pos = 1; pos < path_size; ++pos) {
1001 fixed_transit[pos - 1] = transit_accessor(path[pos - 1], path[pos]);
1002 }
1003 }
1004
1005 if (!ComputeRouteCumulBounds(path, fixed_transit, cumul_offset)) {
1006 return false;
1007 }
1008
1009 // LP Model variables, current_route_cumul_variables_ and lp_slacks.
1010 // Create LP variables for cumuls.
1011 std::vector<int>& lp_cumuls = current_route_cumul_variables_;
1012 lp_cumuls.assign(path_size, -1);
1013 for (int pos = 0; pos < path_size; ++pos) {
1014 const int lp_cumul = solver->CreateNewPositiveVariable();
1015 index_to_cumul_variable_[path[pos]] = lp_cumul;
1016 lp_cumuls[pos] = lp_cumul;
1017 if (!solver->SetVariableBounds(lp_cumul, current_route_min_cumuls_[pos],
1018 current_route_max_cumuls_[pos])) {
1019 return false;
1020 }
1021 const SortedDisjointIntervalList& forbidden =
1022 dimension_->forbidden_intervals()[path[pos]];
1023 if (forbidden.NumIntervals() > 0) {
1024 std::vector<int64_t> starts;
1025 std::vector<int64_t> ends;
1026 for (const ClosedInterval interval :
1027 dimension_->GetAllowedIntervalsInRange(
1028 path[pos], CapAdd(current_route_min_cumuls_[pos], cumul_offset),
1029 CapAdd(current_route_max_cumuls_[pos], cumul_offset))) {
1030 starts.push_back(CapSub(interval.start, cumul_offset));
1031 ends.push_back(CapSub(interval.end, cumul_offset));
1032 }
1033 solver->SetVariableDisjointBounds(lp_cumul, starts, ends);
1034 }
1035 }
1036 // Create LP variables for slacks.
1037 std::vector<int> lp_slacks(path_size - 1, -1);
1038 for (int pos = 0; pos < path_size - 1; ++pos) {
1039 const IntVar* cp_slack = dimension_->SlackVar(path[pos]);
1040 lp_slacks[pos] = solver->CreateNewPositiveVariable();
1041 if (!solver->SetVariableBounds(lp_slacks[pos], cp_slack->Min(),
1042 cp_slack->Max())) {
1043 return false;
1044 }
1045 }
1046
1047 // LP Model constraints and costs.
1048 // Add all path constraints to LP:
1049 // cumul[i] + fixed_transit[i] + slack[i] == cumul[i+1]
1050 // <=> fixed_transit[i] == cumul[i+1] - cumul[i] - slack[i].
1051 for (int pos = 0; pos < path_size - 1; ++pos) {
1052 const int ct =
1053 solver->CreateNewConstraint(fixed_transit[pos], fixed_transit[pos]);
1054 solver->SetCoefficient(ct, lp_cumuls[pos + 1], 1);
1055 solver->SetCoefficient(ct, lp_cumuls[pos], -1);
1056 solver->SetCoefficient(ct, lp_slacks[pos], -1);
1057 }
1058 if (route_cost_offset != nullptr) *route_cost_offset = 0;
1059 if (optimize_costs) {
1060 // Add soft upper bounds.
1061 for (int pos = 0; pos < path_size; ++pos) {
1062 if (!dimension_->HasCumulVarSoftUpperBound(path[pos])) continue;
1063 const int64_t coef =
1064 dimension_->GetCumulVarSoftUpperBoundCoefficient(path[pos]);
1065 if (coef == 0) continue;
1066 int64_t bound = dimension_->GetCumulVarSoftUpperBound(path[pos]);
1067 if (bound < cumul_offset && route_cost_offset != nullptr) {
1068 // Add coef * (cumul_offset - bound) to the cost offset.
1069 *route_cost_offset = CapAdd(*route_cost_offset,
1070 CapProd(CapSub(cumul_offset, bound), coef));
1071 }
1072 bound = std::max<int64_t>(0, CapSub(bound, cumul_offset));
1073 if (current_route_max_cumuls_[pos] <= bound) {
1074 // constraint is never violated.
1075 continue;
1076 }
1077 const int soft_ub_diff = solver->CreateNewPositiveVariable();
1078 solver->SetObjectiveCoefficient(soft_ub_diff, coef);
1079 // cumul - soft_ub_diff <= bound.
1080 const int ct = solver->CreateNewConstraint(
1081 std::numeric_limits<int64_t>::min(), bound);
1082 solver->SetCoefficient(ct, lp_cumuls[pos], 1);
1083 solver->SetCoefficient(ct, soft_ub_diff, -1);
1084 }
1085 // Add soft lower bounds.
1086 for (int pos = 0; pos < path_size; ++pos) {
1087 if (!dimension_->HasCumulVarSoftLowerBound(path[pos])) continue;
1088 const int64_t coef =
1089 dimension_->GetCumulVarSoftLowerBoundCoefficient(path[pos]);
1090 if (coef == 0) continue;
1091 const int64_t bound = std::max<int64_t>(
1092 0, CapSub(dimension_->GetCumulVarSoftLowerBound(path[pos]),
1093 cumul_offset));
1094 if (current_route_min_cumuls_[pos] >= bound) {
1095 // constraint is never violated.
1096 continue;
1097 }
1098 const int soft_lb_diff = solver->CreateNewPositiveVariable();
1099 solver->SetObjectiveCoefficient(soft_lb_diff, coef);
1100 // bound - cumul <= soft_lb_diff
1101 const int ct = solver->CreateNewConstraint(
1102 bound, std::numeric_limits<int64_t>::max());
1103 solver->SetCoefficient(ct, lp_cumuls[pos], 1);
1104 solver->SetCoefficient(ct, soft_lb_diff, 1);
1105 }
1106 }
1107 // Add pickup and delivery limits.
1108 std::vector<int> visited_pairs;
1109 StoreVisitedPickupDeliveryPairsOnRoute(
1110 *dimension_, vehicle, next_accessor, &visited_pairs,
1111 &visited_pickup_delivery_indices_for_pair_);
1112 for (int pair_index : visited_pairs) {
1113 const int64_t pickup_index =
1114 visited_pickup_delivery_indices_for_pair_[pair_index].first;
1115 const int64_t delivery_index =
1116 visited_pickup_delivery_indices_for_pair_[pair_index].second;
1117 visited_pickup_delivery_indices_for_pair_[pair_index] = {-1, -1};
1118
1119 DCHECK_GE(pickup_index, 0);
1120 if (delivery_index < 0) {
1121 // We didn't encounter a delivery for this pickup.
1122 continue;
1123 }
1124
1125 const int64_t limit = dimension_->GetPickupToDeliveryLimitForPair(
1126 pair_index, model->GetPickupIndexPairs(pickup_index)[0].second,
1127 model->GetDeliveryIndexPairs(delivery_index)[0].second);
1128 if (limit < std::numeric_limits<int64_t>::max()) {
1129 // delivery_cumul - pickup_cumul <= limit.
1130 const int ct = solver->CreateNewConstraint(
1131 std::numeric_limits<int64_t>::min(), limit);
1132 solver->SetCoefficient(ct, index_to_cumul_variable_[delivery_index], 1);
1133 solver->SetCoefficient(ct, index_to_cumul_variable_[pickup_index], -1);
1134 }
1135 }
1136
1137 // Add span bound constraint.
1138 const int64_t span_bound = dimension_->GetSpanUpperBoundForVehicle(vehicle);
1139 if (span_bound < std::numeric_limits<int64_t>::max()) {
1140 // end_cumul - start_cumul <= bound
1141 const int ct = solver->CreateNewConstraint(
1142 std::numeric_limits<int64_t>::min(), span_bound);
1143 solver->SetCoefficient(ct, lp_cumuls.back(), 1);
1144 solver->SetCoefficient(ct, lp_cumuls.front(), -1);
1145 }
1146 // Add span cost.
1147 const int64_t span_cost_coef =
1148 dimension_->GetSpanCostCoefficientForVehicle(vehicle);
1149 if (optimize_costs && span_cost_coef > 0) {
1150 solver->SetObjectiveCoefficient(lp_cumuls.back(), span_cost_coef);
1151 solver->SetObjectiveCoefficient(lp_cumuls.front(), -span_cost_coef);
1152 }
1153 // Add soft span cost.
1154 if (optimize_costs && dimension_->HasSoftSpanUpperBounds()) {
1155 SimpleBoundCosts::BoundCost bound_cost =
1156 dimension_->GetSoftSpanUpperBoundForVehicle(vehicle);
1157 if (bound_cost.bound < std::numeric_limits<int64_t>::max() &&
1158 bound_cost.cost > 0) {
1159 const int span_violation = solver->CreateNewPositiveVariable();
1160 // end - start <= bound + span_violation
1161 const int violation = solver->CreateNewConstraint(
1162 std::numeric_limits<int64_t>::min(), bound_cost.bound);
1163 solver->SetCoefficient(violation, lp_cumuls.back(), 1.0);
1164 solver->SetCoefficient(violation, lp_cumuls.front(), -1.0);
1165 solver->SetCoefficient(violation, span_violation, -1.0);
1166 // Add span_violation * cost to objective.
1167 solver->SetObjectiveCoefficient(span_violation, bound_cost.cost);
1168 }
1169 }
1170 // Add global span constraint.
1171 if (optimize_costs && dimension_->global_span_cost_coefficient() > 0) {
1172 // min_start_cumul_ <= cumuls[start]
1173 int ct =
1174 solver->CreateNewConstraint(std::numeric_limits<int64_t>::min(), 0);
1175 solver->SetCoefficient(ct, min_start_cumul_, 1);
1176 solver->SetCoefficient(ct, lp_cumuls.front(), -1);
1177 // max_end_cumul_ >= cumuls[end]
1178 ct = solver->CreateNewConstraint(0, std::numeric_limits<int64_t>::max());
1179 solver->SetCoefficient(ct, max_end_cumul_, 1);
1180 solver->SetCoefficient(ct, lp_cumuls.back(), -1);
1181 }
1182 // Fill transit cost if specified.
1183 if (route_transit_cost != nullptr) {
1184 if (optimize_costs && span_cost_coef > 0) {
1185 const int64_t total_fixed_transit = std::accumulate(
1186 fixed_transit.begin(), fixed_transit.end(), 0, CapAdd);
1187 *route_transit_cost = CapProd(total_fixed_transit, span_cost_coef);
1188 } else {
1189 *route_transit_cost = 0;
1190 }
1191 }
1192 // For every break that must be inside the route, the duration of that break
1193 // must be flowed in the slacks of arcs that can intersect the break.
1194 // This LP modelization is correct but not complete:
1195 // can miss some cases where the breaks cannot fit.
1196 // TODO(user): remove the need for returns in the code below.
1197 current_route_break_variables_.clear();
1198 if (!dimension_->HasBreakConstraints()) return true;
1199 const std::vector<IntervalVar*>& breaks =
1200 dimension_->GetBreakIntervalsOfVehicle(vehicle);
1201 const int num_breaks = breaks.size();
1202 // When there are no breaks, only break distance needs to be modeled,
1203 // and it reduces to a span maximum.
1204 // TODO(user): Also add the case where no breaks can intersect the route.
1205 if (num_breaks == 0) {
1206 int64_t maximum_route_span = std::numeric_limits<int64_t>::max();
1207 for (const auto& distance_duration :
1208 dimension_->GetBreakDistanceDurationOfVehicle(vehicle)) {
1209 maximum_route_span =
1210 std::min(maximum_route_span, distance_duration.first);
1211 }
1212 if (maximum_route_span < std::numeric_limits<int64_t>::max()) {
1213 const int ct = solver->CreateNewConstraint(
1214 std::numeric_limits<int64_t>::min(), maximum_route_span);
1215 solver->SetCoefficient(ct, lp_cumuls.back(), 1);
1216 solver->SetCoefficient(ct, lp_cumuls.front(), -1);
1217 }
1218 return true;
1219 }
1220 // Gather visit information: the visit of node i has [start, end) =
1221 // [cumul[i] - post_travel[i-1], cumul[i] + pre_travel[i]).
1222 // Breaks cannot overlap those visit intervals.
1223 std::vector<int64_t> pre_travel(path_size - 1, 0);
1224 std::vector<int64_t> post_travel(path_size - 1, 0);
1225 {
1226 const int pre_travel_index =
1227 dimension_->GetPreTravelEvaluatorOfVehicle(vehicle);
1228 if (pre_travel_index != -1) {
1229 FillPathEvaluation(path, model->TransitCallback(pre_travel_index),
1230 &pre_travel);
1231 }
1232 const int post_travel_index =
1233 dimension_->GetPostTravelEvaluatorOfVehicle(vehicle);
1234 if (post_travel_index != -1) {
1235 FillPathEvaluation(path, model->TransitCallback(post_travel_index),
1236 &post_travel);
1237 }
1238 }
1239 // If the solver is CPSAT, it will need to represent the times at which
1240 // breaks are scheduled, those variables are used both in the pure breaks
1241 // part and in the break distance part of the model.
1242 // Otherwise, it doesn't need the variables and they are not created.
1243 std::vector<int> lp_break_start;
1244 std::vector<int> lp_break_duration;
1245 std::vector<int> lp_break_end;
1246 if (solver->IsCPSATSolver()) {
1247 lp_break_start.resize(num_breaks, -1);
1248 lp_break_duration.resize(num_breaks, -1);
1249 lp_break_end.resize(num_breaks, -1);
1250 }
1251
1252 std::vector<int> slack_exact_lower_bound_ct(path_size - 1, -1);
1253 std::vector<int> slack_linear_lower_bound_ct(path_size - 1, -1);
1254
1255 const int64_t vehicle_start_min = current_route_min_cumuls_.front();
1256 const int64_t vehicle_start_max = current_route_max_cumuls_.front();
1257 const int64_t vehicle_end_min = current_route_min_cumuls_.back();
1258 const int64_t vehicle_end_max = current_route_max_cumuls_.back();
1259 const int all_break_variables_offset =
1260 vehicle_to_all_break_variables_offset_[vehicle];
1261 for (int br = 0; br < num_breaks; ++br) {
1262 const IntervalVar& break_var = *breaks[br];
1263 if (!break_var.MustBePerformed()) continue;
1264 const int64_t break_start_min = CapSub(break_var.StartMin(), cumul_offset);
1265 const int64_t break_start_max = CapSub(break_var.StartMax(), cumul_offset);
1266 const int64_t break_end_min = CapSub(break_var.EndMin(), cumul_offset);
1267 const int64_t break_end_max = CapSub(break_var.EndMax(), cumul_offset);
1268 const int64_t break_duration_min = break_var.DurationMin();
1269 const int64_t break_duration_max = break_var.DurationMax();
1270 // The CPSAT solver encodes all breaks that can intersect the route,
1271 // the LP solver only encodes the breaks that must intersect the route.
1272 if (solver->IsCPSATSolver()) {
1273 if (break_end_max <= vehicle_start_min ||
1274 vehicle_end_max <= break_start_min) {
1275 all_break_variables_[all_break_variables_offset + 2 * br] = -1;
1276 all_break_variables_[all_break_variables_offset + 2 * br + 1] = -1;
1277 current_route_break_variables_.push_back(-1);
1278 current_route_break_variables_.push_back(-1);
1279 continue;
1280 }
1281 lp_break_start[br] =
1282 solver->AddVariable(break_start_min, break_start_max);
1283 lp_break_end[br] = solver->AddVariable(break_end_min, break_end_max);
1284 lp_break_duration[br] =
1285 solver->AddVariable(break_duration_min, break_duration_max);
1286 // start + duration = end.
1287 solver->AddLinearConstraint(0, 0,
1288 {{lp_break_end[br], 1},
1289 {lp_break_start[br], -1},
1290 {lp_break_duration[br], -1}});
1291 // Record index of variables
1292 all_break_variables_[all_break_variables_offset + 2 * br] =
1293 lp_break_start[br];
1294 all_break_variables_[all_break_variables_offset + 2 * br + 1] =
1295 lp_break_end[br];
1296 current_route_break_variables_.push_back(lp_break_start[br]);
1297 current_route_break_variables_.push_back(lp_break_end[br]);
1298 } else {
1299 if (break_end_min <= vehicle_start_max ||
1300 vehicle_end_min <= break_start_max) {
1301 all_break_variables_[all_break_variables_offset + 2 * br] = -1;
1302 all_break_variables_[all_break_variables_offset + 2 * br + 1] = -1;
1303 current_route_break_variables_.push_back(-1);
1304 current_route_break_variables_.push_back(-1);
1305 continue;
1306 }
1307 }
1308
1309 // Create a constraint for every break, that forces it to be scheduled
1310 // in exactly one place, i.e. one slack or before/after the route.
1311 // sum_i break_in_slack_i == 1.
1312 const int break_in_one_slack_ct = solver->CreateNewConstraint(1, 1);
1313
1314 if (solver->IsCPSATSolver()) {
1315 // Break can be before route.
1316 if (break_end_min <= vehicle_start_max) {
1317 const int ct = solver->AddLinearConstraint(
1318 0, std::numeric_limits<int64_t>::max(),
1319 {{lp_cumuls.front(), 1}, {lp_break_end[br], -1}});
1320 const int break_is_before_route = solver->AddVariable(0, 1);
1321 solver->SetEnforcementLiteral(ct, break_is_before_route);
1322 solver->SetCoefficient(break_in_one_slack_ct, break_is_before_route, 1);
1323 }
1324 // Break can be after route.
1325 if (vehicle_end_min <= break_start_max) {
1326 const int ct = solver->AddLinearConstraint(
1327 0, std::numeric_limits<int64_t>::max(),
1328 {{lp_break_start[br], 1}, {lp_cumuls.back(), -1}});
1329 const int break_is_after_route = solver->AddVariable(0, 1);
1330 solver->SetEnforcementLiteral(ct, break_is_after_route);
1331 solver->SetCoefficient(break_in_one_slack_ct, break_is_after_route, 1);
1332 }
1333 }
1334
1335 // Add the possibility of fitting the break during each slack where it can.
1336 for (int pos = 0; pos < path_size - 1; ++pos) {
1337 // Pass on slacks that cannot start before, cannot end after,
1338 // or are not long enough to contain the break.
1339 const int64_t slack_start_min =
1340 CapAdd(current_route_min_cumuls_[pos], pre_travel[pos]);
1341 if (slack_start_min > break_start_max) break;
1342 const int64_t slack_end_max =
1343 CapSub(current_route_max_cumuls_[pos + 1], post_travel[pos]);
1344 if (break_end_min > slack_end_max) continue;
1345 const int64_t slack_duration_max =
1346 std::min(CapSub(CapSub(current_route_max_cumuls_[pos + 1],
1347 current_route_min_cumuls_[pos]),
1348 fixed_transit[pos]),
1349 dimension_->SlackVar(path[pos])->Max());
1350 if (slack_duration_max < break_duration_min) continue;
1351
1352 // Break can fit into slack: make LP variable, add to break and slack
1353 // constraints.
1354 // Make a linearized slack lower bound (lazily), that represents
1355 // sum_br break_duration_min(br) * break_in_slack(br, pos) <=
1356 // lp_slacks(pos).
1357 const int break_in_slack = solver->AddVariable(0, 1);
1358 solver->SetCoefficient(break_in_one_slack_ct, break_in_slack, 1);
1359 if (slack_linear_lower_bound_ct[pos] == -1) {
1360 slack_linear_lower_bound_ct[pos] = solver->AddLinearConstraint(
1361 std::numeric_limits<int64_t>::min(), 0, {{lp_slacks[pos], -1}});
1362 }
1363 solver->SetCoefficient(slack_linear_lower_bound_ct[pos], break_in_slack,
1364 break_duration_min);
1365 if (solver->IsCPSATSolver()) {
1366 // Exact relation between breaks, slacks and cumul variables.
1367 // Make an exact slack lower bound (lazily), that represents
1368 // sum_br break_duration(br) * break_in_slack(br, pos) <=
1369 // lp_slacks(pos).
1370 const int break_duration_in_slack =
1371 solver->AddVariable(0, slack_duration_max);
1372 solver->AddProductConstraint(break_duration_in_slack,
1373 {break_in_slack, lp_break_duration[br]});
1374 if (slack_exact_lower_bound_ct[pos] == -1) {
1375 slack_exact_lower_bound_ct[pos] = solver->AddLinearConstraint(
1376 std::numeric_limits<int64_t>::min(), 0, {{lp_slacks[pos], -1}});
1377 }
1378 solver->SetCoefficient(slack_exact_lower_bound_ct[pos],
1379 break_duration_in_slack, 1);
1380 // If break_in_slack_i == 1, then
1381 // 1) break_start >= cumul[pos] + pre_travel[pos]
1382 const int break_start_after_current_ct = solver->AddLinearConstraint(
1383 pre_travel[pos], std::numeric_limits<int64_t>::max(),
1384 {{lp_break_start[br], 1}, {lp_cumuls[pos], -1}});
1385 solver->SetEnforcementLiteral(break_start_after_current_ct,
1386 break_in_slack);
1387 // 2) break_end <= cumul[pos+1] - post_travel[pos]
1388 const int break_ends_before_next_ct = solver->AddLinearConstraint(
1389 post_travel[pos], std::numeric_limits<int64_t>::max(),
1390 {{lp_cumuls[pos + 1], 1}, {lp_break_end[br], -1}});
1391 solver->SetEnforcementLiteral(break_ends_before_next_ct,
1392 break_in_slack);
1393 }
1394 }
1395 }
1396
1397 if (!solver->IsCPSATSolver()) return true;
1398 if (!dimension_->GetBreakDistanceDurationOfVehicle(vehicle).empty()) {
1399 // If there is an optional interval, the following model would be wrong.
1400 // TODO(user): support optional intervals.
1401 for (const IntervalVar* interval :
1402 dimension_->GetBreakIntervalsOfVehicle(vehicle)) {
1403 if (!interval->MustBePerformed()) return true;
1404 }
1405 // When this feature is used, breaks are in sorted order.
1406 for (int br = 1; br < num_breaks; ++br) {
1407 solver->AddLinearConstraint(
1408 0, std::numeric_limits<int64_t>::max(),
1409 {{lp_break_end[br - 1], -1}, {lp_break_start[br], 1}});
1410 }
1411 }
1412 for (const auto& distance_duration :
1413 dimension_->GetBreakDistanceDurationOfVehicle(vehicle)) {
1414 const int64_t limit = distance_duration.first;
1415 const int64_t min_break_duration = distance_duration.second;
1416 // Interbreak limit constraint: breaks are interpreted as being in sorted
1417 // order, and the maximum duration between two consecutive
1418 // breaks of duration more than 'min_break_duration' is 'limit'. This
1419 // considers the time until start of route and after end of route to be
1420 // infinite breaks.
1421 // The model for this constraint adds some 'cover_i' variables, such that
1422 // the breaks up to i and the start of route allows to go without a break.
1423 // With s_i the start of break i and e_i its end:
1424 // - the route start covers time from start to start + limit:
1425 // cover_0 = route_start + limit
1426 // - the coverage up to a given break is the largest of the coverage of the
1427 // previous break and if the break is long enough, break end + limit:
1428 // cover_{i+1} = max(cover_i,
1429 // e_i - s_i >= min_break_duration ? e_i + limit : -inf)
1430 // - the coverage of the last break must be at least the route end,
1431 // to ensure the time point route_end-1 is covered:
1432 // cover_{num_breaks} >= route_end
1433 // - similarly, time point s_i-1 must be covered by breaks up to i-1,
1434 // but only if the cover has not reached the route end.
1435 // For instance, a vehicle could have a choice between two days,
1436 // with a potential break on day 1 and a potential break on day 2,
1437 // but the break of day 1 does not have to cover that of day 2!
1438 // cover_{i-1} < route_end => s_i <= cover_{i-1}
1439 // This is sufficient to ensure that the union of the intervals
1440 // (-infinity, route_start], [route_end, +infinity) and all
1441 // [s_i, e_i+limit) where e_i - s_i >= min_break_duration is
1442 // the whole timeline (-infinity, +infinity).
1443 int previous_cover = solver->AddVariable(CapAdd(vehicle_start_min, limit),
1444 CapAdd(vehicle_start_max, limit));
1445 solver->AddLinearConstraint(limit, limit,
1446 {{previous_cover, 1}, {lp_cumuls.front(), -1}});
1447 for (int br = 0; br < num_breaks; ++br) {
1448 if (lp_break_start[br] == -1) continue;
1449 const int64_t break_end_min = CapSub(breaks[br]->EndMin(), cumul_offset);
1450 const int64_t break_end_max = CapSub(breaks[br]->EndMax(), cumul_offset);
1451 // break_is_eligible <=>
1452 // break_end - break_start >= break_minimum_duration.
1453 const int break_is_eligible = solver->AddVariable(0, 1);
1454 const int break_is_not_eligible = solver->AddVariable(0, 1);
1455 {
1456 solver->AddLinearConstraint(
1457 1, 1, {{break_is_eligible, 1}, {break_is_not_eligible, 1}});
1458 const int positive_ct = solver->AddLinearConstraint(
1459 min_break_duration, std::numeric_limits<int64_t>::max(),
1460 {{lp_break_end[br], 1}, {lp_break_start[br], -1}});
1461 solver->SetEnforcementLiteral(positive_ct, break_is_eligible);
1462 const int negative_ct = solver->AddLinearConstraint(
1463 std::numeric_limits<int64_t>::min(), min_break_duration - 1,
1464 {{lp_break_end[br], 1}, {lp_break_start[br], -1}});
1465 solver->SetEnforcementLiteral(negative_ct, break_is_not_eligible);
1466 }
1467 // break_is_eligible => break_cover == break_end + limit.
1468 // break_is_not_eligible => break_cover == vehicle_start_min + limit.
1469 // break_cover's initial domain is the smallest interval that contains the
1470 // union of sets {vehicle_start_min+limit} and
1471 // [break_end_min+limit, break_end_max+limit).
1472 const int break_cover = solver->AddVariable(
1473 CapAdd(std::min(vehicle_start_min, break_end_min), limit),
1474 CapAdd(std::max(vehicle_start_min, break_end_max), limit));
1475 const int limit_cover_ct = solver->AddLinearConstraint(
1476 limit, limit, {{break_cover, 1}, {lp_break_end[br], -1}});
1477 solver->SetEnforcementLiteral(limit_cover_ct, break_is_eligible);
1478 const int empty_cover_ct = solver->AddLinearConstraint(
1479 CapAdd(vehicle_start_min, limit), CapAdd(vehicle_start_min, limit),
1480 {{break_cover, 1}});
1481 solver->SetEnforcementLiteral(empty_cover_ct, break_is_not_eligible);
1482
1483 const int cover =
1484 solver->AddVariable(CapAdd(vehicle_start_min, limit),
1485 std::numeric_limits<int64_t>::max());
1486 solver->AddMaximumConstraint(cover, {previous_cover, break_cover});
1487 // Cover chaining. If route end is not covered, break start must be:
1488 // cover_{i-1} < route_end => s_i <= cover_{i-1}
1489 const int route_end_is_not_covered = solver->AddReifiedLinearConstraint(
1490 1, std::numeric_limits<int64_t>::max(),
1491 {{lp_cumuls.back(), 1}, {previous_cover, -1}});
1492 const int break_start_cover_ct = solver->AddLinearConstraint(
1493 0, std::numeric_limits<int64_t>::max(),
1494 {{previous_cover, 1}, {lp_break_start[br], -1}});
1495 solver->SetEnforcementLiteral(break_start_cover_ct,
1496 route_end_is_not_covered);
1497
1498 previous_cover = cover;
1499 }
1500 solver->AddLinearConstraint(0, std::numeric_limits<int64_t>::max(),
1501 {{previous_cover, 1}, {lp_cumuls.back(), -1}});
1502 }
1503
1504 return true;
1505 }
1506
SetGlobalConstraints(const std::function<int64_t (int64_t)> & next_accessor,int64_t cumul_offset,bool optimize_costs,RoutingLinearSolverWrapper * solver)1507 bool DimensionCumulOptimizerCore::SetGlobalConstraints(
1508 const std::function<int64_t(int64_t)>& next_accessor, int64_t cumul_offset,
1509 bool optimize_costs, RoutingLinearSolverWrapper* solver) {
1510 // Global span cost =
1511 // global_span_cost_coefficient * (max_end_cumul - min_start_cumul).
1512 const int64_t global_span_coeff = dimension_->global_span_cost_coefficient();
1513 if (optimize_costs && global_span_coeff > 0) {
1514 solver->SetObjectiveCoefficient(max_end_cumul_, global_span_coeff);
1515 solver->SetObjectiveCoefficient(min_start_cumul_, -global_span_coeff);
1516 }
1517
1518 // Node precedence constraints, set when both nodes are visited.
1519 for (const RoutingDimension::NodePrecedence& precedence :
1520 dimension_->GetNodePrecedences()) {
1521 const int first_cumul_var = index_to_cumul_variable_[precedence.first_node];
1522 const int second_cumul_var =
1523 index_to_cumul_variable_[precedence.second_node];
1524 if (first_cumul_var < 0 || second_cumul_var < 0) {
1525 // At least one of the nodes is not on any route, skip this precedence
1526 // constraint.
1527 continue;
1528 }
1529 DCHECK_NE(first_cumul_var, second_cumul_var)
1530 << "Dimension " << dimension_->name()
1531 << " has a self-precedence on node " << precedence.first_node << ".";
1532
1533 // cumul[second_node] - cumul[first_node] >= offset.
1534 const int ct = solver->CreateNewConstraint(
1535 precedence.offset, std::numeric_limits<int64_t>::max());
1536 solver->SetCoefficient(ct, second_cumul_var, 1);
1537 solver->SetCoefficient(ct, first_cumul_var, -1);
1538 }
1539
1540 const RoutingModel& model = *dimension_->model();
1541 if (!solver->IsCPSATSolver()) {
1542 // The resource attributes conditional constraints can only be added with
1543 // the CP-SAT MIP solver.
1544 return true;
1545 }
1546
1547 const int num_vehicles = model.vehicles();
1548 const auto& resource_groups = model.GetResourceGroups();
1549 for (int rg_index : model.GetDimensionResourceGroupIndices(dimension_)) {
1550 // Resource domain constraints:
1551 // Every (used) vehicle requiring a resource from this group must be
1552 // assigned to exactly one resource in this group, and each resource must be
1553 // assigned to at most 1 vehicle requiring it.
1554 // For every resource r with Attributes A = resources[r].attributes(dim)
1555 // and every vehicle v, assign(r, v) == 1 -->
1556 // A.start_domain.Min() <= cumul[Start(v)] <= A.start_domain.Max(),
1557 // and
1558 // A.end_domain.Min() <= cumul[End(v)] <= A.end_domain.Max().
1559 const ResourceGroup& resource_group = *resource_groups[rg_index];
1560 DCHECK(!resource_group.GetVehiclesRequiringAResource().empty());
1561
1562 const std::vector<ResourceGroup::Resource>& resources =
1563 resource_group.GetResources();
1564 int num_required_resources = 0;
1565 static const int kNoConstraint = -1;
1566 // Assignment constraints for vehicles: each (used) vehicle must have
1567 // exactly one resource assigned to it.
1568 std::vector<int> vehicle_constraints(model.vehicles(), kNoConstraint);
1569 for (int v : resource_group.GetVehiclesRequiringAResource()) {
1570 if (model.IsEnd(next_accessor(model.Start(v))) &&
1571 !model.IsVehicleUsedWhenEmpty(v)) {
1572 // We don't assign a driver to unused vehicles.
1573 continue;
1574 }
1575 num_required_resources++;
1576 vehicle_constraints[v] = solver->CreateNewConstraint(1, 1);
1577 }
1578 // Assignment constraints for resources: each resource must be assigned to
1579 // at most one (used) vehicle requiring one.
1580 const int num_resources = resources.size();
1581 std::vector<int> resource_constraints(num_resources, kNoConstraint);
1582 int num_available_resources = 0;
1583 for (int r = 0; r < num_resources; r++) {
1584 const ResourceGroup::Attributes& attributes =
1585 resources[r].GetDimensionAttributes(dimension_);
1586 if (attributes.start_domain().Max() < cumul_offset ||
1587 attributes.end_domain().Max() < cumul_offset) {
1588 // This resource's domain has a cumul max lower than the offset, so it's
1589 // not possible to restrict any vehicle start/end to this domain; skip
1590 // it.
1591 continue;
1592 }
1593 num_available_resources++;
1594 resource_constraints[r] = solver->CreateNewConstraint(0, 1);
1595 }
1596
1597 if (num_required_resources > num_available_resources) {
1598 // There aren't enough resources in this group for vehicles requiring one.
1599 return false;
1600 }
1601
1602 std::vector<int>& resource_to_vehicle_assignment_variables =
1603 resource_group_to_resource_to_vehicle_assignment_variables_[rg_index];
1604 resource_to_vehicle_assignment_variables.assign(
1605 num_resources * num_vehicles, -1);
1606 // Create assignment variables, add them to the corresponding constraints,
1607 // and create the reified constraints assign(r, v) == 1 -->
1608 // A(r).start_domain.Min() <= cumul[Start(v)] <= A(r).start_domain.Max(),
1609 // and
1610 // A(r).end_domain.Min() <= cumul[End(v)] <= A(r).end_domain.Max().
1611 for (int r = 0; r < num_resources; r++) {
1612 if (resource_constraints[r] == kNoConstraint) continue;
1613 const ResourceGroup::Attributes& attributes =
1614 resources[r].GetDimensionAttributes(dimension_);
1615 for (int v : resource_group.GetVehiclesRequiringAResource()) {
1616 if (vehicle_constraints[v] == kNoConstraint) continue;
1617
1618 const int assign_r_to_v = solver->AddVariable(0, 1);
1619 resource_to_vehicle_assignment_variables[r * num_vehicles + v] =
1620 assign_r_to_v;
1621 solver->SetCoefficient(vehicle_constraints[v], assign_r_to_v, 1);
1622 solver->SetCoefficient(resource_constraints[r], assign_r_to_v, 1);
1623
1624 const auto& add_domain_constraint =
1625 [&solver, cumul_offset, assign_r_to_v](const Domain& domain,
1626 int cumul_variable) {
1627 if (domain == Domain::AllValues()) {
1628 return;
1629 }
1630 ClosedInterval cumul_bounds;
1631 if (!GetDomainOffsetBounds(domain, cumul_offset, &cumul_bounds)) {
1632 // This domain cannot be assigned to this vehicle.
1633 solver->SetVariableBounds(assign_r_to_v, 0, 0);
1634 return;
1635 }
1636 const int cumul_constraint = solver->AddLinearConstraint(
1637 cumul_bounds.start, cumul_bounds.end, {{cumul_variable, 1}});
1638 solver->SetEnforcementLiteral(cumul_constraint, assign_r_to_v);
1639 };
1640 add_domain_constraint(attributes.start_domain(),
1641 index_to_cumul_variable_[model.Start(v)]);
1642 add_domain_constraint(attributes.end_domain(),
1643 index_to_cumul_variable_[model.End(v)]);
1644 }
1645 }
1646 }
1647 return true;
1648 }
1649
SetValuesFromLP(const std::vector<int> & lp_variables,int64_t offset,RoutingLinearSolverWrapper * solver,std::vector<int64_t> * lp_values) const1650 void DimensionCumulOptimizerCore::SetValuesFromLP(
1651 const std::vector<int>& lp_variables, int64_t offset,
1652 RoutingLinearSolverWrapper* solver, std::vector<int64_t>* lp_values) const {
1653 if (lp_values == nullptr) return;
1654 lp_values->assign(lp_variables.size(), std::numeric_limits<int64_t>::min());
1655 for (int i = 0; i < lp_variables.size(); i++) {
1656 const int lp_var = lp_variables[i];
1657 if (lp_var < 0) continue; // Keep default value, kint64min.
1658 const double lp_value_double = solver->GetValue(lp_var);
1659 const int64_t lp_value_int64 =
1660 (lp_value_double >= std::numeric_limits<int64_t>::max())
1661 ? std::numeric_limits<int64_t>::max()
1662 : MathUtil::FastInt64Round(lp_value_double);
1663 (*lp_values)[i] = CapAdd(lp_value_int64, offset);
1664 }
1665 }
1666
SetResourceIndices(RoutingLinearSolverWrapper * solver,std::vector<std::vector<int>> * resource_indices_per_group) const1667 void DimensionCumulOptimizerCore::SetResourceIndices(
1668 RoutingLinearSolverWrapper* solver,
1669 std::vector<std::vector<int>>* resource_indices_per_group) const {
1670 if (resource_indices_per_group == nullptr ||
1671 resource_group_to_resource_to_vehicle_assignment_variables_.empty()) {
1672 return;
1673 }
1674 const RoutingModel& model = *dimension_->model();
1675 const int num_vehicles = model.vehicles();
1676 DCHECK(!model.GetDimensionResourceGroupIndices(dimension_).empty());
1677 const auto& resource_groups = model.GetResourceGroups();
1678 resource_indices_per_group->resize(resource_groups.size());
1679 for (int rg_index : model.GetDimensionResourceGroupIndices(dimension_)) {
1680 const ResourceGroup& resource_group = *resource_groups[rg_index];
1681 DCHECK(!resource_group.GetVehiclesRequiringAResource().empty());
1682
1683 const int num_resources = resource_group.Size();
1684 std::vector<int>& resource_indices =
1685 resource_indices_per_group->at(rg_index);
1686 resource_indices.assign(num_vehicles, -1);
1687 // Find the resource assigned to each vehicle.
1688 const std::vector<int>& resource_to_vehicle_assignment_variables =
1689 resource_group_to_resource_to_vehicle_assignment_variables_[rg_index];
1690 DCHECK_EQ(resource_to_vehicle_assignment_variables.size(),
1691 num_resources * num_vehicles);
1692 for (int v : resource_group.GetVehiclesRequiringAResource()) {
1693 for (int r = 0; r < num_resources; r++) {
1694 const int assignment_var =
1695 resource_to_vehicle_assignment_variables[r * num_vehicles + v];
1696 if (assignment_var >= 0 && solver->GetValue(assignment_var) == 1) {
1697 // This resource is assigned to this vehicle.
1698 resource_indices[v] = r;
1699 break;
1700 }
1701 }
1702 }
1703 }
1704 }
1705
1706 // GlobalDimensionCumulOptimizer
1707
GlobalDimensionCumulOptimizer(const RoutingDimension * dimension,RoutingSearchParameters::SchedulingSolver solver_type)1708 GlobalDimensionCumulOptimizer::GlobalDimensionCumulOptimizer(
1709 const RoutingDimension* dimension,
1710 RoutingSearchParameters::SchedulingSolver solver_type)
1711 : optimizer_core_(dimension,
1712 /*use_precedence_propagator=*/
1713 !dimension->GetNodePrecedences().empty()) {
1714 switch (solver_type) {
1715 case RoutingSearchParameters::GLOP: {
1716 solver_ = absl::make_unique<RoutingGlopWrapper>(
1717 /*is_relaxation=*/!dimension->model()
1718 ->GetDimensionResourceGroupIndices(dimension)
1719 .empty(),
1720 GetGlopParametersForGlobalLP());
1721 break;
1722 }
1723 case RoutingSearchParameters::CP_SAT: {
1724 solver_ = absl::make_unique<RoutingCPSatWrapper>();
1725 break;
1726 }
1727 default:
1728 LOG(DFATAL) << "Unrecognized solver type: " << solver_type;
1729 }
1730 }
1731
1732 DimensionSchedulingStatus
ComputeCumulCostWithoutFixedTransits(const std::function<int64_t (int64_t)> & next_accessor,int64_t * optimal_cost_without_transits)1733 GlobalDimensionCumulOptimizer::ComputeCumulCostWithoutFixedTransits(
1734 const std::function<int64_t(int64_t)>& next_accessor,
1735 int64_t* optimal_cost_without_transits) {
1736 int64_t cost = 0;
1737 int64_t transit_cost = 0;
1738 DimensionSchedulingStatus status =
1739 optimizer_core_.Optimize(next_accessor, solver_.get(), nullptr, nullptr,
1740 nullptr, &cost, &transit_cost);
1741 if (status != DimensionSchedulingStatus::INFEASIBLE &&
1742 optimal_cost_without_transits != nullptr) {
1743 *optimal_cost_without_transits = CapSub(cost, transit_cost);
1744 }
1745 return status;
1746 }
1747
ComputeCumuls(const std::function<int64_t (int64_t)> & next_accessor,std::vector<int64_t> * optimal_cumuls,std::vector<int64_t> * optimal_breaks,std::vector<std::vector<int>> * optimal_resource_indices)1748 DimensionSchedulingStatus GlobalDimensionCumulOptimizer::ComputeCumuls(
1749 const std::function<int64_t(int64_t)>& next_accessor,
1750 std::vector<int64_t>* optimal_cumuls, std::vector<int64_t>* optimal_breaks,
1751 std::vector<std::vector<int>>* optimal_resource_indices) {
1752 return optimizer_core_.Optimize(next_accessor, solver_.get(), optimal_cumuls,
1753 optimal_breaks, optimal_resource_indices,
1754 nullptr, nullptr);
1755 }
1756
ComputePackedCumuls(const std::function<int64_t (int64_t)> & next_accessor,std::vector<int64_t> * packed_cumuls,std::vector<int64_t> * packed_breaks,std::vector<std::vector<int>> * resource_indices)1757 DimensionSchedulingStatus GlobalDimensionCumulOptimizer::ComputePackedCumuls(
1758 const std::function<int64_t(int64_t)>& next_accessor,
1759 std::vector<int64_t>* packed_cumuls, std::vector<int64_t>* packed_breaks,
1760 std::vector<std::vector<int>>* resource_indices) {
1761 return optimizer_core_.OptimizeAndPack(next_accessor, solver_.get(),
1762 packed_cumuls, packed_breaks,
1763 resource_indices);
1764 }
1765
1766 // ResourceAssignmentOptimizer
1767
ResourceAssignmentOptimizer(const RoutingModel::ResourceGroup * resource_group,LocalDimensionCumulOptimizer * optimizer,LocalDimensionCumulOptimizer * mp_optimizer)1768 ResourceAssignmentOptimizer::ResourceAssignmentOptimizer(
1769 const RoutingModel::ResourceGroup* resource_group,
1770 LocalDimensionCumulOptimizer* optimizer,
1771 LocalDimensionCumulOptimizer* mp_optimizer)
1772 : optimizer_(*optimizer),
1773 mp_optimizer_(*mp_optimizer),
1774 model_(*optimizer->dimension()->model()),
1775 resource_group_(*resource_group) {}
1776
ComputeAssignmentCostsForVehicle(int v,const std::function<int64_t (int64_t)> & next_accessor,bool optimize_vehicle_costs,std::vector<int64_t> * assignment_costs,std::vector<std::vector<int64_t>> * cumul_values,std::vector<std::vector<int64_t>> * break_values)1777 bool ResourceAssignmentOptimizer::ComputeAssignmentCostsForVehicle(
1778 int v, const std::function<int64_t(int64_t)>& next_accessor,
1779 bool optimize_vehicle_costs, std::vector<int64_t>* assignment_costs,
1780 std::vector<std::vector<int64_t>>* cumul_values,
1781 std::vector<std::vector<int64_t>>* break_values) {
1782 DCHECK_NE(assignment_costs, nullptr);
1783 if (!resource_group_.VehicleRequiresAResource(v) ||
1784 (next_accessor(model_.Start(v)) == model_.End(v) &&
1785 !model_.IsVehicleUsedWhenEmpty(v))) {
1786 assignment_costs->clear();
1787 return true;
1788 }
1789 const RoutingDimension& dimension = *optimizer_.dimension();
1790 if (dimension.model()->CheckLimit()) {
1791 // The model's time limit has been reached, stop everything.
1792 return false;
1793 }
1794
1795 const std::vector<ResourceGroup::Resource>& resources =
1796 resource_group_.GetResources();
1797 const int num_resources = resources.size();
1798 std::vector<int> all_resource_indices(num_resources);
1799 std::iota(all_resource_indices.begin(), all_resource_indices.end(), 0);
1800 const bool use_mp_optimizer =
1801 dimension.HasBreakConstraints() &&
1802 !dimension.GetBreakIntervalsOfVehicle(v).empty();
1803 LocalDimensionCumulOptimizer& optimizer =
1804 use_mp_optimizer ? mp_optimizer_ : optimizer_;
1805 std::vector<DimensionSchedulingStatus> statuses =
1806 optimizer.ComputeRouteCumulCostsForResourcesWithoutFixedTransits(
1807 v, next_accessor, resources, all_resource_indices,
1808 optimize_vehicle_costs, assignment_costs, cumul_values, break_values);
1809
1810 if (assignment_costs->empty()) {
1811 // Couldn't assign any resource to this vehicle.
1812 return false;
1813 }
1814 DCHECK_EQ(assignment_costs->size(), num_resources);
1815 DCHECK_EQ(statuses.size(), num_resources);
1816 DCHECK(cumul_values == nullptr || cumul_values->size() == num_resources);
1817 DCHECK(break_values == nullptr || break_values->size() == num_resources);
1818
1819 if (use_mp_optimizer) {
1820 // We already used the mp optimizer, so we don't need to recompute anything.
1821 // If all assignment costs are negative, it means no resource is feasible
1822 // for this vehicle.
1823 return absl::c_any_of(*assignment_costs,
1824 [](int64_t cost) { return cost >= 0; });
1825 }
1826
1827 std::vector<int> mp_optimizer_resource_indices;
1828 for (int r = 0; r < num_resources; r++) {
1829 if (statuses[r] == DimensionSchedulingStatus::RELAXED_OPTIMAL_ONLY) {
1830 mp_optimizer_resource_indices.push_back(r);
1831 }
1832 }
1833
1834 std::vector<int64_t> mp_assignment_costs;
1835 std::vector<std::vector<int64_t>> mp_cumul_values;
1836 std::vector<std::vector<int64_t>> mp_break_values;
1837 mp_optimizer_.ComputeRouteCumulCostsForResourcesWithoutFixedTransits(
1838 v, next_accessor, resources, mp_optimizer_resource_indices,
1839 optimize_vehicle_costs, &mp_assignment_costs,
1840 cumul_values == nullptr ? nullptr : &mp_cumul_values,
1841 break_values == nullptr ? nullptr : &mp_break_values);
1842 if (!mp_optimizer_resource_indices.empty() && mp_assignment_costs.empty()) {
1843 // A timeout was reached during optimization.
1844 return false;
1845 }
1846 DCHECK_EQ(mp_assignment_costs.size(), mp_optimizer_resource_indices.size());
1847 DCHECK(cumul_values == nullptr ||
1848 mp_cumul_values.size() == mp_optimizer_resource_indices.size());
1849 DCHECK(break_values == nullptr ||
1850 mp_break_values.size() == mp_optimizer_resource_indices.size());
1851 for (int i = 0; i < mp_optimizer_resource_indices.size(); i++) {
1852 assignment_costs->at(mp_optimizer_resource_indices[i]) =
1853 mp_assignment_costs[i];
1854 if (cumul_values != nullptr) {
1855 cumul_values->at(mp_optimizer_resource_indices[i])
1856 .swap(mp_cumul_values[i]);
1857 }
1858 if (break_values != nullptr) {
1859 break_values->at(mp_optimizer_resource_indices[i])
1860 .swap(mp_break_values[i]);
1861 }
1862 }
1863 return absl::c_any_of(*assignment_costs,
1864 [](int64_t cost) { return cost >= 0; });
1865 }
1866
ComputeBestAssignmentCost(const std::vector<std::vector<int64_t>> & primary_vehicle_to_resource_assignment_costs,const std::vector<std::vector<int64_t>> & secondary_vehicle_to_resource_assignment_costs,const std::function<bool (int)> & use_primary_for_vehicle,std::vector<int> * resource_indices) const1867 int64_t ResourceAssignmentOptimizer::ComputeBestAssignmentCost(
1868 const std::vector<std::vector<int64_t>>&
1869 primary_vehicle_to_resource_assignment_costs,
1870 const std::vector<std::vector<int64_t>>&
1871 secondary_vehicle_to_resource_assignment_costs,
1872 const std::function<bool(int)>& use_primary_for_vehicle,
1873 std::vector<int>* resource_indices) const {
1874 const int num_vehicles = model_.vehicles();
1875 DCHECK_EQ(primary_vehicle_to_resource_assignment_costs.size(), num_vehicles);
1876 DCHECK_EQ(secondary_vehicle_to_resource_assignment_costs.size(),
1877 num_vehicles);
1878 const int num_resources = resource_group_.Size();
1879
1880 SimpleMinCostFlow flow(
1881 /*reserve_num_nodes*/ 2 + num_vehicles + num_resources,
1882 /*reserve_num_arcs*/ num_vehicles + num_vehicles * num_resources +
1883 num_resources);
1884 const int source_index = num_vehicles + num_resources;
1885 const int sink_index = source_index + 1;
1886 const auto resource_index = [num_vehicles](int r) {
1887 return num_vehicles + r;
1888 };
1889
1890 std::vector<std::vector<ArcIndex>> vehicle_to_resource_arc_index;
1891 if (resource_indices != nullptr) {
1892 vehicle_to_resource_arc_index.resize(
1893 num_vehicles, std::vector<ArcIndex>(num_resources, -1));
1894 }
1895 int num_used_vehicles = 0;
1896 for (int v : resource_group_.GetVehiclesRequiringAResource()) {
1897 DCHECK(use_primary_for_vehicle(v) ||
1898 primary_vehicle_to_resource_assignment_costs[v].empty());
1899 const std::vector<int64_t>& assignment_costs =
1900 use_primary_for_vehicle(v)
1901 ? primary_vehicle_to_resource_assignment_costs[v]
1902 : secondary_vehicle_to_resource_assignment_costs[v];
1903 if (assignment_costs.empty()) {
1904 // We don't need a resource for this vehicle.
1905 continue;
1906 }
1907 DCHECK_EQ(assignment_costs.size(), num_resources);
1908 num_used_vehicles++;
1909 DCHECK_LE(num_used_vehicles, num_resources)
1910 << num_used_vehicles << " used vehicles and only " << num_resources
1911 << " resources available!";
1912
1913 // Add a source->vehicle arc to the flow.
1914 flow.AddArcWithCapacityAndUnitCost(source_index, v, 1, 0);
1915
1916 // Add arcs to the min-cost-flow graph.
1917 bool has_feasible_resource = false;
1918 for (int r = 0; r < num_resources; r++) {
1919 const int64_t assignment_cost = assignment_costs[r];
1920 if (assignment_cost < 0) continue;
1921 const ArcIndex arc_index = flow.AddArcWithCapacityAndUnitCost(
1922 v, resource_index(r), 1, assignment_cost);
1923 if (!vehicle_to_resource_arc_index.empty()) {
1924 vehicle_to_resource_arc_index[v][r] = arc_index;
1925 }
1926 has_feasible_resource = true;
1927 }
1928 DCHECK(has_feasible_resource)
1929 << "No feasible resource for vehicle " << v
1930 << ", should've been caught by ComputeAssignmentCostsForVehicle()";
1931 }
1932
1933 // Add resource->sink arcs to the flow.
1934 for (int r = 0; r < num_resources; r++) {
1935 flow.AddArcWithCapacityAndUnitCost(resource_index(r), sink_index, 1, 0);
1936 }
1937
1938 // Set the flow supply.
1939 flow.SetNodeSupply(source_index, num_used_vehicles);
1940 flow.SetNodeSupply(sink_index, -num_used_vehicles);
1941
1942 // Solve the min-cost flow and return its cost.
1943 if (flow.Solve() != SimpleMinCostFlow::OPTIMAL) {
1944 if (resource_indices != nullptr) resource_indices->clear();
1945 return -1;
1946 }
1947
1948 const int64_t cost = flow.OptimalCost();
1949 if (resource_indices == nullptr) {
1950 return cost;
1951 }
1952
1953 // Fill the resource indices corresponding to the min-cost assignment.
1954 resource_indices->assign(num_vehicles, -1);
1955 for (int v : resource_group_.GetVehiclesRequiringAResource()) {
1956 for (int r = 0; r < num_resources; r++) {
1957 if (vehicle_to_resource_arc_index[v][r] >= 0 &&
1958 flow.Flow(vehicle_to_resource_arc_index[v][r]) > 0) {
1959 resource_indices->at(v) = r;
1960 break;
1961 }
1962 }
1963 }
1964 return cost;
1965 }
1966
1967 } // namespace operations_research
1968