source: trunk/Bonmin/src/CbcBonmin/CbcBonmin.cpp @ 1

Last change on this file since 1 was 1, checked in by andreasw, 13 years ago

imported initial code

  • Property svn:eol-style set to native
  • Property svn:keywords set to "Author Date Id Revision"
File size: 12.9 KB
Line 
1// (C) Copyright Carnegie Mellon University 2006
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// Authors :
6// Pierre Bonami, Carnegie Mellon University,
7//
8// Date : 03/15/2006
9
10#include "BonminConfig.h"
11
12#include "CbcBonmin.hpp"
13#include "CbcModel.hpp"
14#include "CbcBranchActual.hpp"
15#include "CbcCutGenerator.hpp"
16#include "IpoptInterface.hpp"
17#include "OsiClpSolverInterface.hpp"
18
19#include "BonminCbcLpStrategy.hpp"
20#include "BonminCbcNlpStrategy.hpp"
21
22#ifdef COIN_HAS_CPX
23#include "OsiCpxSolverInterface.hpp"
24#endif
25
26
27//MILP machinery
28#include "CglGomory.hpp"
29#include "CglProbing.hpp"
30#include "CglKnapsackCover.hpp"
31#include "CglOddHole.hpp"
32#include "CglClique.hpp"
33#include "CglFlowCover.hpp"
34#include "CglMixedIntegerRounding.hpp"
35#include "CglTwomir.hpp"
36#include "CglPreProcess.hpp"
37#include "CbcCompareUser.hpp"
38#include "CbcCompareActual.hpp"
39#include "CbcBranchUser.hpp"
40
41
42//OA machinery
43#include "IpCbcDummyHeuristic.hpp"
44#include "IpCbcOACutGenerator2.hpp"
45#include "IpCbcOACutGenerator.hpp"
46
47#include "BonminCbcParam.hpp"
48
49#include "OsiAuxInfo.hpp"
50
51
52static CbcModel * currentBranchModel = NULL;
53static IpCbcOACutGenerator2 * currentOA = NULL;
54CbcModel * OAModel = NULL;
55
56#include "CoinSignal.hpp"
57
58extern "C" {
59   static void signal_handler(int whichSignal)
60   {
61      if (currentBranchModel!=NULL) 
62         currentBranchModel->setMaximumNodes(0); // stop at next node
63      if (OAModel!=NULL) 
64         OAModel->setMaximumNodes(0); // stop at next node
65      if (currentOA!=NULL) 
66         currentOA->setMaxLocalSearchTime(0.); // stop at next node
67      return;
68   }
69}
70
71
72/** Constructor.*/
73BonminBB::BonminBB():
74    bestSolution_(NULL),
75    mipStatus_(),
76    bestObj_(1e200),
77    bestBound_(-1e200),
78    continuousRelaxation_(-DBL_MAX),
79    numNodes_(0),
80    mipIterationCount_(0)
81{}
82
83/** Destructor.*/
84BonminBB::~BonminBB()
85{
86  if(bestSolution_) delete [] bestSolution_;
87  bestSolution_ = NULL;
88}
89/** Perform a branch-and-bound on given IpoptInterface using passed parameters.*/
90void
91BonminBB::branchAndBound(IpoptInterface &nlpSolver,
92    const BonminCbcParam &par)
93{
94
95  //Now set-up b&b
96  OsiSolverInterface * si;
97
98
99  nlpSolver.messageHandler()->setLogLevel(par.nlpLogLevel);
100
101  if (par.algo > 0) //OA based
102  {
103    si = new OsiClpSolverInterface;
104    nlpSolver.extractLinearRelaxation(*si);
105    // say bound dubious, does cuts at solution
106    OsiBabSolver * extraStuff = new OsiBabSolver(3);
107    si->setAuxiliaryInfo(extraStuff);
108    delete extraStuff;
109  }
110  else {
111    si = &nlpSolver;
112    nlpSolver.ignoreFailures();
113    OsiBabSolver * extraStuff = new OsiBabSolver(2);
114    si->setAuxiliaryInfo(extraStuff);
115    delete extraStuff;
116  }
117  CbcModel model(*si);
118
119
120
121  if(par.algo==0)//Switch off some feasibility checks and switch on specific nodes info
122  {
123    int specOpt = model.specialOptions();
124    specOpt = 16;
125    model.setSpecialOptions(specOpt);
126    BonminCbcNlpStrategy strat(par.maxFailures, par.maxInfeasible, par.failureBehavior);
127    model.setStrategy(strat);
128  }
129
130  //Setup likely milp cut generator
131  CglGomory miGGen;
132
133  CglProbing probGen;
134  probGen.setUsingObjective(true);
135  probGen.setMaxPass(3);
136  probGen.setMaxProbe(100);
137  probGen.setMaxLook(50);
138
139  CglKnapsackCover knapsackGen;
140  CglMixedIntegerRounding mixedGen;
141
142  //Setup OA generators
143
144  //Resolution of nlp relaxations
145  IpCbcOACutGenerator oaGen(&nlpSolver);
146  oaGen.setMaxDepth(100000);
147  oaGen.setLogLevel(par.oaLogLevel);
148
149
150  //Outer approximation iterations
151  OsiSolverInterface * localSearchSolver=NULL;
152  if(par.milpSubSolver <= 1)/* use cbc */
153  {
154    localSearchSolver = model.solver();
155    //localSearchSolver->messageHandler()->setLogLevel(0);
156  }
157  else if (par.milpSubSolver ==2) /* try to use cplex */
158  {
159#ifdef COIN_HAS_CPX
160    OsiCpxSolverInterface * cpxSolver = new OsiCpxSolverInterface;
161    localSearchSolver = cpxSolver;
162    nlpSolver.extractLinearRelaxation(*localSearchSolver);
163#else
164
165    std::cerr<<"You have set an option to use CPLEX as the milp subsolver in oa decomposition."<<std::endl
166    <<"but apparently CPLEX is not configured to be used in bonmin, see the manual for configuring CPLEX"<<std::endl;
167    throw -1;
168#endif
169
170  }
171  CbcStrategy * strategy = NULL;
172  if(par.milpSubSolver == 1) {
173    strategy = new CbcOaStrategy(par.milpSubSolver_migFreq,
174        par.milpSubSolver_probFreq,
175        par.milpSubSolver_mirFreq,
176        par.milpSubSolver_coverFreq,
177        par.milpSubSolver_minReliability,
178        par.milpSubSolver_numberStrong,
179        par.milpSubSolver_nodeSelection,
180        par.intTol,
181        par.milpLogLevel
182                                );
183  }
184  IpCbcOACutGenerator2 oaDec(&nlpSolver, localSearchSolver, strategy, par.cutoffDecr, par.intTol, 0,1);
185  if(par.algo>0) {
186    oaDec.setLocalSearchNodeLimit(1000000);
187    oaDec.setMaxLocalSearch(100000);
188    oaDec.setMaxLocalSearchPerNode(10000);
189    oaDec. setMaxLocalSearchTime(min(par.maxTime,par.oaDecMaxTime));
190    oaDec.setLogLevel(par.oaLogLevel);
191    oaDec.setLogFrequency(par.oaLogFrequency);
192    oaDec.setSubMilpLogLevel(par.milpLogLevel);
193  }
194  //Setup solver for checking validity of integral solutions
195  IpCbcOACutGenerator2 feasCheck(&nlpSolver, model.solver(),
196      NULL,
197      par.cutoffDecr, par.intTol,
198      0, 0);
199  if(par.algo>0) {
200    feasCheck.setLocalSearchNodeLimit(0);
201    feasCheck.setMaxLocalSearch(0);
202    feasCheck.setMaxLocalSearchPerNode(100000);
203  }
204  IpCbcDummyHeuristic oaHeu(model, &nlpSolver);
205
206  if(par.algo>0) {
207    int numGen = 0;
208    if(par.nlpSolveFrequency != 0)
209      {
210        model.addCutGenerator(&oaGen,par.nlpSolveFrequency,"Outer Approximation Supporting Hyperplanes for NLP optimum");
211        numGen++;
212      }
213    if(par.migFreq != 0)
214      {
215        model.addCutGenerator(&miGGen,par.migFreq,"GMI");
216        numGen++;
217      }
218    if(par.probFreq != 0)
219      {
220        model.addCutGenerator(&probGen,par.probFreq,"Probing");
221        numGen++;
222      }
223    if(par.coverFreq != 0)
224      {
225        model.addCutGenerator(&knapsackGen,par.coverFreq,"covers");
226        numGen++;
227      }
228    if(par.mirFreq != 0)
229      {
230        model.addCutGenerator(&mixedGen,par.mirFreq,"MIR");
231        numGen++;
232      }
233    if(par.oaDecMaxTime>0.)
234      {
235        model.addCutGenerator(&oaDec,1,"Outer Approximation local enumerator");     
236        IpCbcOACutGenerator2 * oaDecCopy = dynamic_cast<IpCbcOACutGenerator2 *>
237          (model.cutGenerators()[numGen]->generator());
238        assert(oaDecCopy);
239        currentOA = oaDecCopy;
240          numGen++;
241      }
242    model.addCutGenerator(&feasCheck,1,"Outer Approximation feasibility checker",false,true);
243    numGen++;
244
245    model.addHeuristic(&oaHeu);
246  }
247
248  //Set true branch-and-bound parameters
249  model.messageHandler()->setLogLevel(par.bbLogLevel);
250  if(par.algo > 0)
251    model.solver()->messageHandler()->setLogLevel(par.lpLogLevel);
252
253
254  //   model.setMaxFailure(par.maxFailures);
255  //   model.setMaxInfeasible(par.maxInfeasible);
256
257  //Pass over user set branching priorities to Cbc
258  {
259    //set priorities, prefered directions...
260    const int * priorities = nlpSolver.getPriorities();
261    const double * upPsCosts = nlpSolver.getUpPsCosts();
262    const double * downPsCosts = nlpSolver.getDownPsCosts();
263    const int * directions = nlpSolver.getBranchingDirections();
264    bool hasPseudo = (upPsCosts!=NULL);
265    model.findIntegers(true,hasPseudo);
266    CbcObject ** simpleIntegerObjects = model.objects();
267    int numberObjects = model.numberObjects();
268    for(int i = 0 ; i < numberObjects ; i++)
269    {
270      int iCol = simpleIntegerObjects[i]->columnNumber();
271      if(priorities)
272        simpleIntegerObjects[i]->setPriority(priorities[iCol]);
273      if(directions)
274        simpleIntegerObjects[i]->setPreferredWay(directions[iCol]);
275      if(upPsCosts) {
276        CbcSimpleIntegerPseudoCost * pscObject =
277          dynamic_cast<CbcSimpleIntegerPseudoCost*> (simpleIntegerObjects[i]);
278        pscObject->setUpPseudoCost(upPsCosts[iCol]);
279        pscObject->setDownPseudoCost(downPsCosts[iCol]);
280      }
281    }
282
283  }
284
285  // Now pass user set Sos constraints (code inspired from CoinSolve.cpp)
286  const TMINLP::SosInfo * sos = nlpSolver.model()->sosConstraints();
287  if(!par.disableSos && sos && sos->num > 0) //we have some sos constraints
288  {
289    const int & numSos = sos->num;
290    CbcObject ** objects = new CbcObject*[numSos];
291    const int * starts = sos->starts;
292    const int * indices = sos->indices;
293    const char * types = sos->types;
294    const double * weights = sos->weights;
295    //verify if model has user set priorities
296    bool hasPriorities = false;
297    const int * varPriorities = nlpSolver.getPriorities();
298    int numberObjects = model.numberObjects();
299    if(varPriorities)
300    {
301      for(int i = 0 ; i < numberObjects ; i++) {
302        if(varPriorities[i]) {
303          hasPriorities = true;
304          break;
305        }
306      }
307    }
308    const int * sosPriorities = sos->priorities;
309    if(sosPriorities)
310    {
311      for(int i = 0 ; i < numSos ; i++) {
312        if(sosPriorities[i]) {
313          hasPriorities = true;
314          break;
315        }
316      }
317    }
318    for(int i = 0 ; i < numSos ; i++)
319    {
320      int start = starts[i];
321      int length = starts[i + 1] - start;
322      objects[i] = new CbcSOS(&model, length, &indices[start],
323          &weights[start], i, types[i]);
324
325      objects[i]->setPriority(10);
326      if(hasPriorities && sosPriorities && sosPriorities[i]) {
327        objects[i]->setPriority(sosPriorities[i]);
328      }
329    }
330    model.addObjects(numSos, objects);
331    for(int i = 0 ; i < numSos ; i++)
332      delete objects[i];
333    delete [] objects;
334  }
335
336  replaceIntegers(model.objects(), model.numberObjects());
337
338  model.setPrintFrequency(par.logInterval);
339
340  model.setDblParam(CbcModel::CbcCutoffIncrement, par.cutoffDecr);
341
342  model.setCutoff(par.cutoff);
343  //  model.setBestObjectiveValue(par.cutoff);
344
345  model.setDblParam(CbcModel::CbcAllowableGap, par.allowableGap);
346  model.setDblParam(CbcModel::CbcAllowableFractionGap, par.allowableFractionGap);
347
348  // Definition of node selection strategy
349  CbcCompareObjective compare0;
350  CbcCompareDepth compare1;
351  CbcCompareUser compare2;
352  if(par.nodeSelection==0) {
353    model.setNodeComparison(compare0);
354  }
355  else if(par.nodeSelection==1) {
356    model.setNodeComparison(compare1);
357  }
358  else if(par.nodeSelection==2) {
359    compare2.setWeight(0.0);
360    model.setNodeComparison(compare2);
361  }
362  else if(par.nodeSelection==3) {
363    model.setNodeComparison(compare2);
364  }
365
366  model.setNumberStrong(par.numberStrong);
367
368  model.setNumberBeforeTrust(par.minReliability);
369
370  model.setNumberPenalties(8);
371
372  model.setDblParam(CbcModel::CbcMaximumSeconds, par.maxTime);
373
374  model.setMaximumNodes(par.maxNodes);
375
376  model.setIntegerTolerance(par.intTol);
377
378
379  // Redundant definition of default branching (as Default == User)
380  CbcBranchUserDecision branch;
381  model.setBranchingMethod(&branch);
382
383  //Get the time and start.
384  model.initialSolve();
385
386  continuousRelaxation_ =model.solver()->getObjValue();
387  if(par.algo == 0)//Set warm start point for Ipopt
388  {
389    const double * colsol = model.solver()->getColSolution();
390    const double * duals = model.solver()->getRowPrice();
391    model.solver()->setColSolution(colsol);
392    model.solver()->setRowPrice(duals);
393  }
394 
395  CoinSighandler_t saveSignal=SIG_DFL;
396  // register signal handler
397  saveSignal = signal(SIGINT,signal_handler);
398
399
400  currentBranchModel = &model;
401 
402  model.branchAndBound();
403  numNodes_ = model.getNodeCount();
404  bestObj_ = model.getObjValue();
405  bestBound_ = model.getBestPossibleObjValue();
406  mipIterationCount_ = model.getIterationCount();
407
408  bool hasFailed = false;
409  if(par.algo==0)//Did we continue branching on a failure
410  {
411    BonminCbcNlpStrategy * nlpStrategy = dynamic_cast<BonminCbcNlpStrategy *>(model.strategy());
412    if(nlpStrategy)
413      hasFailed = nlpStrategy->hasFailed();
414    else
415      throw -1;
416  }
417  else
418    hasFailed = nlpSolver.hasContinuedOnAFailure();
419
420
421  if(hasFailed) {
422    std::cout<<"************************************************************"<<std::endl
423    <<"WARNING : Optimization failed on an NLP during optimization"
424    <<"\n (no optimal value found within tolerances)."<<std::endl
425    <<"Optimization was not stopped because option \n"
426    <<"\"nlp_failure_behavior\" has been set to fathom but"
427    <<" beware that reported solution may not be optimal"<<std::endl
428    <<"************************************************************"<<std::endl;
429  }
430
431  if(model.bestSolution()) {
432    if(bestSolution_)
433      delete [] bestSolution_;
434    bestSolution_ = new double[nlpSolver.getNumCols()];
435    CoinCopyN(model.bestSolution(), nlpSolver.getNumCols(), bestSolution_);
436  }
437  if(!model.status()) {
438    if(bestSolution_)
439      mipStatus_ = FeasibleOptimal;
440    else
441      mipStatus_ = ProvenInfeasible;
442  }
443  else {
444    if(bestSolution_)
445      mipStatus_ = Feasible;
446    else
447      mipStatus_ = NoSolutionKnown;
448  }
449
450
451  if(par.algo > 0)
452    delete si;
453#ifdef COIN_HAS_CPX
454
455  if (par.milpSubSolver ==1)
456    delete localSearchSolver;
457#endif
458  std::cout<<"Finished"<<std::endl;
459  if(strategy)
460    delete strategy;
461
462}
463
464
465/** return the best known lower bound on the objective value*/
466double
467BonminBB::bestBound()
468{
469  if(mipStatus_ == FeasibleOptimal) return bestObj_;
470  else if(mipStatus_ == ProvenInfeasible) return 1e200;
471  else return bestBound_;
472}
Note: See TracBrowser for help on using the repository browser.