In practice, it is very useful to get a good solution reasonably fast. A good bound will greatly reduce the run time and good solutions can satisfy the user on very large problems where a complete search is impossible. Obviously, heuristics are problem dependent although some have more general use. At present there is only one in CBC itself. Hopefully, the number of heuristics will grow. Other hueristics are in the COIN/Cbc/Samples directory. A heuristic tries to obtain a solution to the original problem so it only needs to consider the original rows and does not have to use the current bounds. One to use a greedy heuristic designed for use in the miplib problem fast0507 will be developed later in this section. CBC provides an abstract base class CbcHeuristic and a rounding heuristic in CBC.
This describes how to build a greedy heuristic for a set covering problem. A more general version is in CbcHeuristicGreedy.hpp and CbcHeuristicGreedy.cpp which can be found in the COIN/Cbc/Samples directory, see Chapter 6, More Samples . The heuristic we will code will leave all variables which are at one at this node of the tree to that value and will initially set all others to zero. We then sort all variables in order of their cost divided by the number of entries in rows which are not yet covered. We may randomize that value a bit so that ties will be broken in different ways on different runs of the heuristic. We then choose the best one and set it to one and repeat the exercise. Because this is a set covering problem ≥ we are guaranteed to find a solution (not necessarily a better one though). We could improve the speed by just redoing those affected but in this text we will keep it simple. Also if all elements are 1.0 then we could do faster. The key CbcHeuristic method is int solution(double & solutionValue, double * betterSolution) which returns 0 if no solution found and 1 if found when it fills in the objective value and primal solution. The actual code in CbcHeuristicGreedy.cpp is a little more complicated but this will work and gives the basic idea. For instance the code here assumes all variables are integer. The important bit of data is a copy of the matrix (stored by column) before any cuts have been made. The data used are bounds, objective and the matrix plus two work arrays.
Example 3.4. Data
OsiSolverInterface * solver = model_->solver(); // Get solver from CbcModel
const double * columnLower = solver->getColLower(); // Column Bounds
const double * columnUpper = solver->getColUpper();
const double * rowLower = solver->getRowLower(); // We know we only need lower bounds
const double * solution = solver->getColSolution();
const double * objective = solver->getObjCoefficients(); // In code we also use min/max
double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
double primalTolerance;
solver->getDblParam(OsiPrimalTolerance,primalTolerance);
int numberRows = originalNumberRows_; // This is number of rows when matrix was passed in
// Column copy of matrix (before cuts)
const double * element = matrix_.getElements();
const int * row = matrix_.getIndices();
const CoinBigIndex * columnStart = matrix_.getVectorStarts();
const int * columnLength = matrix_.getVectorLengths();
// Get solution array for heuristic solution
int numberColumns = solver->getNumCols();
double * newSolution = new double [numberColumns];
// And to sum row activities
double * rowActivity = new double[numberRows];
Then we initialize newSolution as rounded down solution.
Example 3.5. initialize newSolution
for (iColumn=0;iColumn<numberColumns;iColumn++) {
CoinBigIndex j;
double value = solution[iColumn];
// Round down integer
if (fabs(floor(value+0.5)-value)<integerTolerance)
value=floor(CoinMax(value+1.0e-3,columnLower[iColumn]));
// make sure clean
value = CoinMin(value,columnUpper[iColumn]);
value = CoinMax(value,columnLower[iColumn]);
newSolution[iColumn]=value;
if (value) {
double cost = direction * objective[iColumn];
newSolutionValue += value*cost;
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int iRow=row[j];
rowActivity[iRow] += value*element[j];
}
}
}
Now some row activities will be below their lower bound so we then find the variable which is cheapest in reducing the sum of infeasibilities. We then repeat. This is a finite process and could be coded to be faster but this is simplest.
Example 3.6. Create feasible new solution
while (true) {
// Get column with best ratio
int bestColumn=-1;
double bestRatio=COIN_DBL_MAX;
for (int iColumn=0;iColumn<numberColumns;iColumn++) {
CoinBigIndex j;
double value = newSolution[iColumn];
double cost = direction * objective[iColumn];
// we could use original upper rather than current
if (value+0.99<columnUpper[iColumn]) {
double sum=0.0; // Compute how much we will reduce infeasibility by
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int iRow=row[j];
double gap = rowLower[iRow]-rowActivity[iRow];
if (gap>1.0e-7) {
sum += CoinMin(element[j],gap);
if (element[j]+rowActivity[iRow]<rowLower[iRow]+1.0e-7) {
sum += element[j];
}
}
if (sum>0.0) {
double ratio = (cost/sum)*(1.0+0.1*CoinDrand48());
if (ratio<bestRatio) {
bestRatio=ratio;
bestColumn=iColumn;
}
}
}
}
if (bestColumn<0)
break; // we have finished
// Increase chosen column
newSolution[bestColumn] += 1.0;
double cost = direction * objective[bestColumn];
newSolutionValue += cost;
for (CoinBigIndex j=columnStart[bestColumn];
j<columnStart[bestColumn]+columnLength[bestColumn];j++) {
int iRow = row[j];
rowActivity[iRow] += element[j];
}
}
We have finished so now we need to see if solution is better and doublecheck we are feasible.
Example 3.7. Check good solution
returnCode=0; // 0 means no good solution
if (newSolutionValue<solutionValue) {
// check feasible
memset(rowActivity,0,numberRows*sizeof(double));
for (iColumn=0;iColumn<numberColumns;iColumn++) {
CoinBigIndex j;
double value = newSolution[iColumn];
if (value) {
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int iRow=row[j];
rowActivity[iRow] += value*element[j];
}
}
}
// check was approximately feasible
bool feasible=true;
for (iRow=0;iRow<numberRows;iRow++) {
if(rowActivity[iRow]<rowLower[iRow]) {
if (rowActivity[iRow]<rowLower[iRow]-10.0*primalTolerance)
feasible = false;
}
}
if (feasible) {
// new solution
memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
solutionValue = newSolutionValue;
// We have good solution
returnCode=1;
}
}