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"
26
27namespace operations_research {
28namespace sat {
29
30using glop::ColIndex;
33using glop::RowIndex;
34
35const 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)) {
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
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
105ColIndex 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
127void 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
178void 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
213void 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
303void 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
362bool 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
405void 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
416double FeasibilityPump::GetLPSolutionValue(IntegerVariable variable) const {
417 return lp_solution_[gtl::FindOrDie(mirror_lp_variable_, variable).value()];
418}
419
420double FeasibilityPump::GetVariableValueAtCpScale(ColIndex var) {
421 return scaler_.UnscaleVariableValue(var, simplex_.GetVariableValue(var));
422}
423
424// ----------------------------------------------------------------
425// -------------------Rounding-------------------------------------
426// ----------------------------------------------------------------
427
428int64 FeasibilityPump::GetIntegerSolutionValue(IntegerVariable variable) const {
429 return integer_solution_[gtl::FindOrDie(mirror_lp_variable_, variable)
430 .value()];
431}
432
433bool 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
451bool 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
460bool 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
503bool 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
547bool 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
675void 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 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
const DenseRow & variable_upper_bounds() const
Definition: lp_data.h:231
void SetVariableType(ColIndex col, VariableType type)
Definition: lp_data.cc:235
const DenseRow & objective_coefficients() const
Definition: lp_data.h:222
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_lower_bounds() const
Definition: lp_data.h:228
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
bool VariableIsPositive(IntegerVariable i)
Definition: integer.h:130
IntegerVariable PositiveVariable(IntegerVariable i)
Definition: integer.h:134
IntType IntTypeAbs(IntType t)
Definition: integer.h:77
double ToDouble(IntegerValue value)
Definition: integer.h:69
std::vector< IntegerVariable > NegationOf(const std::vector< IntegerVariable > &vars)
Definition: integer.cc:27
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