source: stable/0.4/Couenne/src/bound_tightening/CouenneAggrProbing.cpp @ 795

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

merging changes 792-793 trunk -> stable/0.4 (exposing feasible solutions ignored by Cbc)

  • Property svn:keywords set to Author Date Id Revision
File size: 17.5 KB
Line 
1/* $Id: CouenneAggrProbing.cpp 795 2012-01-26 03:19:52Z 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->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      problem->Ub()[indobj] = initCutoff;
196      problem->installCutOff();
197    }
198
199    std::cout << "Iteration " << iter << ", current bound " << currentBound
200              << ", try bound " << tryBound << std::endl;
201
202    /// Now do Branch-and-Bound and see if probing succeeded
203    // Bonmin::Bab bb;
204    // bb.setUsingCouenne(true);
205
206    CouenneBab bb;
207
208    bb(couenne_);
209    if (bb.model().isProvenInfeasible()){
210      /// Problem is infeasible; therefore, probing was successful.
211      currentBound = tryBound;
212      std::cout << "Probing succeeded; brought to finite" << std::endl;
213    }
214    else{
215      /// Problem is not infeasible; we failed
216      std::cout << "Probing failed; still infinity, exit" << std::endl;
217    }   
218    iter++;
219  }
220
221  // Now that we have a finite bound, pick size of the probing interval
222
223  // Override (for testing - will be chosen automatically in final
224  // implementation)
225  intervalSize = 0.1;
226
227  if (intervalSize < COUENNE_AGGR_PROBING_MIN_INTERVAL){
228    intervalSize = COUENNE_AGGR_PROBING_MIN_INTERVAL;
229  }
230 
231  while ((fabs(currentBound) <= COUENNE_AGGR_PROBING_FINITE_BOUND) &&
232         ((CoinCpuTime() - startTime) < maxTime_) &&
233         (failedSteps < maxFailedSteps_) &&
234         (intervalSize >= COUENNE_AGGR_PROBING_MIN_INTERVAL) && 
235         iter < 100){
236
237    /// Set the bound that we want to try
238    if (probeLower){
239      tryBound = currentBound + intervalSize;
240      if (tryBound > initUpper){
241        // It does not make sense to use bounds larger than the initial
242        // ones
243        tryBound = initUpper;
244      }
245      if (lp->isInteger(index)){
246        tryBound = floor(tryBound);
247      }
248      // Relax bounds a little bit
249      lp->setColLower(index, currentBound - COUENNE_AGGR_PROBING_BND_RELAX);
250      problem->Lb()[index] = currentBound - COUENNE_AGGR_PROBING_BND_RELAX;
251      lp->setColUpper(index, tryBound + COUENNE_AGGR_PROBING_BND_RELAX);
252      problem->Ub()[index] = tryBound + COUENNE_AGGR_PROBING_BND_RELAX;
253      if (index < problem->nOrigVars()){
254        nlp->setColLower(index, currentBound - COUENNE_AGGR_PROBING_BND_RELAX);
255        nlp->setColUpper(index, tryBound + COUENNE_AGGR_PROBING_BND_RELAX);
256      }
257    }
258    else{
259      tryBound = currentBound - intervalSize;
260      if (tryBound < initLower){
261        // It does not make sense to use bounds larger than the initial
262        // ones
263        tryBound = initLower;
264      }
265      if (lp->isInteger(index)){
266        tryBound = ceil(tryBound);
267      }
268      // Relax bounds a little bit
269      lp->setColLower(index, tryBound - COUENNE_AGGR_PROBING_BND_RELAX);
270      problem->Lb()[index] = tryBound - COUENNE_AGGR_PROBING_BND_RELAX;
271      lp->setColUpper(index, currentBound + COUENNE_AGGR_PROBING_BND_RELAX);
272      problem->Ub()[index] = currentBound + COUENNE_AGGR_PROBING_BND_RELAX;
273      if (index < problem->nOrigVars()){
274        nlp->setColLower(index, tryBound - COUENNE_AGGR_PROBING_BND_RELAX);
275        nlp->setColUpper(index, currentBound + COUENNE_AGGR_PROBING_BND_RELAX);
276      }
277    }
278
279    lp->resolve();
280    problem->domain()->push(numCols_, lp->getColSolution(),
281                            lp->getColLower(), lp->getColUpper());
282
283    /// Setup Branch-and-Bound limits
284    couenne_->setDoubleParameter(Bonmin::BabSetupBase::MaxTime, 
285                                 CoinMin(maxTime_-(CoinCpuTime()-startTime),
286                                         maxTime_*0.5));
287
288    if (restoreCutoff_){
289      problem->Ub()[indobj] = initCutoff;
290      problem->resetCutOff(initCutoff);
291      problem->installCutOff();
292    }
293
294    std::cout << "Iteration " << iter << ", current bound " << currentBound
295              << ", try bound " << tryBound << std::endl;
296
297    /// Now do Branch-and-Bound and see if probing succeeded
298    // Bonmin::Bab bb;
299    // bb.setUsingCouenne(true);
300
301    CouenneBab bb;
302    bb(couenne_);
303
304    problem->domain()->pop();
305
306    double obj = 0.0;
307    /// Is the search in the current interval complete?
308    bool intervalSearched = (bb.model().isProvenOptimal() || 
309                             bb.model().isProvenInfeasible());
310
311    if ((!intervalSearched) || // If the search is not complete
312        (restoreCutoff_ && // or we don't want to accept new solutions
313         problem->getCutOffSol() && // and we have a new feasible solution
314         problem->checkNLP(problem->getCutOffSol(), obj, true))){
315      /// Try again in a smaller interval
316      if (lp->isInteger(index) && fabs(tryBound-currentBound) < 0.5){
317        /// There is no smaller interval that we can try; bail out
318        failedSteps = maxFailedSteps_;
319      }
320      else{
321        intervalSize /= 2;
322      }
323      failedSteps++;
324      std::cout << "Probing failed; shrinking interval" << std::endl;
325    }
326    else{
327      /// We fully explored the current interval, and there is no
328      /// feasible solution, or there is a solution and we have
329      /// already updated the cutoff. So, we can eliminate the current
330      /// interval. We also double the size of the search interval.
331      if (lp->isInteger(index) && fabs(tryBound-currentBound) < 0.5){
332        /// Make sure we increase by at least one if it is an integer
333        /// variable
334        intervalSize = 1.0;
335      }
336      else{
337        intervalSize *= 2;
338      }
339      currentBound = tryBound;
340      if (lp->isInteger(index)){
341        if (probeLower){
342          currentBound += 1.0;
343        }
344        else {
345          currentBound -= 1.0;
346        }
347      }
348      failedSteps = 0;
349      std::cout << "Probing succeeded; enlarging interval" << std::endl;
350    }
351
352    // Check early termination condition: if we manage to fix the
353    // variable (unlikely), there is nothing more we can do
354    if ((probeLower && fabs(currentBound-initUpper) < COUENNE_EPS) ||
355        (!probeLower && fabs(currentBound-initLower) < COUENNE_EPS)){
356      failedSteps = maxFailedSteps_;
357    }
358
359    // Reset cutoff
360    if (restoreCutoff_){
361      problem->Ub()[indobj] = initCutoff;
362      problem->resetCutOff(initCutoff);
363      problem->installCutOff();
364    }
365
366    problem->domain()->pop();
367
368    iter++;
369  }
370
371  /// Restore initial bounds (we do not want to modify the
372  /// CouenneSetup object: the caller will do that, if needed)
373  lp->setColLower(initLowerLp);
374  lp->setColUpper(initUpperLp);
375  nlp->setColLower(initLowerLp);
376  nlp->setColUpper(initUpperLp);
377  memcpy(problem->Lb(), initLowerLp, numCols_*sizeof(double));
378  memcpy(problem->Ub(), initUpperLp, numCols_*sizeof(double));
379
380  /// Restore parameters and heuristics
381  problem->setCheckAuxBounds(false);
382  //couenne_->nodeComparisonMethod() = initNodeComparison;
383  couenne_->setNodeComparisonMethod(initNodeComparison);
384  couenne_->setIntParameter(Bonmin::BabSetupBase::MaxSolutions, initMaxSol);
385  couenne_->setIntParameter(Bonmin::BabSetupBase::MaxNodes, initMaxNodes);
386  couenne_->setDoubleParameter(Bonmin::BabSetupBase::MaxTime, initMaxTime);
387  couenne_->heuristics() = heuristics;
388
389  /// Restore cutoff
390  if (restoreCutoff_){
391    problem->resetCutOff();
392    problem->setCutOff(initCutoff, initCutoffSol);
393    if (initCutoffSol){
394      delete[] initCutoffSol;
395    }
396  }
397 
398  delete[] initLowerLp;
399  delete[] initUpperLp;
400
401  /// We are done; return best bound found.
402  return currentBound;
403
404}
405
406double CouenneAggrProbing::probeVariable2(int index, bool probeLower){
407  // Does not work! It's impossible to get Maximization problems working...
408  // Adding extra variables doesn't seem to do the trick
409
410  // Useful objects for easy access
411  OsiSolverInterface* lp = couenne_->continuousSolver();
412  CouenneProblem* problem = couenne_->couennePtr()->Problem();
413
414  // Save initial bounds
415  double initUpper = lp->getColUpper()[index];
416  double initLower = lp->getColLower()[index];
417
418  if (initUpper < initLower + COUENNE_EPS){
419    // Variable is fixed, so we can't probe anything
420    return ((probeLower) ? initLower : initUpper);
421  }
422
423  /// Modify the CouenneSetup object to use our options.
424  /// We store the initial values of all parameters that we modify,
425  /// so that we can restore them when we are done.
426  Bonmin::BabSetupBase::NodeComparison initNodeComparison = 
427    couenne_->nodeComparisonMethod();
428  int initMaxNodes = couenne_->getIntParameter(Bonmin::BabSetupBase::MaxNodes);
429  double initMaxTime = couenne_->getDoubleParameter(Bonmin::BabSetupBase::MaxTime);
430  int initMaxSol = couenne_->getIntParameter(Bonmin::BabSetupBase::MaxSolutions);
431  couenne_->setNodeComparisonMethod (Bonmin::BabSetupBase::bestBound);
432  couenne_->setIntParameter(Bonmin::BabSetupBase::MaxNodes, maxNodes_);
433  couenne_->setIntParameter(Bonmin::BabSetupBase::MaxSolutions, COIN_INT_MAX);
434
435  /// First store, then disable all heuristics - we do not need upper
436  /// bounds, plus they probably use a NLP object which we do not know
437  /// how to modify
438  Bonmin::BabSetupBase::HeuristicMethods heuristics = couenne_->heuristics();
439  couenne_->heuristics().clear();
440
441  /// Now, store and modify objective function in the CouenneProblem object
442  double* initLpObj = new double[numCols_];
443  memcpy(initLpObj, lp->getObjCoefficients(), numCols_*sizeof(double));
444  expression* initProbObj = problem->Obj(0)->Body();
445
446  double* newLpObj = new double[numCols_];
447  memset(newLpObj, 0, numCols_*sizeof(double));
448
449  //expression* exprObj = NULL; // PB: unused
450  expression* extraVar = NULL;
451
452  lp->writeLp("before");
453
454  if (probeLower){
455    std::cout << "Probing LOWER" << std::endl;
456    // Set LP objective function
457    newLpObj[index] = 1.0;
458    lp->setObjective(newLpObj);
459
460    lp->writeLp("lower");
461
462    // Set CouenneProblem objective
463    problem->Obj(0)->Body(problem->Variables()[index]);
464    // couenne_->setDoubleParameter(Bonmin::BabSetupBase::Cutoff, COIN_DBL_MAX);
465    // problem->setCutOff(initUpper, lp->getColUpper());
466    // problem->installCutOff();
467
468  }
469  else{
470    // We cannot maximize an objective function in Couenne; so, we
471    // have to introduce an additional variable, equal to the opposite
472    // of the variable that we want to maximize, and minimize that.
473
474    // Add one column and one row to the LP
475    int extraCol = numCols_;
476    lp->setObjective(newLpObj);
477    lp->addCol(0, NULL, NULL, -initUpper, -initLower, 1.0);
478
479    // Create the row x_{extraCol} = -x_{index}
480    int rowIndices[2] = {index, extraCol};
481    double rowElements[2] = {1.0, 1.0};
482    lp->addRow(2, rowIndices, rowElements, 0.0, 0.0);
483    lp->resolve();
484
485    // Create the expression x_{extraCol} = -x_{index} in CouenneProblem
486    extraVar = problem->addVariable(lp->isInteger(index), NULL);
487    // exprObj = new exprOpp(problem->Variables()[index]->clone());
488    // problem->addEQConstraint(extraVar, exprObj);
489    problem->Obj(0)->Body(extraVar);
490
491    // couenne_->setDoubleParameter(Bonmin::BabSetupBase::Cutoff, COIN_DBL_MAX);
492    // problem->setCutOff(-initLower, lp->getColLower());
493    // problem->installCutOff();
494
495    lp->writeLp("upper");
496  }
497
498  couenne_->setNodeComparisonMethod (Bonmin::BabSetupBase::bestBound);
499  couenne_->setIntParameter(Bonmin::BabSetupBase::MaxNodes, maxNodes_);
500  couenne_->setDoubleParameter(Bonmin::BabSetupBase::MaxTime, 
501                               maxTime_);
502
503  // Bonmin::Bab bb;
504  // bb.setUsingCouenne(true);
505
506  CouenneBab bb;
507  bb(couenne_);
508
509  double bestBound = bb.model().getBestPossibleObjValue();
510
511  std::cout << "Obtained bound: " << bb.model().getBestPossibleObjValue() << std::endl;
512
513
514  /// Restore parameters
515  couenne_->setNodeComparisonMethod (initNodeComparison);
516  couenne_->setIntParameter(Bonmin::BabSetupBase::MaxNodes, initMaxNodes);
517  couenne_->setDoubleParameter(Bonmin::BabSetupBase::MaxTime, initMaxTime);
518  couenne_->setIntParameter(Bonmin::BabSetupBase::MaxSolutions, initMaxSol);
519  couenne_->heuristics() = heuristics;
520
521  if (!probeLower){
522    int extra = lp->getNumCols()-1;
523    lp->deleteCols(1, &extra);
524    extra = lp->getNumRows()-1;
525    lp->deleteRows(1, &extra);
526    problem->Variables().pop_back();
527    delete extraVar;
528    // delete exprObj;
529  }
530
531  lp->setObjective(initLpObj);
532  problem->Obj(0)->Body(initProbObj);
533
534  delete[] initLpObj;
535  delete[] newLpObj;
536
537  return ((probeLower) ? bestBound : -bestBound);
538
539}
Note: See TracBrowser for help on using the repository browser.