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

Last change on this file since 1574 was 1574, checked in by lou, 8 years ago

Change to EPL license notice.

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