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

Last change on this file since 1898 was 1898, checked in by stefan, 6 years ago

fixup examples

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