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

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

astyled the devel branch

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