OR-Tools  8.2
feasibility_pump.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 <vector>
18 
21 #include "ortools/sat/cp_model.pb.h"
22 #include "ortools/sat/integer.h"
23 #include "ortools/sat/sat_base.h"
24 #include "ortools/sat/sat_solver.h"
26 
27 namespace operations_research {
28 namespace sat {
29 
30 using glop::ColIndex;
32 using glop::Fractional;
33 using glop::RowIndex;
34 
35 const double FeasibilityPump::kCpEpsilon = 1e-4;
36 
38  : sat_parameters_(*(model->GetOrCreate<SatParameters>())),
39  time_limit_(model->GetOrCreate<TimeLimit>()),
40  integer_trail_(model->GetOrCreate<IntegerTrail>()),
41  trail_(model->GetOrCreate<Trail>()),
42  integer_encoder_(model->GetOrCreate<IntegerEncoder>()),
43  incomplete_solutions_(model->Mutable<SharedIncompleteSolutionManager>()),
44  sat_solver_(model->GetOrCreate<SatSolver>()),
45  domains_(model->GetOrCreate<IntegerDomains>()),
46  mapping_(model->Get<CpModelMapping>()) {
47  // Tweak the default parameters to make the solve incremental.
48  glop::GlopParameters parameters;
49  // Note(user): Primal simplex does better here since we have a limit on
50  // simplex iterations. So dual simplex sometimes fails to find a LP feasible
51  // solution.
52  parameters.set_use_dual_simplex(false);
53  parameters.set_max_number_of_iterations(2000);
54  simplex_.SetParameters(parameters);
55  lp_data_.Clear();
56  integer_lp_.clear();
57 }
58 
60  VLOG(1) << "Feasibility Pump Total number of simplex iterations: "
61  << total_num_simplex_iterations_;
62 }
63 
65  // We still create the mirror variable right away though.
66  for (const IntegerVariable var : ct.vars) {
67  GetOrCreateMirrorVariable(PositiveVariable(var));
68  }
69 
70  integer_lp_.push_back(LinearConstraintInternal());
71  LinearConstraintInternal& new_ct = integer_lp_.back();
72  new_ct.lb = ct.lb;
73  new_ct.ub = ct.ub;
74  const int size = ct.vars.size();
75  CHECK_LE(ct.lb, ct.ub);
76  for (int i = 0; i < size; ++i) {
77  // We only use positive variable inside this class.
78  IntegerVariable var = ct.vars[i];
79  IntegerValue coeff = ct.coeffs[i];
80  if (!VariableIsPositive(var)) {
81  var = NegationOf(var);
82  coeff = -coeff;
83  }
84  new_ct.terms.push_back({GetOrCreateMirrorVariable(var), coeff});
85  }
86  // Important to keep lp_data_ "clean".
87  std::sort(new_ct.terms.begin(), new_ct.terms.end());
88 }
89 
90 void FeasibilityPump::SetObjectiveCoefficient(IntegerVariable ivar,
91  IntegerValue coeff) {
92  objective_is_defined_ = true;
93  const IntegerVariable pos_var =
94  VariableIsPositive(ivar) ? ivar : NegationOf(ivar);
95  if (ivar != pos_var) coeff = -coeff;
96 
97  const auto it = mirror_lp_variable_.find(pos_var);
98  if (it == mirror_lp_variable_.end()) return;
99  const ColIndex col = it->second;
100  integer_objective_.push_back({col, coeff});
101  objective_infinity_norm_ =
102  std::max(objective_infinity_norm_, IntTypeAbs(coeff));
103 }
104 
105 ColIndex FeasibilityPump::GetOrCreateMirrorVariable(
106  IntegerVariable positive_variable) {
107  DCHECK(VariableIsPositive(positive_variable));
108 
109  const auto it = mirror_lp_variable_.find(positive_variable);
110  if (it == mirror_lp_variable_.end()) {
111  const int model_var =
112  mapping_->GetProtoVariableFromIntegerVariable(positive_variable);
113  model_vars_size_ = std::max(model_vars_size_, model_var + 1);
114 
115  const ColIndex col(integer_variables_.size());
116  mirror_lp_variable_[positive_variable] = col;
117  integer_variables_.push_back(positive_variable);
118  var_is_binary_.push_back(false);
119  lp_solution_.push_back(std::numeric_limits<double>::infinity());
120  integer_solution_.push_back(0);
121 
122  return col;
123  }
124  return it->second;
125 }
126 
127 void FeasibilityPump::PrintStats() {
128  if (lp_solution_is_set_) {
129  VLOG(2) << "Fractionality: " << lp_solution_fractionality_;
130  } else {
131  VLOG(2) << "Fractionality: NA";
132  VLOG(2) << "simplex status: " << simplex_.GetProblemStatus();
133  }
134 
135  if (integer_solution_is_set_) {
136  VLOG(2) << "#Infeasible const: " << num_infeasible_constraints_;
137  VLOG(2) << "Infeasibility: " << integer_solution_infeasibility_;
138  } else {
139  VLOG(2) << "Infeasibility: NA";
140  }
141 }
142 
144  if (lp_data_.num_variables() == 0) {
145  InitializeWorkingLP();
146  }
147  UpdateBoundsOfLpVariables();
148  lp_solution_is_set_ = false;
149  integer_solution_is_set_ = false;
150 
151  // Restore the original objective
152  for (ColIndex col(0); col < lp_data_.num_variables(); ++col) {
153  lp_data_.SetObjectiveCoefficient(col, 0.0);
154  }
155  for (const auto& term : integer_objective_) {
156  lp_data_.SetObjectiveCoefficient(term.first, ToDouble(term.second));
157  }
158 
159  mixing_factor_ = 1.0;
160  for (int i = 0; i < max_fp_iterations_; ++i) {
161  if (time_limit_->LimitReached()) break;
162  L1DistanceMinimize();
163  if (!SolveLp()) break;
164  if (lp_solution_is_integer_) break;
165  if (!Round()) break;
166  // We don't end this loop if the integer solutions is feasible in hope to
167  // get better solution.
168  if (integer_solution_is_feasible_) MaybePushToRepo();
169  }
170 
171  if (model_is_unsat_) return false;
172 
173  PrintStats();
174  MaybePushToRepo();
175  return true;
176 }
177 
178 void FeasibilityPump::MaybePushToRepo() {
179  if (incomplete_solutions_ == nullptr) return;
180 
181  std::vector<double> lp_solution(model_vars_size_,
182  std::numeric_limits<double>::infinity());
183  // TODO(user): Consider adding solutions that have low fractionality.
184  if (lp_solution_is_integer_) {
185  // Fill the solution using LP solution values.
186  for (const IntegerVariable positive_var : integer_variables_) {
187  const int model_var =
188  mapping_->GetProtoVariableFromIntegerVariable(positive_var);
189  if (model_var >= 0 && model_var < model_vars_size_) {
190  lp_solution[model_var] = GetLPSolutionValue(positive_var);
191  }
192  }
193  incomplete_solutions_->AddNewSolution(lp_solution);
194  }
195 
196  if (integer_solution_is_feasible_) {
197  // Fill the solution using Integer solution values.
198  for (const IntegerVariable positive_var : integer_variables_) {
199  const int model_var =
200  mapping_->GetProtoVariableFromIntegerVariable(positive_var);
201  if (model_var >= 0 && model_var < model_vars_size_) {
202  lp_solution[model_var] = GetIntegerSolutionValue(positive_var);
203  }
204  }
205  incomplete_solutions_->AddNewSolution(lp_solution);
206  }
207 }
208 
209 // ----------------------------------------------------------------
210 // -------------------LPSolving------------------------------------
211 // ----------------------------------------------------------------
212 
213 void FeasibilityPump::InitializeWorkingLP() {
214  lp_data_.Clear();
215  // Create variables.
216  for (int i = 0; i < integer_variables_.size(); ++i) {
217  CHECK_EQ(ColIndex(i), lp_data_.CreateNewVariable());
218  lp_data_.SetVariableType(ColIndex(i),
220  }
221 
222  // Add constraints.
223  for (const LinearConstraintInternal& ct : integer_lp_) {
224  const ConstraintIndex row = lp_data_.CreateNewConstraint();
225  lp_data_.SetConstraintBounds(row, ToDouble(ct.lb), ToDouble(ct.ub));
226  for (const auto& term : ct.terms) {
227  lp_data_.SetCoefficient(row, term.first, ToDouble(term.second));
228  }
229  }
230 
231  // Add objective.
232  for (const auto& term : integer_objective_) {
233  lp_data_.SetObjectiveCoefficient(term.first, ToDouble(term.second));
234  }
235 
236  const int num_vars = integer_variables_.size();
237  for (int i = 0; i < num_vars; i++) {
238  const IntegerVariable cp_var = integer_variables_[i];
239  const double lb = ToDouble(integer_trail_->LevelZeroLowerBound(cp_var));
240  const double ub = ToDouble(integer_trail_->LevelZeroUpperBound(cp_var));
241  lp_data_.SetVariableBounds(ColIndex(i), lb, ub);
242  }
243 
244  objective_normalization_factor_ = 0.0;
245  glop::ColIndexVector integer_variables;
246  const ColIndex num_cols = lp_data_.num_variables();
247  for (ColIndex col : lp_data_.IntegerVariablesList()) {
248  var_is_binary_[col.value()] = lp_data_.IsVariableBinary(col);
249  if (!var_is_binary_[col.value()]) {
250  integer_variables.push_back(col);
251  }
252 
253  // The aim of this normalization value is to compute a coefficient of the
254  // d_i variables that should be minimized.
255  objective_normalization_factor_ +=
257  }
258  CHECK_GT(lp_data_.IntegerVariablesList().size(), 0);
259  objective_normalization_factor_ =
260  objective_normalization_factor_ / lp_data_.IntegerVariablesList().size();
261 
262  if (!integer_variables.empty()) {
263  // Update the LpProblem with norm variables and constraints.
264  norm_variables_.assign(num_cols, ColIndex(-1));
265  norm_lhs_constraints_.assign(num_cols, RowIndex(-1));
266  norm_rhs_constraints_.assign(num_cols, RowIndex(-1));
267  // For each integer non-binary variable x_i we introduce one new variable
268  // d_i subject to two new constraints:
269  // d_i - x_i >= -round(x'_i)
270  // d_i + x_i >= +round(x'_i)
271  // That's round(x'_i) - d_i <= x_i <= round(x'_i) + d_i, where d_i is an
272  // unbounded non-negative, and x'_i is the value of variable i from the
273  // previous solution obtained during the projection step. Consequently
274  // coefficients of the constraints are set here, but bounds of the
275  // constraints are updated at each iteration of the feasibility pump. Also
276  // coefficients of the objective are set here: x_i's are not present in the
277  // objective (i.e., coefficients set to 0.0), and d_i's are present in the
278  // objective with coefficients set to 1.0.
279  // Note that the treatment of integer non-binary variables is different
280  // from the treatment of binary variables. Binary variables do not impose
281  // any extra variables, nor extra constraints, but their objective
282  // coefficients are changed in the linear projection steps.
283  for (const ColIndex col : integer_variables) {
284  const ColIndex norm_variable = lp_data_.CreateNewVariable();
285  norm_variables_[col] = norm_variable;
286  lp_data_.SetVariableBounds(norm_variable, 0.0, glop::kInfinity);
287  const RowIndex row_a = lp_data_.CreateNewConstraint();
288  norm_lhs_constraints_[col] = row_a;
289  lp_data_.SetCoefficient(row_a, norm_variable, 1.0);
290  lp_data_.SetCoefficient(row_a, col, -1.0);
291  const RowIndex row_b = lp_data_.CreateNewConstraint();
292  norm_rhs_constraints_[col] = row_b;
293  lp_data_.SetCoefficient(row_b, norm_variable, 1.0);
294  lp_data_.SetCoefficient(row_b, col, 1.0);
295  }
296  }
297 
298  scaler_.Scale(&lp_data_);
300  /*detect_integer_constraints=*/false);
301 }
302 
303 void FeasibilityPump::L1DistanceMinimize() {
304  std::vector<double> new_obj_coeffs(lp_data_.num_variables().value(), 0.0);
305 
306  // Set the original subobjective. The coefficients are scaled by mixing factor
307  // and the offset remains at 0 (because it does not affect the solution).
308  const ColIndex num_cols(lp_data_.objective_coefficients().size());
309  for (ColIndex col(0); col < num_cols; ++col) {
310  new_obj_coeffs[col.value()] =
311  mixing_factor_ * lp_data_.objective_coefficients()[col];
312  }
313 
314  // Set the norm subobjective. The coefficients are scaled by 1 - mixing factor
315  // and the offset remains at 0 (because it does not affect the solution).
316  for (const ColIndex col : lp_data_.IntegerVariablesList()) {
317  if (var_is_binary_[col.value()]) {
318  const Fractional objective_coefficient =
319  mixing_factor_ * lp_data_.objective_coefficients()[col] +
320  (1 - mixing_factor_) * objective_normalization_factor_ *
321  (1 - 2 * integer_solution_[col.value()]);
322  new_obj_coeffs[col.value()] = objective_coefficient;
323  } else { // The variable is integer.
324  // Update the bounds of the constraints added in
325  // InitializeIntegerVariables() (see there for more details):
326  // d_i - x_i >= -round(x'_i)
327  // d_i + x_i >= +round(x'_i)
328 
329  // TODO(user): We change both the objective and the bounds, thus
330  // breaking the incrementality. Handle integer variables differently,
331  // e.g., intensify rounding, or use soft fixing from: Fischetti, Lodi,
332  // "Local Branching", Math Program Ser B 98:23-47 (2003).
333  const Fractional objective_coefficient =
334  (1 - mixing_factor_) * objective_normalization_factor_;
335  new_obj_coeffs[norm_variables_[col].value()] = objective_coefficient;
336  // At this point, constraint bounds have already been transformed into
337  // bounds of slack variables. Instead of updating the constraints, we need
338  // to update the slack variables corresponding to them.
339  const ColIndex norm_lhs_slack_variable =
340  lp_data_.GetSlackVariable(norm_lhs_constraints_[col]);
341  const double lhs_scaling_factor =
342  scaler_.VariableScalingFactor(norm_lhs_slack_variable);
343  lp_data_.SetVariableBounds(
344  norm_lhs_slack_variable, -glop::kInfinity,
345  lhs_scaling_factor * integer_solution_[col.value()]);
346  const ColIndex norm_rhs_slack_variable =
347  lp_data_.GetSlackVariable(norm_rhs_constraints_[col]);
348  const double rhs_scaling_factor =
349  scaler_.VariableScalingFactor(norm_rhs_slack_variable);
350  lp_data_.SetVariableBounds(
351  norm_rhs_slack_variable, -glop::kInfinity,
352  -rhs_scaling_factor * integer_solution_[col.value()]);
353  }
354  }
355  for (ColIndex col(0); col < lp_data_.num_variables(); ++col) {
356  lp_data_.SetObjectiveCoefficient(col, new_obj_coeffs[col.value()]);
357  }
358  // TODO(user): Tune this or expose as parameter.
359  mixing_factor_ *= 0.8;
360 }
361 
362 bool FeasibilityPump::SolveLp() {
363  const int num_vars = integer_variables_.size();
364  VLOG(3) << "LP relaxation: " << lp_data_.GetDimensionString() << ".";
365 
366  const auto status = simplex_.Solve(lp_data_, time_limit_);
367  total_num_simplex_iterations_ += simplex_.GetNumberOfIterations();
368  if (!status.ok()) {
369  VLOG(1) << "The LP solver encountered an error: " << status.error_message();
370  simplex_.ClearStateForNextSolve();
371  return false;
372  }
373 
374  // TODO(user): This shouldn't really happen except if the problem is UNSAT.
375  // But we can't just rely on a potentially imprecise LP to close the problem.
376  // The rest of the solver should do that with exact precision.
377  VLOG(3) << "simplex status: " << simplex_.GetProblemStatus();
379  return false;
380  }
381 
382  lp_solution_fractionality_ = 0.0;
387  lp_solution_is_set_ = true;
388  for (int i = 0; i < num_vars; i++) {
389  const double value = GetVariableValueAtCpScale(ColIndex(i));
390  lp_solution_[i] = value;
391  lp_solution_fractionality_ = std::max(
392  lp_solution_fractionality_, std::abs(value - std::round(value)));
393  }
394 
395  // Compute the objective value.
396  lp_objective_ = 0;
397  for (const auto& term : integer_objective_) {
398  lp_objective_ += lp_solution_[term.first.value()] * term.second.value();
399  }
400  lp_solution_is_integer_ = lp_solution_fractionality_ < kCpEpsilon;
401  }
402  return true;
403 }
404 
405 void FeasibilityPump::UpdateBoundsOfLpVariables() {
406  const int num_vars = integer_variables_.size();
407  for (int i = 0; i < num_vars; i++) {
408  const IntegerVariable cp_var = integer_variables_[i];
409  const double lb = ToDouble(integer_trail_->LevelZeroLowerBound(cp_var));
410  const double ub = ToDouble(integer_trail_->LevelZeroUpperBound(cp_var));
411  const double factor = scaler_.VariableScalingFactor(ColIndex(i));
412  lp_data_.SetVariableBounds(ColIndex(i), lb * factor, ub * factor);
413  }
414 }
415 
416 double FeasibilityPump::GetLPSolutionValue(IntegerVariable variable) const {
417  return lp_solution_[gtl::FindOrDie(mirror_lp_variable_, variable).value()];
418 }
419 
420 double FeasibilityPump::GetVariableValueAtCpScale(ColIndex var) {
421  return scaler_.UnscaleVariableValue(var, simplex_.GetVariableValue(var));
422 }
423 
424 // ----------------------------------------------------------------
425 // -------------------Rounding-------------------------------------
426 // ----------------------------------------------------------------
427 
428 int64 FeasibilityPump::GetIntegerSolutionValue(IntegerVariable variable) const {
429  return integer_solution_[gtl::FindOrDie(mirror_lp_variable_, variable)
430  .value()];
431 }
432 
433 bool FeasibilityPump::Round() {
434  bool rounding_successful = true;
435  if (sat_parameters_.fp_rounding() == SatParameters::NEAREST_INTEGER) {
436  rounding_successful = NearestIntegerRounding();
437  } else if (sat_parameters_.fp_rounding() == SatParameters::LOCK_BASED) {
438  rounding_successful = LockBasedRounding();
439  } else if (sat_parameters_.fp_rounding() ==
440  SatParameters::ACTIVE_LOCK_BASED) {
441  rounding_successful = ActiveLockBasedRounding();
442  } else if (sat_parameters_.fp_rounding() ==
443  SatParameters::PROPAGATION_ASSISTED) {
444  rounding_successful = PropagationRounding();
445  }
446  if (!rounding_successful) return false;
447  FillIntegerSolutionStats();
448  return true;
449 }
450 
451 bool FeasibilityPump::NearestIntegerRounding() {
452  if (!lp_solution_is_set_) return false;
453  for (int i = 0; i < lp_solution_.size(); ++i) {
454  integer_solution_[i] = static_cast<int64>(std::round(lp_solution_[i]));
455  }
456  integer_solution_is_set_ = true;
457  return true;
458 }
459 
460 bool FeasibilityPump::LockBasedRounding() {
461  if (!lp_solution_is_set_) return false;
462  const int num_vars = integer_variables_.size();
463 
464  // We compute the number of locks based on variable coefficient in constraints
465  // and constraint bounds. This doesn't change over time so we cache it.
466  if (var_up_locks_.empty()) {
467  var_up_locks_.resize(num_vars, 0);
468  var_down_locks_.resize(num_vars, 0);
469  for (int i = 0; i < num_vars; ++i) {
470  for (const auto entry : lp_data_.GetSparseColumn(ColIndex(i))) {
471  ColIndex slack = lp_data_.GetSlackVariable(entry.row());
472  const bool constraint_upper_bounded =
473  lp_data_.variable_lower_bounds()[slack] > -glop::kInfinity;
474 
475  const bool constraint_lower_bounded =
476  lp_data_.variable_upper_bounds()[slack] < glop::kInfinity;
477 
478  if (entry.coefficient() > 0) {
479  var_up_locks_[i] += constraint_upper_bounded;
480  var_down_locks_[i] += constraint_lower_bounded;
481  } else {
482  var_up_locks_[i] += constraint_lower_bounded;
483  var_down_locks_[i] += constraint_upper_bounded;
484  }
485  }
486  }
487  }
488 
489  for (int i = 0; i < lp_solution_.size(); ++i) {
490  if (std::abs(lp_solution_[i] - std::round(lp_solution_[i])) < 0.1 ||
491  var_up_locks_[i] == var_down_locks_[i]) {
492  integer_solution_[i] = static_cast<int64>(std::round(lp_solution_[i]));
493  } else if (var_up_locks_[i] > var_down_locks_[i]) {
494  integer_solution_[i] = static_cast<int64>(std::floor(lp_solution_[i]));
495  } else {
496  integer_solution_[i] = static_cast<int64>(std::ceil(lp_solution_[i]));
497  }
498  }
499  integer_solution_is_set_ = true;
500  return true;
501 }
502 
503 bool FeasibilityPump::ActiveLockBasedRounding() {
504  if (!lp_solution_is_set_) return false;
505  const int num_vars = integer_variables_.size();
506 
507  // We compute the number of locks based on variable coefficient in constraints
508  // and constraint bounds of active constraints. We consider the bound of the
509  // constraint that is tight for the current lp solution.
510  for (int i = 0; i < num_vars; ++i) {
511  if (std::abs(lp_solution_[i] - std::round(lp_solution_[i])) < 0.1) {
512  integer_solution_[i] = static_cast<int64>(std::round(lp_solution_[i]));
513  }
514 
515  int up_locks = 0;
516  int down_locks = 0;
517  for (const auto entry : lp_data_.GetSparseColumn(ColIndex(i))) {
518  const ConstraintStatus row_status =
519  simplex_.GetConstraintStatus(entry.row());
520  if (row_status == ConstraintStatus::AT_LOWER_BOUND) {
521  if (entry.coefficient() > 0) {
522  down_locks++;
523  } else {
524  up_locks++;
525  }
526  } else if (row_status == ConstraintStatus::AT_UPPER_BOUND) {
527  if (entry.coefficient() > 0) {
528  up_locks++;
529  } else {
530  down_locks++;
531  }
532  }
533  }
534  if (up_locks == down_locks) {
535  integer_solution_[i] = static_cast<int64>(std::round(lp_solution_[i]));
536  } else if (up_locks > down_locks) {
537  integer_solution_[i] = static_cast<int64>(std::floor(lp_solution_[i]));
538  } else {
539  integer_solution_[i] = static_cast<int64>(std::ceil(lp_solution_[i]));
540  }
541  }
542 
543  integer_solution_is_set_ = true;
544  return true;
545 }
546 
547 bool FeasibilityPump::PropagationRounding() {
548  if (!lp_solution_is_set_) return false;
549  sat_solver_->ResetToLevelZero();
550 
551  // Compute an order in which we will fix variables and do the propagation.
552  std::vector<int> rounding_order;
553  {
554  std::vector<std::pair<double, int>> binary_fractionality_vars;
555  std::vector<std::pair<double, int>> general_fractionality_vars;
556  for (int i = 0; i < lp_solution_.size(); ++i) {
557  const double fractionality =
558  std::abs(std::round(lp_solution_[i]) - lp_solution_[i]);
559  if (var_is_binary_[i]) {
560  binary_fractionality_vars.push_back({fractionality, i});
561  } else {
562  general_fractionality_vars.push_back({fractionality, i});
563  }
564  }
565  std::sort(binary_fractionality_vars.begin(),
566  binary_fractionality_vars.end());
567  std::sort(general_fractionality_vars.begin(),
568  general_fractionality_vars.end());
569 
570  for (int i = 0; i < binary_fractionality_vars.size(); ++i) {
571  rounding_order.push_back(binary_fractionality_vars[i].second);
572  }
573  for (int i = 0; i < general_fractionality_vars.size(); ++i) {
574  rounding_order.push_back(general_fractionality_vars[i].second);
575  }
576  }
577 
578  for (const int var_index : rounding_order) {
579  if (time_limit_->LimitReached()) return false;
580  // Get the bounds of the variable.
581  const IntegerVariable var = integer_variables_[var_index];
582  const Domain& domain = (*domains_)[var];
583 
584  const IntegerValue lb = integer_trail_->LowerBound(var);
585  const IntegerValue ub = integer_trail_->UpperBound(var);
586  if (lb == ub) {
587  integer_solution_[var_index] = lb.value();
588  continue;
589  }
590 
591  const int64 rounded_value =
592  static_cast<int64>(std::round(lp_solution_[var_index]));
593  const int64 floor_value =
594  static_cast<int64>(std::floor(lp_solution_[var_index]));
595  const int64 ceil_value =
596  static_cast<int64>(std::ceil(lp_solution_[var_index]));
597 
598  const bool floor_is_in_domain =
599  (domain.Contains(floor_value) && lb.value() <= floor_value);
600  const bool ceil_is_in_domain =
601  (domain.Contains(ceil_value) && ub.value() >= ceil_value);
602  if (domain.IsEmpty()) {
603  integer_solution_[var_index] = rounded_value;
604  model_is_unsat_ = true;
605  return false;
606  }
607 
608  if (ceil_value < lb.value()) {
609  integer_solution_[var_index] = lb.value();
610  } else if (floor_value > ub.value()) {
611  integer_solution_[var_index] = ub.value();
612  } else if (ceil_is_in_domain && floor_is_in_domain) {
613  DCHECK(domain.Contains(rounded_value));
614  integer_solution_[var_index] = rounded_value;
615  } else if (ceil_is_in_domain) {
616  integer_solution_[var_index] = ceil_value;
617  } else if (floor_is_in_domain) {
618  integer_solution_[var_index] = floor_value;
619  } else {
620  const std::pair<IntegerLiteral, IntegerLiteral> values_in_domain =
621  integer_encoder_->Canonicalize(
622  IntegerLiteral::GreaterOrEqual(var, IntegerValue(rounded_value)));
623  const int64 lower_value = values_in_domain.first.bound.value();
624  const int64 higher_value = -values_in_domain.second.bound.value();
625  const int64 distance_from_lower_value =
626  std::abs(lower_value - rounded_value);
627  const int64 distance_from_higher_value =
628  std::abs(higher_value - rounded_value);
629 
630  integer_solution_[var_index] =
631  (distance_from_lower_value < distance_from_higher_value)
632  ? lower_value
633  : higher_value;
634  }
635 
636  CHECK(domain.Contains(integer_solution_[var_index]));
637  CHECK_GE(integer_solution_[var_index], lb);
638  CHECK_LE(integer_solution_[var_index], ub);
639 
640  // Propagate the value.
641  //
642  // When we want to fix the variable at its lb or ub, we do not create an
643  // equality literal to minimize the number of new literal we create. This
644  // is because creating an "== value" literal will implicitly also create
645  // a ">= value" and a "<= value" literals.
646  Literal to_enqueue;
647  const IntegerValue value(integer_solution_[var_index]);
648  if (value == lb) {
649  to_enqueue = integer_encoder_->GetOrCreateAssociatedLiteral(
651  } else if (value == ub) {
652  to_enqueue = integer_encoder_->GetOrCreateAssociatedLiteral(
654  } else {
655  to_enqueue =
657  }
658 
659  if (!sat_solver_->FinishPropagation()) {
660  model_is_unsat_ = true;
661  return false;
662  }
663  sat_solver_->EnqueueDecisionAndBacktrackOnConflict(to_enqueue);
664 
665  if (sat_solver_->IsModelUnsat()) {
666  model_is_unsat_ = true;
667  return false;
668  }
669  }
670  sat_solver_->ResetToLevelZero();
671  integer_solution_is_set_ = true;
672  return true;
673 }
674 
675 void FeasibilityPump::FillIntegerSolutionStats() {
676  // Compute the objective value.
677  integer_solution_objective_ = 0;
678  for (const auto& term : integer_objective_) {
679  integer_solution_objective_ +=
680  integer_solution_[term.first.value()] * term.second.value();
681  }
682 
683  integer_solution_is_feasible_ = true;
684  num_infeasible_constraints_ = 0;
685  integer_solution_infeasibility_ = 0;
686  for (RowIndex i(0); i < integer_lp_.size(); ++i) {
687  int64 activity = 0;
688  for (const auto& term : integer_lp_[i].terms) {
689  const int64 prod =
690  CapProd(integer_solution_[term.first.value()], term.second.value());
691  if (prod <= kint64min || prod >= kint64max) {
692  activity = prod;
693  break;
694  }
695  activity = CapAdd(activity, prod);
696  if (activity <= kint64min || activity >= kint64max) break;
697  }
698  if (activity > integer_lp_[i].ub || activity < integer_lp_[i].lb) {
699  integer_solution_is_feasible_ = false;
700  num_infeasible_constraints_++;
701  const int64 ub_infeasibility = activity > integer_lp_[i].ub.value()
702  ? activity - integer_lp_[i].ub.value()
703  : 0;
704  const int64 lb_infeasibility = activity < integer_lp_[i].lb.value()
705  ? integer_lp_[i].lb.value() - activity
706  : 0;
707  integer_solution_infeasibility_ =
708  std::max(integer_solution_infeasibility_,
709  std::max(ub_infeasibility, lb_infeasibility));
710  }
711  }
712 }
713 
714 } // namespace sat
715 } // namespace operations_research
int64 max
Definition: alldiff_cst.cc:139
#define CHECK(condition)
Definition: base/logging.h:495
#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(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 push_back(const value_type &x)
A simple class to enforce both an elapsed time limit and a deterministic time limit in the same threa...
Definition: time_limit.h:105
bool LimitReached()
Returns true when the external limit is true, or the deterministic time is over the deterministic lim...
Definition: time_limit.h:532
void SetVariableBounds(ColIndex col, Fractional lower_bound, Fractional upper_bound)
Definition: lp_data.cc:248
void SetCoefficient(RowIndex row, ColIndex col, Fractional value)
Definition: lp_data.cc:316
ColIndex GetSlackVariable(RowIndex row) const
Definition: lp_data.cc:753
const DenseRow & variable_lower_bounds() const
Definition: lp_data.h:228
const DenseRow & objective_coefficients() const
Definition: lp_data.h:222
const std::vector< ColIndex > & IntegerVariablesList() const
Definition: lp_data.cc:279
Fractional GetObjectiveCoefficientForMinimizationVersion(ColIndex col) const
Definition: lp_data.cc:418
void SetConstraintBounds(RowIndex row, Fractional lower_bound, Fractional upper_bound)
Definition: lp_data.cc:308
void SetVariableType(ColIndex col, VariableType type)
Definition: lp_data.cc:235
void AddSlackVariablesWhereNecessary(bool detect_integer_constraints)
Definition: lp_data.cc:695
void SetObjectiveCoefficient(ColIndex col, Fractional value)
Definition: lp_data.cc:325
bool IsVariableBinary(ColIndex col) const
Definition: lp_data.cc:299
const DenseRow & variable_upper_bounds() const
Definition: lp_data.h:231
std::string GetDimensionString() const
Definition: lp_data.cc:424
const SparseColumn & GetSparseColumn(ColIndex col) const
Definition: lp_data.cc:408
Fractional VariableScalingFactor(ColIndex col) const
Fractional UnscaleVariableValue(ColIndex col, Fractional value) const
Fractional GetVariableValue(ColIndex col) const
ABSL_MUST_USE_RESULT Status Solve(const LinearProgram &lp, TimeLimit *time_limit)
ConstraintStatus GetConstraintStatus(RowIndex row) const
void SetParameters(const GlopParameters &parameters)
void assign(IntType size, const T &v)
Definition: lp_types.h:274
int GetProtoVariableFromIntegerVariable(IntegerVariable var) const
double GetLPSolutionValue(IntegerVariable variable) const
int64 GetIntegerSolutionValue(IntegerVariable variable) const
void AddLinearConstraint(const LinearConstraint &ct)
void SetObjectiveCoefficient(IntegerVariable ivar, IntegerValue coeff)
Literal GetOrCreateLiteralAssociatedToEquality(IntegerVariable var, IntegerValue value)
Definition: integer.cc:248
std::pair< IntegerLiteral, IntegerLiteral > Canonicalize(IntegerLiteral i_lit) const
Definition: integer.cc:184
Literal GetOrCreateAssociatedLiteral(IntegerLiteral i_lit)
Definition: integer.cc:202
IntegerValue UpperBound(IntegerVariable i) const
Definition: integer.h:1304
IntegerValue LevelZeroUpperBound(IntegerVariable var) const
Definition: integer.h:1355
IntegerValue LevelZeroLowerBound(IntegerVariable var) const
Definition: integer.h:1350
IntegerValue LowerBound(IntegerVariable i) const
Definition: integer.h:1300
Class that owns everything related to a particular optimization model.
Definition: sat/model.h:38
int EnqueueDecisionAndBacktrackOnConflict(Literal true_literal)
Definition: sat_solver.cc:861
void AddNewSolution(const std::vector< double > &lp_solution)
SatParameters parameters
const Constraint * ct
int64 value
IntVar * var
Definition: expr_array.cc:1858
GRBmodel * model
static const int64 kint64max
int64_t int64
ColIndex col
Definition: markowitz.cc:176
RowIndex row
Definition: markowitz.cc:175
const Collection::value_type::second_type & FindOrDie(const Collection &collection, const typename Collection::value_type::first_type &key)
Definition: map_util.h:176
std::vector< ColIndex > ColIndexVector
Definition: lp_types.h:308
const double kInfinity
Definition: lp_types.h:83
IntType IntTypeAbs(IntType t)
Definition: integer.h:77
IntegerVariable PositiveVariable(IntegerVariable i)
Definition: integer.h:134
std::vector< IntegerVariable > NegationOf(const std::vector< IntegerVariable > &vars)
Definition: integer.cc:27
bool VariableIsPositive(IntegerVariable i)
Definition: integer.h:130
double ToDouble(IntegerValue value)
Definition: integer.h:69
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
int64 CapAdd(int64 x, int64 y)
int64 CapProd(int64 x, int64 y)
static IntegerLiteral LowerOrEqual(IntegerVariable i, IntegerValue bound)
Definition: integer.h:1270
static IntegerLiteral GreaterOrEqual(IntegerVariable i, IntegerValue bound)
Definition: integer.h:1264