/home/coin/SVN-release/OS-2.2.0/Couenne/src/bound_tightening/aggressiveBT.cpp

Go to the documentation of this file.
00001 /* $Id: aggressiveBT.cpp 230 2009-07-18 11:42:59Z pbelotti $
00002  *
00003  * Name:    aggressiveBT.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: aggressive bound tightening -- fake bounds in variables to
00006  *          exclude parts of the solution space through fathoming on
00007  *          bounds/infeasibility
00008  *
00009  * (C) Carnegie-Mellon University, 2007-09.
00010  * This file is licensed under the Common Public License (CPL)
00011  */
00012 
00013 #include "BonTNLPSolver.hpp"
00014 #include "BonNlpHeuristic.hpp"
00015 #include "CoinHelperFunctions.hpp"
00016 #include "BonCouenneInfo.hpp"
00017 
00018 #include "CouenneCutGenerator.hpp"
00019 #include "CouenneProblem.hpp"
00020 
00021 #define MAX_ABT_ITER           1  // max # aggressive BT iterations
00022 #define THRES_ABT_IMPROVED     0  // only continue ABT if at least these bounds have improved
00023 #define THRES_ABT_ORIG      1000  // only do ABT on originals if they are more than this 
00024 
00025 static double distanceToBound (int n, const double* xOrig,
00026                                const double* lower, const double* upper) {
00027 
00028   double Xdist = 0.;
00029 
00030   for (int i=0, j=n; j--; i++) {
00031 
00032     CouNumber 
00033       diff, 
00034       xO = xOrig [i];
00035 
00036     if      ((diff = lower [i] - xO) > 0) Xdist += diff;
00037     else if ((diff = xO - upper [i]) > 0) Xdist += diff;
00038   }
00039 
00040   return Xdist;
00041 }
00042 
00043 
00044 // Aggressive Bound Tightening: for each variable, fake new bounds
00045 // [l,b] or [b,u] and apply bound tightening. If the interval is
00046 // fathomed on bounds or on infeasibility, the complementary bound
00047 // interval is a valid tightening.
00048 
00049 bool CouenneProblem::aggressiveBT (Bonmin::OsiTMINLPInterface *nlp,
00050                                    t_chg_bounds *chg_bds, 
00051                                    Bonmin::BabInfo * babInfo) const {
00052 
00053   Jnlst () -> Printf (J_ITERSUMMARY, J_BOUNDTIGHTENING, "Aggressive FBBT\n");
00054 
00055   Bonmin::CouenneInfo* couInfo =
00056     dynamic_cast <Bonmin::CouenneInfo *> (babInfo);
00057 
00058   int  ncols  = nVars ();
00059   bool retval = false;
00060 
00061   CouNumber
00062     *olb = new CouNumber [ncols],
00063     *oub = new CouNumber [ncols];
00064 
00065   // save current bounds
00066   CoinCopyN (Lb (), ncols, olb);
00067   CoinCopyN (Ub (), ncols, oub);
00068 
00069   // Find the solution that is closest to the current bounds
00070   // TODO: Also check obj value
00071   SmartPtr<const Bonmin::CouenneInfo::NlpSolution> closestSol;
00072   double dist = 1e50;
00073 
00074   const std::list<SmartPtr<const Bonmin::CouenneInfo::NlpSolution> >& solList =
00075     couInfo->NlpSolutions();
00076 
00077   for (std::list<SmartPtr<const Bonmin::CouenneInfo::NlpSolution> >::const_iterator 
00078          i = solList.begin();
00079        i != solList.end(); i++) {
00080     assert(nOrigVars_ == (*i)->nVars());
00081     const double thisDist = distanceToBound(nOrigVars_, (*i)->solution(), olb, oub);
00082     if (thisDist < dist) {
00083       closestSol = *i;
00084       dist = thisDist;
00085     }
00086   }
00087 
00088   jnlst_ -> Printf (J_VECTOR, J_BOUNDTIGHTENING, "best dist = %e\n", dist);
00089 
00090   // If this solution is not sufficiently inside the bounds, we solve the NLP now
00091   if (dist > 0.1) { // TODO: Find tolerance
00092 
00093     // find integer initial point /////////////////////////
00094 
00095     int nvars = nVars ();
00096 
00097     double *lower = new double [nvars];
00098     double *upper = new double [nvars];
00099 
00100     CoinFillN (lower, nvars, -COUENNE_INFINITY);
00101     CoinFillN (upper, nvars,  COUENNE_INFINITY);
00102 
00103     CoinCopyN (nlp -> getColLower (), nOrigVars_, lower);
00104     CoinCopyN (nlp -> getColUpper (), nOrigVars_, upper);
00105 
00106     double *Y = new double [nvars];
00107 
00108     CoinFillN (Y,    nvars,      0.);
00109     CoinCopyN (X (), nOrigVars_, Y);
00110 
00111     if (getIntegerCandidate (nlp -> getColSolution (), Y, lower, upper) < 0) {
00112 
00113       jnlst_ -> Printf (J_ITERSUMMARY, J_BOUNDTIGHTENING, "TODO: find NLP point in ABT failed\n");
00114       retval = true;
00115 
00116     } else {
00117 
00119 
00120       nlp -> setColLower    (lower);
00121       nlp -> setColUpper    (upper);
00122       nlp -> setColSolution (Y);
00123 
00124       try {nlp -> initialSolve ();}
00125       catch (Bonmin::TNLPSolver::UnsolvedError *E) {}
00126     
00127       delete [] Y;
00128       delete [] lower;
00129       delete [] upper;
00130 
00131       if (nlp->isProvenOptimal()) {
00132         closestSol = new Bonmin::CouenneInfo::NlpSolution 
00133           (nOrigVars_, nlp->getColSolution(), nlp->getObjValue());
00134         couInfo->addSolution(closestSol);
00135         dist = 0.;
00136       }
00137       else {
00138         jnlst_ -> Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING, "TODO: NLP solve in ABT failed\n");
00139         retval = true;
00140       }
00141     }
00142   }
00143 
00144   if (!retval && (dist < 1e10)) {
00145 
00146     retval = true;
00147 
00148     // X is now the NLP solution, but in a low-dimensional space. We
00149     // have to get the corresponding point in higher dimensional space
00150     // through getAuxs()
00151 
00152     double *X = new double [ncols];
00153     CoinCopyN (closestSol->solution(), nOrigVars_, X);
00154     getAuxs (X);
00155 
00156     // create a new, fictitious, bound bookkeeping structure
00157     t_chg_bounds *f_chg = new t_chg_bounds [ncols];
00158 
00159     if (Jnlst()->ProduceOutput(J_ITERSUMMARY, J_BOUNDTIGHTENING)) {
00160       //    CouNumber cutoff = getCutOff ();
00161       int       objind = Obj (0) -> Body  () -> Index ();
00162       for (int i=0; i<nOrigVars_; i++)
00163         Jnlst()->Printf(J_MOREVECTOR, J_BOUNDTIGHTENING,
00164                         "   %2d %+20g [%+20g %+20g]\n",
00165                         i, X [i], Lb (i), Ub (i));
00166       Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING,
00167                       "-------------\nAggressive BT. Current bound = %g, cutoff = %g, %d vars\n", 
00168                       Lb (objind), getCutOff (), ncols);
00169     }
00170 
00171     int improved, second, iter = 0;
00172 
00173     // Repeatedly fake tightening bounds on both sides of every variable
00174     // to concentrate around current NLP point.
00175     //
00176     // MAX_ABT_ITER is the maximum # of outer cycles. Each call to
00177     // fake_tighten in turn has an iterative algorithm for a
00178     // derivative-free, uni-dimensional optimization problem on a
00179     // monotone function.
00180 
00181     do {
00182 
00183       improved = 0;
00184 
00185       // scan all variables
00186       for (int i=0; i<ncols; i++) {
00187 
00188         if (CoinCpuTime () > maxCpuTime_)
00189           break;
00190 
00191         int index = evalOrder (i);
00192 
00193         if (Var (index) -> Multiplicity () <= 0) 
00194           continue;
00195 
00196         // AW: We only want to do the loop that temporarily changes
00197         // bounds around the NLP solution only for those points from the
00198         // NLP solution (no auxiliary vars)?
00199 
00200         // PBe: if we do want that, index should be initialized as i, as
00201         // evalOrder gives a variable index out of an array index.
00202 
00203         // PBe: That makes a lot of sense when problems are really
00204         // big. Instances arki000[24].nl spend a lot of time here
00205 
00206         if ((nOrigVars_ < THRES_ABT_ORIG) || (index < nOrigVars_)) {
00207 
00208           // if (index == objind) continue; // don't do it on objective function
00209 
00210           Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING,
00211                           "------------- tighten left x%d\n", index);
00212 
00213           // tighten on left
00214           if ((X [index] >= Lb (index) + COUENNE_EPS)
00215               && ((improved = fake_tighten (0, index, X, olb, oub, chg_bds, f_chg)) < 0)) {
00216             retval = false;
00217             break;
00218           }
00219 
00220           second = 0;
00221 
00222           Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING,
00223                           "------------- tighten right x%d\n", index);
00224 
00225           // tighten on right
00226           if ((X [index] <= Ub (index) - COUENNE_EPS)
00227               && ((second = fake_tighten (1, index, X, olb, oub, chg_bds, f_chg) < 0))) {
00228             retval = false;
00229             break;
00230           }
00231 
00232           improved += second;
00233         }
00234       }
00235     } while (retval && (improved > THRES_ABT_IMPROVED) && (iter++ < MAX_ABT_ITER));
00236 
00237     // store new valid bounds, or restore old ones if none changed
00238     CoinCopyN (olb, ncols, Lb ());
00239     CoinCopyN (oub, ncols, Ub ());
00240 
00241     if (Jnlst()->ProduceOutput(J_ITERSUMMARY, J_BOUNDTIGHTENING)) {
00242       Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING,"------------------\n");
00243 
00244       if (!retval) Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING,
00245                                    "Couenne infeasible node from aggressive BT\n");
00246 
00247       int objind = Obj (0) -> Body  () -> Index ();
00248 
00249       Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING,
00250                       "-------------\ndone Aggressive BT. Current bound = %g, cutoff = %g, %d vars\n", 
00251                       Lb (objind), getCutOff (), ncols);
00252 
00253       if (Jnlst()->ProduceOutput(J_DETAILED, J_BOUNDTIGHTENING))
00254         for (int i=0; i<nOrigVars_; i++)
00255           printf("   x%02d [%+20g %+20g]  | %+20g\n",
00256                  i, Lb (i), Ub (i), X [i]);
00257 
00258       if (Jnlst()->ProduceOutput(J_MOREDETAILED, J_BOUNDTIGHTENING))
00259         for (int i=nOrigVars_; i<ncols; i++)
00260           printf ("   w%02d [%+20g %+20g]  | %+20g\n", i, Lb (i), Ub (i), X [i]);
00261     }
00262 
00263     delete [] X;
00264     delete [] f_chg;
00265     
00266   } else
00267 
00268     if ((dist > 1e10) && !retval)
00269       jnlst_ -> Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING, "TODO: Don't have point for ABT\n");
00270 
00271   delete [] olb;
00272   delete [] oub;
00273 
00274   return retval;// && btCore (psi, cs, chg_bds, babInfo, true); // !!!
00275   //return retval && btCore (psi, cs, chg_bds, babInfo, true);
00276 }

Generated on Thu Aug 5 03:02:56 2010 by  doxygen 1.4.7