source: trunk/Couenne/src/bound_tightening/aggressiveBT.cpp @ 131

Last change on this file since 131 was 131, checked in by pbelotti, 13 years ago

looser checks on bound tightening -- should prevent good upper/lower bound on objective from fathoming a node

File size: 7.9 KB
Line 
1/*
2 * Name:    aggressiveBT.cpp
3 * Author:  Pietro Belotti
4 * Purpose: aggressive bound tightening -- fake bounds in variables to
5 *          exclude parts of the solution space through fathoming on
6 *          bounds/infeasibility
7 *
8 * (C) Carnegie-Mellon University, 2007-08.
9 * This file is licensed under the Common Public License (CPL)
10 */
11
12#include "CoinHelperFunctions.hpp"
13#include "BonCouenneInfo.hpp"
14
15#include "CouenneCutGenerator.hpp"
16#include "CouenneProblem.hpp"
17
18#define MAX_ABT_ITER           1  // max # aggressive BT iterations
19#define THRES_ABT_IMPROVED     0  // only continue ABT if at least these bounds have improved
20#define THRES_ABT_ORIG      1000  // only do ABT on originals if they are more than this
21
22static double distanceToBound(int n, const double* xOrig,
23                              const double* lower, const double* upper)
24{
25  double Xdist = 0.;
26  for (int i=0; i<n; i++) {
27    if (lower[i] > xOrig[i]) {
28      Xdist += lower[i] - xOrig[i];
29    }
30    if (upper[i] < xOrig[i]) {
31      Xdist += xOrig[i] - upper[i];
32    }
33  }
34  return Xdist;
35}
36
37// Aggressive Bound Tightening: for each variable, fake new bounds
38// [l,b] or [b,u] and apply bound tightening. If the interval is
39// fathomed on bounds or on infeasibility, the complementary bound
40// interval is a valid tightening.
41
42bool CouenneProblem::aggressiveBT (Bonmin::OsiTMINLPInterface *nlp,
43                                   t_chg_bounds *chg_bds, 
44                                   Bonmin::BabInfo * babInfo) const {
45
46  Jnlst()->Printf (J_ITERSUMMARY, J_BOUNDTIGHTENING, "Aggressive FBBT\n");
47
48  Bonmin::CouenneInfo* couInfo =
49    dynamic_cast<Bonmin::CouenneInfo*>(babInfo);
50
51  int  ncols  = nVars ();
52  bool retval = true;
53
54  CouNumber
55    *olb = new CouNumber [ncols],
56    *oub = new CouNumber [ncols];
57
58  // save current bounds
59  CoinCopyN (Lb (), ncols, olb);
60  CoinCopyN (Ub (), ncols, oub);
61
62  // Find the solution that is closest to the current bounds
63  // TODO: Also check obj value
64  SmartPtr<const Bonmin::CouenneInfo::NlpSolution> closestSol;
65  double dist = 1e50;
66
67  const std::list<SmartPtr<const Bonmin::CouenneInfo::NlpSolution> >& solList =
68    couInfo->NlpSolutions();
69
70  for (std::list<SmartPtr<const Bonmin::CouenneInfo::NlpSolution> >::const_iterator
71         i = solList.begin();
72       i != solList.end(); i++) {
73    assert(nOrigVars_ == (*i)->nVars());
74    const double thisDist = distanceToBound(nOrigVars_, (*i)->solution(), olb, oub);
75    if (thisDist < dist) {
76      closestSol = *i;
77      dist = thisDist;
78    }
79  }
80
81  jnlst_ -> Printf(J_VECTOR, J_BOUNDTIGHTENING, "best dist = %e\n", dist);
82
83  // If this solution is not sufficiently inside the bounds, we solve the NLP now
84  if (dist > 0.1) { // TODO: Find tolerance
85
86    // find integer initial point /////////////////////////
87
88    int nvars = nVars ();
89
90    double *lower = new double [nvars];
91    double *upper = new double [nvars];
92
93    CoinFillN (lower, nvars, -COUENNE_INFINITY);
94    CoinFillN (upper, nvars,  COUENNE_INFINITY);
95
96    CoinCopyN (nlp -> getColLower (), nOrigVars_, lower);
97    CoinCopyN (nlp -> getColUpper (), nOrigVars_, upper);
98
99    double *Y = new double [nvars];
100    CoinFillN (Y, nvars, 0.);
101    CoinCopyN (X (), nOrigVars_, Y);
102
103    if (getIntegerCandidate (nlp -> getColSolution (), Y, lower, upper) < 0) {
104      jnlst_ -> Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING, "TODO: find NLP point in ABT failed\n");
105      delete [] Y;
106      delete [] lower;
107      delete [] upper;
108      return true;
109    }
110
111    //////////////////////////////////////////////////////
112
113    nlp -> setColLower    (lower);
114    nlp -> setColUpper    (upper);
115    nlp -> setColSolution (Y);
116
117    nlp -> initialSolve ();
118
119    delete [] Y;
120    delete [] lower;
121    delete [] upper;
122
123    if (nlp->isProvenOptimal()) {
124      closestSol = new Bonmin::CouenneInfo::NlpSolution
125        (nOrigVars_, nlp->getColSolution(), nlp->getObjValue());
126      couInfo->addSolution(closestSol);
127      dist = 0.;
128    }
129    else {
130      jnlst_ -> Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING, "TODO: NLP solve in ABT failed\n");
131      return true;
132    }
133  }
134
135  if (dist>1e10) {
136    jnlst_ -> Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING, "TODO: Don't have point for ABT\n");
137    return true;
138  }
139
140  // X is now the NLP solution, but in a low-dimensional space. We
141  // have to get the corresponding point in higher dimensional space
142  // through getAuxs()
143
144  double *X = new double [ncols];
145  CoinCopyN (closestSol->solution(), nOrigVars_, X);
146  getAuxs (X);
147
148  // create a new, fictitious, bound bookkeeping structure
149  t_chg_bounds *f_chg = new t_chg_bounds [ncols];
150
151  if (Jnlst()->ProduceOutput(J_ITERSUMMARY, J_BOUNDTIGHTENING)) {
152    //    CouNumber cutoff = getCutOff ();
153    int       objind = Obj (0) -> Body  () -> Index ();
154    for (int i=0; i<nOrigVars_; i++)
155      Jnlst()->Printf(J_MOREVECTOR, J_BOUNDTIGHTENING,
156                      "   %2d %+20g [%+20g %+20g]\n",
157                      i, X [i], Lb (i), Ub (i));
158    Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING,
159                    "-------------\nAggressive BT. Current bound = %g, cutoff = %g, %d vars\n", 
160                    Lb (objind), getCutOff (), ncols);
161  }
162
163  int improved, second, iter = 0;
164
165  // Repeatedly fake tightening bounds on both sides of every variable
166  // to concentrate around current NLP point.
167  //
168  // MAX_ABT_ITER is the maximum # of outer cycles. Each call to
169  // fake_tighten in turn has an iterative algorithm for a
170  // derivative-free, uni-dimensional optimization problem on a
171  // monotone function.
172
173  do {
174
175    improved = 0;
176
177    // scan all variables
178    for (int i=0; i<ncols; i++) {
179
180      if (CoinCpuTime () > maxCpuTime_)
181        break;
182
183      int index = evalOrder (i);
184
185      if (Var (index) -> Multiplicity () <= 0) 
186        continue;
187
188      // AW: We only want to do the loop that temporarily changes
189      // bounds around the NLP solution only for those points from the
190      // NLP solution (no auxiliary vars)?
191
192      // PBe: if we do want that, index should be initialized as i, as
193      // evalOrder gives a variable index out of an array index.
194
195      // PBe: That makes a lot of sense when problems are really
196      // big. Instances arki000[24].nl spend a lot of time here
197
198      if ((nOrigVars_ < THRES_ABT_ORIG) || (index < nOrigVars_)) {
199
200        // if (index == objind) continue; // don't do it on objective function
201
202        Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING,
203                        "------------- tighten left x%d\n", index);
204
205        // tighten on left
206        if ((X [index] >= Lb (index) + COUENNE_EPS)
207            && ((improved = fake_tighten (0, index, X, olb, oub, chg_bds, f_chg)) < 0)) {
208          retval = false;
209          break;
210        }
211
212        second = 0;
213
214        Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING,
215                        "------------- tighten right x%d\n", index);
216
217        // tighten on right
218        if ((X [index] <= Ub (index) - COUENNE_EPS)
219            && ((second = fake_tighten (1, index, X, olb, oub, chg_bds, f_chg) < 0))) {
220          retval = false;
221          break;
222        }
223
224        improved += second;
225      }
226    }
227  } while (retval && (improved > THRES_ABT_IMPROVED) && (iter++ < MAX_ABT_ITER));
228
229  // store new valid bounds, or restore old ones if none changed
230  CoinCopyN (olb, ncols, Lb ());
231  CoinCopyN (oub, ncols, Ub ());
232
233  if (Jnlst()->ProduceOutput(J_ITERSUMMARY, J_BOUNDTIGHTENING)) {
234    Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING,"------------------\n");
235
236    if (!retval) Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING,
237                                 "Couenne infeasible node from aggressive BT\n");
238
239    int objind = Obj (0) -> Body  () -> Index ();
240
241    Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING,
242                    "-------------\ndone Aggressive BT. Current bound = %g, cutoff = %g, %d vars\n", 
243                    Lb (objind), getCutOff (), ncols);
244
245    if (Jnlst()->ProduceOutput(J_DETAILED, J_BOUNDTIGHTENING))
246      for (int i=0; i<nOrigVars_; i++)
247        printf("   x%02d [%+20g %+20g]  | %+20g\n",
248               i, Lb (i), Ub (i), X [i]);
249
250    if (Jnlst()->ProduceOutput(J_MOREDETAILED, J_BOUNDTIGHTENING))
251      for (int i=nOrigVars_; i<ncols; i++)
252        printf ("   w%02d [%+20g %+20g]  | %+20g\n", i, Lb (i), Ub (i), X [i]);
253  }
254
255  delete [] X;
256  delete [] f_chg;
257  delete [] olb;
258  delete [] oub;
259
260  return retval;// && btCore (psi, cs, chg_bds, babInfo, true); // !!!
261  //return retval && btCore (psi, cs, chg_bds, babInfo, true);
262}
Note: See TracBrowser for help on using the repository browser.