coin-or / Clp

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

Trouble getting updated dual solution #167

Closed spoorendonk closed 3 years ago

spoorendonk commented 3 years ago

I have a simple column generation problem where I need to add the columns twice (with resolve() twice) to get the correct dual variables. Am I missing a setting somewhere?

I use OsiClpSolverInterface, but get the same result with ClpSimplex directly.

My examples is:

In iteration 0, to obtain feasibility, I start with a phase 1 and have only artificial variables like

\Problem name: ClpDefaultName

Minimize
obj: x0 + x1
Subject To
cons0:  x0 = 1
cons1:  x1 = 1
cons2:  x0 + x1 <= 2
cons2_low: x0 + x1 >= 1
Bounds
 0 <= x0 <= 1
 0 <= x1 <= 1
End

with duals [1, 1, 0] after an initialSolve().

Now I price in variables similar to x0 and x1 but with cost 0 to price out the artificial ones, so iteration 1 is

\Problem name: ClpDefaultName

Minimize
obj: x0 + x1
Subject To
cons0:  x0 + x2 = 1
cons1:  x1 = 1
cons2:  x0 + x1 + x2 <= 2
cons2_low: x0 + x1 + x2 >= 1
Bounds
 0 <= x0 <= 1
 0 <= x1 <= 1
 0 <= x2 <= 1
End

then after resolve() i get duals [1, 1, 0] again. The problem is proven optimal and the objevtive value is 1, so the dual solution is off?

If a I price in the same variable again, ie.,

\Problem name: ClpDefaultName

Minimize
obj: x0 + x1
Subject To
cons0:  x0 + x2 + x3 = 1
cons1:  x1 = 1
cons2:  x0 + x1 + x2 + x3 <= 2
cons2_low: x0 + x1 + x2 + x3 >= 1
Bounds
 0 <= x0 <= 1
 0 <= x1 <= 1
 0 <= x2 <= 1
 0 <= x3 <= 1
End

and resolve() I get the correct duals [0, 1, 0]. The same problem occurs when pricing in the next variable (and pricing out x1).

My code is

  OsiClpSolverInterface clpSolver;

  int numCols = 2;
  int numRows = 3;
  CoinBigIndex start[3] = {0, 2, 4};
  int index[4] = {0, 2, 1, 2};
  double value[4] = {1, 1, 1, 1};
  double collb[2] = {0, 0};
  double colub[2] = {1, 1};
  double obj[2] = {1, 1};
  double rowlb[3] = {1, 1, 1};
  double rowub[3] = {1, 1, 2};
  clpSolver.loadProblem(numCols, numRows, start, index, value, collb, colub, obj, rowlb, rowub);
  auto filename = "LP_" + std::to_string(0);
  clpSolver.writeLp(filename.c_str());
  clpSolver.initialSolve();

  // dump duals
  std::cout << "iteration 0 - Duals: ";
  for(auto i=0; i< clpSolver.getNumRows(); i++){
    std::cout <<  clpSolver.getRowPrice()[i] << " ";
  }
  std::cout << std::endl;

  int rows[2] = {0, 2};
  double elements[2] = {1, 1};
  clpSolver.addCol(2, rows, elements, 0, 1, 0);
  filename = "LP_" + std::to_string(1);
    clpSolver.writeLp(filename.c_str());
  clpSolver.resolve();

  // dump duals
  std::cout << "iteration 1 - Duals: ";
  for(auto i=0; i< clpSolver.getNumRows(); i++){
    std::cout <<  clpSolver.getRowPrice()[i] << " ";
  }
  std::cout << std::endl;

  rows[0] = 0;
  rows[1] = 2;
  elements[0] = 1;
  elements[1] = 1;
  clpSolver.addCol(2, rows, elements, 0, 1, 0);
  filename = "LP_" + std::to_string(2);
    clpSolver.writeLp(filename.c_str());
  clpSolver.resolve();

  // dump duals
  std::cout << "iteration 2  - Duals: ";
  for(auto i=0; i< clpSolver.getNumRows(); i++){
    std::cout <<  clpSolver.getRowPrice()[i] << " ";
  }
  std::cout << std::endl;

  rows[0] = 1;
  rows[1] = 2;
  elements[0] = 1;
  elements[1] = 1;
  clpSolver.addCol(2, rows, elements, 0, 1, 0);
  filename = "LP_" + std::to_string(3);
    clpSolver.writeLp(filename.c_str());
  clpSolver.resolve();

  // dump duals
  std::cout << "iteration 3  - Duals: ";
  for(auto i=0; i< clpSolver.getNumRows(); i++){
    std::cout <<  clpSolver.getRowPrice()[i] << " ";
  }
  std::cout << std::endl;

  rows[0] = 1;
  rows[1] = 2;
  elements[0] = 1;
  elements[1] = 1;
  clpSolver.addCol(2, rows, elements, 0, 1, 0);
  filename = "LP_" + std::to_string(4);
    clpSolver.writeLp(filename.c_str());
  clpSolver.resolve();

  // dump duals
  std::cout << "iteration 4  - Duals: ";
  for(auto i=0; i< clpSolver.getNumRows(); i++){
    std::cout <<  clpSolver.getRowPrice()[i] << " ";
  }
  std::cout << std::endl;
