21 #include "ortools/sat/cp_model.pb.h"
35 const double FeasibilityPump::kCpEpsilon = 1e-4;
38 : sat_parameters_(*(
model->GetOrCreate<SatParameters>())),
60 VLOG(1) <<
"Feasibility Pump Total number of simplex iterations: "
61 << total_num_simplex_iterations_;
66 for (
const IntegerVariable
var :
ct.vars) {
70 integer_lp_.
push_back(LinearConstraintInternal());
71 LinearConstraintInternal& new_ct = integer_lp_.
back();
74 const int size =
ct.vars.size();
76 for (
int i = 0; i < size; ++i) {
78 IntegerVariable
var =
ct.vars[i];
79 IntegerValue coeff =
ct.coeffs[i];
84 new_ct.terms.push_back({GetOrCreateMirrorVariable(
var), coeff});
87 std::sort(new_ct.terms.begin(), new_ct.terms.end());
92 objective_is_defined_ =
true;
93 const IntegerVariable pos_var =
95 if (ivar != pos_var) coeff = -coeff;
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_ =
105 ColIndex FeasibilityPump::GetOrCreateMirrorVariable(
106 IntegerVariable positive_variable) {
109 const auto it = mirror_lp_variable_.find(positive_variable);
110 if (it == mirror_lp_variable_.end()) {
111 const int model_var =
113 model_vars_size_ =
std::max(model_vars_size_, model_var + 1);
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);
127 void FeasibilityPump::PrintStats() {
128 if (lp_solution_is_set_) {
129 VLOG(2) <<
"Fractionality: " << lp_solution_fractionality_;
131 VLOG(2) <<
"Fractionality: NA";
135 if (integer_solution_is_set_) {
136 VLOG(2) <<
"#Infeasible const: " << num_infeasible_constraints_;
137 VLOG(2) <<
"Infeasibility: " << integer_solution_infeasibility_;
139 VLOG(2) <<
"Infeasibility: NA";
145 InitializeWorkingLP();
147 UpdateBoundsOfLpVariables();
148 lp_solution_is_set_ =
false;
149 integer_solution_is_set_ =
false;
155 for (
const auto& term : integer_objective_) {
159 mixing_factor_ = 1.0;
160 for (
int i = 0; i < max_fp_iterations_; ++i) {
162 L1DistanceMinimize();
163 if (!SolveLp())
break;
164 if (lp_solution_is_integer_)
break;
168 if (integer_solution_is_feasible_) MaybePushToRepo();
171 if (model_is_unsat_)
return false;
178 void FeasibilityPump::MaybePushToRepo() {
179 if (incomplete_solutions_ ==
nullptr)
return;
181 std::vector<double> lp_solution(model_vars_size_,
182 std::numeric_limits<double>::infinity());
184 if (lp_solution_is_integer_) {
186 for (
const IntegerVariable positive_var : integer_variables_) {
187 const int model_var =
189 if (model_var >= 0 && model_var < model_vars_size_) {
196 if (integer_solution_is_feasible_) {
198 for (
const IntegerVariable positive_var : integer_variables_) {
199 const int model_var =
201 if (model_var >= 0 && model_var < model_vars_size_) {
213 void FeasibilityPump::InitializeWorkingLP() {
216 for (
int i = 0; i < integer_variables_.size(); ++i) {
223 for (
const LinearConstraintInternal&
ct : integer_lp_) {
226 for (
const auto& term :
ct.terms) {
232 for (
const auto& term : integer_objective_) {
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];
244 objective_normalization_factor_ = 0.0;
249 if (!var_is_binary_[
col.value()]) {
250 integer_variables.push_back(
col);
255 objective_normalization_factor_ +=
259 objective_normalization_factor_ =
262 if (!integer_variables.empty()) {
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));
283 for (
const ColIndex
col : integer_variables) {
285 norm_variables_[
col] = norm_variable;
288 norm_lhs_constraints_[
col] = row_a;
292 norm_rhs_constraints_[
col] = row_b;
298 scaler_.
Scale(&lp_data_);
303 void FeasibilityPump::L1DistanceMinimize() {
304 std::vector<double> new_obj_coeffs(lp_data_.
num_variables().value(), 0.0);
309 for (ColIndex
col(0);
col < num_cols; ++
col) {
310 new_obj_coeffs[
col.value()] =
317 if (var_is_binary_[
col.value()]) {
320 (1 - mixing_factor_) * objective_normalization_factor_ *
321 (1 - 2 * integer_solution_[
col.value()]);
322 new_obj_coeffs[
col.value()] = objective_coefficient;
334 (1 - mixing_factor_) * objective_normalization_factor_;
335 new_obj_coeffs[norm_variables_[
col].value()] = objective_coefficient;
339 const ColIndex norm_lhs_slack_variable =
341 const double lhs_scaling_factor =
345 lhs_scaling_factor * integer_solution_[
col.value()]);
346 const ColIndex norm_rhs_slack_variable =
348 const double rhs_scaling_factor =
352 -rhs_scaling_factor * integer_solution_[
col.value()]);
359 mixing_factor_ *= 0.8;
362 bool FeasibilityPump::SolveLp() {
363 const int num_vars = integer_variables_.size();
366 const auto status = simplex_.
Solve(lp_data_, time_limit_);
369 VLOG(1) <<
"The LP solver encountered an error: " << status.error_message();
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)));
397 for (
const auto& term : integer_objective_) {
398 lp_objective_ += lp_solution_[term.first.value()] * term.second.value();
400 lp_solution_is_integer_ = lp_solution_fractionality_ < kCpEpsilon;
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];
417 return lp_solution_[
gtl::FindOrDie(mirror_lp_variable_, variable).value()];
420 double FeasibilityPump::GetVariableValueAtCpScale(ColIndex
var) {
429 return integer_solution_[
gtl::FindOrDie(mirror_lp_variable_, variable)
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();
446 if (!rounding_successful)
return false;
447 FillIntegerSolutionStats();
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]));
456 integer_solution_is_set_ =
true;
460 bool FeasibilityPump::LockBasedRounding() {
461 if (!lp_solution_is_set_)
return false;
462 const int num_vars = integer_variables_.size();
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) {
472 const bool constraint_upper_bounded =
475 const bool constraint_lower_bounded =
478 if (entry.coefficient() > 0) {
479 var_up_locks_[i] += constraint_upper_bounded;
480 var_down_locks_[i] += constraint_lower_bounded;
482 var_up_locks_[i] += constraint_lower_bounded;
483 var_down_locks_[i] += constraint_upper_bounded;
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]));
496 integer_solution_[i] =
static_cast<int64>(std::ceil(lp_solution_[i]));
499 integer_solution_is_set_ =
true;
503 bool FeasibilityPump::ActiveLockBasedRounding() {
504 if (!lp_solution_is_set_)
return false;
505 const int num_vars = integer_variables_.size();
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]));
520 if (row_status == ConstraintStatus::AT_LOWER_BOUND) {
521 if (entry.coefficient() > 0) {
526 }
else if (row_status == ConstraintStatus::AT_UPPER_BOUND) {
527 if (entry.coefficient() > 0) {
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]));
539 integer_solution_[i] =
static_cast<int64>(std::ceil(lp_solution_[i]));
543 integer_solution_is_set_ =
true;
547 bool FeasibilityPump::PropagationRounding() {
548 if (!lp_solution_is_set_)
return false;
552 std::vector<int> rounding_order;
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});
562 general_fractionality_vars.push_back({fractionality, i});
565 std::sort(binary_fractionality_vars.begin(),
566 binary_fractionality_vars.end());
567 std::sort(general_fractionality_vars.begin(),
568 general_fractionality_vars.end());
570 for (
int i = 0; i < binary_fractionality_vars.size(); ++i) {
571 rounding_order.push_back(binary_fractionality_vars[i].second);
573 for (
int i = 0; i < general_fractionality_vars.size(); ++i) {
574 rounding_order.push_back(general_fractionality_vars[i].second);
578 for (
const int var_index : rounding_order) {
581 const IntegerVariable
var = integer_variables_[var_index];
582 const Domain& domain = (*domains_)[
var];
587 integer_solution_[var_index] = lb.value();
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]));
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;
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;
620 const std::pair<IntegerLiteral, IntegerLiteral> values_in_domain =
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);
630 integer_solution_[var_index] =
631 (distance_from_lower_value < distance_from_higher_value)
636 CHECK(domain.Contains(integer_solution_[var_index]));
637 CHECK_GE(integer_solution_[var_index], lb);
638 CHECK_LE(integer_solution_[var_index], ub);
647 const IntegerValue
value(integer_solution_[var_index]);
651 }
else if (
value == ub) {
660 model_is_unsat_ =
true;
666 model_is_unsat_ =
true;
671 integer_solution_is_set_ =
true;
675 void FeasibilityPump::FillIntegerSolutionStats() {
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();
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) {
688 for (
const auto& term : integer_lp_[i].terms) {
690 CapProd(integer_solution_[term.first.value()], term.second.value());
691 if (prod <= kint64min || prod >=
kint64max) {
695 activity =
CapAdd(activity, prod);
696 if (activity <= kint64min || activity >=
kint64max)
break;
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()
704 const int64 lb_infeasibility = activity < integer_lp_[i].lb.value()
705 ? integer_lp_[i].lb.value() - activity
707 integer_solution_infeasibility_ =
708 std::max(integer_solution_infeasibility_,
709 std::max(ub_infeasibility, lb_infeasibility));
#define CHECK_EQ(val1, val2)
#define CHECK_GE(val1, val2)
#define CHECK_GT(val1, val2)
#define DCHECK(condition)
#define CHECK_LE(val1, val2)
#define VLOG(verboselevel)
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...
bool LimitReached()
Returns true when the external limit is true, or the deterministic time is over the deterministic lim...
void SetVariableBounds(ColIndex col, Fractional lower_bound, Fractional upper_bound)
void SetCoefficient(RowIndex row, ColIndex col, Fractional value)
ColIndex GetSlackVariable(RowIndex row) const
const DenseRow & variable_lower_bounds() const
const DenseRow & objective_coefficients() const
const std::vector< ColIndex > & IntegerVariablesList() const
Fractional GetObjectiveCoefficientForMinimizationVersion(ColIndex col) const
void SetConstraintBounds(RowIndex row, Fractional lower_bound, Fractional upper_bound)
ColIndex CreateNewVariable()
void SetVariableType(ColIndex col, VariableType type)
void AddSlackVariablesWhereNecessary(bool detect_integer_constraints)
void SetObjectiveCoefficient(ColIndex col, Fractional value)
bool IsVariableBinary(ColIndex col) const
const DenseRow & variable_upper_bounds() const
ColIndex num_variables() const
RowIndex CreateNewConstraint()
std::string GetDimensionString() const
const SparseColumn & GetSparseColumn(ColIndex col) const
void Scale(LinearProgram *lp)
Fractional VariableScalingFactor(ColIndex col) const
Fractional UnscaleVariableValue(ColIndex col, Fractional value) const
Fractional GetVariableValue(ColIndex col) const
int64 GetNumberOfIterations() const
ABSL_MUST_USE_RESULT Status Solve(const LinearProgram &lp, TimeLimit *time_limit)
ProblemStatus GetProblemStatus() const
ConstraintStatus GetConstraintStatus(RowIndex row) const
void ClearStateForNextSolve()
void SetParameters(const GlopParameters ¶meters)
void assign(IntType size, const T &v)
int GetProtoVariableFromIntegerVariable(IntegerVariable var) const
FeasibilityPump(Model *model)
glop::RowIndex ConstraintIndex
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)
std::pair< IntegerLiteral, IntegerLiteral > Canonicalize(IntegerLiteral i_lit) const
Literal GetOrCreateAssociatedLiteral(IntegerLiteral i_lit)
IntegerValue UpperBound(IntegerVariable i) const
IntegerValue LevelZeroUpperBound(IntegerVariable var) const
IntegerValue LevelZeroLowerBound(IntegerVariable var) const
IntegerValue LowerBound(IntegerVariable i) const
Class that owns everything related to a particular optimization model.
int EnqueueDecisionAndBacktrackOnConflict(Literal true_literal)
bool IsModelUnsat() const
void AddNewSolution(const std::vector< double > &lp_solution)
static const int64 kint64max
const Collection::value_type::second_type & FindOrDie(const Collection &collection, const typename Collection::value_type::first_type &key)
std::vector< ColIndex > ColIndexVector
IntType IntTypeAbs(IntType t)
IntegerVariable PositiveVariable(IntegerVariable i)
std::vector< IntegerVariable > NegationOf(const std::vector< IntegerVariable > &vars)
bool VariableIsPositive(IntegerVariable i)
double ToDouble(IntegerValue value)
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)
static IntegerLiteral GreaterOrEqual(IntegerVariable i, IntegerValue bound)