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

Last change on this file since 2496 was 2469, checked in by unxusr, 8 months ago

formatting

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