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

Last change on this file since 945 was 945, checked in by stefan, 8 years ago

generateCuts in Cgl >= 0.58 will not be const anymore

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