source: stable/0.1/Bonmin/src/CbcBonmin/CbcBonmin.cpp @ 470

Last change on this file since 470 was 470, checked in by pbonami, 12 years ago

Fix a stupid bug in retrieving optimal solution

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