Help Accessing the Dual Solutions for the Initial Solve of LP Relaxation

[spkelle2]

I have a subclass of a ClpEventHandler that I have set up to record dual solutions whenever CBC is solving an LP relaxation with dual simplex. I'm having good luck printing out the dual solutions (i.e. dualColumnSolution() and dualRowSolution()) while any heuristics are running or after resolving a node once cuts have been added, but I can't seem to get the solution to print while doing the initial solve. To clarify, I do get a solution, but all dual variable values print out to 0 with exception of possibly the last iteration. I would expect that since I am performing the initial solve with dual simplex that I would get nonzero values for these printouts for at least some of the variables. Could someone help me understand why this is the case/what I am misunderstanding and if there is anything else I can do to find the dual solution values during initial solve? If it helps, I have a copy of my code and its stdout below.

call arguments: -import /Users/sean/coin-or/build-debug/Data/Sample/1.2/p0033.mps -preprocess off -presolve off -solve

code:

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <OsiCbcSolverInterface.hpp>
#include "OsiClpSolverInterface.hpp"
#include <CbcModel.hpp>
#include <CglPreProcess.hpp>
#include <CbcSolver.hpp>
#include <Cbc_C_Interface.h>
#include "CbcEventHandler.hpp"
#include "ClpEventHandler.hpp"

static int callBack(CbcModel *model, int whereFrom)
{
  return 0;
}

class MyEventHandler3 : public ClpEventHandler {

public:
  /**@name Overrides */
  //@{
  virtual int event(Event whichEvent);
  //@}

  /**@name Constructors, destructor etc*/
  //@{
  /** Default constructor. */
  MyEventHandler3();
  /// Constructor with pointer to model (redundant as setEventHandler does)
  MyEventHandler3(ClpSimplex *model);
  /** Destructor */
  virtual ~MyEventHandler3();
  /** The copy constructor. */
  MyEventHandler3(const MyEventHandler3 &rhs);
  /// Assignment
  MyEventHandler3 &operator=(const MyEventHandler3 &rhs);
  /// Clone
  virtual ClpEventHandler *clone() const;
  //@}

protected:
  // data goes here
};
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
MyEventHandler3::MyEventHandler3()
    : ClpEventHandler()
{
}

//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
MyEventHandler3::MyEventHandler3(const MyEventHandler3 &rhs)
    : ClpEventHandler(rhs)
{
}

// Constructor with pointer to model
MyEventHandler3::MyEventHandler3(ClpSimplex *model)
    : ClpEventHandler(model)
{
}

//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
MyEventHandler3::~MyEventHandler3()
{
}

