source: html/ch03s02.html @ 2261

Last change on this file since 2261 was 561, checked in by ladanyi, 12 years ago

Cbc static pages moved

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.6 KB
Line 
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. &#10;  Selecting a Node in the Search Tree&#10;  "/><link rel="previous" href="ch03.html" title="Chapter 3. &#10;  Selecting a Node in the Search Tree&#10;  "/><link rel="next" href="ch04.html" title="Chapter 4. &#10;  Branching&#10; "/></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. &#10;More Samples&#10;">Chapter 6, <i>
21More 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 &#8805; 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 &amp; 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_-&gt;solver(); // Get solver from CbcModel
46  const double * columnLower = solver-&gt;getColLower(); // Column Bounds
47  const double * columnUpper = solver-&gt;getColUpper();
48  const double * rowLower = solver-&gt;getRowLower(); // We know we only need lower bounds
49  const double * solution = solver-&gt;getColSolution();
50  const double * objective = solver-&gt;getObjCoefficients(); // In code we also use min/max
51  double integerTolerance = model_-&gt;getDblParam(CbcModel::CbcIntegerTolerance);
52  double primalTolerance;
53  solver-&gt;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-&gt;getNumCols();
63  double * newSolution = new double [numberColumns];
64  // And to sum row activities
65  double * rowActivity = new double[numberRows];
66     
67  </pre></div><p>
68Then 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&lt;numberColumns;iColumn++) {
72    CoinBigIndex j;
73    double value = solution[iColumn];
74    // Round down integer
75    if (fabs(floor(value+0.5)-value)&lt;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&lt;columnStart[iColumn]+columnLength[iColumn];j++) {
86        int iRow=row[j];
87        rowActivity[iRow] += value*element[j];
88      }
89    }
90  }
91     
92  </pre></div><p>
93Now some row activities will be below their lower bound so
94we then find the variable which is cheapest in reducing the sum of
95infeasibilities.  We then repeat.  This is a finite process and could be coded
96to 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&lt;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&lt;columnUpper[iColumn]) {
109        double sum=0.0; // Compute how much we will reduce infeasibility by
110        for (j=columnStart[iColumn];
111             j&lt;columnStart[iColumn]+columnLength[iColumn];j++) {
112          int iRow=row[j];
113          double gap = rowLower[iRow]-rowActivity[iRow];
114          if (gap&gt;1.0e-7) {
115            sum += CoinMin(element[j],gap);
116          if (element[j]+rowActivity[iRow]&lt;rowLower[iRow]+1.0e-7) {
117            sum += element[j];
118          }
119        }
120        if (sum&gt;0.0) {
121          double ratio = (cost/sum)*(1.0+0.1*CoinDrand48());
122          if (ratio&lt;bestRatio) {
123            bestRatio=ratio;
124            bestColumn=iColumn;
125          }
126        }
127      }
128    }
129    if (bestColumn&lt;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&lt;columnStart[bestColumn]+columnLength[bestColumn];j++) {
137      int iRow = row[j];
138      rowActivity[iRow] += element[j];
139    }
140  }
141     
142  </pre></div><p>
143We have finished so now we need to see if solution is better and doublecheck
144we 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&lt;solutionValue) {
149    // check feasible
150    memset(rowActivity,0,numberRows*sizeof(double));
151    for (iColumn=0;iColumn&lt;numberColumns;iColumn++) {
152      CoinBigIndex j;
153      double value = newSolution[iColumn];
154      if (value) {
155        for (j=columnStart[iColumn];
156             j&lt;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&lt;numberRows;iRow++) {
165      if(rowActivity[iRow]&lt;rowLower[iRow]) {
166        if (rowActivity[iRow]&lt;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>
Note: See TracBrowser for help on using the repository browser.