| 1 | <?xml version="1.0" encoding="ISO-8859-1" standalone="no"?> |
|---|
| 2 | <html xmlns="http://www.w3.org/1999/xhtml"><head><title>CbcHeuristic - Heuristic Methods</title><meta name="generator" content="DocBook XSL Stylesheets V1.61.2"/><link rel="home" href="index.html" title="CBC User Guide"/><link rel="up" href="ch03.html" title="Chapter 3. Selecting a Node in the Search Tree "/><link rel="previous" href="ch03.html" title="Chapter 3. Selecting a Node in the Search Tree "/><link rel="next" href="ch04.html" title="Chapter 4. Branching "/></head><body><div class="navheader"><table width="100%" summary="Navigation header"><tr><th colspan="3" align="center">CbcHeuristic - Heuristic Methods</th></tr><tr><td width="20%" align="left"><a accesskey="p" href="ch03.html">Prev</a> </td><th width="60%" align="center">Chapter 3. |
|---|
| 3 | Selecting a Node in the Search Tree |
|---|
| 4 | </th><td width="20%" align="right"> <a accesskey="n" href="ch04.html">Next</a></td></tr></table><hr/></div><div class="section" lang="en"><div class="titlepage"><div><div><h2 class="title" style="clear: both"><a id="heuristics"/>CbcHeuristic - Heuristic Methods</h2></div></div><div/></div><p> |
|---|
| 5 | In practice, it is very useful to get a good solution reasonably fast. |
|---|
| 6 | A good bound will greatly reduce the run time and good solutions can satisfy the user |
|---|
| 7 | on very large problems where a complete search is impossible. Obviously, heuristics are |
|---|
| 8 | problem dependent although some have more general use. |
|---|
| 9 | At present there is only one in CBC itself. Hopefully, the number of heuristics will grow. |
|---|
| 10 | Other hueristics are in the <tt class="filename">COIN/Cbc/Samples</tt> |
|---|
| 11 | directory. A heuristic tries to obtain a solution to the original |
|---|
| 12 | problem so it only needs to consider the original rows and does not have to use the |
|---|
| 13 | current bounds. |
|---|
| 14 | One to use a greedy heuristic designed for use in the miplib problem |
|---|
| 15 | fast0507 will be developed later in this section. |
|---|
| 16 | CBC provides an abstract base class <tt class="classname">CbcHeuristic</tt> and a rounding heuristic in CBC. |
|---|
| 17 | </p><p> |
|---|
| 18 | This describes how to build a greedy heuristic for a set covering problem. |
|---|
| 19 | A more general version is in <tt class="filename">CbcHeuristicGreedy.hpp</tt> and |
|---|
| 20 | <tt class="filename">CbcHeuristicGreedy.cpp</tt> which can be found in the <tt class="filename">COIN/Cbc/Samples</tt> directory, see <a href="ch06.html" title="Chapter 6. More Samples ">Chapter 6, <i> |
|---|
| 21 | More Samples |
|---|
| 22 | </i></a>. |
|---|
| 23 | |
|---|
| 24 | The heuristic we will code will leave all variables which are at one at this node of the |
|---|
| 25 | tree to that value and will |
|---|
| 26 | initially set all others to zero. We then sort all variables in order of their cost |
|---|
| 27 | divided by the number of entries in rows which are not yet covered. We may randomize that |
|---|
| 28 | value a bit so that ties will be broken in different ways on different runs of the heuristic. |
|---|
| 29 | We then choose the best one and set it to one and repeat the exercise. Because this is |
|---|
| 30 | a set covering problem ≥ we are guaranteed to find a solution (not necessarily a |
|---|
| 31 | better one though). We could |
|---|
| 32 | improve the speed by just redoing those affected but in this text we will keep it simple. |
|---|
| 33 | Also if all elements are 1.0 then we could do faster. |
|---|
| 34 | The key CbcHeuristic method is <tt class="function">int solution(double & solutionValue, |
|---|
| 35 | double * betterSolution)</tt> |
|---|
| 36 | which returns 0 if no solution found and 1 if found when it fills in the objective value |
|---|
| 37 | and primal solution. The actual code in <tt class="filename">CbcHeuristicGreedy.cpp</tt> is |
|---|
| 38 | a little more complicated but this will work and gives the basic idea. For instance |
|---|
| 39 | the code here assumes all variables are integer. |
|---|
| 40 | The important bit of data is a copy of the matrix (stored by column) |
|---|
| 41 | before any cuts have been made. The data used are bounds, objective and the matrix |
|---|
| 42 | plus two work arrays. |
|---|
| 43 | </p><div class="example"><a id="id2984921"/><p class="title"><b>Example 3.4. Data</b></p><pre class="programlisting"> |
|---|
| 44 | |
|---|
| 45 | OsiSolverInterface * solver = model_->solver(); // Get solver from CbcModel |
|---|
| 46 | const double * columnLower = solver->getColLower(); // Column Bounds |
|---|
| 47 | const double * columnUpper = solver->getColUpper(); |
|---|
| 48 | const double * rowLower = solver->getRowLower(); // We know we only need lower bounds |
|---|
| 49 | const double * solution = solver->getColSolution(); |
|---|
| 50 | const double * objective = solver->getObjCoefficients(); // In code we also use min/max |
|---|
| 51 | double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance); |
|---|
| 52 | double primalTolerance; |
|---|
| 53 | solver->getDblParam(OsiPrimalTolerance,primalTolerance); |
|---|
| 54 | int numberRows = originalNumberRows_; // This is number of rows when matrix was passed in |
|---|
| 55 | // Column copy of matrix (before cuts) |
|---|
| 56 | const double * element = matrix_.getElements(); |
|---|
| 57 | const int * row = matrix_.getIndices(); |
|---|
| 58 | const CoinBigIndex * columnStart = matrix_.getVectorStarts(); |
|---|
| 59 | const int * columnLength = matrix_.getVectorLengths(); |
|---|
| 60 | |
|---|
| 61 | // Get solution array for heuristic solution |
|---|
| 62 | int numberColumns = solver->getNumCols(); |
|---|
| 63 | double * newSolution = new double [numberColumns]; |
|---|
| 64 | // And to sum row activities |
|---|
| 65 | double * rowActivity = new double[numberRows]; |
|---|
| 66 | |
|---|
| 67 | </pre></div><p> |
|---|
| 68 | Then we initialize newSolution as rounded down solution. |
|---|
| 69 | </p><div class="example"><a id="id2984946"/><p class="title"><b>Example 3.5. initialize newSolution</b></p><pre class="programlisting"> |
|---|
| 70 | |
|---|
| 71 | for (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]; |
|---|
| 85 | j<columnStart[iColumn]+columnLength[iColumn];j++) { |
|---|
| 86 | int iRow=row[j]; |
|---|
| 87 | rowActivity[iRow] += value*element[j]; |
|---|
| 88 | } |
|---|
| 89 | } |
|---|
| 90 | } |
|---|
| 91 | |
|---|
| 92 | </pre></div><p> |
|---|
| 93 | Now some row activities will be below their lower bound so |
|---|
| 94 | we then find the variable which is cheapest in reducing the sum of |
|---|
| 95 | infeasibilities. We then repeat. This is a finite process and could be coded |
|---|
| 96 | to be faster but this is simplest. |
|---|
| 97 | </p><div class="example"><a id="id2984998"/><p class="title"><b>Example 3.6. Create feasible new solution</b></p><pre class="programlisting"> |
|---|
| 98 | |
|---|
| 99 | while (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 | } |
|---|
| 120 | if (sum>0.0) { |
|---|
| 121 | double ratio = (cost/sum)*(1.0+0.1*CoinDrand48()); |
|---|
| 122 | if (ratio<bestRatio) { |
|---|
| 123 | bestRatio=ratio; |
|---|
| 124 | bestColumn=iColumn; |
|---|
| 125 | } |
|---|
| 126 | } |
|---|
| 127 | } |
|---|
| 128 | } |
|---|
| 129 | if (bestColumn<0) |
|---|
| 130 | break; // we have finished |
|---|
| 131 | // Increase chosen column |
|---|
| 132 | newSolution[bestColumn] += 1.0; |
|---|
| 133 | double cost = direction * objective[bestColumn]; |
|---|
| 134 | newSolutionValue += cost; |
|---|
| 135 | for (CoinBigIndex j=columnStart[bestColumn]; |
|---|
| 136 | j<columnStart[bestColumn]+columnLength[bestColumn];j++) { |
|---|
| 137 | int iRow = row[j]; |
|---|
| 138 | rowActivity[iRow] += element[j]; |
|---|
| 139 | } |
|---|
| 140 | } |
|---|
| 141 | |
|---|
| 142 | </pre></div><p> |
|---|
| 143 | We have finished so now we need to see if solution is better and doublecheck |
|---|
| 144 | we are feasible. |
|---|
| 145 | </p><div class="example"><a id="id2985096"/><p class="title"><b>Example 3.7. Check good solution</b></p><pre class="programlisting"> |
|---|
| 146 | |
|---|
| 147 | returnCode=0; // 0 means no good solution |
|---|
| 148 | if (newSolutionValue<solutionValue) { |
|---|
| 149 | // check feasible |
|---|
| 150 | memset(rowActivity,0,numberRows*sizeof(double)); |
|---|
| 151 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
|---|
| 152 | CoinBigIndex j; |
|---|
| 153 | double value = newSolution[iColumn]; |
|---|
| 154 | if (value) { |
|---|
| 155 | for (j=columnStart[iColumn]; |
|---|
| 156 | j<columnStart[iColumn]+columnLength[iColumn];j++) { |
|---|
| 157 | int iRow=row[j]; |
|---|
| 158 | rowActivity[iRow] += value*element[j]; |
|---|
| 159 | } |
|---|
| 160 | } |
|---|
| 161 | } |
|---|
| 162 | // check was approximately feasible |
|---|
| 163 | bool feasible=true; |
|---|
| 164 | for (iRow=0;iRow<numberRows;iRow++) { |
|---|
| 165 | if(rowActivity[iRow]<rowLower[iRow]) { |
|---|
| 166 | if (rowActivity[iRow]<rowLower[iRow]-10.0*primalTolerance) |
|---|
| 167 | feasible = false; |
|---|
| 168 | } |
|---|
| 169 | } |
|---|
| 170 | if (feasible) { |
|---|
| 171 | // new solution |
|---|
| 172 | memcpy(betterSolution,newSolution,numberColumns*sizeof(double)); |
|---|
| 173 | solutionValue = newSolutionValue; |
|---|
| 174 | // We have good solution |
|---|
| 175 | returnCode=1; |
|---|
| 176 | } |
|---|
| 177 | } |
|---|
| 178 | |
|---|
| 179 | </pre></div></div><div class="navfooter"><hr/><table width="100%" summary="Navigation footer"><tr><td width="40%" align="left"><a accesskey="p" href="ch03.html">Prev</a> </td><td width="20%" align="center"><a accesskey="u" href="ch03.html">Up</a></td><td width="40%" align="right"> <a accesskey="n" href="ch04.html">Next</a></td></tr><tr><td width="40%" align="left" valign="top">Chapter 3. |
|---|
| 180 | Selecting a Node in the Search Tree |
|---|
| 181 | </td><td width="20%" align="center"><a accesskey="h" href="index.html">Home</a></td><td width="40%" align="right" valign="top"> Chapter 4. |
|---|
| 182 | Branching |
|---|
| 183 | </td></tr></table></div></body></html> |
|---|