//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
MyEventHandler3 &
MyEventHandler3::operator=(const MyEventHandler3 &rhs)
{
  if (this != &rhs) {
    ClpEventHandler::operator=(rhs);
  }
  return *this;
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
ClpEventHandler *MyEventHandler3::clone() const
{
  return new MyEventHandler3(*this);
}

int MyEventHandler3::event(Event whichEvent)
{
  // After each dual iteration
  if (whichEvent == endOfIteration && model_->algorithm() < 0){
    std::cout << "Iteration " << this->model_->numberIterations() << " dual solutions \n";
    // column dual variables
    std::cout << "s = [";
    for (int col_idx = 0; col_idx < this->model_->numberColumns(); col_idx++) {
      std::cout << this->model_->dualColumnSolution()[col_idx] << " ";
    }
    std::cout << "]\n";
    // row dual variables
    std::cout << "y = [";
    for (int row_idx = 0; row_idx < this->model_->numberRows(); row_idx++) {
      std::cout << this->model_->dualRowSolution()[row_idx] << " ";
    }
    std::cout << "]\n";
    std::cout << "\n";
    return -1;
  } else {
    return -1;
  }
}

int main(int argc, char **argv)
{
  if (argc < 2) {
    fprintf(stderr, "enter instance name\n");
    exit(EXIT_FAILURE);
  }

  CbcSolverUsefulData cbcData;
  cbcData.noPrinting_ = false;

  // set up lp solver
  OsiClpSolverInterface lp;
  lp.readMps(argv[1]);
  ClpSimplex *lp_model = lp.getModelPtr();
  MyEventHandler3 eventHandler(lp_model);
  lp_model->passInEventHandler(&eventHandler);

  // set up branch and cut solver
  CbcModel model(lp);
  model.setMaximumNodes(1);

  CbcMain0(model, cbcData);
  CbcMain1(argc - 1, (const char **)(argv + 1), model, callBack, cbcData);

  return 0;
}

stdout (through initial solve - I assume I don't have an issue reading the mps file since it otherwise solves):

/Users/sean/coin-or/build-debug/Cbc/2.10/examples/nik_simple -import /Users/sean/coin-or/build-debug/Data/Sample/1.2/p0033.mps -preprocess off -presolve off -solve
Coin6001E Unable to open mps input file -import.mps
Coin0008I  read with -1 errors
Welcome to the CBC MILP Solver 
Version: 2.10 
Build Date: Jul 14 2022 

command line - -import /Users/sean/coin-or/build-debug/Data/Sample/1.2/p0033.mps -preprocess off -presolve off -solve (default strategy 1)
At line 15 NAME          P0033
At line 16 ROWS
At line 34 COLUMNS
At line 109 RHS
At line 118 BOUNDS
At line 152 ENDATA
Problem P0033 has 16 rows, 33 columns and 98 elements
Coin0008I P0033 read with 0 errors
Option for preprocess changed from sos to off
Option for presolve changed from on to off
Iteration 1 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 2 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 3 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 4 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 5 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 6 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 7 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 8 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 9 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 10 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 11 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 12 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 13 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 14 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 15 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 16 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 17 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 18 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 19 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Iteration 20 dual solutions 
s = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
y = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Continuous objective value is 2520.57 - 0.00 seconds
Cutoff increment increased from 1e-05 to 0.9999

[jjhforrest]

Because of scaling and other issues, the external dual and reduced cost regions are not current at each iteration. You will need to access model->djRegion(0/1) and if model->rowScale() is null then use them, if not null scale by rowScale() and columnScale().

[spkelle2]

That is super helpful. Thank you for that call out and quick reply @jjhforrest. Let me give that a try.

[spkelle2]

For anyone curious, here's the updated version of my code printing out dual values on the primal constraints and variables for the root relaxation. I added a couple of attributes to the event handler so that I could filter for only reporting duals when the problem was in the space of the original formulation. Again, I needed to turn off presolve and preprocessing for this to work (i.e. add -preprocess off -presolve off to your call args).

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <OsiCbcSolverInterface.hpp>
#include "OsiClpSolverInterface.hpp"
#include <CbcModel.hpp>
#include <CglPreProcess.hpp>
#include <CbcSolver.hpp>
#include <Cbc_C_Interface.h>
#include "CbcEventHandler.hpp"
#include "ClpEventHandler.hpp"

static int callBack(CbcModel *model, int whereFrom)
{
  return 0;
}

class MyEventHandler3 : public ClpEventHandler {

public:
  /**@name Overrides */
  //@{
  virtual int event(Event whichEvent);
  //@}

  /**@name Constructors, destructor etc*/
  //@{
  /** Default constructor. */
  MyEventHandler3();
  /// Constructor with pointer to model (redundant as setEventHandler does)
  MyEventHandler3(ClpSimplex *model);
  /** Destructor */
  virtual ~MyEventHandler3();
  /** The copy constructor. */
  MyEventHandler3(const MyEventHandler3 &rhs);
  /// Assignment
  MyEventHandler3 &operator=(const MyEventHandler3 &rhs);
  /// Clone
  virtual ClpEventHandler *clone() const;
  //@}

  int rows_;
  int cols_;
  bool called_before_;

protected:
  // data goes here
};
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
MyEventHandler3::MyEventHandler3()
    : ClpEventHandler()
{
  rows_ = 0;
  cols_ = 0;
  called_before_ = false;
}

//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
MyEventHandler3::MyEventHandler3(const MyEventHandler3 &rhs)
    : ClpEventHandler(rhs)
{
  rows_ = rhs.rows_;
  cols_ = rhs.cols_;
  called_before_ = rhs.called_before_;
}

// Constructor with pointer to model
MyEventHandler3::MyEventHandler3(ClpSimplex *model)
    : ClpEventHandler(model)
{
  rows_ = 0;
  cols_ = 0;
  called_before_ = false;
}

//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
MyEventHandler3::~MyEventHandler3()
{
}

//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
MyEventHandler3 &
MyEventHandler3::operator=(const MyEventHandler3 &rhs)
{
  if (this != &rhs) {
    ClpEventHandler::operator=(rhs);
  }
  return *this;
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
ClpEventHandler *MyEventHandler3::clone() const
{
  return new MyEventHandler3(*this);
}

int MyEventHandler3::event(Event whichEvent)
{
  // get the size of the problem if we haven't before
  if (!called_before_){
    rows_ = this->model_->numberRows();
    cols_ = this->model_->numberColumns();
    called_before_ = true;
  }
  // After each dual iteration
  if (whichEvent == endOfIteration && model_->algorithm() < 0){
    // if working on the full problem, get the dual solutions
    if (this->model_->numberColumns() == cols_ && this->model_->numberRows() == rows_) {
      std::cout << "Iteration " << this->model_->numberIterations() << " dual solutions \n";
      // column dual variables
      std::cout << "s = [ ";
      for (int col_idx = 0; col_idx < this->model_->numberColumns(); col_idx++) {
        if (this->model_->rowScale()) {
          std::cout << this->model_->djRegion(1)[col_idx] /
                       this->model_->columnScale()[col_idx] << " ";
        }
        else {
          std::cout << this->model_->djRegion(1)[col_idx] << " ";
        }
      }
      std::cout << "]\n";
      // row dual variables
      std::cout << "y = [ ";
      for (int row_idx = 0; row_idx < this->model_->numberRows(); row_idx++) {
        if (this->model_->rowScale()) {
          std::cout << this->model_->djRegion(0)[row_idx] *
                       this->model_->rowScale()[row_idx] << " ";
        }
        else {
          std::cout << this->model_->djRegion(0)[row_idx] << " ";
        }
      }
      std::cout << "]\n";
      std::cout << "\n";
    }
  }
  return -1; // carry on
}

int main(int argc, char **argv)
{
  if (argc < 2) {
    fprintf(stderr, "enter instance name\n");
    exit(EXIT_FAILURE);
  }

  CbcSolverUsefulData cbcData;
  cbcData.noPrinting_ = false;

  // set up lp solver
  OsiClpSolverInterface lp;
  lp.readMps(argv[1]);
  ClpSimplex *lp_model = lp.getModelPtr();
  MyEventHandler3 eventHandler(lp_model);
  lp_model->passInEventHandler(&eventHandler);

  // set up branch and cut solver
  CbcModel model(lp);
  model.setMaximumNodes(1);

  CbcMain0(model, cbcData);
  CbcMain1(argc - 1, (const char **)(argv + 1), model, callBack, cbcData);

  return 0;
}