source: branches/devel/Bonmin/src/CbcBonmin/BonCbc.cpp @ 58

Last change on this file since 58 was 58, checked in by pbonami, 13 years ago

Finish renaming

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