ERGO-Code / HiGHS

Linear optimization software
MIT License
913 stars 170 forks source link

Improve logging when dual ratio test fails #1484

Closed NPC closed 10 months ago

NPC commented 10 months ago

Hello,

I've encountered this issue that I can't find an explanation for, nor can I identify the cause by looking at the log, can you please help?

Here's the LP MPS file, running on Windows 11 from command line (same result called from the DLL via C# interface): LP-problem-2023-11-03.zip

Run via highs.exe --parallel on LP-problem.mps, output:

Running HiGHS 1.6.0: Copyright (c) 2023 HiGHS under MIT licence terms
LP   LP-problem has 18247 rows; 57099 cols; 107808 nonzeros
Presolving model
9111 rows, 30394 cols, 53144 nonzeros
8939 rows, 29003 cols, 50965 nonzeros
Presolve : Reductions: rows 8939(-9308); columns 29003(-28096); elements 50965(-56843)
Solving the presolved LP
Using EKK parallel dual simplex solver - PAMI with concurrency of 8
  Iteration        Objective     Infeasibilities num(sum)
          0    -6.3260410194e+02 Ph1: 1(1); Du: 1(632.604) 0s
Model   status      : Not Set
HiGHS run time      :          0.14

I tried parallel off, presolve off, versions 1.5.3 and 1.4.0, no difference.

As I understand, “Not Set” means that the algorithm ended unexpectedly, before it could report on any errors or infeasibilities, is this correct? Could the values Ph1: 1(1); Du: 1(632.604) in “Infeasibilities” column provide some valuable insight?

Simpler problems continue to solve without issues, as well as complex but different problems, these's something off about this one.

I'll be grateful for any help.

jajhall commented 10 months ago

Your problem is badly scaled. You have very large costs

5055 values satisfy 10^( 11) <= v < 10^( 12)

upper bounds on variables

12 values satisfy 10^( 10) <= v < 10^( 11)

and RHS values for equality constraints

4 values satisfy 10^( 10) <= v < 10^( 11)

The HiGHS dual simplex solver fails due to the excessive costs. I see that you're running HiGHS from the command line, so the logging output that you quote is all that you're told. If calling HiGHS from an application, you'd get an error code returned from Highs::run()

With the option log_dev_level=1, set in a file (say) Options.set, and running

bin/highs --options_file=Options.set LP-problem-2023-11-03.mps

HiGHS gives the following logging information

Cost perturbation for LP-problem-2023-11-03 Initially have 29003 nonzero costs (100%) with min / average / max = 26 / 6.23951e+12 / 2.00573e+14 Large so set max_abs_cost = sqrt(sqrt(max_abs_cost)) = 3763.29 Perturbation column base = 0.00188165 Perturbation row base = 1e-12 Iteration Objective Infeasibilities num(sum) DuPh1 0 -6.3260410194e+02 Ph1: 1(1); Du: 1(632.604) No reason DuPh1 1 0.0000000000e+00 Ph1: 0(0) Possibly optimal DuPh2 1 2.6075283087e+22 Pr: 4928(2.86593e+11) No reason DuPh2 267 2.6236821231e+22 Pr: 4740(6.43333e+10) Synthetic clock DuPh2 526 2.6762934176e+22 Pr: 4580(4.59647e+10) Synthetic clock DualChuzC: No change in loop 2 so return error DualChuzC: workCount = 43; selectTheta=4.60303e+13; remainTheta=4.60303e+13 DualChuzC: workDataNorm = 7.39874; workDualNorm = 3.70191e+15 DuPh2 669 2.7326666760e+22 Pr: 4538(4.17281e+10) Choose column failure DualChuzC: No change in loop 2 so return error DualChuzC: workCount = 43; selectTheta=4.60303e+13; remainTheta=4.60303e+13 DualChuzC: workDataNorm = 7.39874; workDualNorm = 3.70191e+15 dual-phase-2-not-solved HEkkDual::solve return of HighsStatus::Error Simplex iterations: DuPh1 1; DuPh2 668; Total 669 solveLpSimplex return of HighsStatus::Error callSolveLp return of HighsStatus::Error Model status : Not Set HiGHS run time : 0.11

This is after your constraint matrix has been scaled - which is necessary as (originally) there is a very large range of values

Matrix sparsity of dimension 107808 with 107808 nonzeros (100%) in [ 1.77e-09, 1] 73763 values satisfy 10^( 0) <= v < 10^( 1) 33123 values satisfy 10^( -4) <= v < 10^( -3) 21 values satisfy 10^( -6) <= v < 10^( -5) 103 values satisfy 10^( -7) <= v < 10^( -6) 798 values satisfy 10^( -8) <= v < 10^( -7)

Unfortunately, scaling up the small values increases some of the costs used by the solver by a couple of orders of magnitude. This means that, using 64 bit doubles, there's not enough numerical precision for the dual simplex algorithm to run correctly. Hence the failure.

jajhall commented 10 months ago

Essentially what's happening is that you're asking HiGHS to get dual feasibility to within a tolerance of 1e-7 using values of order 1e14, a relative difference that's 1e5 greater than machine precision of 1e-16.

To you I suggest scaling your costs down by 1e6

To me, I suggest a more informative logging message when this error occurs

NPC commented 10 months ago

@jajhall Thank you, this is very helpful! I'll look into cost scaling, this is something we haven't considered (or run into issues with) before. I actually call HiGHS from C#, but it didn't return an error other than kNotset.

jajhall commented 10 months ago

kNotset is a model status, indicating that the model status hasn't been set - because an error occurred.

However, there should be a run status error visible to you.

I'm not a C# user, and the C# interface was written externally. However (in src/interfaces/highs_csharp_api.cs) I see

public HighsStatus run()

{ return (HighsStatus)HighsLpSolver.Highs_run(this.highs); }

which is calling the C API method

HighsInt Highs_run(void highs) { return (HighsInt)((Highs)highs)->run(); }

which is calling the native C++ Highs class method

HighsStatus Highs::run()

So, the enum class member HighsStatus::kError returned from C++ should be cast to an integer and get passed down (eventually!) to C#. Do you see this value?

jajhall commented 10 months ago

Closed by #1485