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 
29 namespace operations_research {
30 namespace 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();
49  Synchronize();
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 
105 void 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 
182 bool 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 
305 double 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 
392 namespace {
393 
394 void 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 
900 void 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
const SharedResponseManager & shared_response() const
Definition: cp_model_lns.h:144
const std::vector< std::vector< int > > & VarToConstraint() const
Definition: cp_model_lns.h:113
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
Neighborhood RelaxGivenVariables(const CpSolverResponse &initial_solution, const std::vector< int > &relaxed_variables) const
std::vector< int > GetActiveIntervals(const CpSolverResponse &initial_solution) const
const std::vector< int > & ActiveVariables() const
Definition: cp_model_lns.h:106
Neighborhood RemoveMarkedConstraints(const std::vector< int > &constraints_to_remove) const
const std::vector< std::vector< int > > & ConstraintToVar() const
Definition: cp_model_lns.h:110
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)
bool RefIsPositive(int ref)
std::vector< int > UsedIntervals(const ConstraintProto &ct)
void FillDomainInProto(const Domain &domain, ProtoWithDomain *proto)
Domain ReadDomainFromProto(const ProtoWithDomain &proto)
Neighborhood GenerateSchedulingNeighborhoodForRelaxation(const absl::Span< const int > intervals_to_relax, const CpSolverResponse &initial_solution, const NeighborhoodGeneratorHelper &helper)
RINSNeighborhood GetRINSNeighborhood(const SharedResponseManager *response_manager, const SharedRelaxationSolutionRepository *relaxation_solutions, const SharedLPSolutionRepository *lp_solutions, SharedIncompleteSolutionManager *incomplete_solutions, absl::BitGenRef random)
Definition: rins.cc:99
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