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> |
---|