[Coin-lpsolver] Primal search algorithm for nonlinear objective

[rmj at itu.dk]

Hi,

I'm considering a specific problem on the form

min c(x) s.t. Ax <= b, x >= 0.

The cost function c(x) is nonlinear, concave and there are no derivatives
available. For the considered problem, however, optimal solutions can be
assumed to be corner points of the polyhedron Ax <= b.  A prototype
implemented in Octave shows that it may be a good idea to consider a
primal simplex-like improvement search method where a pivot is taken if it
leads to a reduction in c(x). That is, in each iteration, consider a set
of pivot column candidates and their associated pivot rows (found using
usual ratio test to maintain x >= 0) and perform the pivot that reduces
c(x) the most.

The prototype, however, is too slow for the real problem and the question
is whether it is possible to implement the algorithm using the highly
optimized pivot routines of Clp. I’m new to Clp and so far I have
considered the following approaches:

1)	Use the Simplex primal algorithm with a pivot algorithm derived from
ClpPrimalColumnDantzig. I’m unsure, however, if this will work out. There
are many assumptions in ClpSimplexPrimal that do not fit with my
algorithm. For instance, the pivot algorithm in ClpSimplexPrimal only
picks a pivot column, the pivot row is computed from a ratio test. But in
an improvement search, it will be necessary to consider the delta cost
from a candidate set of complete pivots, so picking the column and row
cannot be separated in two functions. In addition, there is no linear
cost, so whether a pivot reduces it is irrelevant. Making changes in
ClpSimplexPrimal to accommodate all this seems very hard unless you know
the Clp implementation in detail.

2)	For the OSI interface, Clp implements a pivot function that simply
performs a pivot on an ingoing and outgoing variable. I wonder if it could
be used to implement the algorithm (without knowing the details of Clp),
basically performing a series of pivots reducing c(x). One issue, though,
is that the function not only requires sequenceIn_, pivotRow_, and
directionIn_ and Out to be set but also e.g. the  lower_ array (lower_ is
instantiated by Primal simplex but set to null before primal returns). Is
the pivot function intended for a more specific use? If so, is there a
simple way to get to a state where it can be used as required by my
algorithm?  The code example below gives a segment fault due to lower_ ==
Null, but is there a simple way to extend it to make it work fine for a
series of pivots (given that the pivots are legal)?

3)	Use one of the nonlinear optimization algorithms that already exist
such as ClpSimplex::nonlinearSLP or  ClpSimplex:: reducedGradient. But it
seems that these functions are mainly for quadratic problems or at least
require a convex objective. Moreover, I have not found a way to load a
general cost function, so I have assumed that it is not possible. Is that
correct?

Thanks a lot for any help on these questions and hints on which approach
to use.  My impression is that it’s not too hard. I just don’t have enough
knowledge about the inner workings of Clp to do it reliably.

Best regards,
Rune M. Jensen

...
  // Create a new model
  ClpSimplex model;

  //  optimal solution: x* = (1,1)
  //
  //  minimize -1 x0 - 1 x1
  //  s.t       1 x0 + 2 x1 <= 3
  //            2 x0 + 1 x1 <= 3
  //              x0        >= 0
  //              x1        >= 0

  int n_cols = 2;
  double *objective    = new double[n_cols];//the objective coefficients
  double *col_lb       = new double[n_cols];//the column lower bounds
  double *col_ub       = new double[n_cols];//the column upper bounds


  //Define the variable lower/upper bounds.
  // x0 >= 0   =>  0 <= x0 <= infinity
  // x1 >= 0   =>  0 <= x1 <= infinity
  col_lb[0] = 0.0;
  col_lb[1] = 0.0;
  col_ub[0] = OsiClpInfinity;
  col_ub[1] = OsiClpInfinity;

  int n_rows = 2;
  double *row_lb = new double[n_rows]; //the row lower bounds
  double *row_ub = new double[n_rows]; //the row upper bounds

  //Define the constraint matrix.
  CoinPackedMatrix *matrix =  new CoinPackedMatrix(false,0,0);
  matrix->setDimensions(0, n_cols);

  //1 x0 + 2 x1 <= 3  =>  -infinity <= 1 x0 + 2 x2 <= 3
  CoinPackedVector row1;
  row1.insert(0, 1.0);
  row1.insert(1, 2.0);
  row_lb[0] = -1.0 * OsiClpInfinity;
  row_ub[0] = 3.0;
  matrix->appendRow(row1);

  //2 x0 + 1 x1 <= 3  =>  -infinity <= 2 x0 + 1 x1 <= 3
  CoinPackedVector row2;
  row2.insert(0, 2.0);
  row2.insert(1, 1.0);
  row_lb[1] = -1.0 * OsiClpInfinity;
  row_ub[1] = 3.0;
  matrix->appendRow(row2);

  //Define the objective coefficients.
  //minimize -1 x0 - 1 x1
  objective[0] = -1.0;
  objective[1] = -1.0;

  model.setLogLevel(4);

  //load the problem to OSI
  model.loadProblem(*matrix, col_lb, col_ub, objective, row_lb, row_ub);

  model.setSequenceIn(0);  // this is the first pivot taken by primal
  model.setPivotRow(1);    // simplex
  model.setDirectionIn(1);
  model.setDirectionOut(-1);
  model.scaling(0);
  model.pivot();
...