source: trunk/Cbc/examples/qmip2.cpp

Last change on this file was 2469, checked in by unxusr, 3 months ago

formatting

  • Property svn:keywords set to Author Date Id Revision
File size: 6.3 KB
RevLine 
[1574]1// $Id: qmip2.cpp 2469 2019-01-06 23:17:46Z forrest $
[394]2// Copyright (C) 2006, International Business Machines
3// Corporation and others.  All Rights Reserved.
[1574]4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
[394]6#include <cassert>
7
[1898]8#include "CoinPragma.hpp"
9
[394]10// For Branch and bound
11#include "OsiClpSolverInterface.hpp"
[2469]12#include "CbcModel.hpp"
[394]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"
[2469]28class CbcFeasibilityNoStrong : public CbcFeasibilityBase {
[394]29public:
[2469]30  // Default Constructor
31  CbcFeasibilityNoStrong() {}
[394]32
[706]33  virtual ~CbcFeasibilityNoStrong() {}
[2469]34  // Copy constructor
35  CbcFeasibilityNoStrong(const CbcFeasibilityNoStrong &rhs) {}
[394]36
[2469]37  // Assignment operator
38  CbcFeasibilityNoStrong &operator=(const CbcFeasibilityNoStrong &rhs)
39  {
40    return *this;
41  }
42
[394]43  /// Clone
[2469]44  virtual CbcFeasibilityBase *clone() const
45  {
46    return new CbcFeasibilityNoStrong();
47  }
[394]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  */
[2469]58  virtual int feasible(CbcModel *model, int mode)
59  {
60    return mode;
61  }
[394]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
[2469]85int main(int argc, const char *argv[])
[394]86{
87
88  // Define a Solver which inherits from OsiClpsolverInterface -> OsiSolverInterface
[2469]89
[394]90  OsiClpSolverInterface solver1;
91
[2469]92  int nX = 4;
[394]93
[2469]94  int nY = (nX * (nX - 1) / 2);
[394]95  // All columns
[2469]96  double *obj = new double[nX + nY];
97  double *clo = new double[nX + nY];
98  double *cup = new double[nX + nY];
[394]99  int i;
[2469]100  for (i = 0; i < nX; i++) {
101    obj[i] = -(i + 1);
102    clo[i] = 0.0;
103    cup[i] = 1.0;
[394]104  }
[2469]105  for (i = nX; i < nX + nY; i++) {
[394]106    obj[i] = 0.0;
[2469]107    clo[i] = 0.0;
108    cup[i] = 1.0;
[394]109  }
110  // Just ordinary rows
[2469]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;
[394]117  }
118  // and first row
[2469]119  rup[0] = nX / 2.0;
[394]120  // Matrix
[2469]121  int nEl = nX + nX * nX;
122  int *row = new int[nEl];
123  int *col = new int[nEl];
124  double *el = new double[nEl];
[394]125  // X
[2469]126  nEl = 0;
[394]127  // May need scale to make plausible
128  double scaleFactor = 1.0;
[2469]129  for (i = 0; i < nX; i++) {
130    row[nEl] = 0;
131    col[nEl] = i;
132    el[nEl++] = 1.0;
[394]133    // and diagonal
[2469]134    row[nEl] = i + 1;
135    col[nEl] = i;
136    double value = CoinDrand48() * scaleFactor;
[394]137    // make reasonable (so multiples of 0.000001)
138    value *= 1.0e6;
[2469]139    int iValue = (int)(value + 1.0);
[394]140    value = iValue;
141    value *= 1.0e-6;
[2469]142    el[nEl++] = value;
[394]143  }
144  // Y
145  nY = nX;
146  // And stored cuts
147  CglStored stored;
[2469]148  double cutEls[3] = { 1.0, 1.0, -1.0 };
[394]149  int cutIndices[3];
[2469]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;
[394]155      // add cut
[2469]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;
[394]160      // multiply to make ones with most negative objective violated
161      // make reasonable (so multiples of 0.000001)
[2469]162      value *= 1.0e6 + 1.0e6 * j;
163      int iValue = (int)(value + 1.0);
[394]164      value = iValue;
165      value *= 1.0e-6;
[2469]166      el[nEl++] = value;
[394]167      // and other
[2469]168      if (i != j) {
169        row[nEl] = j + 1;
170        col[nEl] = nY;
171        el[nEl++] = value;
[394]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);
[2469]179  delete[] obj;
180  delete[] clo;
181  delete[] cup;
182  delete[] rlo;
183  delete[] rup;
184  delete[] row;
185  delete[] col;
186  delete[] el;
[394]187  // Integers
[2469]188  for (i = 0; i < nX; i++)
[394]189    solver1.setInteger(i);
190  // Reduce printout
[2469]191  solver1.setHintParam(OsiDoReducePrint, true, OsiHintTry);
192  // This clones solver
[394]193  CbcModel model(solver1);
[640]194  // Add stored cuts (making sure at all depths)
[2469]195  model.addCutGenerator(&stored, 1, "Stored", true, false, false, -100, 1, -1);
[394]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.
[640]200      c) a solution found by strong branching will be ignored.
201      d) don't recompute a solution once found
[394]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);
[640]212  // Say don't recompute solution d)
213  model.setSpecialOptions(4);
[2469]214
[394]215  double time1 = CoinCpuTime();
216
217  // Do complete search
[2469]218
[394]219  model.branchAndBound();
220
[2469]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;
[394]226
227  // Print solution if finished - we can't get names from Osi!
228
[2469]229  if (!model.status() && model.getMinimizationObjValue() < 1.0e50) {
[394]230    int numberColumns = model.solver()->getNumCols();
[2469]231
[394]232    //const double * solution = model.bestSolution();
[2469]233    const double *solution = model.solver()->getColSolution();
234
[394]235    int iColumn;
[2469]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);
[394]240    }
241  }
242  return 0;
[2469]243}
Note: See TracBrowser for help on using the repository browser.