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

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

merge r945 from trunk: allow compiling against Cgl 0.58

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