source: trunk/Couenne/src/bound_tightening/CouenneAggrProbing.cpp @ 793

Last change on this file since 793 was 793, checked in by pbelotti, 9 years ago

Created CouenneBab? from Bonmin::Bab. Putting CouenneRecordBestSol? to good use: upper bound now saved onto .sol file using real best saved solution (should fix a lot of "infeasible" instances that were pruned by good bounds ignored by Cbc).

  • Property svn:keywords set to Author Date Id Revision
File size: 17.5 KB
Line 
1/* $Id: CouenneAggrProbing.cpp 793 2012-01-26 03:07:16Z pbelotti $
2 *
3 * Name:    CouenneAggrProbing.cpp
4 * Author:  Giacomo Nannicini
5 * Purpose: Aggressive probing
6 *
7 * (C) Giacomo Nannicini, 2010.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11#include "CouenneAggrProbing.hpp"
12#include "CouenneProblemElem.hpp"
13#include "CouenneExprVar.hpp"
14#include "CouenneExprOpp.hpp"
15//#include "BonCbc.hpp"
16#include "CouenneBab.hpp"
17#include "CouenneCutGenerator.hpp"
18#include <string>
19
20#define COUENNE_AGGR_PROBING_FINITE_BOUND 1.0e+10
21#define COUENNE_AGGR_PROBING_MIN_INTERVAL 1.0e-2
22#define COUENNE_AGGR_PROBING_BND_RELAX COUENNE_EPS
23
24using namespace Couenne;
25
26CouenneAggrProbing::CouenneAggrProbing(CouenneSetup *setup,
27                                       const Ipopt::SmartPtr<Ipopt::OptionsList> options)
28{
29  couenne_ = setup;
30  numCols_ = couenne_->couennePtr()->Problem()->nVars();
31  maxTime_ = COIN_DBL_MAX;
32  maxFailedSteps_ = 10;
33  maxNodes_ = 1000;
34  initCutoff_ = COUENNE_INFINITY;
35  restoreCutoff_ = false;
36
37}
38
39CouenneAggrProbing::CouenneAggrProbing(const CouenneAggrProbing &rhs){
40  couenne_ = new CouenneSetup(*rhs.couenne_);
41  numCols_ = rhs.numCols_;
42  maxTime_ = rhs.maxTime_;
43  maxFailedSteps_ = rhs.maxFailedSteps_;
44  maxNodes_ = rhs.maxNodes_;
45  initCutoff_ = rhs.initCutoff_;
46  restoreCutoff_ = rhs.restoreCutoff_;
47}
48
49CouenneAggrProbing::~CouenneAggrProbing(){
50}
51
52void CouenneAggrProbing::registerOptions(Ipopt::SmartPtr <Bonmin::RegisteredOptions> roptions) {
53  // Nothing for the moment, but will be added later as needed
54}
55
56void CouenneAggrProbing::generateCuts(const OsiSolverInterface& solver,
57                                      OsiCuts& cuts, 
58                                      const CglTreeInfo info) const {
59  // Empty for the moment (cannot be used automatically in Branch-and-Bound)
60}
61
62double CouenneAggrProbing::getMaxTime() const {
63  return maxTime_;
64}
65
66void CouenneAggrProbing::setMaxTime(double value){
67  maxTime_ = value;
68}
69
70int CouenneAggrProbing::getMaxFailedSteps() const {
71  return maxFailedSteps_;
72}
73
74void CouenneAggrProbing::setMaxFailedSteps(int value){
75  maxFailedSteps_ = value;
76}
77
78int CouenneAggrProbing::getMaxNodes() const {
79  return maxNodes_;
80}
81
82void CouenneAggrProbing::setMaxNodes(int value){
83  maxNodes_ = value;
84}
85
86bool CouenneAggrProbing::getRestoreCutoff() const {
87  return restoreCutoff_;
88}
89
90void CouenneAggrProbing::setRestoreCutoff(bool value){
91  restoreCutoff_ = value;
92}
93
94double CouenneAggrProbing::probeVariable(int index, bool probeLower){
95
96  // Useful objects for easy access
97  OsiSolverInterface* nlp = couenne_->nonlinearSolver();
98  OsiSolverInterface* lp = couenne_->continuousSolver();
99  CouenneProblem* problem = couenne_->couennePtr()->Problem();
100
101  // Save initial bounds
102  double initUpper = lp->getColUpper()[index];
103  double initLower = lp->getColLower()[index];
104
105  double* initLowerLp = new double[numCols_];
106  double* initUpperLp = new double[numCols_];
107
108  memcpy(initLowerLp, lp->getColLower(), numCols_*sizeof(double));
109  memcpy(initUpperLp, lp->getColUpper(), numCols_*sizeof(double));
110
111  if (initUpper < initLower + COUENNE_EPS){
112    // Variable is fixed, so we can't tighten
113    return ((probeLower) ? initLower : initUpper);
114  }
115
116  // Index of the aux variable representing the objective function
117  int indobj = problem->Obj(0)->Body()->Index();
118
119  // Initial cutoff value
120  double initCutoff = problem->getCutOff (); //Ub()[indobj];
121
122  double* initCutoffSol = NULL; 
123
124  if (restoreCutoff_ && problem->getCutOff() < COUENNE_INFINITY){
125    initCutoffSol = new double[numCols_];
126    memcpy(initCutoffSol, problem->getCutOffSol(), numCols_*sizeof(double));
127  }
128
129  // Save parameters
130  Bonmin::BabSetupBase::NodeComparison initNodeComparison = 
131    couenne_->nodeComparisonMethod();
132  int initMaxNodes = couenne_->getIntParameter(Bonmin::BabSetupBase::MaxNodes);
133  double initMaxTime = couenne_->getDoubleParameter(Bonmin::BabSetupBase::MaxTime);
134  int initMaxSol = couenne_->getIntParameter(Bonmin::BabSetupBase::MaxSolutions);
135  couenne_->setNodeComparisonMethod(Bonmin::BabSetupBase::bestBound);
136  //couenne_->nodeComparisonMethod() = Bonmin::BabSetupBase::bestBound;
137  couenne_->setIntParameter(Bonmin::BabSetupBase::MaxNodes, maxNodes_);
138  couenne_->setIntParameter(Bonmin::BabSetupBase::MaxSolutions, COIN_INT_MAX);
139  problem->setCheckAuxBounds(true);
140
141  /// First store, then disable all heuristics.
142  Bonmin::BabSetupBase::HeuristicMethods heuristics = couenne_->heuristics();
143  couenne_->heuristics().clear();
144
145  double currentBound = (probeLower) ? initLower : initUpper;
146  double startTime = CoinCpuTime();
147  int failedSteps = 0;
148  double intervalSize = 0.0;
149  double tryBound = 0.0;
150
151  int iter = 0;
152
153  if (probeLower)
154    std::cout << "Probing lower on var " << index << std::endl;
155  else
156    std::cout << "Probing upper on var " << index << std::endl;
157
158  if ((fabs(currentBound) > COUENNE_AGGR_PROBING_FINITE_BOUND) &&
159      ((probeLower && initUpper > -COUENNE_AGGR_PROBING_FINITE_BOUND) ||
160       (!probeLower && initLower < COUENNE_AGGR_PROBING_FINITE_BOUND))){
161    // The bound is too large to apply the standard probing method;
162    // try to reduce it to a finite value. We only do this if we want
163    // to probe a variable on an infinite (or close to) bound, and the
164    // other bound of the variable is sufficiently far away
165    if (probeLower){
166      tryBound = -COUENNE_AGGR_PROBING_FINITE_BOUND;
167      lp->setColLower(index, currentBound);
168      problem->Lb()[index] = currentBound;
169      lp->setColUpper(index, tryBound);
170      problem->Ub()[index] = tryBound;
171      if (index < problem->nOrigVars()){
172        nlp->setColLower(index, currentBound);
173        nlp->setColUpper(index, tryBound);
174      }
175    }
176    else{
177      tryBound = COUENNE_AGGR_PROBING_FINITE_BOUND;
178      lp->setColLower(index, tryBound);
179      problem->Lb()[index] = tryBound;
180      lp->setColUpper(index, currentBound);
181      problem->Ub()[index] = currentBound;
182      if (index < problem->nOrigVars()){
183        nlp->setColLower(index, tryBound);
184        nlp->setColUpper(index, currentBound);
185      }
186    }
187
188    /// Setup Branch-and-Bound limits
189    couenne_->setDoubleParameter(Bonmin::BabSetupBase::MaxTime, 
190                                 CoinMin(maxTime_-(CoinCpuTime()-startTime),
191                                         maxTime_*0.5));
192
193    if (restoreCutoff_){
194      problem->resetCutOff(initCutoff);
195      if (indobj >= 0)
196        problem->Ub(indobj) = initCutoff;
197      problem->installCutOff();
198    }
199
200    std::cout << "Iteration " << iter << ", current bound " << currentBound
201              << ", try bound " << tryBound << std::endl;
202
203    /// Now do Branch-and-Bound and see if probing succeeded
204    // Bonmin::Bab bb;
205    // bb.setUsingCouenne(true);
206
207    CouenneBab bb;
208
209    bb(couenne_);
210    if (bb.model().isProvenInfeasible()){
211      /// Problem is infeasible; therefore, probing was successful.
212      currentBound = tryBound;
213      std::cout << "Probing succeeded; brought to finite" << std::endl;
214    }
215    else{
216      /// Problem is not infeasible; we failed
217      std::cout << "Probing failed; still infinity, exit" << std::endl;
218    }   
219    iter++;
220  }
221
222  // Now that we have a finite bound, pick size of the probing interval
223
224  // Override (for testing - will be chosen automatically in final
225  // implementation)
226  intervalSize = 0.1;
227
228  if (intervalSize < COUENNE_AGGR_PROBING_MIN_INTERVAL){
229    intervalSize = COUENNE_AGGR_PROBING_MIN_INTERVAL;
230  }
231 
232  while ((fabs(currentBound) <= COUENNE_AGGR_PROBING_FINITE_BOUND) &&
233         ((CoinCpuTime() - startTime) < maxTime_) &&
234         (failedSteps < maxFailedSteps_) &&
235         (intervalSize >= COUENNE_AGGR_PROBING_MIN_INTERVAL) && 
236         iter < 100){
237
238    /// Set the bound that we want to try
239    if (probeLower){
240      tryBound = currentBound + intervalSize;
241      if (tryBound > initUpper){
242        // It does not make sense to use bounds larger than the initial
243        // ones
244        tryBound = initUpper;
245      }
246      if (lp->isInteger(index)){
247        tryBound = floor(tryBound);
248      }
249      // Relax bounds a little bit
250      lp->setColLower(index, currentBound - COUENNE_AGGR_PROBING_BND_RELAX);
251      problem->Lb()[index] = currentBound - COUENNE_AGGR_PROBING_BND_RELAX;
252      lp->setColUpper(index, tryBound + COUENNE_AGGR_PROBING_BND_RELAX);
253      problem->Ub()[index] = tryBound + COUENNE_AGGR_PROBING_BND_RELAX;
254      if (index < problem->nOrigVars()){
255        nlp->setColLower(index, currentBound - COUENNE_AGGR_PROBING_BND_RELAX);
256        nlp->setColUpper(index, tryBound + COUENNE_AGGR_PROBING_BND_RELAX);
257      }
258    }
259    else{
260      tryBound = currentBound - intervalSize;
261      if (tryBound < initLower){
262        // It does not make sense to use bounds larger than the initial
263        // ones
264        tryBound = initLower;
265      }
266      if (lp->isInteger(index)){
267        tryBound = ceil(tryBound);
268      }
269      // Relax bounds a little bit
270      lp->setColLower(index, tryBound - COUENNE_AGGR_PROBING_BND_RELAX);
271      problem->Lb()[index] = tryBound - COUENNE_AGGR_PROBING_BND_RELAX;
272      lp->setColUpper(index, currentBound + COUENNE_AGGR_PROBING_BND_RELAX);
273      problem->Ub()[index] = currentBound + COUENNE_AGGR_PROBING_BND_RELAX;
274      if (index < problem->nOrigVars()){
275        nlp->setColLower(index, tryBound - COUENNE_AGGR_PROBING_BND_RELAX);
276        nlp->setColUpper(index, currentBound + COUENNE_AGGR_PROBING_BND_RELAX);
277      }
278    }
279
280    lp->resolve();
281    problem->domain()->push(numCols_, lp->getColSolution(),
282                            lp->getColLower(), lp->getColUpper());
283
284    /// Setup Branch-and-Bound limits
285    couenne_->setDoubleParameter(Bonmin::BabSetupBase::MaxTime, 
286                                 CoinMin(maxTime_-(CoinCpuTime()-startTime),
287                                         maxTime_*0.5));
288
289    if (restoreCutoff_){
290      if (indobj >= 0)
291        problem->Ub(indobj) = initCutoff;
292      problem->resetCutOff(initCutoff);
293      problem->installCutOff();
294    }
295
296    std::cout << "Iteration " << iter << ", current bound " << currentBound
297              << ", try bound " << tryBound << std::endl;
298
299    /// Now do Branch-and-Bound and see if probing succeeded
300    // Bonmin::Bab bb;
301    // bb.setUsingCouenne(true);
302
303    CouenneBab bb;
304    bb(couenne_);
305
306    problem->domain()->pop();
307
308    double obj = 0.0;
309    /// Is the search in the current interval complete?
310    bool intervalSearched = (bb.model().isProvenOptimal() || 
311                             bb.model().isProvenInfeasible());
312
313    if ((!intervalSearched) || // If the search is not complete
314        (restoreCutoff_ && // or we don't want to accept new solutions
315         problem->getCutOffSol() && // and we have a new feasible solution
316         problem->checkNLP(problem->getCutOffSol(), obj, true))){
317      /// Try again in a smaller interval
318      if (lp->isInteger(index) && fabs(tryBound-currentBound) < 0.5){
319        /// There is no smaller interval that we can try; bail out
320        failedSteps = maxFailedSteps_;
321      }
322      else{
323        intervalSize /= 2;
324      }
325      failedSteps++;
326      std::cout << "Probing failed; shrinking interval" << std::endl;
327    }
328    else{
329      /// We fully explored the current interval, and there is no
330      /// feasible solution, or there is a solution and we have
331      /// already updated the cutoff. So, we can eliminate the current
332      /// interval. We also double the size of the search interval.
333      if (lp->isInteger(index) && fabs(tryBound-currentBound) < 0.5){
334        /// Make sure we increase by at least one if it is an integer
335        /// variable
336        intervalSize = 1.0;
337      }
338      else{
339        intervalSize *= 2;
340      }
341      currentBound = tryBound;
342      if (lp->isInteger(index)){
343        if (probeLower){
344          currentBound += 1.0;
345        }
346        else {
347          currentBound -= 1.0;
348        }
349      }
350      failedSteps = 0;
351      std::cout << "Probing succeeded; enlarging interval" << std::endl;
352    }
353
354    // Check early termination condition: if we manage to fix the
355    // variable (unlikely), there is nothing more we can do
356    if ((probeLower && fabs(currentBound-initUpper) < COUENNE_EPS) ||
357        (!probeLower && fabs(currentBound-initLower) < COUENNE_EPS)){
358      failedSteps = maxFailedSteps_;
359    }
360
361    // Reset cutoff
362    if (restoreCutoff_){
363      if (indobj >= 0)
364        problem->Ub(indobj) = initCutoff;
365      problem->resetCutOff(initCutoff);
366      problem->installCutOff();
367    }
368
369    problem->domain()->pop();
370
371    iter++;
372  }
373
374  /// Restore initial bounds (we do not want to modify the
375  /// CouenneSetup object: the caller will do that, if needed)
376  lp->setColLower(initLowerLp);
377  lp->setColUpper(initUpperLp);
378  nlp->setColLower(initLowerLp);
379  nlp->setColUpper(initUpperLp);
380  memcpy(problem->Lb(), initLowerLp, numCols_*sizeof(double));
381  memcpy(problem->Ub(), initUpperLp, numCols_*sizeof(double));
382
383  /// Restore parameters and heuristics
384  problem->setCheckAuxBounds(false);
385  //couenne_->nodeComparisonMethod() = initNodeComparison;
386  couenne_->setNodeComparisonMethod(initNodeComparison);
387  couenne_->setIntParameter(Bonmin::BabSetupBase::MaxSolutions, initMaxSol);
388  couenne_->setIntParameter(Bonmin::BabSetupBase::MaxNodes, initMaxNodes);
389  couenne_->setDoubleParameter(Bonmin::BabSetupBase::MaxTime, initMaxTime);
390  couenne_->heuristics() = heuristics;
391
392  /// Restore cutoff
393  if (restoreCutoff_){
394    problem->resetCutOff();
395    problem->setCutOff(initCutoff, initCutoffSol);
396    if (initCutoffSol){
397      delete[] initCutoffSol;
398    }
399  }
400 
401  delete[] initLowerLp;
402  delete[] initUpperLp;
403
404  /// We are done; return best bound found.
405  return currentBound;
406
407}
408
409double CouenneAggrProbing::probeVariable2(int index, bool probeLower){
410  // Does not work! It's impossible to get Maximization problems working...
411  // Adding extra variables doesn't seem to do the trick
412
413  // Useful objects for easy access
414  OsiSolverInterface* lp = couenne_->continuousSolver();
415  CouenneProblem* problem = couenne_->couennePtr()->Problem();
416
417  // Save initial bounds
418  double initUpper = lp->getColUpper()[index];
419  double initLower = lp->getColLower()[index];
420
421  if (initUpper < initLower + COUENNE_EPS){
422    // Variable is fixed, so we can't probe anything
423    return ((probeLower) ? initLower : initUpper);
424  }
425
426  /// Modify the CouenneSetup object to use our options.
427  /// We store the initial values of all parameters that we modify,
428  /// so that we can restore them when we are done.
429  Bonmin::BabSetupBase::NodeComparison initNodeComparison = 
430    couenne_->nodeComparisonMethod();
431  int initMaxNodes = couenne_->getIntParameter(Bonmin::BabSetupBase::MaxNodes);
432  double initMaxTime = couenne_->getDoubleParameter(Bonmin::BabSetupBase::MaxTime);
433  int initMaxSol = couenne_->getIntParameter(Bonmin::BabSetupBase::MaxSolutions);
434  couenne_->setNodeComparisonMethod (Bonmin::BabSetupBase::bestBound);
435  couenne_->setIntParameter(Bonmin::BabSetupBase::MaxNodes, maxNodes_);
436  couenne_->setIntParameter(Bonmin::BabSetupBase::MaxSolutions, COIN_INT_MAX);
437
438  /// First store, then disable all heuristics - we do not need upper
439  /// bounds, plus they probably use a NLP object which we do not know
440  /// how to modify
441  Bonmin::BabSetupBase::HeuristicMethods heuristics = couenne_->heuristics();
442  couenne_->heuristics().clear();
443
444  /// Now, store and modify objective function in the CouenneProblem object
445  double* initLpObj = new double[numCols_];
446  memcpy(initLpObj, lp->getObjCoefficients(), numCols_*sizeof(double));
447  expression* initProbObj = problem->Obj(0)->Body();
448
449  double* newLpObj = new double[numCols_];
450  memset(newLpObj, 0, numCols_*sizeof(double));
451
452  //expression* exprObj = NULL; // PB: unused
453  expression* extraVar = NULL;
454
455  lp->writeLp("before");
456
457  if (probeLower){
458    std::cout << "Probing LOWER" << std::endl;
459    // Set LP objective function
460    newLpObj[index] = 1.0;
461    lp->setObjective(newLpObj);
462
463    lp->writeLp("lower");
464
465    // Set CouenneProblem objective
466    problem->Obj(0)->Body(problem->Variables()[index]);
467    // couenne_->setDoubleParameter(Bonmin::BabSetupBase::Cutoff, COIN_DBL_MAX);
468    // problem->setCutOff(initUpper, lp->getColUpper());
469    // problem->installCutOff();
470
471  }
472  else{
473    // We cannot maximize an objective function in Couenne; so, we
474    // have to introduce an additional variable, equal to the opposite
475    // of the variable that we want to maximize, and minimize that.
476
477    // Add one column and one row to the LP
478    int extraCol = numCols_;
479    lp->setObjective(newLpObj);
480    lp->addCol(0, NULL, NULL, -initUpper, -initLower, 1.0);
481
482    // Create the row x_{extraCol} = -x_{index}
483    int rowIndices[2] = {index, extraCol};
484    double rowElements[2] = {1.0, 1.0};
485    lp->addRow(2, rowIndices, rowElements, 0.0, 0.0);
486    lp->resolve();
487
488    // Create the expression x_{extraCol} = -x_{index} in CouenneProblem
489    extraVar = problem->addVariable(lp->isInteger(index), NULL);
490    // exprObj = new exprOpp(problem->Variables()[index]->clone());
491    // problem->addEQConstraint(extraVar, exprObj);
492    problem->Obj(0)->Body(extraVar);
493
494    // couenne_->setDoubleParameter(Bonmin::BabSetupBase::Cutoff, COIN_DBL_MAX);
495    // problem->setCutOff(-initLower, lp->getColLower());
496    // problem->installCutOff();
497
498    lp->writeLp("upper");
499  }
500
501  couenne_->setNodeComparisonMethod (Bonmin::BabSetupBase::bestBound);
502  couenne_->setIntParameter(Bonmin::BabSetupBase::MaxNodes, maxNodes_);
503  couenne_->setDoubleParameter(Bonmin::BabSetupBase::MaxTime, 
504                               maxTime_);
505
506  // Bonmin::Bab bb;
507  // bb.setUsingCouenne(true);
508
509  CouenneBab bb;
510  bb(couenne_);
511
512  double bestBound = bb.model().getBestPossibleObjValue();
513
514  std::cout << "Obtained bound: " << bb.model().getBestPossibleObjValue() << std::endl;
515
516
517  /// Restore parameters
518  couenne_->setNodeComparisonMethod (initNodeComparison);
519  couenne_->setIntParameter(Bonmin::BabSetupBase::MaxNodes, initMaxNodes);
520  couenne_->setDoubleParameter(Bonmin::BabSetupBase::MaxTime, initMaxTime);
521  couenne_->setIntParameter(Bonmin::BabSetupBase::MaxSolutions, initMaxSol);
522  couenne_->heuristics() = heuristics;
523
524  if (!probeLower){
525    int extra = lp->getNumCols()-1;
526    lp->deleteCols(1, &extra);
527    extra = lp->getNumRows()-1;
528    lp->deleteRows(1, &extra);
529    problem->Variables().pop_back();
530    delete extraVar;
531    // delete exprObj;
532  }
533
534  lp->setObjective(initLpObj);
535  problem->Obj(0)->Body(initProbObj);
536
537  delete[] initLpObj;
538  delete[] newLpObj;
539
540  return ((probeLower) ? bestBound : -bestBound);
541
542}
Note: See TracBrowser for help on using the repository browser.