coin-or / Clp

COIN-OR Linear Programming Solver
Other
396 stars 82 forks source link

CLP produces wrong results #217

Open hongkai-dai opened 2 years ago

hongkai-dai commented 2 years ago

First I would like to thank the authors of CLP for providing this open-source solver.

I ran into a problem when using CLP, and I think it can produce incorrect solution under some situations.

Here is my problem: I want to find the point with the smallest infinity-norm to the point(0.99, 1.99), where the point lies within the polytope with four vertices

(eps, eps)
(1, 1)
(eps, 2)
(1, 2)

where eps is a small number. Since the point (0.99, 1.99) is within this polytope, I would expect that the program finds the infinity norm equal to 0.

I create a simple LP to do this.

minimize z
- z + 5.5511151231257827e-17 * x(0) + x(1) + 5.5511151231257827e-17 * x(2) + x(3))<= 0.99
- z + x(1) + 2 * x(2) + 2 * x(3)) <= 1.99
0.99 <= z + 5.5511151231257827e-17 * x(0) + x(1) + 5.5511151231257827e-17 * x(2) + x(3)) 
1.99 <= z + x(1) + 2 * x(2) + 2 * x(3))
x(0) + x(1) + x(2) + x(3) == 1
0 <= x(0) <= 1
0 <= x(1) <= 1
0 <= x(2) <= 1
0 <= x(3) <= 1

Here x is the weights for each vertex, such that the closest point is written as ∑ᵢ vertex[i] * x[i], and z represents the infinity norm from the closest point to (0.99, 1.99).

Here is my program in c++

#include "ClpSimplex.hpp" 
#include <limits>
#include <iostream>

void solve(double z_lower, double eps) {
  ClpSimplex model;
  model.setLogLevel(0);
  const int numcols = 5;
  const int numrows = 5;
  const double kInf = std::numeric_limits<double>::infinity();
  double collb[5] = {0, 0, 0, 0, z_lower};

  double colub[5] = {1, 1, 1, 1, kInf};
  int start[6] = {0, 3, 8, 13, 18, 22};
  int index[22] = {0, 2, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3};
  double val[22] = {eps, eps, 1, 1, 1, 1, 1, 1, eps, 2, eps, 2, 1, 1, 2, 1, 2, 1, -1, -1, 1, 1};
  double rowlb[5] = {-kInf, -kInf, 0.99, 1.99, 1};
  double rowub[5] = {0.99, 1.99, kInf, kInf, 1};
  double obj[5] = {0, 0, 0, 0, 1};
  model.loadProblem(numcols, numrows, start, index, val, collb, colub, obj, rowlb, rowub);
  model.primal();
  std::cout << "solution: ";
  for (int i = 0; i < numcols; ++i) {
    std::cout << *(model.getColSolution() + i) << " ";
  }
  std::cout << "\n";
}

int main() {
  const double kInf = std::numeric_limits<double>::infinity();
  for (double z_lower : {0., -kInf}) {
    for (double eps : {5.5511151231257827e-17, 0.}) {
      std::cout << "collb[4]=" << z_lower << ", eps=" << eps << "\n";
      solve(z_lower, eps);
    }
  }
}

In this code, the variables are x[0], x[1], x[2], x[3], z.

I tested four settings, based on

  1. Whether to add the constraint z>=0 explicitly (by specifying collb[4]). I think the other constraints would already imply z>=0, so I didn't expect setting collb[4] = 0 would make a difference.
  2. Whether eps is a tiny number, or exactly 0.

Here is my findings:

  1. When eps is small (eps = 5.5511151231257827e-17), and I don't set collb[4]=0, but leave collb[4]=-inf, then CLP finds incorrect result, it thinks the smallest infinity norm is 0.49.
  2. If either eps is 0 or collb[4] = 0, then CLP finds the correct result, that the smallest infinity norm is 0.

Is this due to some numerical error? Why setting collb[4] = 0 would change the result?

jjhforrest commented 2 years ago

A scaling issue - 5.55...e-17 is below standard zero tolerance.

If you add model.scaling(); before primal all 4 runs give same result.

On 22/12/2021 06:55, Hongkai Dai wrote:

First I would like to thank the authors of CLP for providing this open-source solver.

I ran into a problem when using CLP, and I think it can produce incorrect solution under some situations.

