source: branches/gh-pages/heuristics.md @ 2521

Last change on this file since 2521 was 2521, checked in by stefan, 2 months ago

migrate Cbc documentation to gh-pages

  • converted docbook to markdown
  • source now in gh-pages branch, so remove Cbc/doc
  • some cleanup
  • update link in README.md
  • Property svn:eol-style set to native
File size: 7.0 KB
Line 
1## Getting Good Bounds in CBC
2
3# CbcHeuristic - Heuristic Methods
4
5In practice, it is very useful to get a good solution reasonably fast.
6Any MIP-feasible solution produces an upper bound, and a good bound will
7greatly reduce the run time. Good solutions can satisfy the user on very
8large problems where a complete search is impossible. Obviously,
9heuristics are problem dependent, although some do have more general
10use. At present there is only one heuristic in CBC itself,
11`CbcRounding`. Hopefully, the number will grow. Other heuristics are in
12the `Cbc/examples` directory. A heuristic tries to obtain a solution
13to the original problem so it only needs to consider the original rows
14and does not have to use the current bounds. CBC provides an abstract
15base class `CbcHeuristic` and a rounding heuristic in CBC.
16
17This chapter describes how to build a greedy heuristic for a set
18covering problem, e.g., the miplib problem fast0507. A more general (and
19efficient) version of the heuristic is in `CbcHeuristicGreedy.hpp` and
20`CbcHeuristicGreedy.cpp` located in the `Cbc/examples` directory.
21
22The greedy heuristic will leave all variables taking value one at this
23node of the tree at value one, and will initially set all other variable
24to value zero. All variables are then sorted in order of their cost
25divided by the number of entries in rows which are not yet covered. (We
26may randomize that value a bit so that ties will be broken in different
27ways on different runs of the heuristic.) The best one is choosen, and
28set to one. The process is repeated. Because this is a set covering
29problem (i.e., all constraints are GE), the heuristic is guaranteed to
30find a solution (but not necessarily an improved solution). The speed of
31the heuristic could be improved by just redoing those affected, but for
32illustrative purposes we will keep it simple.(The speed could also be
33improved if all elements are 1.0).
34
35The key `CbcHeuristic` method is `int solution(double & solutionValue,
36double * betterSolution)`. The `solution()` method returns 0 if no
37solution found, and returns 1 if a solution is found, in which case it
38fills in the objective value and primal solution. The code in
39`CbcHeuristicGreedy.cpp` is a little more complicated than this
40following example. For instance, the code here assumes all variables are
41integer. The important bit of data is a copy of the matrix (stored by
42column) before any cuts have been made. The data used are bounds,
43objective and the matrix plus two work arrays.
44
45```
46OsiSolverInterface * solver = model_->solver(); // Get solver from CbcModel
47const double * columnLower = solver->getColLower(); // Column Bounds
48const double * columnUpper = solver->getColUpper();
49const double * rowLower = solver->getRowLower(); // We know we only need lower bounds
50const double * solution = solver->getColSolution();
51const double * objective = solver->getObjCoefficients(); // In code we also use min/max
52double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
53double primalTolerance;
54solver->getDblParam(OsiPrimalTolerance,primalTolerance);
55int numberRows = originalNumberRows_; // This is number of rows when matrix was passed in
56// Column copy of matrix (before cuts)
57const double * element = matrix_.getElements();
58const int * row = matrix_.getIndices();
59const CoinBigIndex * columnStart = matrix_.getVectorStarts();
60const int * columnLength = matrix_.getVectorLengths();
61
62// Get solution array for heuristic solution
63int numberColumns = solver->getNumCols();
64double * newSolution = new double [numberColumns];
65// And to sum row activities
66double * rowActivity = new double[numberRows];
67```
68
69The `newSolution` is then initialized to the rounded down solution:
70```
71for (iColumn=0;iColumn<numberColumns;iColumn++) {
72  CoinBigIndex j;
73  double value = solution[iColumn];
74  // Round down integer
75  if (fabs(floor(value+0.5)-value)<integerTolerance)
76    value=floor(CoinMax(value+1.0e-3,columnLower[iColumn]));
77  // make sure clean
78  value = CoinMin(value,columnUpper[iColumn]);
79  value = CoinMax(value,columnLower[iColumn]);
80  newSolution[iColumn]=value;
81  if (value) {
82    double cost = direction * objective[iColumn];
83    newSolutionValue += value*cost;
84    for (j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
85      int iRow=row[j];
86      rowActivity[iRow] += value*element[j];
87    }
88  }
89}
90```
91
92At this point some row activities may be below their lower bound. To
93correct this infeasibility, the variable which is cheapest in reducing
94the sum of infeasibilities is found and updated, and the process
95repeats. This is a finite process. (Theimplementation could be faster,
96but is kept simple for illustrative purposes.)
97
98```
99while (true) {
100  // Get column with best ratio
101  int bestColumn=-1;
102  double bestRatio=COIN_DBL_MAX;
103  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
104    CoinBigIndex j;
105    double value = newSolution[iColumn];
106    double cost = direction * objective[iColumn];
107    // we could use original upper rather than current
108    if (value+0.99<columnUpper[iColumn]) {
109      double sum=0.0; // Compute how much we will reduce infeasibility by
110      for (j=columnStart[iColumn];
111           j<columnStart[iColumn]+columnLength[iColumn];j++) {
112        int iRow=row[j];
113        double gap = rowLower[iRow]-rowActivity[iRow];
114        if (gap>1.0e-7)
115          sum += CoinMin(element[j],gap);
116        if (element[j]+rowActivity[iRow]<rowLower[iRow]+1.0e-7)
117          sum += element[j];
118      }
119      if (sum>0.0) {
120        double ratio = (cost/sum)*(1.0+0.1*CoinDrand48());
121        if (ratio<bestRatio) {
122          bestRatio=ratio;
123          bestColumn=iColumn;
124        }
125      }
126    }
127  }
128  if (bestColumn<0)
129    break; // we have finished
130  // Increase chosen column
131  newSolution[bestColumn] += 1.0;
132  double cost = direction * objective[bestColumn];
133  newSolutionValue += cost;
134  for (CoinBigIndex j=columnStart[bestColumn];
135    j<columnStart[bestColumn]+columnLength[bestColumn];j++) {
136    int iRow = row[j];
137    rowActivity[iRow] += element[j];
138  }
139}
140```
141
142A solution value of `newSolution` is compared to the best solution
143value. If `newSolution` is an improvement, its feasibility is validated.
144```
145returnCode=0; // 0 means no good solution
146if (newSolutionValue<solutionValue) { // minimization
147  // check feasible
148  memset(rowActivity,0,numberRows*sizeof(double));
149  for (iColumn=0;iColumn<numberColumns;iColumn++) {
150    CoinBigIndex j;
151    double value = newSolution[iColumn];
152    if (value) {
153      for (j=columnStart[iColumn];
154           j<columnStart[iColumn]+columnLength[iColumn];j++) {
155        int iRow=row[j];
156        rowActivity[iRow] += value*element[j];
157      }
158    }
159  }
160  // check was approximately feasible
161  bool feasible=true;
162  for (iRow=0;iRow<numberRows;iRow++) {
163    if(rowActivity[iRow]<rowLower[iRow])
164      if (rowActivity[iRow]<rowLower[iRow]-10.0*primalTolerance)
165        feasible = false;
166  }
167  if (feasible) {
168    // new solution
169    memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
170    solutionValue = newSolutionValue;
171    // We have good solution
172    returnCode=1;
173  }
174}
175```
Note: See TracBrowser for help on using the repository browser.