source: html/trunk/Cbc/ch03s02.html @ 554

Last change on this file since 554 was 554, checked in by rlh, 16 years ago

initial import of Cbc documentation

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.5 KB
Line 
1<html><head><meta http-equiv="Content-Type" content="text/html; charset=ISO-8859-1"><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. 
2  Other Classes and examples
3  "><link rel="previous" href="ch03.html" title="Chapter 3. 
4  Other Classes and examples
5  "><link rel="next" href="ch03s03.html" title="Branching"></head><body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"><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. 
6  Other Classes and examples
7  </th><td width="20%" align="right"> <a accesskey="n" href="ch03s03.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 name="heuristics"></a>CbcHeuristic - heuristic methods</h2></div></div><div></div></div><p>
8  For practical use it is very useful to be able to get a good solution reasonably fast.
9  A good bound will greatly reduce the run time and good solutions can satisfy the user
10  on very large problems where a complete search is impossible.  Heuristics are obviously
11  problem dependant although some have more general use.  I hope to increase the number.
12  At present there is only one in Cbc itself although there are others in the Samples
13  directory.  The heuristic is trying to obtain a solution to the original
14  problem so it need only consider original rows and does not have to use the
15  current bounds.
16  One to use a greedy heuristic designed for use in the miplib problem
17  fast0507 will be developed later in this section.
18  Cbc provides a abstract base class CbcHeuristic and a rounding heuristic in Cbc.
19  </p><p>
20  This describes how to build a greedy heuristic for a set covering problem.
21   A more general version is in <tt class="filename">CbcHeuristicGreedy.hpp</tt> and
22   <tt class="filename">CbcHeuristicGreedy.cpp</tt>
23  (this code can be found in the CBC Samples directory, see
24  <a href="ch05.html" title="Chapter 5. 
25More Samples
26">Chapter 5, <i>
27More Samples
28</i></a>).
29
30  The heuristic we will code will leave all variables which are at one at this node of the
31  tree to that value and will
32  initially set all others to zero.  We then sort all variables in order of their cost
33  divided by the number of entries in rows which are not yet covered.  We may randomize that
34  value a bit so that ties will be broken in different ways on different runs of the heuristic.
35  We then choose the best one and set it to one and repeat the exercise. Because this is
36  a set covering problem &#8805; we are guaranteed to find a solution (not necessarily a
37  better one though).  We could
38  improve the speed by just redoing those affected but in this text we will keep it simple.
39  Also if all elements are 1.0 then we could do faster.
40  The key CbcHeuristic method is <tt class="function">int solution(double &amp; solutionValue,
41                                              double * betterSolution)</tt>
42  which returns 0 if no solution found and 1 if found when it fills in the objective value
43  and primal solution.  The actual code in <tt class="filename">CbcHeuristicGreedy.cpp</tt> is
44  a little more complicated but this will work and gives the basic idea.  For instance
45  the code here assumes all variables are integer.
46  The important bit of data is a copy of the matrix (stored by column)
47  before any cuts have been made.  The data used are bounds, objective and the matrix
48  plus two work arrays.
49  </p><div class="example"><a name="id2903012"></a><p class="title"><b>Example 3.4. Data</b></p><pre class="programlisting">
50   
51  OsiSolverInterface * solver = model_-&gt;solver(); // Get solver from CbcModel
52  const double * columnLower = solver-&gt;getColLower(); // Column Bounds
53  const double * columnUpper = solver-&gt;getColUpper();
54  const double * rowLower = solver-&gt;getRowLower(); // We know we only need lower bounds
55  const double * solution = solver-&gt;getColSolution();
56  const double * objective = solver-&gt;getObjCoefficients(); // In code we also use min/max
57  double integerTolerance = model_-&gt;getDblParam(CbcModel::CbcIntegerTolerance);
58  double primalTolerance;
59  solver-&gt;getDblParam(OsiPrimalTolerance,primalTolerance);
60  int numberRows = originalNumberRows_; // This is number of rows when matrix was passed in
61  // Column copy of matrix (before cuts)
62  const double * element = matrix_.getElements();
63  const int * row = matrix_.getIndices();
64  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
65  const int * columnLength = matrix_.getVectorLengths();
66
67  // Get solution array for heuristic solution
68  int numberColumns = solver-&gt;getNumCols();
69  double * newSolution = new double [numberColumns];
70  // And to sum row activities
71  double * rowActivity = new double[numberRows];
72     
73  </pre></div><p>
74Then we initialize newSolution as rounded down solution.
75</p><div class="example"><a name="id2903038"></a><p class="title"><b>Example 3.5. initialize newSolution</b></p><pre class="programlisting">
76   
77  for (iColumn=0;iColumn&lt;numberColumns;iColumn++) {
78    CoinBigIndex j;
79    double value = solution[iColumn];
80    // Round down integer
81    if (fabs(floor(value+0.5)-value)&lt;integerTolerance)
82      value=floor(CoinMax(value+1.0e-3,columnLower[iColumn]));
83    // make sure clean
84    value = CoinMin(value,columnUpper[iColumn]);
85    value = CoinMax(value,columnLower[iColumn]);
86    newSolution[iColumn]=value;
87    if (value) {
88      double cost = direction * objective[iColumn];
89      newSolutionValue += value*cost;
90      for (j=columnStart[iColumn];
91           j&lt;columnStart[iColumn]+columnLength[iColumn];j++) {
92        int iRow=row[j];
93        rowActivity[iRow] += value*element[j];
94      }
95    }
96  }
97     
98  </pre></div><p>
99Now some row activities will be below their lower bound so
100we then find the variable which is cheapest in reducing the sum of
101infeasibilities.  We then repeat.  This is a finite process and could be coded
102to be faster but this is simplest.
103  </p><div class="example"><a name="id2903094"></a><p class="title"><b>Example 3.6. Create feasible new solution</b></p><pre class="programlisting">
104   
105  while (true) {
106    // Get column with best ratio
107    int bestColumn=-1;
108    double bestRatio=COIN_DBL_MAX;
109    for (int iColumn=0;iColumn&lt;numberColumns;iColumn++) {
110      CoinBigIndex j;
111      double value = newSolution[iColumn];
112      double cost = direction * objective[iColumn];
113      // we could use original upper rather than current
114      if (value+0.99&lt;columnUpper[iColumn]) {
115        double sum=0.0; // Compute how much we will reduce infeasibility by
116        for (j=columnStart[iColumn];
117             j&lt;columnStart[iColumn]+columnLength[iColumn];j++) {
118          int iRow=row[j];
119          double gap = rowLower[iRow]-rowActivity[iRow];
120          if (gap&gt;1.0e-7) {
121            sum += CoinMin(element[j],gap);
122          if (element[j]+rowActivity[iRow]&lt;rowLower[iRow]+1.0e-7) {
123            sum += element[j];
124          }
125        }
126        if (sum&gt;0.0) {
127          double ratio = (cost/sum)*(1.0+0.1*CoinDrand48());
128          if (ratio&lt;bestRatio) {
129            bestRatio=ratio;
130            bestColumn=iColumn;
131          }
132        }
133      }
134    }
135    if (bestColumn&lt;0)
136      break; // we have finished
137    // Increase chosen column
138    newSolution[bestColumn] += 1.0;
139    double cost = direction * objective[bestColumn];
140    newSolutionValue += cost;
141    for (CoinBigIndex j=columnStart[bestColumn];
142         j&lt;columnStart[bestColumn]+columnLength[bestColumn];j++) {
143      int iRow = row[j];
144      rowActivity[iRow] += element[j];
145    }
146  }
147     
148  </pre></div><p>
149We have finished so now we need to see if solution is better and doublecheck
150we are feasible.
151  </p><div class="example"><a name="id2903189"></a><p class="title"><b>Example 3.7. Check good solution</b></p><pre class="programlisting">
152   
153  returnCode=0; // 0 means no good solution
154  if (newSolutionValue&lt;solutionValue) {
155    // check feasible
156    memset(rowActivity,0,numberRows*sizeof(double));
157    for (iColumn=0;iColumn&lt;numberColumns;iColumn++) {
158      CoinBigIndex j;
159      double value = newSolution[iColumn];
160      if (value) {
161        for (j=columnStart[iColumn];
162             j&lt;columnStart[iColumn]+columnLength[iColumn];j++) {
163          int iRow=row[j];
164          rowActivity[iRow] += value*element[j];
165        }
166      }
167    }
168    // check was approximately feasible
169    bool feasible=true;
170    for (iRow=0;iRow&lt;numberRows;iRow++) {
171      if(rowActivity[iRow]&lt;rowLower[iRow]) {
172        if (rowActivity[iRow]&lt;rowLower[iRow]-10.0*primalTolerance)
173          feasible = false;
174      }
175    }
176    if (feasible) {
177      // new solution
178      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
179      solutionValue = newSolutionValue;
180      // We have good solution
181      returnCode=1;
182    }
183  }
184     
185  </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="ch03s03.html">Next</a></td></tr><tr><td width="40%" align="left" valign="top">Chapter 3. 
186  Other Classes and examples
187   </td><td width="20%" align="center"><a accesskey="h" href="index.html">Home</a></td><td width="40%" align="right" valign="top"> Branching</td></tr></table></div></body></html>
Note: See TracBrowser for help on using the repository browser.