source: stable/0.1/Couenne/src/bound_tightening/aggressiveBT.cpp @ 134

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

merged trunk revision with Ipopt exception handling

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