source: html/ch04s02.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: 12.1 KB
Line 
1<?xml version="1.0" encoding="ISO-8859-1" standalone="no"?>
2<html xmlns="http://www.w3.org/1999/xhtml"><head><title>Advanced use of solver</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="ch04.html" title="Chapter 4. &#10;  Branching&#10; "/><link rel="previous" href="ch04.html" title="Chapter 4. &#10;  Branching&#10; "/><link rel="next" href="ch05.html" title="Chapter 5. &#10;  Building and Modifying a Model&#10;  "/></head><body><div class="navheader"><table width="100%" summary="Navigation header"><tr><th colspan="3" align="center">Advanced use of solver</th></tr><tr><td width="20%" align="left"><a accesskey="p" href="ch04.html">Prev</a> </td><th width="60%" align="center">Chapter 4. 
3  Branching
4 </th><td width="20%" align="right"> <a accesskey="n" href="ch05.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="solver"/>Advanced use of solver</h2></div></div><div/></div><p>
5  Coin Branch and Cut uses a generic OsiSolverInterface and its <tt class="function">resolve</tt> capability.
6  This does not give much flexibility so advanced users can inherit from the interface
7  of choice.  This describes such a solver for a long thin problem e.g. fast0507 again.
8  As with all these examples it is not guaranteed that this is the fastest way to solve
9  any of these problems - they are to illustrate techniques.
10   The full code is in <tt class="filename">CbcSolver2.hpp</tt> and
11   <tt class="filename">CbcSolver2.cpp</tt>
12  (this code can be found in the CBC Samples directory, see
13  <a href="ch06.html" title="Chapter 6. &#10;More Samples&#10;">Chapter 6, <i>
14More Samples
15</i></a>).
16  <tt class="function">initialSolve</tt> is called a few times so although we will not gain much
17  this is a simpler place to start.  The example derives from OsiClpSolverInterface and the code
18  is:
19  </p><div class="example"><a id="id2985361"/><p class="title"><b>Example 4.3. initialSolve</b></p><pre class="programlisting">
20   
21  // modelPtr_ is of type ClpSimplex *
22  modelPtr_-&gt;setLogLevel(1); // switch on a bit of printout
23  modelPtr_-&gt;scaling(0); // We don't want scaling for fast0507
24  setBasis(basis_,modelPtr_); // Put basis into ClpSimplex
25  // Do long thin by sprint
26  ClpSolve options;
27  options.setSolveType(ClpSolve::usePrimalorSprint);
28  options.setPresolveType(ClpSolve::presolveOff);
29  options.setSpecialOption(1,3,15); // Do 15 sprint iterations
30  modelPtr_-&gt;initialSolve(options); // solve problem
31  basis_ = getBasis(modelPtr_); // save basis
32  modelPtr_-&gt;setLogLevel(0); // switch off printout
33     
34  </pre></div><p>
35The <tt class="function">resolve</tt> method is more complicated.  The main pieces of data are
36a counter count_ which is incremented each solve and an int array node_ which stores the last time
37a variable was active in a solution.  For the first few times normal dual is called and
38node_ array is updated.
39</p><div class="example"><a id="id2985410"/><p class="title"><b>Example 4.4. First few solves</b></p><pre class="programlisting">
40   
41  if (count_&lt;10) {
42    OsiClpSolverInterface::resolve(); // Normal resolve
43    if (modelPtr_-&gt;status()==0) {
44      count_++; // feasible - save any nonzero or basic
45      const double * solution = modelPtr_-&gt;primalColumnSolution();
46      for (int i=0;i&lt;numberColumns;i++) {
47        if (solution[i]&gt;1.0e-6||modelPtr_-&gt;getStatus(i)==ClpSimplex::basic) {
48          node_[i]=CoinMax(count_,node_[i]);
49          howMany_[i]++;
50        }
51      }
52    } else {
53      printf("infeasible early on\n");
54    }
55  }
56     
57  </pre></div><p>
58After the first few solves we only use those which took part in a solution in the last so many
59solves.  As fast0507 is a set covering problem we can also take out any rows which are
60already covered.
61  </p><div class="example"><a id="id2985437"/><p class="title"><b>Example 4.5. Create small problem</b></p><pre class="programlisting">
62   
63    int * whichRow = new int[numberRows]; // Array to say which rows used
64    int * whichColumn = new int [numberColumns]; // Array to say which columns used
65    int i;
66    const double * lower = modelPtr_-&gt;columnLower();
67    const double * upper = modelPtr_-&gt;columnUpper();
68    setBasis(basis_,modelPtr_); // Set basis
69    int nNewCol=0; // Number of columns in small model
70    // Column copy of matrix
71    const double * element = modelPtr_-&gt;matrix()-&gt;getElements();
72    const int * row = modelPtr_-&gt;matrix()-&gt;getIndices();
73    const CoinBigIndex * columnStart = modelPtr_-&gt;matrix()-&gt;getVectorStarts();
74    const int * columnLength = modelPtr_-&gt;matrix()-&gt;getVectorLengths();
75   
76    int * rowActivity = new int[numberRows]; // Number of columns with entries in each row
77    memset(rowActivity,0,numberRows*sizeof(int));
78    int * rowActivity2 = new int[numberRows]; // Lower bound on row activity for each row
79    memset(rowActivity2,0,numberRows*sizeof(int));
80    char * mark = (char *) modelPtr_-&gt;dualColumnSolution(); // Get some space to mark columns
81    memset(mark,0,numberColumns);
82    for (i=0;i&lt;numberColumns;i++) {
83      bool choose = (node_[i]&gt;count_-memory_&amp;&amp;node_[i]&gt;0); // Choose if used recently
84      // Take if used recently or active in some sense
85      if ((choose&amp;&amp;upper[i])
86          ||(modelPtr_-&gt;getStatus(i)!=ClpSimplex::atLowerBound&amp;&amp;
87             modelPtr_-&gt;getStatus(i)!=ClpSimplex::isFixed)
88          ||lower[i]&gt;0.0) {
89        mark[i]=1; // mark as used
90        whichColumn[nNewCol++]=i; // add to list
91        CoinBigIndex j;
92        double value = upper[i];
93        if (value) {
94          for (j=columnStart[i];
95               j&lt;columnStart[i]+columnLength[i];j++) {
96            int iRow=row[j];
97            assert (element[j]==1.0);
98            rowActivity[iRow] ++; // This variable can cover this row
99          }
100          if (lower[i]&gt;0.0) {
101            for (j=columnStart[i];
102                 j&lt;columnStart[i]+columnLength[i];j++) {
103              int iRow=row[j];
104              rowActivity2[iRow] ++; // This row redundant
105            }
106          }
107        }
108      }
109    }
110    int nOK=0; // Use to count rows which can be covered
111    int nNewRow=0; // Use to make list of rows needed
112    for (i=0;i&lt;numberRows;i++) {
113      if (rowActivity[i])
114        nOK++;
115      if (!rowActivity2[i])
116        whichRow[nNewRow++]=i; // not satisfied
117      else
118        modelPtr_-&gt;setRowStatus(i,ClpSimplex::basic); // make slack basic
119    }
120    if (nOK&lt;numberRows) {
121      // The variables we have do not cover rows - see if we can find any that do
122      for (i=0;i&lt;numberColumns;i++) {
123        if (!mark[i]&amp;&amp;upper[i]) {
124          CoinBigIndex j;
125          int good=0;
126          for (j=columnStart[i];
127               j&lt;columnStart[i]+columnLength[i];j++) {
128            int iRow=row[j];
129            if (!rowActivity[iRow]) {
130              rowActivity[iRow] ++;
131              good++;
132            }
133          }
134          if (good) {
135            nOK+=good; // This covers - put in list
136            whichColumn[nNewCol++]=i;
137          }
138        }
139      }
140    }
141    delete [] rowActivity;
142    delete [] rowActivity2;
143    if (nOK&lt;numberRows) {
144      // By inspection the problem is infeasible - no need to solve
145      modelPtr_-&gt;setProblemStatus(1);
146      delete [] whichRow;
147      delete [] whichColumn;
148      printf("infeasible by inspection\n");
149      return;
150    }
151    // Now make up a small model with the right rows and columns
152    ClpSimplex *  temp = new ClpSimplex(modelPtr_,nNewRow,whichRow,nNewCol,whichColumn);
153     
154  </pre></div><p>
155If the variables cover the rows then we know that the problem is feasible (We are not using
156cuts).  If the rows
157were E rows then this might not be the case and we would have to do more work.  When we have solved
158then we see if there are any negative reduced costs and if there are then we have to go to the
159full problem and use primal to clean up.
160  </p><div class="example"><a id="id2985478"/><p class="title"><b>Example 4.6. Check optimal solution</b></p><pre class="programlisting">
161   
162    temp-&gt;setDualObjectiveLimit(1.0e50); // Switch off dual cutoff as problem is restricted
163    temp-&gt;dual(); // solve
164    double * solution = modelPtr_-&gt;primalColumnSolution(); // put back solution
165    const double * solution2 = temp-&gt;primalColumnSolution();
166    memset(solution,0,numberColumns*sizeof(double));
167    for (i=0;i&lt;nNewCol;i++) {
168      int iColumn = whichColumn[i];
169      solution[iColumn]=solution2[i];
170      modelPtr_-&gt;setStatus(iColumn,temp-&gt;getStatus(i));
171    }
172    double * rowSolution = modelPtr_-&gt;primalRowSolution();
173    const double * rowSolution2 = temp-&gt;primalRowSolution();
174    double * dual = modelPtr_-&gt;dualRowSolution();
175    const double * dual2 = temp-&gt;dualRowSolution();
176    memset(dual,0,numberRows*sizeof(double));
177    for (i=0;i&lt;nNewRow;i++) {
178      int iRow=whichRow[i];
179      modelPtr_-&gt;setRowStatus(iRow,temp-&gt;getRowStatus(i));
180      rowSolution[iRow]=rowSolution2[i];
181      dual[iRow]=dual2[i];
182    }
183    // See if optimal
184    double * dj = modelPtr_-&gt;dualColumnSolution();
185    // get reduced cost for large problem
186    // this assumes minimization
187    memcpy(dj,modelPtr_-&gt;objective(),numberColumns*sizeof(double));
188    modelPtr_-&gt;transposeTimes(-1.0,dual,dj);
189    modelPtr_-&gt;setObjectiveValue(temp-&gt;objectiveValue());
190    modelPtr_-&gt;setProblemStatus(0);
191    int nBad=0;
192     
193    for (i=0;i&lt;numberColumns;i++) {
194      if (modelPtr_-&gt;getStatus(i)==ClpSimplex::atLowerBound
195          &amp;&amp;upper[i]&gt;lower[i]&amp;&amp;dj[i]&lt;-1.0e-5)
196        nBad++;
197    }
198    // If necessary claen up with primal (and save some statistics)
199    if (nBad) {
200      timesBad_++;
201      modelPtr_-&gt;primal(1);
202      iterationsBad_ += modelPtr_-&gt;numberIterations();
203    }
204     
205  </pre></div><p>
206We then update node_ array as for the first few solves.  To give some idea of the effect of this
207tactic fast0507 has 63,009 variables and the small problem never has more than 4,000 variables.
208In just over ten percent of solves did we have to resolve and then the average number of iterations
209on full problem was less than 20.
210To give another example - again only for illustrative purposes it is possible to do quadratic
211mip.  In this case we make <tt class="function">resolve</tt> the same as
212<tt class="function">initialSolve</tt>.
213   The full code is in <tt class="filename">ClpQuadInterface.hpp</tt> and
214   <tt class="filename">ClpQuadInterface.cpp</tt>
215  (this code can be found in the CBC Samples directory, see
216  <a href="ch06.html" title="Chapter 6. &#10;More Samples&#10;">Chapter 6, <i>
217More Samples
218</i></a>).
219  </p><div class="example"><a id="id2985520"/><p class="title"><b>Example 4.7. Solve a quadratic mip</b></p><pre class="programlisting">
220   
221  // save cutoff
222  double cutoff = modelPtr_-&gt;dualObjectiveLimit();
223  modelPtr_-&gt;setDualObjectiveLimit(1.0e50);
224  modelPtr_-&gt;scaling(0);
225  modelPtr_-&gt;setLogLevel(0);
226  // solve with no objective to get feasible solution
227  setBasis(basis_,modelPtr_);
228  modelPtr_-&gt;dual();
229  basis_ = getBasis(modelPtr_);
230  modelPtr_-&gt;setDualObjectiveLimit(cutoff);
231  if (modelPtr_-&gt;problemStatus())
232    return; // problem was infeasible
233  // Now pass in quadratic objective
234  ClpObjective * saveObjective  = modelPtr_-&gt;objectiveAsObject();
235  modelPtr_-&gt;setObjectivePointer(quadraticObjective_);
236  modelPtr_-&gt;primal();
237  modelPtr_-&gt;setDualObjectiveLimit(cutoff);
238  if (modelPtr_-&gt;objectiveValue()&gt;cutoff)
239    modelPtr_-&gt;setProblemStatus(1);
240  modelPtr_-&gt;setObjectivePointer(saveObjective);
241     
242  </pre></div></div><div class="navfooter"><hr/><table width="100%" summary="Navigation footer"><tr><td width="40%" align="left"><a accesskey="p" href="ch04.html">Prev</a> </td><td width="20%" align="center"><a accesskey="u" href="ch04.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 4. 
243  Branching
244  </td><td width="20%" align="center"><a accesskey="h" href="index.html">Home</a></td><td width="40%" align="right" valign="top"> Chapter 5. 
245  Building and Modifying a Model
246  </td></tr></table></div></body></html>
Note: See TracBrowser for help on using the repository browser.