[557] | 1 | <?xml version="1.0" encoding="ISO-8859-1" standalone="no"?> |
---|
| 2 | <html xmlns="http://www.w3.org/1999/xhtml"><head><title>Chapter 7. |
---|
| 3 | Advanced Solver Uses |
---|
| 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="ch06.html" title="Chapter 6. Cutting planes"/><link rel="next" href="ch07s02.html" title="Quadratic MIP"/></head><body><div class="navheader"><table width="100%" summary="Navigation header"><tr><th colspan="3" align="center">Chapter 7. |
---|
| 5 | Advanced Solver Uses |
---|
| 6 | </th></tr><tr><td width="20%" align="left"><a accesskey="p" href="ch06.html">Prev</a> </td><th width="60%" align="center"> </th><td width="20%" align="right"> <a accesskey="n" href="ch07s02.html">Next</a></td></tr></table><hr/></div><div class="chapter" lang="en"><div class="titlepage"><div><div><h2 class="title"><a id="SolverChap"/>Chapter 7. |
---|
| 7 | Advanced Solver Uses |
---|
| 8 | </h2></div></div><div/></div><div class="toc"><p><b>Table of Contents</b></p><dl><dt><a href="ch07.html#solver">Creating a Solver via Inheritance</a></dt><dt><a href="ch07s02.html">Quadratic MIP</a></dt></dl></div><div class="section" lang="en"><div class="titlepage"><div><div><h2 class="title" style="clear: both"><a id="solver"/>Creating a Solver via Inheritance</h2></div></div><div/></div><p> |
---|
| 9 | CBC uses a generic <tt class="classname">OsiSolverInterface</tt> and its <tt class="function">resolve</tt> capability. |
---|
| 10 | This does not give much flexibility so advanced users can inherit from their interface |
---|
| 11 | of choice. This section illustrates how to implement a specialized solver for a long thin problem, e.g., fast0507 again. As with the other examples in the Guide, the sample code is not guaranteed to be the fastest way to solve the problem. The main purpose of the example is to illustrate techniques. The full source is in <tt class="filename">CbcSolver2.hpp</tt> and <tt class="filename">CbcSolver2.cpp</tt> located in the CBC Samples directory, see |
---|
| 12 | <a href="ch08.html" title="Chapter 8. More Samples ">Chapter 8, <i> |
---|
| 13 | More Samples |
---|
| 14 | </i></a>. |
---|
| 15 | </p><p> |
---|
| 16 | The method <tt class="function">initialSolve</tt> is called a few times in CBC, and provides a convenient starting point. The <tt class="varname">modelPtr_</tt> derives from <tt class="classname">OsiClpSolverInterface</tt>. |
---|
| 17 | </p><div class="example"><a id="initialSolve"/><p class="title"><b>Example 7.1. <tt class="function">initialSolve()</tt></b></p><pre class="programlisting"> |
---|
| 18 | |
---|
| 19 | // modelPtr_ is of type ClpSimplex * |
---|
| 20 | modelPtr_->setLogLevel(1); // switch on a bit of printout |
---|
| 21 | modelPtr_->scaling(0); // We don't want scaling for fast0507 |
---|
| 22 | setBasis(basis_,modelPtr_); // Put basis into ClpSimplex |
---|
| 23 | // Do long thin by sprint |
---|
| 24 | ClpSolve options; |
---|
| 25 | options.setSolveType(ClpSolve::usePrimalorSprint); |
---|
| 26 | options.setPresolveType(ClpSolve::presolveOff); |
---|
| 27 | options.setSpecialOption(1,3,15); // Do 15 sprint iterations |
---|
| 28 | modelPtr_->initialSolve(options); // solve problem |
---|
| 29 | basis_ = getBasis(modelPtr_); // save basis |
---|
| 30 | modelPtr_->setLogLevel(0); // switch off printout |
---|
| 31 | |
---|
| 32 | </pre></div><p> |
---|
| 33 | The <tt class="function">resolve()</tt> method is more complicated than <tt class="function">initialSolve()</tt>. The main pieces of data are a counter <tt class="varname">count_</tt> which is incremented each solve and an integer array <tt class="varname">node_</tt> which stores the last time |
---|
| 34 | a variable was active in a solution. For the first few solves, the normal Dual Simplex is called and |
---|
| 35 | <tt class="varname">node_</tt> array is updated. |
---|
[559] | 36 | </p><div class="example"><a id="id3000834"/><p class="title"><b>Example 7.2. First Few Solves</b></p><pre class="programlisting"> |
---|
[557] | 37 | |
---|
| 38 | if (count_<10) { |
---|
| 39 | OsiClpSolverInterface::resolve(); // Normal resolve |
---|
| 40 | if (modelPtr_->status()==0) { |
---|
| 41 | count_++; // feasible - save any nonzero or basic |
---|
| 42 | const double * solution = modelPtr_->primalColumnSolution(); |
---|
| 43 | for (int i=0;i<numberColumns;i++) { |
---|
| 44 | if (solution[i]>1.0e-6||modelPtr_->getStatus(i)==ClpSimplex::basic) { |
---|
| 45 | node_[i]=CoinMax(count_,node_[i]); |
---|
| 46 | howMany_[i]++; |
---|
| 47 | } |
---|
| 48 | } |
---|
| 49 | } else { |
---|
| 50 | printf("infeasible early on\n"); |
---|
| 51 | } |
---|
| 52 | } |
---|
| 53 | |
---|
| 54 | </pre></div><p> |
---|
| 55 | After the first few solves, only those variables which took part in a solution in the last so many |
---|
| 56 | solves are used. As fast0507 is a set covering problem, any rows which are already covered can be taken out. |
---|
[559] | 57 | </p><div class="example"><a id="id3000862"/><p class="title"><b>Example 7.3. Create Small Sub-Problem</b></p><pre class="programlisting"> |
---|
[557] | 58 | |
---|
| 59 | int * whichRow = new int[numberRows]; // Array to say which rows used |
---|
| 60 | int * whichColumn = new int [numberColumns]; // Array to say which columns used |
---|
| 61 | int i; |
---|
| 62 | const double * lower = modelPtr_->columnLower(); |
---|
| 63 | const double * upper = modelPtr_->columnUpper(); |
---|
| 64 | setBasis(basis_,modelPtr_); // Set basis |
---|
| 65 | int nNewCol=0; // Number of columns in small model |
---|
| 66 | // Column copy of matrix |
---|
| 67 | const double * element = modelPtr_->matrix()->getElements(); |
---|
| 68 | const int * row = modelPtr_->matrix()->getIndices(); |
---|
| 69 | const CoinBigIndex * columnStart = modelPtr_->matrix()->getVectorStarts(); |
---|
| 70 | const int * columnLength = modelPtr_->matrix()->getVectorLengths(); |
---|
| 71 | |
---|
| 72 | int * rowActivity = new int[numberRows]; // Number of columns with entries in each row |
---|
| 73 | memset(rowActivity,0,numberRows*sizeof(int)); |
---|
| 74 | int * rowActivity2 = new int[numberRows]; // Lower bound on row activity for each row |
---|
| 75 | memset(rowActivity2,0,numberRows*sizeof(int)); |
---|
| 76 | char * mark = (char *) modelPtr_->dualColumnSolution(); // Get some space to mark columns |
---|
| 77 | memset(mark,0,numberColumns); |
---|
| 78 | for (i=0;i<numberColumns;i++) { |
---|
| 79 | bool choose = (node_[i]>count_-memory_&&node_[i]>0); // Choose if used recently |
---|
| 80 | // Take if used recently or active in some sense |
---|
| 81 | if ((choose&&upper[i]) |
---|
| 82 | ||(modelPtr_->getStatus(i)!=ClpSimplex::atLowerBound&& |
---|
| 83 | modelPtr_->getStatus(i)!=ClpSimplex::isFixed) |
---|
| 84 | ||lower[i]>0.0) { |
---|
| 85 | mark[i]=1; // mark as used |
---|
| 86 | whichColumn[nNewCol++]=i; // add to list |
---|
| 87 | CoinBigIndex j; |
---|
| 88 | double value = upper[i]; |
---|
| 89 | if (value) { |
---|
| 90 | for (j=columnStart[i]; |
---|
| 91 | j<columnStart[i]+columnLength[i];j++) { |
---|
| 92 | int iRow=row[j]; |
---|
| 93 | assert (element[j]==1.0); |
---|
| 94 | rowActivity[iRow] ++; // This variable can cover this row |
---|
| 95 | } |
---|
| 96 | if (lower[i]>0.0) { |
---|
| 97 | for (j=columnStart[i]; |
---|
| 98 | j<columnStart[i]+columnLength[i];j++) { |
---|
| 99 | int iRow=row[j]; |
---|
| 100 | rowActivity2[iRow] ++; // This row redundant |
---|
| 101 | } |
---|
| 102 | } |
---|
| 103 | } |
---|
| 104 | } |
---|
| 105 | } |
---|
| 106 | int nOK=0; // Use to count rows which can be covered |
---|
| 107 | int nNewRow=0; // Use to make list of rows needed |
---|
| 108 | for (i=0;i<numberRows;i++) { |
---|
| 109 | if (rowActivity[i]) |
---|
| 110 | nOK++; |
---|
| 111 | if (!rowActivity2[i]) |
---|
| 112 | whichRow[nNewRow++]=i; // not satisfied |
---|
| 113 | else |
---|
| 114 | modelPtr_->setRowStatus(i,ClpSimplex::basic); // make slack basic |
---|
| 115 | } |
---|
| 116 | if (nOK<numberRows) { |
---|
| 117 | // The variables we have do not cover rows - see if we can find any that do |
---|
| 118 | for (i=0;i<numberColumns;i++) { |
---|
| 119 | if (!mark[i]&&upper[i]) { |
---|
| 120 | CoinBigIndex j; |
---|
| 121 | int good=0; |
---|
| 122 | for (j=columnStart[i]; |
---|
| 123 | j<columnStart[i]+columnLength[i];j++) { |
---|
| 124 | int iRow=row[j]; |
---|
| 125 | if (!rowActivity[iRow]) { |
---|
| 126 | rowActivity[iRow] ++; |
---|
| 127 | good++; |
---|
| 128 | } |
---|
| 129 | } |
---|
| 130 | if (good) { |
---|
| 131 | nOK+=good; // This covers - put in list |
---|
| 132 | whichColumn[nNewCol++]=i; |
---|
| 133 | } |
---|
| 134 | } |
---|
| 135 | } |
---|
| 136 | } |
---|
| 137 | delete [] rowActivity; |
---|
| 138 | delete [] rowActivity2; |
---|
| 139 | if (nOK<numberRows) { |
---|
| 140 | // By inspection the problem is infeasible - no need to solve |
---|
| 141 | modelPtr_->setProblemStatus(1); |
---|
| 142 | delete [] whichRow; |
---|
| 143 | delete [] whichColumn; |
---|
| 144 | printf("infeasible by inspection\n"); |
---|
| 145 | return; |
---|
| 146 | } |
---|
| 147 | // Now make up a small model with the right rows and columns |
---|
| 148 | ClpSimplex * temp = new ClpSimplex(modelPtr_,nNewRow,whichRow,nNewCol,whichColumn); |
---|
| 149 | |
---|
| 150 | </pre></div><p> |
---|
| 151 | If the variables cover the rows, then the problem is feasible (no cuts are being used). (If the rows |
---|
| 152 | were equality constraints, then this might not be the case. More work would be needed.) After the solution to the subproblem, the reduced costs of the full problem are checked. If the reduced cost of any variable not in the subproblem is negative, the code goes back to the full problem and cleans up with Primal Simplex. |
---|
[559] | 153 | </p><div class="example"><a id="id3000908"/><p class="title"><b>Example 7.4. Check Optimal Solution</b></p><pre class="programlisting"> |
---|
[557] | 154 | |
---|
| 155 | temp->setDualObjectiveLimit(1.0e50); // Switch off dual cutoff as problem is restricted |
---|
| 156 | temp->dual(); // solve |
---|
| 157 | double * solution = modelPtr_->primalColumnSolution(); // put back solution |
---|
| 158 | const double * solution2 = temp->primalColumnSolution(); |
---|
| 159 | memset(solution,0,numberColumns*sizeof(double)); |
---|
| 160 | for (i=0;i<nNewCol;i++) { |
---|
| 161 | int iColumn = whichColumn[i]; |
---|
| 162 | solution[iColumn]=solution2[i]; |
---|
| 163 | modelPtr_->setStatus(iColumn,temp->getStatus(i)); |
---|
| 164 | } |
---|
| 165 | double * rowSolution = modelPtr_->primalRowSolution(); |
---|
| 166 | const double * rowSolution2 = temp->primalRowSolution(); |
---|
| 167 | double * dual = modelPtr_->dualRowSolution(); |
---|
| 168 | const double * dual2 = temp->dualRowSolution(); |
---|
| 169 | memset(dual,0,numberRows*sizeof(double)); |
---|
| 170 | for (i=0;i<nNewRow;i++) { |
---|
| 171 | int iRow=whichRow[i]; |
---|
| 172 | modelPtr_->setRowStatus(iRow,temp->getRowStatus(i)); |
---|
| 173 | rowSolution[iRow]=rowSolution2[i]; |
---|
| 174 | dual[iRow]=dual2[i]; |
---|
| 175 | } |
---|
| 176 | // See if optimal |
---|
| 177 | double * dj = modelPtr_->dualColumnSolution(); |
---|
| 178 | // get reduced cost for large problem |
---|
| 179 | // this assumes minimization |
---|
| 180 | memcpy(dj,modelPtr_->objective(),numberColumns*sizeof(double)); |
---|
| 181 | modelPtr_->transposeTimes(-1.0,dual,dj); |
---|
| 182 | modelPtr_->setObjectiveValue(temp->objectiveValue()); |
---|
| 183 | modelPtr_->setProblemStatus(0); |
---|
| 184 | int nBad=0; |
---|
| 185 | |
---|
| 186 | for (i=0;i<numberColumns;i++) { |
---|
| 187 | if (modelPtr_->getStatus(i)==ClpSimplex::atLowerBound |
---|
| 188 | &&upper[i]>lower[i]&&dj[i]<-1.0e-5) |
---|
| 189 | nBad++; |
---|
| 190 | } |
---|
| 191 | // If necessary clean up with primal (and save some statistics) |
---|
| 192 | if (nBad) { |
---|
| 193 | timesBad_++; |
---|
| 194 | modelPtr_->primal(1); |
---|
| 195 | iterationsBad_ += modelPtr_->numberIterations(); |
---|
| 196 | } |
---|
| 197 | |
---|
| 198 | </pre></div><p> |
---|
| 199 | The array <tt class="varname">node_</tt> is updated, as for the first few solves. To give some idea of the effect of this tactic, the problem fast0507 has 63,009 variables but the small problem never has more than 4,000 variables. In only about ten percent of solves was it necessary to resolve, and then the average number of iterations |
---|
| 200 | on full problem was less than 20. |
---|
| 201 | </p></div></div><div class="navfooter"><hr/><table width="100%" summary="Navigation footer"><tr><td width="40%" align="left"><a accesskey="p" href="ch06.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="ch07s02.html">Next</a></td></tr><tr><td width="40%" align="left" valign="top">Chapter 6. Cutting planes </td><td width="20%" align="center"><a accesskey="h" href="index.html">Home</a></td><td width="40%" align="right" valign="top"> Quadratic MIP</td></tr></table></div></body></html> |
---|