Skip to content

Commit 1571759

Browse files
jajhalljajhall
authored
Merge pull request #2709 from ERGO-Code/fix-2643
Avoid primal simplex when unscaled problem has primal infeasibilities but is dual feasible: use dual simplex to clean up
2 parents a568d00 + 6d7e3a7 commit 1571759

File tree

10 files changed

+252
-19
lines changed

10 files changed

+252
-19
lines changed

‎check/TestMipSolver.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -519,7 +519,7 @@ TEST_CASE("MIP-infeasible-start", "[highs_test_mip_solver]") {
519519

520520
// Stefan's example
521521
std::string filename;
522-
filename = std::string(HIGHS_DIR) + "/check/instances/infeasible.mps";
522+
filename = std::string(HIGHS_DIR) + "/check/instances/infeasible-mip1.mps";
523523

524524
highs.readModel(filename);
525525
sol.col_value = {75, 0, 275, 300, 300, 0, 0, 0, 50, 0, 0,

‎check/instances/infeasible-mip0.mps‎

Lines changed: 140 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,140 @@
1+
NAME issue-2643 11-row 9-col 3-conts 6-integer infeasible MIP
2+
ROWS
3+
N OBJ
4+
L c1
5+
L c8
6+
G c11
7+
G c15
8+
G c16
9+
L c18
10+
L c19
11+
G c20
12+
L c21
13+
L c30
14+
L c_b
15+
COLUMNS
16+
MARK 'MARKER' 'INTORG'
17+
a c1 -3.285000000000e+01
18+
a c8 1.643000000000e+01
19+
a c11 -9.010000000000e+01
20+
a c15 1.044000000000e+01
21+
a c16 5.839000000000e+01
22+
a c18 4.890000000000e+00
23+
a c19 -8.793000000000e+01
24+
a c30 4.052000000000e+01
25+
a c_b 3.425000000000e+02
26+
MARK 'MARKER' 'INTEND'
27+
b c1 -7.836000000000e+01
28+
b c8 -5.884000000000e+01
29+
b c11 6.600000000000e+00
30+
b c15 3.756000000000e+01
31+
b c16 -7.458000000000e+01
32+
b c18 -2.943000000000e+01
33+
b c19 3.706000000000e+01
34+
b c21 3.627000000000e+01
35+
b c30 8.378000000000e+01
36+
b c_b -3.432200000000e+02
37+
b OBJ -2.388000000000e+01
38+
MARK 'MARKER' 'INTORG'
39+
c c1 8.447000000000e+01
40+
c c8 7.200000000000e-01
41+
c c11 -5.095000000000e+01
42+
c c15 2.361000000000e+01
43+
c c16 -4.415000000000e+01
44+
c c18 4.096000000000e+01
45+
c c19 6.282000000000e+01
46+
c c21 9.855000000000e+01
47+
c c30 -5.153000000000e+01
48+
c c_b -4.321800000000e+02
49+
c OBJ 3.222700000000e+02
50+
MARK 'MARKER' 'INTEND'
51+
MARK 'MARKER' 'INTORG'
52+
d c1 -6.273000000000e+01
53+
d c8 -3.610000000000e+00
54+
d c11 7.613000000000e+01
55+
d c15 4.711000000000e+01
56+
d c16 6.418000000000e+01
57+
d c18 8.969000000000e+01
58+
d c19 7.961000000000e+01
59+
d c30 -6.288000000000e+01
60+
d c_b 1.897600000000e+02
61+
d OBJ 1.437100000000e+02
62+
MARK 'MARKER' 'INTEND'
63+
e c1 3.142000000000e+01
64+
e c8 1.298000000000e+01
65+
e c11 -4.320000000000e+01
66+
e c15 -7.787000000000e+01
67+
e c16 -5.097000000000e+01
68+
e c18 5.004000000000e+01
69+
e c19 -2.433000000000e+01
70+
e c30 1.527000000000e+01
71+
e c_b -4.748400000000e+02
72+
e OBJ 1.394300000000e+02
73+
f c1 -1.334000000000e+01
74+
f c8 7.861000000000e+01
75+
f c11 5.318000000000e+01
76+
f c15 3.386000000000e+01
77+
f c16 -8.273000000000e+01
78+
f c18 7.153000000000e+01
79+
f c19 5.374000000000e+01
80+
f c30 -7.809000000000e+01
81+
f c_b 2.149000000000e+01
82+
MARK 'MARKER' 'INTORG'
83+
g c1 -9.900000000000e+00
84+
g c11 3.000000000000e-02
85+
g c15 -4.859000000000e+01
86+
g c16 4.934000000000e+01
87+
g c18 -6.935000000000e+01
88+
g c19 1.416000000000e+01
89+
g c30 7.693000000000e+01
90+
g c_b -4.488500000000e+02
91+
g OBJ 3.465700000000e+02
92+
MARK 'MARKER' 'INTEND'
93+
MARK 'MARKER' 'INTORG'
94+
h c1 -3.243000000000e+01
95+
h c8 -6.963000000000e+01
96+
h c11 8.690000000000e+01
97+
h c15 1.032000000000e+01
98+
h c16 -1.539000000000e+01
99+
h c18 2.638000000000e+01
100+
h c19 9.907000000000e+01
101+
h c20 -6.612000000000e+01
102+
h c21 9.098000000000e+01
103+
h c30 3.970000000000e+01
104+
h c_b 1.005500000000e+02
105+
MARK 'MARKER' 'INTEND'
106+
MARK 'MARKER' 'INTORG'
107+
i c1 -3.877000000000e+01
108+
i c8 5.683000000000e+01
109+
i c11 -7.974000000000e+01
110+
i c15 -1.480000000000e+01
111+
i c16 -6.928000000000e+01
112+
i c18 5.720000000000e+01
113+
i c19 4.524000000000e+01
114+
i c30 -6.992000000000e+01
115+
i c_b -6.559000000000e+01
116+
i OBJ -1.916800000000e+02
117+
MARK 'MARKER' 'INTEND'
118+
RHS
119+
RHS c1 -1.709860000000e+03
120+
RHS c8 -7.222500000000e+02
121+
RHS c11 1.751420000000e+03
122+
RHS c15 -1.926560000000e+03
123+
RHS c16 -1.379760000000e+03
124+
RHS c18 1.821870000000e+03
125+
RHS c19 -1.226870000000e+03
126+
RHS c30 2.301000000000e+01
127+
RHS c_b -6.806710000000e+04
128+
BOUNDS
129+
LO BND a -1.000000000000e+02
130+
UP BND b 1.000000000000e+02
131+
LO BND c -1.000000000000e+02
132+
LO BND d -1.000000000000e+02
133+
UP BND d 1.000000000000e+02
134+
LO BND f -1.000000000000e+02
135+
LO BND g -1.000000000000e+02
136+
UP BND g 1.000000000000e+02
137+
LO BND h -1.000000000000e+02
138+
LO BND i -1.000000000000e+02
139+
UP BND i 1.000000000000e+02
140+
ENDATA

check/instances/infeasible.mps renamed to check/instances/infeasible-mip1.mps

File renamed without changes.

‎highs/lp_data/HighsSolution.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -78,6 +78,7 @@ void getLpKktFailures(const HighsOptions& options, const HighsLp& lp,
7878
HighsPrimalDualErrors& primal_dual_errors,
7979
const bool get_residuals = false);
8080

81+
// Inner getKktFailures
8182
void getKktFailures(const HighsOptions& options, const bool is_qp,
8283
const HighsLp& lp, const std::vector<double>& gradient,
8384
const HighsSolution& solution, HighsInfo& highs_info,

‎highs/simplex/HApp.h‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -440,6 +440,19 @@ inline HighsStatus solveLpSimplex(HighsLpSolverObject& solver_object) {
440440
(int)ekk_instance.iteration_count_);
441441
}
442442
} else {
443+
// There are unscaled primal infeasibilities, so use dual
444+
// simplex
445+
assert(num_unscaled_primal_infeasibilities > 0);
446+
if (options.simplex_strategy != kSimplexStrategyDual)
447+
highsLogDev(
448+
options.log_options, HighsLogType::kInfo,
449+
"Forcing change from %s to %s\n",
450+
ekk_instance.simplexStrategyToString(options.simplex_strategy)
451+
.c_str(),
452+
ekk_instance.simplexStrategyToString(kSimplexStrategyDual)
453+
.c_str());
454+
455+
options.simplex_strategy = kSimplexStrategyDual;
443456
// Using dual simplex, so force Devex if starting from an advanced
444457
// basis with no steepest edge weights
445458
if ((status.has_basis || basis.valid) &&

‎highs/simplex/HEkk.cpp‎

Lines changed: 42 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -192,6 +192,7 @@ void HEkk::clearEkkDataInfo() {
192192
info.backtracking_basis_.clear();
193193
info.backtracking_basis_costs_shifted_ = false;
194194
info.backtracking_basis_costs_perturbed_ = false;
195+
info.backtracking_basis_bounds_shifted_ = false;
195196
info.backtracking_basis_bounds_perturbed_ = false;
196197
info.backtracking_basis_workShift_.clear();
197198
info.backtracking_basis_workLowerShift_.clear();
@@ -221,6 +222,7 @@ void HEkk::clearEkkDataInfo() {
221222
info.allow_bound_perturbation = true;
222223
info.costs_shifted = false;
223224
info.costs_perturbed = false;
225+
info.bounds_shifted = false;
224226
info.bounds_perturbed = false;
225227

226228
info.num_primal_infeasibilities = kHighsIllegalInfeasibilityCount;
@@ -1052,8 +1054,8 @@ HighsStatus HEkk::solve(const bool force_phase2) {
10521054
algorithm_name = "primal";
10531055
reportSimplexPhaseIterations(options_->log_options, iteration_count_, info_,
10541056
true);
1055-
highsLogUser(options_->log_options, HighsLogType::kInfo,
1056-
"Using EKK primal simplex solver\n");
1057+
highsLogUser(options_->log_options, HighsLogType::kInfo, "Using %s\n",
1058+
simplexStrategyToString(kSimplexStrategyPrimal).c_str());
10571059
HEkkPrimal primal_solver(*this);
10581060
call_status = primal_solver.solve(force_phase2);
10591061
assert(called_return_from_solve_);
@@ -1066,17 +1068,17 @@ HighsStatus HEkk::solve(const bool force_phase2) {
10661068
// Solve, depending on the particular strategy
10671069
if (simplex_strategy == kSimplexStrategyDualTasks) {
10681070
highsLogUser(options_->log_options, HighsLogType::kInfo,
1069-
"Using EKK parallel dual simplex solver - SIP with "
1070-
"concurrency of %d\n",
1071+
"Using %s with concurrency of %d\n",
1072+
simplexStrategyToString(kSimplexStrategyDualTasks).c_str(),
10711073
int(info_.num_concurrency));
10721074
} else if (simplex_strategy == kSimplexStrategyDualMulti) {
10731075
highsLogUser(options_->log_options, HighsLogType::kInfo,
1074-
"Using EKK parallel dual simplex solver - PAMI with "
1075-
"concurrency of %d\n",
1076+
"Using %s with concurrency of %d\n",
1077+
simplexStrategyToString(kSimplexStrategyDualMulti).c_str(),
10761078
int(info_.num_concurrency));
10771079
} else {
1078-
highsLogUser(options_->log_options, HighsLogType::kInfo,
1079-
"Using EKK dual simplex solver - serial\n");
1080+
highsLogUser(options_->log_options, HighsLogType::kInfo, "Using %s\n",
1081+
simplexStrategyToString(kSimplexStrategyDual).c_str());
10801082
}
10811083
HEkkDual dual_solver(*this);
10821084
call_status = dual_solver.solve(force_phase2);
@@ -1101,7 +1103,7 @@ HighsStatus HEkk::solve(const bool force_phase2) {
11011103
if (return_status == HighsStatus::kError)
11021104
return returnFromEkkSolve(return_status);
11031105
highsLogDev(options_->log_options, HighsLogType::kInfo,
1104-
"EKK %s simplex solver returns %" HIGHSINT_FORMAT
1106+
"%s simplex solver returns %" HIGHSINT_FORMAT
11051107
" primal and %" HIGHSINT_FORMAT
11061108
" dual infeasibilities: "
11071109
"Status %s\n",
@@ -1877,6 +1879,8 @@ bool HEkk::getBacktrackingBasis() {
18771879
basis_ = info_.backtracking_basis_;
18781880
info_.costs_shifted = (info_.backtracking_basis_costs_shifted_ != 0);
18791881
info_.costs_perturbed = (info_.backtracking_basis_costs_perturbed_ != 0);
1882+
info_.bounds_shifted = (info_.backtracking_basis_bounds_shifted_ != 0);
1883+
info_.bounds_perturbed = (info_.backtracking_basis_bounds_perturbed_ != 0);
18801884
info_.workShift_ = info_.backtracking_basis_workShift_;
18811885
const HighsInt num_tot = lp_.num_col_ + lp_.num_row_;
18821886
for (HighsInt iVar = 0; iVar < num_tot; iVar++)
@@ -1901,6 +1905,7 @@ void HEkk::putBacktrackingBasis(
19011905
info_.backtracking_basis_.basicIndex_ = basicIndex_before_compute_factor;
19021906
info_.backtracking_basis_costs_shifted_ = info_.costs_shifted;
19031907
info_.backtracking_basis_costs_perturbed_ = info_.costs_perturbed;
1908+
info_.backtracking_basis_bounds_shifted_ = info_.bounds_shifted;
19041909
info_.backtracking_basis_bounds_perturbed_ = info_.bounds_perturbed;
19051910
info_.backtracking_basis_workShift_ = info_.workShift_;
19061911
const HighsInt num_tot = lp_.num_col_ + lp_.num_row_;
@@ -2291,6 +2296,24 @@ bool HEkk::lpFactorRowCompatible(const HighsInt expectedNumRow) const {
22912296
return consistent_num_row;
22922297
}
22932298

2299+
std::string HEkk::simplexStrategyToString(
2300+
const HighsInt simplex_strategy) const {
2301+
assert(kSimplexStrategyMin <= simplex_strategy &&
2302+
simplex_strategy <= kSimplexStrategyMax);
2303+
if (simplex_strategy == kSimplexStrategyChoose)
2304+
return "choose simplex solver";
2305+
if (simplex_strategy == kSimplexStrategyDual) return "dual simplex solver";
2306+
if (simplex_strategy == kSimplexStrategyDualPlain)
2307+
return "serial dual simplex solver";
2308+
if (simplex_strategy == kSimplexStrategyDualTasks)
2309+
return "parallel dual simplex solver - SIP";
2310+
if (simplex_strategy == kSimplexStrategyDualMulti)
2311+
return "parallel dual simplex solver - PAMI";
2312+
if (simplex_strategy == kSimplexStrategyPrimal)
2313+
return "primal simplex solver";
2314+
return "Unknown";
2315+
}
2316+
22942317
void HEkk::zeroBasicDuals() {
22952318
for (HighsInt iRow = 0; iRow < lp_.num_row_; iRow++)
22962319
info_.workDual_[basis_.basicIndex_[iRow]] = 0;
@@ -2544,6 +2567,7 @@ void HEkk::initialiseBound(const SimplexAlgorithm algorithm,
25442567
const HighsInt solve_phase, const bool perturb) {
25452568
initialiseLpColBound();
25462569
initialiseLpRowBound();
2570+
info_.bounds_shifted = false;
25472571
info_.bounds_perturbed = false;
25482572
// Primal simplex bounds are either from the LP or perturbed
25492573
if (algorithm == SimplexAlgorithm::kPrimal) {
@@ -3523,7 +3547,8 @@ HighsStatus HEkk::returnFromSolve(const HighsStatus return_status) {
35233547
case HighsModelStatus::kInfeasible: {
35243548
// Primal infeasibility has been identified in primal phase 1,
35253549
// or proved in dual phase 2. There should be no primal
3526-
// perturbations
3550+
// perturbations or shifts
3551+
assert(!info_.bounds_shifted);
35273552
assert(!info_.bounds_perturbed);
35283553
if (exit_algorithm_ == SimplexAlgorithm::kPrimal) {
35293554
// Reset the simplex costs and recompute duals after primal
@@ -3540,6 +3565,7 @@ HighsStatus HEkk::returnFromSolve(const HighsStatus return_status) {
35403565
// Dual simplex has identified dual infeasibility in phase
35413566
// 1. There should be no dual perturbations
35423567
assert(exit_algorithm_ == SimplexAlgorithm::kDual);
3568+
assert(!info_.costs_shifted);
35433569
assert(!info_.costs_perturbed);
35443570
// Reset the simplex bounds and recompute primals
35453571
initialiseBound(SimplexAlgorithm::kDual, kSolvePhase2);
@@ -3553,7 +3579,10 @@ HighsStatus HEkk::returnFromSolve(const HighsStatus return_status) {
35533579
// Primal simplex has identified unboundedness in phase 2. There
35543580
// should be no primal or dual perturbations
35553581
assert(exit_algorithm_ == SimplexAlgorithm::kPrimal);
3556-
assert(!info_.costs_perturbed && !info_.bounds_perturbed);
3582+
assert(!info_.costs_shifted);
3583+
assert(!info_.costs_perturbed);
3584+
assert(!info_.bounds_shifted);
3585+
assert(!info_.bounds_perturbed);
35573586
computeSimplexInfeasible();
35583587
// Primal solution should be feasible
35593588
assert(info_.num_primal_infeasibilities == 0);
@@ -3583,7 +3612,7 @@ HighsStatus HEkk::returnFromSolve(const HighsStatus return_status) {
35833612
default: {
35843613
highsLogDev(
35853614
options_->log_options, HighsLogType::kError,
3586-
"EKK %s simplex solver returns status %s\n",
3615+
"%s simplex solver returns status %s\n",
35873616
exit_algorithm_ == SimplexAlgorithm::kPrimal ? "primal" : "dual",
35883617
utilModelStatusToString(model_status_).c_str());
35893618
return HighsStatus::kError;
@@ -3603,6 +3632,7 @@ HighsStatus HEkk::returnFromSolve(const HighsStatus return_status) {
36033632
} else {
36043633
return_dual_solution_status_ = kSolutionStatusInfeasible;
36053634
}
3635+
assert(debugNoShiftsOrPerturbations());
36063636
computePrimalObjectiveValue();
36073637
if (!options_->log_dev_level) {
36083638
const bool force = true;

‎highs/simplex/HEkk.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -139,6 +139,7 @@ class HEkk {
139139
void initialisePartitionedRowwiseMatrix();
140140
bool lpFactorRowCompatible() const;
141141
bool lpFactorRowCompatible(const HighsInt expectedNumRow) const;
142+
std::string simplexStrategyToString(const HighsInt simplex_strategy) const;
142143

143144
void zeroBasicDuals();
144145

@@ -408,7 +409,7 @@ class HEkk {
408409
HighsDebugStatus debugSimplexDualInfeasible(const std::string message,
409410
const bool force_report = false);
410411
HighsDebugStatus debugComputeDual(const bool initialise = false) const;
411-
412+
bool debugNoShiftsOrPerturbations() const;
412413
friend class HEkkPrimal;
413414
friend class HEkkDual;
414415
friend class HEkkDualRow;

‎highs/simplex/HEkkDebug.cpp‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1693,3 +1693,30 @@ HighsDebugStatus HEkk::debugSimplexDualInfeasible(const std::string message,
16931693
info.sum_dual_infeasibilities);
16941694
return HighsDebugStatus::kOk;
16951695
}
1696+
1697+
bool HEkk::debugNoShiftsOrPerturbations() const {
1698+
switch (model_status_) {
1699+
case HighsModelStatus::kOptimal: {
1700+
if (this->info_.costs_shifted || this->info_.costs_perturbed ||
1701+
this->info_.bounds_shifted || this->info_.bounds_perturbed) {
1702+
return false;
1703+
}
1704+
break;
1705+
}
1706+
case HighsModelStatus::kInfeasible: {
1707+
if (this->info_.bounds_shifted || this->info_.bounds_perturbed) {
1708+
return false;
1709+
}
1710+
break;
1711+
}
1712+
case HighsModelStatus::kUnboundedOrInfeasible: {
1713+
if (this->info_.costs_shifted || this->info_.costs_perturbed) {
1714+
return false;
1715+
}
1716+
break;
1717+
}
1718+
default:
1719+
return true;
1720+
}
1721+
return true;
1722+
}

0 commit comments

Comments
 (0)