source: stable/0.2/Couenne/src/bound_tightening/aggressiveBT.cpp @ 159

Last change on this file since 159 was 159, checked in by pbelotti, 11 years ago

created new stable branch 0.2 from trunk (rev. 157)

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