OR-Tools  8.2
cp_model_lns.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 <limits>
17#include <numeric>
18#include <vector>
19
20#include "ortools/sat/cp_model.pb.h"
23#include "ortools/sat/integer.h"
25#include "ortools/sat/rins.h"
28
29namespace operations_research {
30namespace sat {
31
33 CpModelProto const* model_proto, SatParameters const* parameters,
34 SharedResponseManager* shared_response, SharedTimeLimit* shared_time_limit,
35 SharedBoundsManager* shared_bounds)
36 : SubSolver(""),
37 parameters_(*parameters),
38 model_proto_(*model_proto),
39 shared_time_limit_(shared_time_limit),
40 shared_bounds_(shared_bounds),
41 shared_response_(shared_response) {
42 CHECK(shared_response_ != nullptr);
43 if (shared_bounds_ != nullptr) {
44 shared_bounds_id_ = shared_bounds_->RegisterNewId();
45 }
46 *model_proto_with_only_variables_.mutable_variables() =
47 model_proto_.variables();
48 RecomputeHelperData();
50}
51
53 absl::MutexLock mutex_lock(&mutex_);
54 if (shared_bounds_ != nullptr) {
55 std::vector<int> model_variables;
56 std::vector<int64> new_lower_bounds;
57 std::vector<int64> new_upper_bounds;
58 shared_bounds_->GetChangedBounds(shared_bounds_id_, &model_variables,
59 &new_lower_bounds, &new_upper_bounds);
60
61 for (int i = 0; i < model_variables.size(); ++i) {
62 const int var = model_variables[i];
63 const int64 new_lb = new_lower_bounds[i];
64 const int64 new_ub = new_upper_bounds[i];
65 if (VLOG_IS_ON(3)) {
66 const auto& domain =
67 model_proto_with_only_variables_.variables(var).domain();
68 const int64 old_lb = domain.Get(0);
69 const int64 old_ub = domain.Get(domain.size() - 1);
70 VLOG(3) << "Variable: " << var << " old domain: [" << old_lb << ", "
71 << old_ub << "] new domain: [" << new_lb << ", " << new_ub
72 << "]";
73 }
74 const Domain old_domain =
75 ReadDomainFromProto(model_proto_with_only_variables_.variables(var));
76 const Domain new_domain =
77 old_domain.IntersectionWith(Domain(new_lb, new_ub));
78 if (new_domain.IsEmpty()) {
79 // This can mean two things:
80 // 1/ This variable is a normal one and the problem is UNSAT or
81 // 2/ This variable is optional, and its associated literal must be
82 // set to false.
83 //
84 // Currently, we wait for any full solver to pick the crossing bounds
85 // and do the correct stuff on their own. We do not want to have empty
86 // domain in the proto as this would means INFEASIBLE. So we just ignore
87 // such bounds here.
88 //
89 // TODO(user): We could set the optional literal to false directly in
90 // the bound sharing manager. We do have to be careful that all the
91 // different solvers have the same optionality definition though.
92 continue;
93 }
95 new_domain, model_proto_with_only_variables_.mutable_variables(var));
96 }
97
98 // Only trigger the computation if needed.
99 if (!model_variables.empty()) {
100 RecomputeHelperData();
101 }
102 }
103}
104
105void NeighborhoodGeneratorHelper::RecomputeHelperData() {
106 // Recompute all the data in case new variables have been fixed.
107 //
108 // TODO(user): Ideally we should ignore trivially true/false constraint, but
109 // this will duplicate already existing code :-( we should probably still do
110 // at least enforcement literal and clauses? We could maybe run a light
111 // presolve?
112 var_to_constraint_.assign(model_proto_.variables_size(), {});
113 constraint_to_var_.assign(model_proto_.constraints_size(), {});
114 const auto register_var = [&](int var, int ct_index) {
116 if (IsConstant(var)) return;
117 var_to_constraint_[var].push_back(ct_index);
118 constraint_to_var_[ct_index].push_back(var);
119 CHECK_GE(var, 0);
120 CHECK_LT(var, model_proto_.variables_size());
121 };
122
123 for (int ct_index = 0; ct_index < model_proto_.constraints_size();
124 ++ct_index) {
125 for (const int var : UsedVariables(model_proto_.constraints(ct_index))) {
126 register_var(var, ct_index);
127 }
128
129 // We replace intervals by their underlying integer variables.
130 if (parameters_.lns_expand_intervals_in_constraint_graph()) {
131 for (const int interval :
132 UsedIntervals(model_proto_.constraints(ct_index))) {
133 for (const int var :
134 UsedVariables(model_proto_.constraints(interval))) {
135 register_var(var, ct_index);
136 }
137 }
138 }
139 }
140
141 type_to_constraints_.clear();
142 const int num_constraints = model_proto_.constraints_size();
143 for (int c = 0; c < num_constraints; ++c) {
144 const int type = model_proto_.constraints(c).constraint_case();
145 if (type >= type_to_constraints_.size()) {
146 type_to_constraints_.resize(type + 1);
147 }
148 type_to_constraints_[type].push_back(c);
149 }
150
151 active_variables_.clear();
152 active_variables_set_.assign(model_proto_.variables_size(), false);
153
154 if (parameters_.lns_focus_on_decision_variables()) {
155 for (const auto& search_strategy : model_proto_.search_strategy()) {
156 for (const int var : search_strategy.variables()) {
157 const int pos_var = PositiveRef(var);
158 if (!active_variables_set_[pos_var] && !IsConstant(pos_var)) {
159 active_variables_set_[pos_var] = true;
160 active_variables_.push_back(pos_var);
161 }
162 }
163 }
164
165 // Revert to no focus if active_variables_ is empty().
166 if (!active_variables_.empty()) return;
167 }
168
169 // Add all non-constant variables.
170 for (int i = 0; i < model_proto_.variables_size(); ++i) {
171 if (!IsConstant(i)) {
172 active_variables_.push_back(i);
173 active_variables_set_[i] = true;
174 }
175 }
176}
177
179 return active_variables_set_[var];
180}
181
182bool NeighborhoodGeneratorHelper::IsConstant(int var) const {
183 return model_proto_with_only_variables_.variables(var).domain_size() == 2 &&
184 model_proto_with_only_variables_.variables(var).domain(0) ==
185 model_proto_with_only_variables_.variables(var).domain(1);
186}
187
189 Neighborhood neighborhood;
190 neighborhood.is_reduced = false;
191 neighborhood.is_generated = true;
192 neighborhood.cp_model = model_proto_;
193 *neighborhood.cp_model.mutable_variables() =
194 model_proto_with_only_variables_.variables();
195 return neighborhood;
196}
197
199 const CpSolverResponse& initial_solution) const {
200 std::vector<int> active_intervals;
201 for (const int i : TypeToConstraints(ConstraintProto::kInterval)) {
202 const ConstraintProto& interval_ct = ModelProto().constraints(i);
203 // We only look at intervals that are performed in the solution. The
204 // unperformed intervals should be automatically freed during the generation
205 // phase.
206 if (interval_ct.enforcement_literal().size() == 1) {
207 const int enforcement_ref = interval_ct.enforcement_literal(0);
208 const int enforcement_var = PositiveRef(enforcement_ref);
209 const int value = initial_solution.solution(enforcement_var);
210 if (RefIsPositive(enforcement_ref) == (value == 0)) {
211 continue;
212 }
213 }
214
215 // We filter out fixed intervals. Because of presolve, if there is an
216 // enforcement literal, it cannot be fixed.
217 if (interval_ct.enforcement_literal().empty() &&
218 IsConstant(PositiveRef(interval_ct.interval().start())) &&
219 IsConstant(PositiveRef(interval_ct.interval().size())) &&
220 IsConstant(PositiveRef(interval_ct.interval().end()))) {
221 continue;
222 }
223
224 active_intervals.push_back(i);
225 }
226 return active_intervals;
227}
228
230 const CpSolverResponse& initial_solution,
231 const std::vector<int>& variables_to_fix) const {
232 // TODO(user,user): Do not include constraint with all fixed variables to
233 // save memory and speed-up LNS presolving.
234 Neighborhood neighborhood = FullNeighborhood();
235
236 // Set the current solution as a hint.
237 neighborhood.cp_model.clear_solution_hint();
238 for (int var = 0; var < neighborhood.cp_model.variables_size(); ++var) {
239 neighborhood.cp_model.mutable_solution_hint()->add_vars(var);
240 neighborhood.cp_model.mutable_solution_hint()->add_values(
241 initial_solution.solution(var));
242 }
243
244 neighborhood.is_reduced = !variables_to_fix.empty();
245 if (!neighborhood.is_reduced) return neighborhood;
246 CHECK_EQ(initial_solution.solution_size(),
247 neighborhood.cp_model.variables_size());
248 for (const int var : variables_to_fix) {
249 neighborhood.cp_model.mutable_variables(var)->clear_domain();
250 neighborhood.cp_model.mutable_variables(var)->add_domain(
251 initial_solution.solution(var));
252 neighborhood.cp_model.mutable_variables(var)->add_domain(
253 initial_solution.solution(var));
254 }
255
256 // TODO(user): force better objective? Note that this is already done when the
257 // hint above is successfully loaded (i.e. if it passes the presolve
258 // correctly) since the solver will try to find better solution than the
259 // current one.
260 return neighborhood;
261}
262
264 const std::vector<int>& constraints_to_remove) const {
265 // TODO(user,user): Do not include constraint with all fixed variables to
266 // save memory and speed-up LNS presolving.
267 Neighborhood neighborhood = FullNeighborhood();
268
269 if (constraints_to_remove.empty()) return neighborhood;
270 neighborhood.is_reduced = false;
271 for (const int constraint : constraints_to_remove) {
272 neighborhood.cp_model.mutable_constraints(constraint)->Clear();
273 }
274
275 return neighborhood;
276}
277
279 const CpSolverResponse& initial_solution,
280 const std::vector<int>& relaxed_variables) const {
281 std::vector<bool> relaxed_variables_set(model_proto_.variables_size(), false);
282 for (const int var : relaxed_variables) relaxed_variables_set[var] = true;
283 std::vector<int> fixed_variables;
284 for (const int i : active_variables_) {
285 if (!relaxed_variables_set[i]) {
286 fixed_variables.push_back(i);
287 }
288 }
289 return FixGivenVariables(initial_solution, fixed_variables);
290}
291
293 const CpSolverResponse& initial_solution) const {
294 std::vector<int> fixed_variables;
295 for (const int i : active_variables_) {
296 fixed_variables.push_back(i);
297 }
298 return FixGivenVariables(initial_solution, fixed_variables);
299}
300
303}
304
305double NeighborhoodGenerator::GetUCBScore(int64 total_num_calls) const {
306 absl::MutexLock mutex_lock(&mutex_);
307 DCHECK_GE(total_num_calls, num_calls_);
308 if (num_calls_ <= 10) return std::numeric_limits<double>::infinity();
309 return current_average_ + sqrt((2 * log(total_num_calls)) / num_calls_);
310}
311
313 absl::MutexLock mutex_lock(&mutex_);
314
315 // To make the whole update process deterministic, we currently sort the
316 // SolveData.
317 std::sort(solve_data_.begin(), solve_data_.end());
318
319 // This will be used to update the difficulty of this neighborhood.
320 int num_fully_solved_in_batch = 0;
321 int num_not_fully_solved_in_batch = 0;
322
323 for (const SolveData& data : solve_data_) {
325 ++num_calls_;
326
327 // INFEASIBLE or OPTIMAL means that we "fully solved" the local problem.
328 // If we didn't, then we cannot be sure that there is no improving solution
329 // in that neighborhood.
330 if (data.status == CpSolverStatus::INFEASIBLE ||
331 data.status == CpSolverStatus::OPTIMAL) {
332 ++num_fully_solved_calls_;
333 ++num_fully_solved_in_batch;
334 } else {
335 ++num_not_fully_solved_in_batch;
336 }
337
338 // It seems to make more sense to compare the new objective to the base
339 // solution objective, not the best one. However this causes issue in the
340 // logic below because on some problems the neighborhood can always lead
341 // to a better "new objective" if the base solution wasn't the best one.
342 //
343 // This might not be a final solution, but it does work ok for now.
344 const IntegerValue best_objective_improvement =
346 ? IntegerValue(CapSub(data.new_objective_bound.value(),
347 data.initial_best_objective_bound.value()))
348 : IntegerValue(CapSub(data.initial_best_objective.value(),
349 data.new_objective.value()));
350 if (best_objective_improvement > 0) {
351 num_consecutive_non_improving_calls_ = 0;
352 } else {
353 ++num_consecutive_non_improving_calls_;
354 }
355
356 // TODO(user): Weight more recent data.
357 // degrade the current average to forget old learnings.
358 const double gain_per_time_unit =
359 std::max(0.0, static_cast<double>(best_objective_improvement.value())) /
360 (1.0 + data.deterministic_time);
361 if (num_calls_ <= 100) {
362 current_average_ += (gain_per_time_unit - current_average_) / num_calls_;
363 } else {
364 current_average_ = 0.9 * current_average_ + 0.1 * gain_per_time_unit;
365 }
366
367 deterministic_time_ += data.deterministic_time;
368 }
369
370 // Update the difficulty.
371 difficulty_.Update(/*num_decreases=*/num_not_fully_solved_in_batch,
372 /*num_increases=*/num_fully_solved_in_batch);
373
374 // Bump the time limit if we saw no better solution in the last few calls.
375 // This means that as the search progress, we likely spend more and more time
376 // trying to solve individual neighborhood.
377 //
378 // TODO(user): experiment with resetting the time limit if a solution is
379 // found.
380 if (num_consecutive_non_improving_calls_ > 50) {
381 num_consecutive_non_improving_calls_ = 0;
382 deterministic_limit_ *= 1.02;
383
384 // We do not want the limit to go to high. Intuitively, the goal is to try
385 // out a lot of neighborhoods, not just spend a lot of time on a few.
386 deterministic_limit_ = std::min(60.0, deterministic_limit_);
387 }
388
389 solve_data_.clear();
390}
391
392namespace {
393
394void GetRandomSubset(double relative_size, std::vector<int>* base,
395 absl::BitGenRef random) {
396 // TODO(user): we could generate this more efficiently than using random
397 // shuffle.
398 std::shuffle(base->begin(), base->end(), random);
399 const int target_size = std::round(relative_size * base->size());
400 base->resize(target_size);
401}
402
403} // namespace
404
406 const CpSolverResponse& initial_solution, double difficulty,
407 absl::BitGenRef random) {
408 std::vector<int> fixed_variables = helper_.ActiveVariables();
409 GetRandomSubset(1.0 - difficulty, &fixed_variables, random);
410 return helper_.FixGivenVariables(initial_solution, fixed_variables);
411}
412
414 const CpSolverResponse& initial_solution, double difficulty,
415 absl::BitGenRef random) {
416 std::vector<int> active_constraints;
417 for (int ct = 0; ct < helper_.ModelProto().constraints_size(); ++ct) {
418 if (helper_.ModelProto().constraints(ct).constraint_case() ==
419 ConstraintProto::CONSTRAINT_NOT_SET) {
420 continue;
421 }
422 active_constraints.push_back(ct);
423 }
424
425 const int num_active_vars = helper_.ActiveVariables().size();
426 const int num_model_vars = helper_.ModelProto().variables_size();
427 const int target_size = std::ceil(difficulty * num_active_vars);
428 const int num_constraints = helper_.ConstraintToVar().size();
429 if (num_constraints == 0 || target_size == num_active_vars) {
430 return helper_.FullNeighborhood();
431 }
432 CHECK_GT(target_size, 0);
433
434 std::shuffle(active_constraints.begin(), active_constraints.end(), random);
435
436 std::vector<bool> visited_variables_set(num_model_vars, false);
437 std::vector<int> relaxed_variables;
438
439 for (const int constraint_index : active_constraints) {
440 CHECK_LT(constraint_index, num_constraints);
441 for (const int var : helper_.ConstraintToVar()[constraint_index]) {
442 if (visited_variables_set[var]) continue;
443 visited_variables_set[var] = true;
444 if (helper_.IsActive(var)) {
445 relaxed_variables.push_back(var);
446 if (relaxed_variables.size() == target_size) break;
447 }
448 }
449 if (relaxed_variables.size() == target_size) break;
450 }
451 return helper_.RelaxGivenVariables(initial_solution, relaxed_variables);
452}
453
455 const CpSolverResponse& initial_solution, double difficulty,
456 absl::BitGenRef random) {
457 const int num_active_vars = helper_.ActiveVariables().size();
458 const int num_model_vars = helper_.ModelProto().variables_size();
459 const int target_size = std::ceil(difficulty * num_active_vars);
460 if (target_size == num_active_vars) {
461 return helper_.FullNeighborhood();
462 }
463 CHECK_GT(target_size, 0) << difficulty << " " << num_active_vars;
464
465 std::vector<bool> visited_variables_set(num_model_vars, false);
466 std::vector<int> relaxed_variables;
467 std::vector<int> visited_variables;
468
469 const int first_var =
470 helper_.ActiveVariables()[absl::Uniform<int>(random, 0, num_active_vars)];
471 visited_variables_set[first_var] = true;
472 visited_variables.push_back(first_var);
473 relaxed_variables.push_back(first_var);
474
475 std::vector<int> random_variables;
476 for (int i = 0; i < visited_variables.size(); ++i) {
477 random_variables.clear();
478 // Collect all the variables that appears in the same constraints as
479 // visited_variables[i].
480 for (const int ct : helper_.VarToConstraint()[visited_variables[i]]) {
481 for (const int var : helper_.ConstraintToVar()[ct]) {
482 if (visited_variables_set[var]) continue;
483 visited_variables_set[var] = true;
484 random_variables.push_back(var);
485 }
486 }
487 // We always randomize to change the partial subgraph explored afterwards.
488 std::shuffle(random_variables.begin(), random_variables.end(), random);
489 for (const int var : random_variables) {
490 if (relaxed_variables.size() < target_size) {
491 visited_variables.push_back(var);
492 if (helper_.IsActive(var)) {
493 relaxed_variables.push_back(var);
494 }
495 } else {
496 break;
497 }
498 }
499 if (relaxed_variables.size() >= target_size) break;
500 }
501
502 return helper_.RelaxGivenVariables(initial_solution, relaxed_variables);
503}
504
506 const CpSolverResponse& initial_solution, double difficulty,
507 absl::BitGenRef random) {
508 const int num_active_vars = helper_.ActiveVariables().size();
509 const int num_model_vars = helper_.ModelProto().variables_size();
510 const int target_size = std::ceil(difficulty * num_active_vars);
511 const int num_constraints = helper_.ConstraintToVar().size();
512 if (num_constraints == 0 || target_size == num_active_vars) {
513 return helper_.FullNeighborhood();
514 }
515 CHECK_GT(target_size, 0);
516
517 std::vector<bool> visited_variables_set(num_model_vars, false);
518 std::vector<int> relaxed_variables;
519 std::vector<bool> added_constraints(num_constraints, false);
520 std::vector<int> next_constraints;
521
522 // Start by a random constraint.
523 next_constraints.push_back(absl::Uniform<int>(random, 0, num_constraints));
524 added_constraints[next_constraints.back()] = true;
525
526 std::vector<int> random_variables;
527 while (relaxed_variables.size() < target_size) {
528 // Stop if we have a full connected component.
529 if (next_constraints.empty()) break;
530
531 // Pick a random unprocessed constraint.
532 const int i = absl::Uniform<int>(random, 0, next_constraints.size());
533 const int constraint_index = next_constraints[i];
534 std::swap(next_constraints[i], next_constraints.back());
535 next_constraints.pop_back();
536
537 // Add all the variable of this constraint and increase the set of next
538 // possible constraints.
539 CHECK_LT(constraint_index, num_constraints);
540 random_variables = helper_.ConstraintToVar()[constraint_index];
541 std::shuffle(random_variables.begin(), random_variables.end(), random);
542 for (const int var : random_variables) {
543 if (visited_variables_set[var]) continue;
544 visited_variables_set[var] = true;
545 if (helper_.IsActive(var)) {
546 relaxed_variables.push_back(var);
547 }
548 if (relaxed_variables.size() == target_size) break;
549
550 for (const int ct : helper_.VarToConstraint()[var]) {
551 if (added_constraints[ct]) continue;
552 added_constraints[ct] = true;
553 next_constraints.push_back(ct);
554 }
555 }
556 }
557 return helper_.RelaxGivenVariables(initial_solution, relaxed_variables);
558}
559
561 const absl::Span<const int> intervals_to_relax,
562 const CpSolverResponse& initial_solution,
563 const NeighborhoodGeneratorHelper& helper) {
564 Neighborhood neighborhood = helper.FullNeighborhood();
565 neighborhood.is_reduced =
566 (intervals_to_relax.size() <
567 helper.TypeToConstraints(ConstraintProto::kInterval).size());
568
569 // We will extend the set with some interval that we cannot fix.
570 std::set<int> ignored_intervals(intervals_to_relax.begin(),
571 intervals_to_relax.end());
572
573 // Fix the presence/absence of non-relaxed intervals.
574 for (const int i : helper.TypeToConstraints(ConstraintProto::kInterval)) {
575 if (ignored_intervals.count(i)) continue;
576
577 const ConstraintProto& interval_ct = neighborhood.cp_model.constraints(i);
578 if (interval_ct.enforcement_literal().empty()) continue;
579
580 CHECK_EQ(interval_ct.enforcement_literal().size(), 1);
581 const int enforcement_ref = interval_ct.enforcement_literal(0);
582 const int enforcement_var = PositiveRef(enforcement_ref);
583 const int value = initial_solution.solution(enforcement_var);
584
585 // If the interval is not enforced, we just relax it. If it belongs to an
586 // exactly one constraint, and the enforced interval is not relaxed, then
587 // propagation will force this interval to stay not enforced. Otherwise,
588 // LNS will be able to change which interval will be enforced among all
589 // alternatives.
590 if (RefIsPositive(enforcement_ref) == (value == 0)) {
591 ignored_intervals.insert(i);
592 continue;
593 }
594
595 // Fix the value.
596 neighborhood.cp_model.mutable_variables(enforcement_var)->clear_domain();
597 neighborhood.cp_model.mutable_variables(enforcement_var)->add_domain(value);
598 neighborhood.cp_model.mutable_variables(enforcement_var)->add_domain(value);
599 }
600
601 for (const int c : helper.TypeToConstraints(ConstraintProto::kNoOverlap)) {
602 // Sort all non-relaxed intervals of this constraint by current start
603 // time.
604 std::vector<std::pair<int64, int>> start_interval_pairs;
605 for (const int i :
606 neighborhood.cp_model.constraints(c).no_overlap().intervals()) {
607 if (ignored_intervals.count(i)) continue;
608 const ConstraintProto& interval_ct = neighborhood.cp_model.constraints(i);
609
610 // TODO(user): we ignore size zero for now.
611 const int size_var = interval_ct.interval().size();
612 if (initial_solution.solution(size_var) == 0) continue;
613
614 const int start_var = interval_ct.interval().start();
615 const int64 start_value = initial_solution.solution(start_var);
616 start_interval_pairs.push_back({start_value, i});
617 }
618 std::sort(start_interval_pairs.begin(), start_interval_pairs.end());
619
620 // Add precedence between the remaining intervals, forcing their order.
621 for (int i = 0; i + 1 < start_interval_pairs.size(); ++i) {
622 const int before_var =
623 neighborhood.cp_model.constraints(start_interval_pairs[i].second)
624 .interval()
625 .end();
626 const int after_var =
627 neighborhood.cp_model.constraints(start_interval_pairs[i + 1].second)
628 .interval()
629 .start();
630 CHECK_LE(initial_solution.solution(before_var),
631 initial_solution.solution(after_var));
632
633 LinearConstraintProto* linear =
634 neighborhood.cp_model.add_constraints()->mutable_linear();
635 linear->add_domain(kint64min);
636 linear->add_domain(0);
637 linear->add_vars(before_var);
638 linear->add_coeffs(1);
639 linear->add_vars(after_var);
640 linear->add_coeffs(-1);
641 }
642 }
643
644 // Set the current solution as a hint.
645 //
646 // TODO(user): Move to common function?
647 neighborhood.cp_model.clear_solution_hint();
648 for (int var = 0; var < neighborhood.cp_model.variables_size(); ++var) {
649 neighborhood.cp_model.mutable_solution_hint()->add_vars(var);
650 neighborhood.cp_model.mutable_solution_hint()->add_values(
651 initial_solution.solution(var));
652 }
653 neighborhood.is_generated = true;
654
655 return neighborhood;
656}
657
659 const CpSolverResponse& initial_solution, double difficulty,
660 absl::BitGenRef random) {
661 std::vector<int> intervals_to_relax =
662 helper_.GetActiveIntervals(initial_solution);
663 GetRandomSubset(difficulty, &intervals_to_relax, random);
664
665 return GenerateSchedulingNeighborhoodForRelaxation(intervals_to_relax,
666 initial_solution, helper_);
667}
668
670 const CpSolverResponse& initial_solution, double difficulty,
671 absl::BitGenRef random) {
672 std::vector<std::pair<int64, int>> start_interval_pairs;
673 for (const int i : helper_.GetActiveIntervals(initial_solution)) {
674 const ConstraintProto& interval_ct = helper_.ModelProto().constraints(i);
675 const int start_var = interval_ct.interval().start();
676 const int64 start_value = initial_solution.solution(start_var);
677 start_interval_pairs.push_back({start_value, i});
678 }
679 std::sort(start_interval_pairs.begin(), start_interval_pairs.end());
680 const int relaxed_size = std::floor(difficulty * start_interval_pairs.size());
681
682 std::uniform_int_distribution<int> random_var(
683 0, start_interval_pairs.size() - relaxed_size - 1);
684 const int random_start_index = random_var(random);
685 std::vector<int> intervals_to_relax;
686 // TODO(user,user): Consider relaxing more than one time window intervals.
687 // This seems to help with Giza models.
688 for (int i = random_start_index; i < relaxed_size; ++i) {
689 intervals_to_relax.push_back(start_interval_pairs[i].second);
690 }
691 return GenerateSchedulingNeighborhoodForRelaxation(intervals_to_relax,
692 initial_solution, helper_);
693}
694
696 if (incomplete_solutions_ != nullptr) {
697 return incomplete_solutions_->HasNewSolution();
698 }
699
700 if (response_manager_ != nullptr) {
701 if (response_manager_->SolutionsRepository().NumSolutions() == 0) {
702 return false;
703 }
704 }
705
706 // At least one relaxation solution should be available to generate a
707 // neighborhood.
708 if (lp_solutions_ != nullptr && lp_solutions_->NumSolutions() > 0) {
709 return true;
710 }
711
712 if (relaxation_solutions_ != nullptr &&
713 relaxation_solutions_->NumSolutions() > 0) {
714 return true;
715 }
716 return false;
717}
718
720 const CpSolverResponse& initial_solution, double difficulty,
721 absl::BitGenRef random) {
722 Neighborhood neighborhood = helper_.FullNeighborhood();
723 neighborhood.is_generated = false;
724
725 const bool lp_solution_available =
726 (lp_solutions_ != nullptr && lp_solutions_->NumSolutions() > 0);
727
728 const bool relaxation_solution_available =
729 (relaxation_solutions_ != nullptr &&
730 relaxation_solutions_->NumSolutions() > 0);
731
732 const bool incomplete_solution_available =
733 (incomplete_solutions_ != nullptr &&
734 incomplete_solutions_->HasNewSolution());
735
736 if (!lp_solution_available && !relaxation_solution_available &&
737 !incomplete_solution_available) {
738 return neighborhood;
739 }
740
741 RINSNeighborhood rins_neighborhood;
742 // Randomly select the type of relaxation if both lp and relaxation solutions
743 // are available.
744 // TODO(user): Tune the probability value for this.
745 std::bernoulli_distribution random_bool(0.5);
746 const bool use_lp_relaxation =
747 (lp_solution_available && relaxation_solution_available)
748 ? random_bool(random)
749 : lp_solution_available;
750 if (use_lp_relaxation) {
751 rins_neighborhood =
752 GetRINSNeighborhood(response_manager_,
753 /*relaxation_solutions=*/nullptr, lp_solutions_,
754 incomplete_solutions_, random);
755 neighborhood.source_info =
756 incomplete_solution_available ? "incomplete" : "lp";
757 } else {
758 CHECK(relaxation_solution_available || incomplete_solution_available);
759 rins_neighborhood = GetRINSNeighborhood(
760 response_manager_, relaxation_solutions_,
761 /*lp_solutions=*/nullptr, incomplete_solutions_, random);
762 neighborhood.source_info =
763 incomplete_solution_available ? "incomplete" : "relaxation";
764 }
765
766 if (rins_neighborhood.fixed_vars.empty() &&
767 rins_neighborhood.reduced_domain_vars.empty()) {
768 return neighborhood;
769 }
770
771 // Fix the variables in the local model.
772 for (const std::pair</*model_var*/ int, /*value*/ int64> fixed_var :
773 rins_neighborhood.fixed_vars) {
774 const int var = fixed_var.first;
775 const int64 value = fixed_var.second;
776 if (var >= neighborhood.cp_model.variables_size()) continue;
777 if (!helper_.IsActive(var)) continue;
778
779 const Domain domain =
780 ReadDomainFromProto(neighborhood.cp_model.variables(var));
781 if (!domain.Contains(value)) {
782 // TODO(user): Instead of aborting, pick the closest point in the domain?
783 return neighborhood;
784 }
785
786 neighborhood.cp_model.mutable_variables(var)->clear_domain();
787 neighborhood.cp_model.mutable_variables(var)->add_domain(value);
788 neighborhood.cp_model.mutable_variables(var)->add_domain(value);
789 neighborhood.is_reduced = true;
790 }
791
792 for (const std::pair</*model_var*/ int, /*domain*/ std::pair<int64, int64>>
793 reduced_var : rins_neighborhood.reduced_domain_vars) {
794 const int var = reduced_var.first;
795 const int64 lb = reduced_var.second.first;
796 const int64 ub = reduced_var.second.second;
797 if (var >= neighborhood.cp_model.variables_size()) continue;
798 if (!helper_.IsActive(var)) continue;
799 Domain domain = ReadDomainFromProto(neighborhood.cp_model.variables(var));
800 domain = domain.IntersectionWith(Domain(lb, ub));
801 if (domain.IsEmpty()) {
802 // TODO(user): Instead of aborting, pick the closest point in the domain?
803 return neighborhood;
804 }
805 FillDomainInProto(domain, neighborhood.cp_model.mutable_variables(var));
806 neighborhood.is_reduced = true;
807 }
808 neighborhood.is_generated = true;
809 return neighborhood;
810}
811
813 const CpSolverResponse& initial_solution, double difficulty,
814 absl::BitGenRef random) {
815 std::vector<int> removable_constraints;
816 const int num_constraints = helper_.ModelProto().constraints_size();
817 removable_constraints.reserve(num_constraints);
818 for (int c = 0; c < num_constraints; ++c) {
819 // Removing intervals is not easy because other constraint might require
820 // them, so for now, we don't remove them.
821 if (helper_.ModelProto().constraints(c).constraint_case() ==
822 ConstraintProto::kInterval) {
823 continue;
824 }
825 removable_constraints.push_back(c);
826 }
827
828 const int target_size =
829 std::round((1.0 - difficulty) * removable_constraints.size());
830
831 const int random_start_index =
832 absl::Uniform<int>(random, 0, removable_constraints.size());
833 std::vector<int> removed_constraints;
834 removed_constraints.reserve(target_size);
835 int c = random_start_index;
836 while (removed_constraints.size() < target_size) {
837 removed_constraints.push_back(removable_constraints[c]);
838 ++c;
839 if (c == removable_constraints.size()) {
840 c = 0;
841 }
842 }
843
844 return helper_.RemoveMarkedConstraints(removed_constraints);
845}
846
849 NeighborhoodGeneratorHelper const* helper, const std::string& name)
850 : NeighborhoodGenerator(name, helper) {
851 std::vector<int> removable_constraints;
852 const int num_constraints = helper_.ModelProto().constraints_size();
853 constraint_weights_.reserve(num_constraints);
854 // TODO(user): Experiment with different starting weights.
855 for (int c = 0; c < num_constraints; ++c) {
856 switch (helper_.ModelProto().constraints(c).constraint_case()) {
857 case ConstraintProto::kCumulative:
858 case ConstraintProto::kAllDiff:
859 case ConstraintProto::kElement:
860 case ConstraintProto::kRoutes:
861 case ConstraintProto::kCircuit:
862 constraint_weights_.push_back(3.0);
863 num_removable_constraints_++;
864 break;
865 case ConstraintProto::kBoolOr:
866 case ConstraintProto::kBoolAnd:
867 case ConstraintProto::kBoolXor:
868 case ConstraintProto::kIntProd:
869 case ConstraintProto::kIntDiv:
870 case ConstraintProto::kIntMod:
871 case ConstraintProto::kIntMax:
872 case ConstraintProto::kLinMax:
873 case ConstraintProto::kIntMin:
874 case ConstraintProto::kLinMin:
875 case ConstraintProto::kNoOverlap:
876 case ConstraintProto::kNoOverlap2D:
877 constraint_weights_.push_back(2.0);
878 num_removable_constraints_++;
879 break;
880 case ConstraintProto::kLinear:
881 case ConstraintProto::kTable:
882 case ConstraintProto::kAutomaton:
883 case ConstraintProto::kInverse:
884 case ConstraintProto::kReservoir:
885 case ConstraintProto::kAtMostOne:
886 case ConstraintProto::kExactlyOne:
887 constraint_weights_.push_back(1.0);
888 num_removable_constraints_++;
889 break;
890 case ConstraintProto::CONSTRAINT_NOT_SET:
891 case ConstraintProto::kInterval:
892 // Removing intervals is not easy because other constraint might require
893 // them, so for now, we don't remove them.
894 constraint_weights_.push_back(0.0);
895 break;
896 }
897 }
898}
899
900void WeightedRandomRelaxationNeighborhoodGenerator::
901 AdditionalProcessingOnSynchronize(const SolveData& solve_data) {
902 const IntegerValue best_objective_improvement =
903 solve_data.new_objective_bound - solve_data.initial_best_objective_bound;
904
905 const std::vector<int>& removed_constraints =
906 removed_constraints_[solve_data.neighborhood_id];
907
908 // Heuristic: We change the weights of the removed constraints if the
909 // neighborhood is solved (status is OPTIMAL or INFEASIBLE) or we observe an
910 // improvement in objective bounds. Otherwise we assume that the
911 // difficulty/time wasn't right for us to record feedbacks.
912 //
913 // If the objective bounds are improved, we bump up the weights. If the
914 // objective bounds are worse and the problem status is OPTIMAL, we bump down
915 // the weights. Otherwise if the new objective bounds are same as current
916 // bounds (which happens a lot on some instances), we do not update the
917 // weights as we do not have a clear signal whether the constraints removed
918 // were good choices or not.
919 // TODO(user): We can improve this heuristic with more experiments.
920 if (best_objective_improvement > 0) {
921 // Bump up the weights of all removed constraints.
922 for (int c : removed_constraints) {
923 if (constraint_weights_[c] <= 90.0) {
924 constraint_weights_[c] += 10.0;
925 } else {
926 constraint_weights_[c] = 100.0;
927 }
928 }
929 } else if (solve_data.status == CpSolverStatus::OPTIMAL &&
930 best_objective_improvement < 0) {
931 // Bump down the weights of all removed constraints.
932 for (int c : removed_constraints) {
933 if (constraint_weights_[c] > 0.5) {
934 constraint_weights_[c] -= 0.5;
935 }
936 }
937 }
938 removed_constraints_.erase(solve_data.neighborhood_id);
939}
940
942 const CpSolverResponse& initial_solution, double difficulty,
943 absl::BitGenRef random) {
944 const int target_size =
945 std::round((1.0 - difficulty) * num_removable_constraints_);
946
947 std::vector<int> removed_constraints;
948
949 // Generate a random number between (0,1) = u[i] and use score[i] =
950 // u[i]^(1/w[i]) and then select top k items with largest scores.
951 // Reference: https://utopia.duth.gr/~pefraimi/research/data/2007EncOfAlg.pdf
952 std::vector<std::pair<double, int>> constraint_removal_scores;
953 std::uniform_real_distribution<double> random_var(0.0, 1.0);
954 for (int c = 0; c < constraint_weights_.size(); ++c) {
955 if (constraint_weights_[c] <= 0) continue;
956 const double u = random_var(random);
957 const double score = std::pow(u, (1 / constraint_weights_[c]));
958 constraint_removal_scores.push_back({score, c});
959 }
960 std::sort(constraint_removal_scores.rbegin(),
961 constraint_removal_scores.rend());
962 for (int i = 0; i < target_size; ++i) {
963 removed_constraints.push_back(constraint_removal_scores[i].second);
964 }
965
966 Neighborhood result = helper_.RemoveMarkedConstraints(removed_constraints);
967 absl::MutexLock mutex_lock(&mutex_);
968 result.id = next_available_id_;
969 next_available_id_++;
970 removed_constraints_.insert({result.id, removed_constraints});
971 return result;
972}
973
974} // namespace sat
975} // namespace operations_research
int64 min
Definition: alldiff_cst.cc:138
int64 max
Definition: alldiff_cst.cc:139
#define CHECK(condition)
Definition: base/logging.h:495
#define CHECK_LT(val1, val2)
Definition: base/logging.h:700
#define CHECK_EQ(val1, val2)
Definition: base/logging.h:697
#define CHECK_GE(val1, val2)
Definition: base/logging.h:701
#define CHECK_GT(val1, val2)
Definition: base/logging.h:702
#define DCHECK_GE(val1, val2)
Definition: base/logging.h:889
#define DCHECK(condition)
Definition: base/logging.h:884
#define CHECK_LE(val1, val2)
Definition: base/logging.h:699
#define VLOG(verboselevel)
Definition: base/logging.h:978
void Update(int num_decreases, int num_increases)
We call domain any subset of Int64 = [kint64min, kint64max].
Domain IntersectionWith(const Domain &domain) const
Returns the intersection of D and domain.
bool IsEmpty() const
Returns true if this is the empty set.
bool Contains(int64 value) const
Returns true iff value is in Domain.
Neighborhood Generate(const CpSolverResponse &initial_solution, double difficulty, absl::BitGenRef random) final
Neighborhood Generate(const CpSolverResponse &initial_solution, double difficulty, absl::BitGenRef random) final
NeighborhoodGeneratorHelper(CpModelProto const *model_proto, SatParameters const *parameters, SharedResponseManager *shared_response, SharedTimeLimit *shared_time_limit=nullptr, SharedBoundsManager *shared_bounds=nullptr)
Definition: cp_model_lns.cc:32
Neighborhood FixAllVariables(const CpSolverResponse &initial_solution) const
Neighborhood FixGivenVariables(const CpSolverResponse &initial_solution, const std::vector< int > &variables_to_fix) const
const absl::Span< const int > TypeToConstraints(ConstraintProto::ConstraintCase type) const
Definition: cp_model_lns.h:118
const std::vector< int > & ActiveVariables() const
Definition: cp_model_lns.h:106
Neighborhood RelaxGivenVariables(const CpSolverResponse &initial_solution, const std::vector< int > &relaxed_variables) const
std::vector< int > GetActiveIntervals(const CpSolverResponse &initial_solution) const
const SharedResponseManager & shared_response() const
Definition: cp_model_lns.h:144
const std::vector< std::vector< int > > & ConstraintToVar() const
Definition: cp_model_lns.h:110
Neighborhood RemoveMarkedConstraints(const std::vector< int > &constraints_to_remove) const
const std::vector< std::vector< int > > & VarToConstraint() const
Definition: cp_model_lns.h:113
double GetUCBScore(int64 total_num_calls) const
virtual void AdditionalProcessingOnSynchronize(const SolveData &solve_data)
Definition: cp_model_lns.h:325
const NeighborhoodGeneratorHelper & helper_
Definition: cp_model_lns.h:328
Neighborhood Generate(const CpSolverResponse &initial_solution, double difficulty, absl::BitGenRef random) final
Neighborhood Generate(const CpSolverResponse &initial_solution, double difficulty, absl::BitGenRef random) final
Neighborhood Generate(const CpSolverResponse &initial_solution, double difficulty, absl::BitGenRef random) final
void GetChangedBounds(int id, std::vector< int > *variables, std::vector< int64 > *new_lower_bounds, std::vector< int64 > *new_upper_bounds)
const SharedSolutionRepository< int64 > & SolutionsRepository() const
Neighborhood Generate(const CpSolverResponse &initial_solution, double difficulty, absl::BitGenRef random) final
Neighborhood Generate(const CpSolverResponse &initial_solution, double difficulty, absl::BitGenRef random) final
Neighborhood Generate(const CpSolverResponse &initial_solution, double difficulty, absl::BitGenRef random) final
WeightedRandomRelaxationNeighborhoodGenerator(NeighborhoodGeneratorHelper const *helper, const std::string &name)
Neighborhood Generate(const CpSolverResponse &initial_solution, double difficulty, absl::BitGenRef random) final
SatParameters parameters
CpModelProto const * model_proto
const std::string name
const Constraint * ct
int64 value
IntVar * var
Definition: expr_array.cc:1858
int64_t int64
static const int64 kint64min
std::vector< int > UsedVariables(const ConstraintProto &ct)
std::vector< int > UsedIntervals(const ConstraintProto &ct)
Domain ReadDomainFromProto(const ProtoWithDomain &proto)
Neighborhood GenerateSchedulingNeighborhoodForRelaxation(const absl::Span< const int > intervals_to_relax, const CpSolverResponse &initial_solution, const NeighborhoodGeneratorHelper &helper)
bool RefIsPositive(int ref)
RINSNeighborhood GetRINSNeighborhood(const SharedResponseManager *response_manager, const SharedRelaxationSolutionRepository *relaxation_solutions, const SharedLPSolutionRepository *lp_solutions, SharedIncompleteSolutionManager *incomplete_solutions, absl::BitGenRef random)
Definition: rins.cc:99
void FillDomainInProto(const Domain &domain, ProtoWithDomain *proto)
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
int64 CapSub(int64 x, int64 y)
IntervalVar * interval
Definition: resource.cc:98
std::vector< std::pair< int, std::pair< int64, int64 > > > reduced_domain_vars
Definition: rins.h:60
std::vector< std::pair< int, int64 > > fixed_vars
Definition: rins.h:58
#define VLOG_IS_ON(verboselevel)
Definition: vlog_is_on.h:41