source: html/ch04.html @ 2246

Last change on this file since 2246 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: 10.5 KB
Line 
1<?xml version="1.0" encoding="ISO-8859-1" standalone="no"?>
2<html xmlns="http://www.w3.org/1999/xhtml"><head><title>Chapter 4. 
3  Getting Good Bounds in CBC
4  </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="index.html" title="CBC User Guide"/><link rel="previous" href="ch03.html" title="Chapter 3. &#10;  Selecting the Next Node in the Search Tree&#10;  "/><link rel="next" href="ch05.html" title="Chapter 5. &#10;  Branching&#10; "/></head><body><div class="navheader"><table width="100%" summary="Navigation header"><tr><th colspan="3" align="center">Chapter 4. 
5  Getting Good Bounds in CBC
6  </th></tr><tr><td width="20%" align="left"><a accesskey="p" href="ch03.html">Prev</a> </td><th width="60%" align="center"> </th><td width="20%" align="right"> <a accesskey="n" href="ch05.html">Next</a></td></tr></table><hr/></div><div class="chapter" lang="en"><div class="titlepage"><div><div><h2 class="title"><a id="hueristicChap"/>Chapter 4. 
7  Getting Good Bounds in CBC
8  </h2></div></div><div/></div><div class="toc"><p><b>Table of Contents</b></p><dl><dt><a href="ch04.html#heuristics">CbcHeuristic - Heuristic Methods</a></dt></dl></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>
9  In practice, it is very useful to get a good solution reasonably fast. Any MIP-feasible solution produces an upper bound, and a good bound will greatly reduce the run time. Good solutions can satisfy the user
10  on very large problems where a complete search is impossible.  Obviously, heuristics are
11  problem dependent, although some do have more general use.
12  At present there is only one heuristic in CBC itself, <tt class="classname">CbcRounding</tt>. Hopefully, the number will grow. Other heuristics are in the <tt class="filename">COIN/Cbc/Samples</tt>
13  directory.  A heuristic tries to obtain a solution to the original
14  problem so it only needs to consider the original rows and does not have to use the
15  current bounds. CBC provides an abstract base class <tt class="classname">CbcHeuristic</tt> and a rounding heuristic in CBC.
16  </p><p>
17  This chapter describes how to build a greedy heuristic for a set covering problem, e.g., the miplib problem fast0507. A more general (and efficient) version of the heuristic is in <tt class="filename">CbcHeuristicGreedy.hpp</tt> and <tt class="filename">CbcHeuristicGreedy.cpp</tt> located in the <tt class="filename">COIN/Cbc/Samples</tt> directory, see <a href="ch08.html" title="Chapter 8. &#10;More Samples&#10;">Chapter 8, <i>
18More Samples
19</i></a>.
20</p><p>
21  The greedy heuristic will leave all variables taking value one at this node of the
22  tree at value one, and will initially set all other variables to value zero. 
23  All variables are then sorted in order of their cost
24  divided by the number of entries in rows which are not yet covered. (We may randomize that
25  value a bit so that ties will be broken in different ways on different runs of the heuristic.)
26  The best one is choosen, and set to one. The process is repeated. Because this is
27  a set covering problem (i.e., all constraints are &#8805;), the heuristic is guaranteed to find a solution (but not necessarily an improved solution). The speed of the heuristic could be improved by just redoing those affected, but for illustrative purposes we will keep it simple. (The speed could also be improved if all elements are 1.0).
28</p><p>
29  The key <tt class="classname">CbcHeuristic</tt> method is <tt class="function">int solution(double &amp; solutionValue,
30                                              double * betterSolution)</tt>.
31  The <tt class="function">solution()</tt> method returns 0 if no solution found, and returns 1 if a solution is found, in which case it fills in the objective value and primal solution.  The code in <tt class="filename">CbcHeuristicGreedy.cpp</tt> is a little more complicated than this following example. 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.
32  </p><div class="example"><a id="id3000043"/><p class="title"><b>Example 4.1. Data</b></p><pre class="programlisting">
33   
34  OsiSolverInterface * solver = model_-&gt;solver(); // Get solver from CbcModel
35  const double * columnLower = solver-&gt;getColLower(); // Column Bounds
36  const double * columnUpper = solver-&gt;getColUpper();
37  const double * rowLower = solver-&gt;getRowLower(); // We know we only need lower bounds
38  const double * solution = solver-&gt;getColSolution();
39  const double * objective = solver-&gt;getObjCoefficients(); // In code we also use min/max
40  double integerTolerance = model_-&gt;getDblParam(CbcModel::CbcIntegerTolerance);
41  double primalTolerance;
42  solver-&gt;getDblParam(OsiPrimalTolerance,primalTolerance);
43  int numberRows = originalNumberRows_; // This is number of rows when matrix was passed in
44  // Column copy of matrix (before cuts)
45  const double * element = matrix_.getElements();
46  const int * row = matrix_.getIndices();
47  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
48  const int * columnLength = matrix_.getVectorLengths();
49
50  // Get solution array for heuristic solution
51  int numberColumns = solver-&gt;getNumCols();
52  double * newSolution = new double [numberColumns];
53  // And to sum row activities
54  double * rowActivity = new double[numberRows];
55     
56  </pre></div><p>
57The <tt class="varname">newSolution</tt> is then initialized to the rounded down solution.
58</p><div class="example"><a id="id3000072"/><p class="title"><b>Example 4.2. Initialize <tt class="varname">newSolution</tt></b></p><pre class="programlisting">
59   
60  for (iColumn=0;iColumn&lt;numberColumns;iColumn++) {
61    CoinBigIndex j;
62    double value = solution[iColumn];
63    // Round down integer
64    if (fabs(floor(value+0.5)-value)&lt;integerTolerance)
65      value=floor(CoinMax(value+1.0e-3,columnLower[iColumn]));
66    // make sure clean
67    value = CoinMin(value,columnUpper[iColumn]);
68    value = CoinMax(value,columnLower[iColumn]);
69    newSolution[iColumn]=value;
70    if (value) {
71      double cost = objective[iColumn];
72      newSolutionValue += value*cost;
73      for (j=columnStart[iColumn];
74           j&lt;columnStart[iColumn]+columnLength[iColumn];j++) {
75        int iRow=row[j];
76        rowActivity[iRow] += value*element[j];
77      }
78    }
79  }
80     
81  </pre></div><p>
82
83
84At this point some row activities are below their lower bound. To correct the infeasibility, the variable which is cheapest in reducing the sum of infeasibilities is found and updated, and the process repeats.  This is a finite process. (The implementation could be faster, but is kept simple for illustrative purposes.)
85  </p><div class="example"><a id="id3000153"/><p class="title"><b>Example 4.3. Create Feasible <tt class="varname">newSolution</tt> from Initial <tt class="varname">newSolution</tt></b></p><pre class="programlisting">
86   
87  while (true) {
88    // Get column with best ratio
89    int bestColumn=-1;
90    double bestRatio=COIN_DBL_MAX;
91    for (int iColumn=0;iColumn&lt;numberColumns;iColumn++) {
92      CoinBigIndex j;
93      double value = newSolution[iColumn];
94      double cost = direction * objective[iColumn];
95      // we could use original upper rather than current
96      if (value+0.99&lt;columnUpper[iColumn]) {
97        double sum=0.0; // Compute how much we will reduce infeasibility by
98        for (j=columnStart[iColumn];
99             j&lt;columnStart[iColumn]+columnLength[iColumn];j++) {
100          int iRow=row[j];
101          double gap = rowLower[iRow]-rowActivity[iRow];
102          if (gap&gt;1.0e-7) {
103            sum += CoinMin(element[j],gap);
104          if (element[j]+rowActivity[iRow]&lt;rowLower[iRow]+1.0e-7) {
105            sum += element[j];
106          }
107        }
108        if (sum&gt;0.0) {
109          double ratio = (cost/sum)*(1.0+0.1*CoinDrand48());
110          if (ratio&lt;bestRatio) {
111            bestRatio=ratio;
112            bestColumn=iColumn;
113          }
114        }
115      }
116    }
117    if (bestColumn&lt;0)
118      break; // we have finished
119    // Increase chosen column
120    newSolution[bestColumn] += 1.0;
121    double cost = direction * objective[bestColumn];
122    newSolutionValue += cost;
123    for (CoinBigIndex j=columnStart[bestColumn];
124         j&lt;columnStart[bestColumn]+columnLength[bestColumn];j++) {
125      int iRow = row[j];
126      rowActivity[iRow] += element[j];
127    }
128  }
129     
130  </pre></div><p>
131A solution value of <tt class="varname">newSolution</tt> is compared to the best solution value. If <tt class="varname">newSolution</tt> is an improvement, its feasibility is validated. We expect <tt class="varname">newSolution</tt> to be feasible, and are trapping for unexpected numerical errors.
132  </p><div class="example"><a id="id3000197"/><p class="title"><b>Example 4.4. Check Solution Quality of <tt class="varname">newSolution</tt></b></p><pre class="programlisting">
133   
134  returnCode=0; // 0 means no good solution
135  if (newSolutionValue&lt;solutionValue) { // minimization
136    // check feasible
137    memset(rowActivity,0,numberRows*sizeof(double));
138    for (iColumn=0;iColumn&lt;numberColumns;iColumn++) {
139      CoinBigIndex j;
140      double value = newSolution[iColumn];
141      if (value) {
142        for (j=columnStart[iColumn];
143             j&lt;columnStart[iColumn]+columnLength[iColumn];j++) {
144          int iRow=row[j];
145          rowActivity[iRow] += value*element[j];
146        }
147      }
148    }
149    // check was approximately feasible
150    bool feasible=true;
151    for (iRow=0;iRow&lt;numberRows;iRow++) {
152      if(rowActivity[iRow]&lt;rowLower[iRow]) {
153        if (rowActivity[iRow]&lt;rowLower[iRow]-10.0*primalTolerance)
154          feasible = false;
155      }
156    }
157    if (feasible) {
158      // new solution
159      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
160      solutionValue = newSolutionValue;
161      // We have good solution
162      returnCode=1;
163    }
164  }
165     
166  </pre></div></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="index.html">Up</a></td><td width="40%" align="right"> <a accesskey="n" href="ch05.html">Next</a></td></tr><tr><td width="40%" align="left" valign="top">Chapter 3. 
167  Selecting the Next Node in the Search Tree
168   </td><td width="20%" align="center"><a accesskey="h" href="index.html">Home</a></td><td width="40%" align="right" valign="top"> Chapter 5. 
169  Branching
170 </td></tr></table></div></body></html>
Note: See TracBrowser for help on using the repository browser.