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

Last change on this file since 307 was 307, checked in by pbelotti, 12 years ago

introducing semi-auxiliaries, i.e., variables defined with an inequality constraint. Instead of w := f(x), they can be defined as w >= f(x) or w <= f(x). Saves on number of auxiliary variables and constraints in the linearization. NOTE: Not totally stable yet, optima can be wrong.

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