Here is my problem: I want to find the point with the smallest infinity-norm to the point(0.99, 1.99), where the point lies within the polytope with four vertices

|(eps, eps) (1, 1) (eps, 2) (1, 2) |

where |eps| is a small number. Since the point (0.99, 1.99) is within this polytope, I would expect that the program finds the infinity norm equal to 0.

I create a simple LP to do this.

|minimize z - z + 5.5511151231257827e-17 x(0) + x(1) + 5.5511151231257827e-17 x(2) + x(3))<= 0.99 - z + x(1) + 2 * x(2) + 2

  • x(3)) <= 1.99 0.99 <= z + 5.5511151231257827e-17 x(0) + x(1) + 5.5511151231257827e-17 x(2) + x(3)) 1.99 <= z + x(1) + 2 * x(2) + 2
  • x(3)) x(0) + x(1) + x(2) + x(3) == 1 0 <= x(0) <= 1 0 <= x(1) <= 1 0 <= x(2) <= 1 0 <= x(3) <= 1 |

Here |x| is the weights for each vertex, such that the closest point is written as ∑ᵢ vertex[i] * x[i], and |z| represents the infinity norm from the closest point to (0.99, 1.99).

Here is my program in c++

include "ClpSimplex.hpp"

include

include

void solve(double z_lower,double eps) { ClpSimplex model; model.setLogLevel(0); const int numcols =5; const int numrows =5; const double kInf = std::numeric_limits::infinity(); double collb[5] = {0,0,0,0, z_lower};

double colub[5] = {1,1,1,1,kInf}; int start[6] = {0,3,8,13,18,22}; int index[22] = {0,2,4,0,1,2,3,4,0,1,2,3,4,0,1,2,3,4,0,1,2,3}; double val[22] = {eps, eps,1,1,1,1,1,1, eps,2, eps,2,1,1,2,1,2,1, -1, -1,1,1}; double rowlb[5] = {-kInf, -kInf,0.99,1.99,1}; double rowub[5] = {0.99,1.99,kInf,kInf,1}; double obj[5] = {0,0,0,0,1}; model.loadProblem(numcols, numrows, start,index, val, collb, colub, obj, rowlb, rowub); model.primal(); std::cout <<"solution: "; for (int i =0; i < numcols; ++i) { std::cout << *(model.getColSolution() + i) <<" "; } std::cout <<"\n"; }

int main() { int option =0; double x4_lb; double eps; const double kInf = std::numeric_limits::infinity(); for (double z_lower : {0., -kInf}) { for (double eps : {5.5511151231257827e-17,0.}) { std::cout <<"collb[4]=" << z_lower <<", eps=" << eps <<"\n"; solve(z_lower, eps); } } }

In this code, the variables are x[0], x[1], x[2], x[3], z.

I tested four settings, based on

  1. Whether to add the constraint z>=0 explicitly (by specifying collb[4]). I think the other constraints would already imply z>=0, so I didn't expect setting collb[4] = 0 would make a difference.
  2. Whether |eps| is a tiny number, or exactly 0.

Here is my findings:

  1. When eps is small (eps = 5.5511151231257827e-17), and I don't set collb[4]=0, but leave collb[4]=-inf, then CLP finds incorrect result, it thinks the smallest infinity norm is 0.49.
  2. If either eps is 0 or collb[4] = 0, then CLP finds the correct result, that the smallest infinity norm is 0.

Is this due to some numerical error? Why setting collb[4] = 0 would change the result?

— Reply to this email directly, view it on GitHub https://github.com/coin-or/Clp/issues/217, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABWJYHHKOHF5OD6SNN7TDK3USFY6TANCNFSM5KR246YQ. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you are subscribed to this thread.Message ID: @.***>

hongkai-dai commented 2 years ago

Thanks @jjhforrest for the answer, calling model.scaling() solves my problem. Much appreciated!

hongkai-dai commented 2 years ago

I have a question about calling model.scaling(0). The documentation says 0 means scaling is off. Hence I would expect that calling model.scaling(0) is the same as not calling model.scaling(0). But when I add model.scaling(0) all 4 cases produce the right results. What should be the expected behavior of setting scaling mode to 0?

hongkai-dai commented 2 years ago

And setting model.scaling(2) produces the wrong result (x = (0, 0.5, 0.5, 0), z = 0.49).