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

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

probing now applied even in absence of NLP solution

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