jjhforrest commented 3 years ago

Simon,

The reason you are getting those results is because after iteration 1 - see below -

On 07/12/2020 13:12, Simon Spoorendonk wrote:

I have a simple column generation problem where I need to add the columns twice (with |resolve()| twice) to get the correct dual variables. Am I missing a setting somewhere?

I use OsiClpSolverInterface, but get the same result with ClpSimplex directly.

My examples is:

In iteration 0, to obtain feasibility, I start with a phase 1 and have only artificial variables like

|\Problem name: ClpDefaultName Minimize obj: x0 + x1 Subject To cons0: x0 = 1 cons1: x1 = 1 cons2: x0 + x1 <= 2 cons2_low: x0 + x1 >= 1 Bounds 0 <= x0 <= 1 0 <= x1 <= 1 End |

with duals |[1, 1, 0] | after an |initialSolve()|.

Now I price in variables similar to |x0| and |x1| but with cost 0 to price out the artificial ones, so iteration 1 is

|\Problem name: ClpDefaultName Minimize obj: x0 + x1 Subject To cons0: x0 + x2 = 1 cons1: x1 = 1 cons2: x0 + x1 + x2 <= 2 cons2_low: x0 + x1 + x2 >= 1 Bounds 0 <= x0 <= 1 0 <= x1 <= 1 0 <= x2 <= 1 End |

then after |resolve()| i get duals |[1, 1, 0] | again. The problem is proven optimal and the objevtive value is 1, so the dual solution is off?

On the resolve x2 flips from 0 to 1 so - the primal solution is 0, 1, 1 but the reduced costs are 0, 0, -1 as x2 is at its upper bound. With the implied row x2 <= 1 the arithmetic works out.

If upper bounds are large then it works as expected.

John Forrest

spoorendonk commented 3 years ago

Thanks.

Forgive me if my LP theory is a little rusty. Is it possible to force a resolve/pivot such that all reduced costs are >= 0? Which I guess would lead to the dual solution I am looking for?

jjhforrest commented 3 years ago

Not if any upper bounds are active.

So for instance if you modify LP_1.lp to give x2 an upper bound of 0.9, then the optimal solution of 1.1 has to have x2 at upper bound with reduced cost of -1. As I wrote in my previous note you would need to have upper bounds on all variables as infinity.

On 07/12/2020 18:40, Simon Spoorendonk wrote:

Thanks.

Forgive me if my LP theory is a little rusty. Is it possible to force a resolve/pivot such that all reduced costs are >= 0? Which I guess would lead to the dual solution I am looking for?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/coin-or/Clp/issues/167#issuecomment-740103569, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABWJYHFMXI53IDLPJWCBHP3STUORRANCNFSM4UQOCYOA.

spoorendonk commented 3 years ago

Allright. Infinity it is. Thanks again.