source: trunk/Cbc/examples/qmip2.cpp @ 1173

Last change on this file since 1173 was 1173, checked in by forrest, 10 years ago

add $id

  • Property svn:keywords set to Id
File size: 6.1 KB
Line 
1/* $Id: qmip2.cpp 1173 2009-06-04 09:44:10Z forrest $ */
2// Copyright (C) 2006, International Business Machines
3// Corporation and others.  All Rights Reserved.
4#if defined(_MSC_VER)
5// Turn off compiler warning about long names
6#  pragma warning(disable:4786)
7#endif
8
9#include <cassert>
10
11// For Branch and bound
12#include "OsiClpSolverInterface.hpp"
13#include "CbcModel.hpp"
14#include "CbcCutGenerator.hpp"
15#include "CoinHelperFunctions.hpp"
16#include "CbcStrategy.hpp"
17
18// Need stored cuts
19
20#include "CglStored.hpp"
21
22// For saying about solution validity
23#include "OsiAuxInfo.hpp"
24
25
26// Time
27#include "CoinTime.hpp"
28// Class to disallow strong branching solutions
29#include "CbcFeasibilityBase.hpp"
30class CbcFeasibilityNoStrong : public CbcFeasibilityBase{
31public:
32  // Default Constructor
33  CbcFeasibilityNoStrong () {}
34
35  virtual ~CbcFeasibilityNoStrong() {}
36  // Copy constructor
37  CbcFeasibilityNoStrong ( const CbcFeasibilityNoStrong &rhs) {}
38   
39  // Assignment operator
40  CbcFeasibilityNoStrong & operator=( const CbcFeasibilityNoStrong& rhs)
41  { return * this;}
42
43  /// Clone
44  virtual CbcFeasibilityBase * clone() const
45  { return new CbcFeasibilityNoStrong();}
46
47  /**
48     On input mode:
49     0 - called after a solve but before any cuts
50     -1 - called after strong branching
51     Returns :
52     0 - no opinion
53     -1 pretend infeasible
54     1 pretend integer solution
55  */
56  virtual int feasible(CbcModel * model, int mode) 
57  {return mode;}
58};
59
60
61/************************************************************************
62
63This main program solves the following 0-1 problem:
64
65min -x0 - 2x1 - 3x2 - 4x3
66
67subject to
68 
69x0 + x1 + x2 + x3 <= 2
70
71and quadratic constraints with positive random numbers
72
73It does it creating extra yij variables and constraints xi + xj -1 <= yij
74and putting quadratic elements on y
75
76The extra constraints are treated as stored cuts.
77
78This is to show how to keep branching even if we have a solution
79
80************************************************************************/
81
82int main (int argc, const char *argv[])
83{
84
85  // Define a Solver which inherits from OsiClpsolverInterface -> OsiSolverInterface
86 
87  OsiClpSolverInterface solver1;
88
89  int nX=4;
90
91  int nY = (nX * (nX-1)/2);
92  // All columns
93  double * obj = new double [nX+nY];
94  double * clo = new double[nX+nY];
95  double * cup = new double[nX+nY];
96  int i;
97  for (i=0;i<nX;i++) {
98    obj[i] = -(i+1);
99    clo[i]=0.0;
100    cup[i]=1.0;
101  }
102  for (i=nX;i<nX+nY;i++) {
103    obj[i] = 0.0;
104    clo[i]=0.0;
105    cup[i]=1.0;
106  }
107  // Just ordinary rows
108  int nRow = 1+nX;
109  double * rlo = new double[nRow];
110  double * rup = new double[nRow];
111  for (i=0;i<nRow;i++) {
112    rlo[i]=-COIN_DBL_MAX;
113    rup[i]=1.0;
114  }
115  // and first row
116  rup[0]=nX/2.0;
117  // Matrix
118  int nEl = nX+nX*nX;
119  int * row = new int[nEl];
120  int * col = new int[nEl];
121  double * el = new double[nEl];
122  // X
123  nEl=0;
124  // May need scale to make plausible
125  double scaleFactor = 1.0;
126  for (i=0;i<nX;i++) {
127    row[nEl]=0;
128    col[nEl]=i;
129    el[nEl++]=1.0;
130    // and diagonal
131    row[nEl]=i+1;
132    col[nEl]=i;
133    double value = CoinDrand48()*scaleFactor;
134    // make reasonable (so multiples of 0.000001)
135    value *= 1.0e6;
136    int iValue = (int) (value+1.0);
137    value = iValue;
138    value *= 1.0e-6;
139    el[nEl++]=value;
140  }
141  // Y
142  nY = nX;
143  // And stored cuts
144  CglStored stored;
145  double cutEls[3]={1.0,1.0,-1.0};
146  int cutIndices[3];
147  for (i=0;i<nX;i++) {
148    cutIndices[0]=i;
149    for (int j=i+1;j<nX;j++) {
150      cutIndices[1]=j;
151      cutIndices[2]=nY;
152      // add cut
153      stored.addCut(-COIN_DBL_MAX,1.0,3,cutIndices,cutEls);
154      row[nEl]=i+1;
155      col[nEl]=nY;
156      double value = CoinDrand48()*scaleFactor;
157      // multiply to make ones with most negative objective violated
158      // make reasonable (so multiples of 0.000001)
159      value *= 1.0e6+1.0e6*j;
160      int iValue = (int) (value+1.0);
161      value = iValue;
162      value *= 1.0e-6;
163      el[nEl++]=value;
164      // and other
165      if (i!=j) {
166        row[nEl]=j+1;
167        col[nEl]=nY;
168        el[nEl++]=value;
169      }
170      nY++;
171    }
172  }
173  // Create model
174  CoinPackedMatrix matrix(true, row, col, el, nEl);
175  solver1.loadProblem(matrix, clo, cup, obj, rlo, rup);
176  delete [] obj;
177  delete [] clo;
178  delete [] cup;
179  delete [] rlo;
180  delete [] rup;
181  delete [] row;
182  delete [] col;
183  delete [] el;
184  // Integers
185  for (i=0;i<nX;i++)
186    solver1.setInteger(i);
187  // Reduce printout
188  solver1.setHintParam(OsiDoReducePrint,true,OsiHintTry);
189  // This clones solver
190  CbcModel model(solver1);
191  // Add stored cuts (making sure at all depths)
192  model.addCutGenerator(&stored,1,"Stored",true,false,false,-100,1,-1);
193  /*  You need the next few lines -
194      a) so that cut generator will always be called again if it generated cuts
195      b) it is known that matrix is not enough to define problem so do cuts even
196         if it looks integer feasible at continuous optimum.
197      c) a solution found by strong branching will be ignored.
198      d) don't recompute a solution once found
199  */
200  // Make sure cut generator called correctly (a)
201  model.cutGenerator(0)->setMustCallAgain(true);
202  // Say cuts needed at continuous (b)
203  OsiBabSolver oddCuts;
204  oddCuts.setSolverType(4);
205  model.passInSolverCharacteristics(&oddCuts);
206  // Say no to all solutions by strong branching (c)
207  CbcFeasibilityNoStrong noStrong;
208  model.setProblemFeasibility(noStrong);
209  // Say don't recompute solution d)
210  model.setSpecialOptions(4);
211 
212  double time1 = CoinCpuTime();
213
214  // Do complete search
215 
216  model.branchAndBound();
217
218  std::cout<<argv[1]<<" took "<<CoinCpuTime()-time1<<" seconds, "
219           <<model.getNodeCount()<<" nodes with objective "
220           <<model.getObjValue()
221           <<(!model.status() ? " Finished" : " Not finished")
222           <<std::endl;
223
224  // Print solution if finished - we can't get names from Osi!
225
226  if (!model.status()&&model.getMinimizationObjValue()<1.0e50) {
227    int numberColumns = model.solver()->getNumCols();
228   
229    //const double * solution = model.bestSolution();
230    const double * solution = model.solver()->getColSolution();
231   
232    int iColumn;
233    for (iColumn=0;iColumn<numberColumns;iColumn++) {
234      double value=solution[iColumn];
235      if (fabs(value)>1.0e-7&&model.solver()->isInteger(iColumn)) 
236        printf("Column %d has value %g\n",iColumn,value);
237    }
238  }
239  return 0;
240}   
Note: See TracBrowser for help on using the repository browser.