source: stable/0.4/Couenne/src/bound_tightening/aggressiveBT.cpp @ 903

Last change on this file since 903 was 903, checked in by pbelotti, 8 years ago

merged change 902 on consistent lb/ub passed to Ipopt

  • Property svn:keywords set to Author Date Id Revision
File size: 9.9 KB
Line 
1/* $Id: aggressiveBT.cpp 903 2012-09-04 03:58:38Z pbelotti $
2 *
3 * Name:    aggressiveBT.cpp
4 * Author:  Pietro Belotti
5 * Purpose: probing -- fake bounds in variables to exclude parts of
6 *          the solution space through fathoming on
7 *          bounds/infeasibility
8 *
9 * (C) Carnegie-Mellon University, 2007-11.
10 * This file is licensed under the Eclipse Public License (EPL)
11 */
12
13#include "CouenneCutGenerator.hpp"
14
15#include "BonTNLPSolver.hpp"
16#include "BonNlpHeuristic.hpp"
17#include "CoinHelperFunctions.hpp"
18#include "BonCouenneInfo.hpp"
19
20#include "CouenneProblem.hpp"
21#include "CouenneExprVar.hpp"
22#include "CouenneProblemElem.hpp"
23
24using namespace Ipopt;
25using namespace Couenne;
26
27namespace Bonmin {
28  class OsiTMINLPInterface;
29  class BabInfo;
30  class TNLPSolver;
31}
32
33#define MAX_ABT_ITER           4  // max # aggressive BT iterations
34#define THRES_ABT_IMPROVED     0  // only continue ABT if at least these bounds have improved
35#define THRES_ABT_ORIG       100  // only do ABT on auxiliaries if they are less originals than this
36
37static double distanceToBound (register int n, 
38                               register const double* xOrig,
39                               register const double* lower, 
40                               register const double* upper,
41                               register double cutoff_distance) { // stop if distance is above this
42
43  register double Xdist = 0.;
44
45  for (; n--; ++upper, ++xOrig) {
46
47    register CouNumber diff = *lower++ - *xOrig;
48
49    if      ( diff                    > 0.) {if ((Xdist += diff) > cutoff_distance) break;}
50    else if ((diff = *xOrig - *upper) > 0.) {if ((Xdist += diff) > cutoff_distance) break;}
51  }
52
53  return Xdist;
54}
55
56
57// Probing: for each variable, fake new bounds [l,b] or [b,u], with
58// given b, and apply bound tightening. If the interval is fathomed on
59// bounds or on infeasibility, the complementary bound interval is a
60// valid tightening.
61
62bool CouenneProblem::aggressiveBT (Bonmin::OsiTMINLPInterface *nlp,
63                                   t_chg_bounds *chg_bds, 
64                                   const CglTreeInfo &info,
65                                   Bonmin::BabInfo * babInfo) const {
66
67  if (info.level <= 0 && !(info.inTree))  {
68    jnlst_ -> Printf (J_ERROR, J_COUENNE, "Probing: ");
69    fflush (stdout);
70  }
71
72  CouenneInfo* couInfo =
73    dynamic_cast <CouenneInfo *> (babInfo);
74
75  int  ncols  = nVars ();
76  bool retval = false;
77
78  CouNumber
79    *olb = CoinCopyOfArray (Lb (), ncols),
80    *oub = CoinCopyOfArray (Ub (), ncols);
81
82  // Find the solution that is closest to the current bounds
83  // TODO: Also check obj value
84  SmartPtr<const CouenneInfo::NlpSolution> closestSol;
85  double dist = 1.e50;
86
87  if (couInfo) {
88
89    const std::list<SmartPtr<const CouenneInfo::NlpSolution> >& solList =
90      couInfo->NlpSolutions();
91
92    for (std::list<SmartPtr<const CouenneInfo::NlpSolution> >::const_iterator
93           i = solList.begin();
94         i  != solList.end(); ++i) {
95
96      assert (nOrigVars_ - ndefined_ == (*i)->nVars());
97
98      const double thisDist = distanceToBound (nOrigVars_ - ndefined_, (*i)->solution(), olb, oub, dist);
99
100      if (thisDist < dist) {
101        closestSol = *i;
102        dist = thisDist;
103      }
104    }
105  }
106
107  jnlst_ -> Printf (J_VECTOR, J_BOUNDTIGHTENING, "best dist = %e\n", dist);
108
109  bool haveNLPsol = false;
110
111  // If this solution is not sufficiently inside the bounds, we solve the NLP now
112  if (dist > 0.1) { // TODO: Find tolerance
113
114    // find integer initial point /////////////////////////
115
116    int nvars = nVars ();
117
118    double *lower = new double [nvars];
119    double *upper = new double [nvars];
120
121    CoinFillN (lower, nvars, -COUENNE_INFINITY);
122    CoinFillN (upper, nvars,  COUENNE_INFINITY);
123
124    CoinCopyN (nlp -> getColLower (), nOrigVars_ - ndefined_, lower);
125    CoinCopyN (nlp -> getColUpper (), nOrigVars_ - ndefined_, upper);
126
127    double *Y = new double [nvars];
128
129    CoinZeroN (Y,    nvars);
130    CoinCopyN (X (), nOrigVars_ - ndefined_, Y);
131
132    if (getIntegerCandidate (nlp -> getColSolution () ? nlp -> getColSolution () : Y, Y, lower, upper) < 0) {
133
134      jnlst_ -> Printf (J_ITERSUMMARY, J_BOUNDTIGHTENING, "TODO: find NLP point in ABT failed\n");
135      retval = true;
136
137    } else {
138
139      //////////////////////////////////////////////////////
140
141      nlp -> setColLower    (lower);
142      nlp -> setColUpper    (upper);
143      nlp -> setColSolution (Y);
144
145      try {
146
147        nlp -> options () -> SetNumericValue ("max_cpu_time", CoinMax (0., getMaxCpuTime () - CoinCpuTime ()));
148        nlp -> initialSolve ();
149      }
150
151      catch (Bonmin::TNLPSolver::UnsolvedError *E) {}
152   
153      delete [] Y;
154      delete [] lower;
155      delete [] upper;
156
157      if (nlp -> isProvenOptimal ()) {
158
159        if (couInfo) {
160          closestSol = new CouenneInfo::NlpSolution
161            (nOrigVars_ - ndefined_, nlp->getColSolution(), nlp->getObjValue());
162          couInfo->addSolution(closestSol);
163          dist = 0.;
164          haveNLPsol = true;     
165        }
166      }
167      else {
168        jnlst_ -> Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING, "TODO: NLP solve in ABT failed\n");
169        retval = true;
170      }
171    }
172  }
173
174  int nTotImproved = 0;
175
176  // Probing can also run on an LP point.
177
178  //if (true || (retval && (dist < 1e10))) {
179
180  {
181    retval = true;
182
183    // if the NLP solver succeeded, X is now the NLP solution, but in
184    // a low-dimensional space. We have to get the corresponding point
185    // in higher dimensional space through getAuxs()
186
187    double *X = NULL;
188
189    if (haveNLPsol) {
190
191      X = new double [ncols];
192      CoinZeroN (X, nVars ());
193      CoinCopyN (closestSol -> solution(), nOrigVars_ - ndefined_, X);
194      getAuxs (X);
195    } else X = domain () -> x ();
196
197    // create a new, fictitious, bound bookkeeping structure
198    t_chg_bounds *f_chg = new t_chg_bounds [ncols];
199
200    if (Jnlst()->ProduceOutput(J_ITERSUMMARY, J_BOUNDTIGHTENING)) {
201      //    CouNumber cutoff = getCutOff ();
202      int       objind = Obj (0) -> Body  () -> Index ();
203      for (int i=0; i<nOrigVars_ - ndefined_; i++)
204        Jnlst()->Printf(J_MOREVECTOR, J_BOUNDTIGHTENING,
205                        "   %2d %+20g [%+20g %+20g]\n",
206                        i, X [i], Lb (i), Ub (i));
207      if (objind >= 0)
208        Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING,
209                        "-------------\nAggressive BT. Current bound = %g, cutoff = %g, %d vars\n", 
210                        Lb (objind), getCutOff (), ncols);
211    }
212
213    int improved, second, iter = 0;
214
215    // Repeatedly fake tightening bounds on both sides of every variable
216    // to concentrate around current NLP point.
217    //
218    // MAX_ABT_ITER is the maximum # of outer cycles. Each call to
219    // fake_tighten in turn has an iterative algorithm for a
220    // derivative-free, uni-dimensional optimization problem on a
221    // monotone function.
222
223    bool maxTimeReached = false;
224
225    do {
226
227      improved = 0;
228
229      // scan all variables
230      for (int i=0; i<ncols; i++) {
231
232        if (CoinCpuTime () > maxCpuTime_) {
233          maxTimeReached = true; // avoids getrusage after loop...
234          break;
235        }
236
237        int index = evalOrder (i);
238
239        if ((Var (index) -> Multiplicity () <= 0) ||
240            (fabs (Lb (index) - Ub (index)) < COUENNE_EPS) ||
241            ((nOrigVars_ >= THRES_ABT_ORIG) &&
242             (index >= nOrigVars_)))
243          continue;
244
245        // AW: We only want to do the loop that temporarily changes
246        // bounds around the NLP solution only for those points from the
247        // NLP solution (no auxiliary vars)?
248
249        // PBe: if we do want that, index should be initialized as i, as
250        // evalOrder gives a variable index out of an array index.
251
252        // PBe: That makes a lot of sense when problems are really
253        // big. Instances arki000[24].nl spend a lot of time here
254
255        { //if ((nOrigVars_ < THRES_ABT_ORIG) || (index < nOrigVars_)) {
256
257          // if (index == objind) continue; // don't do it on objective function
258
259          improved = 0;
260
261          if ((variables_ [index] -> sign () != expression::AUX_GEQ) &&
262              (X [index] >= Lb (index) + COUENNE_EPS)) {
263
264            Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING,
265                            "------------- tighten left %d-th = x%d @%g [%g,%g]\n", 
266                            i, index, X [index], olb [index], oub [index]);
267
268            // tighten on left
269            if ((improved = fake_tighten (0, index, X, olb, oub, chg_bds, f_chg)) < 0) {
270
271              retval = false;
272              break;
273            }
274          }
275
276          second = 0;
277
278          if (retval && (variables_ [index] -> sign () != expression::AUX_LEQ) &&
279              (X [index] <= Ub (index) - COUENNE_EPS)) {
280
281            Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING,
282                            "------------- tighten right %d-th = x%d @%g [%g,%g]\n", 
283                            i, index, X [index], olb [index], oub [index]);
284
285            // tighten on right
286            if ((second = fake_tighten (1, index, X, olb, oub, chg_bds, f_chg) < 0)) {
287              retval = false;
288              break;
289            }
290          }
291
292          improved += second;
293          nTotImproved += improved;
294        }
295      }
296
297    } while (!maxTimeReached && retval && (improved > THRES_ABT_IMPROVED) && (iter++ < MAX_ABT_ITER));
298
299    // store new valid bounds, or restore old ones if none changed
300    CoinCopyN (olb, ncols, Lb ());
301    CoinCopyN (oub, ncols, Ub ());
302
303    if (Jnlst()->ProduceOutput(J_ITERSUMMARY, J_BOUNDTIGHTENING)) {
304
305      Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING,"------------------\n");
306
307      if (!retval) Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING,
308                                   "Couenne infeasible node from aggressive BT\n");
309
310      int objind = Obj (0) -> Body  () -> Index ();
311
312      if (objind >= 0)
313        Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING,
314                        "-------------\ndone Aggressive BT. Current bound = %g, cutoff = %g, %d vars\n", 
315                        Lb (objind), getCutOff (), ncols);
316
317      if (Jnlst()->ProduceOutput(J_DETAILED, J_BOUNDTIGHTENING))
318        for (int i=0; i<nOrigVars_; i++)
319          printf("   x%02d [%+20g %+20g]  | %+20g\n",
320                 i, Lb (i), Ub (i), X [i]);
321
322      if (Jnlst()->ProduceOutput(J_MOREDETAILED, J_BOUNDTIGHTENING))
323        for (int i=nOrigVars_; i<ncols; i++)
324          printf ("   w%02d [%+20g %+20g]  | %+20g\n", i, Lb (i), Ub (i), X [i]);
325    }
326
327    if (haveNLPsol)
328      delete [] X;
329    delete [] f_chg;   
330  }
331  // else
332  // if ((dist > 1e10) && !retval)
333  //   jnlst_ -> Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING, "TODO: Don't have point for ABT\n");
334
335  delete [] olb;
336  delete [] oub;
337
338  if (info.level <= 0 && !(info.inTree))  {
339    if (!retval) jnlst_ -> Printf (J_ERROR, J_COUENNE, "infeasible\n");
340    else         jnlst_ -> Printf (J_ERROR, J_COUENNE, "%d improved bounds\n", nTotImproved);
341  }
342
343  return retval; // && btCore (psi, cs, chg_bds, babInfo, true); // !!!
344  //return retval && btCore (psi, cs, chg_bds, babInfo, true);
345}
Note: See TracBrowser for help on using the repository browser.