/home/coin/SVN-release/OS-2.1.1/Couenne/src/bound_tightening/boundTightening.cpp

Go to the documentation of this file.
00001 /* $Id: boundTightening.cpp 151 2009-06-16 03:31:06Z pbelotti $
00002  *
00003  * Name:    boundTightening.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: tighten bounds prior to convexification cuts
00006  *
00007  * (C) Carnegie-Mellon University, 2006. 
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 #include "CouenneCutGenerator.hpp"
00012 #include "CouenneProblem.hpp"
00013 #include "BonBabInfos.hpp"
00014 #include "BonCbc.hpp"
00015 
00016 // max # bound tightening iterations
00017 #define MAX_BT_ITER 3
00018 #define THRES_IMPROVED 0
00019 
00020 
00021 // core of the bound tightening procedure
00022 
00023 bool CouenneProblem::btCore (t_chg_bounds *chg_bds) const {
00024 
00025   // Bound propagation and implied bounds ////////////////////
00026 
00027   int   ntightened = 0,
00028       nbwtightened = 0,
00029       niter = 0;
00030 
00031   bool first = true;
00032 
00033   installCutOff ();
00034 
00035   // check if bt cuts the optimal solution -- now and after bound tightening
00036   bool contains_optimum = false;
00037 
00038   if (optimum_ != NULL) {
00039     contains_optimum = true;
00040     for (int i=0; i < nVars (); i++)
00041       if ((optimum_ [i] < Lb (i) * (1 - COUENNE_EPS) - COUENNE_EPS) ||
00042           (optimum_ [i] > Ub (i) * (1 + COUENNE_EPS) + COUENNE_EPS)) {
00043         /*printf ("won't check BT: %d [%g,%g] (%g) -- %g\n", 
00044                 i, Lb (i), Ub (i), optimum_ [i],
00045                 CoinMax (- optimum_ [i] + (Lb (i) * (1 - COUENNE_EPS) - COUENNE_EPS),
00046                 optimum_ [i] - (Ub (i) * (1 + COUENNE_EPS) + COUENNE_EPS)));*/
00047         contains_optimum = false;
00048         break;
00049       }
00050   }
00051 
00052   do {
00053 
00054     if (CoinCpuTime () > maxCpuTime_)
00055       break;
00056 
00057     // propagate bounds to auxiliary variables
00058 
00059     //    if ((nbwtightened > 0) || (ntightened > 0))
00060     ntightened = ((nbwtightened > 0) || first) ? 
00061       tightenBounds (chg_bds) : 0;
00062 
00063     // implied bounds. Call also at the beginning, as some common
00064     // expression may have non-propagated bounds
00065 
00066     // if last call didn't signal infeasibility
00067     nbwtightened = ((ntightened > 0) || ((ntightened==0) && first)) ? impliedBounds (chg_bds) : 0;
00068 
00069     if (first)
00070       first = false;
00071 
00072     if ((ntightened < 0) || (nbwtightened < 0)) {
00073       Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING, "infeasible BT\n");
00074       return false;
00075     }
00076 
00077     // continue if EITHER procedures gave (positive) results, as
00078     // expression structure is not a tree.
00079 
00080     if (contains_optimum) {
00081       for (int i=0; i<nVars (); i++)
00082         if ((optimum_[i] < Lb(i) - COUENNE_EPS * (1. + CoinMin (fabs(optimum_ [i]), fabs (Lb(i))))) ||
00083             (optimum_[i] > Ub(i) + COUENNE_EPS * (1. + CoinMin (fabs(optimum_ [i]), fabs (Ub(i)))))) {
00084           printf ("bound tightening FAIL: %d [%e,%e] (%e) -- %e\n", 
00085                   i, Lb (i), Ub (i), optimum_ [i],
00086                   CoinMax (- optimum_ [i] + Lb (i),
00087                            optimum_ [i] - Ub (i)));
00088           contains_optimum = false;
00089         }
00090     }
00091 
00092   } while (((ntightened > 0) || (nbwtightened > 0)) && 
00093            (ntightened + nbwtightened > THRES_IMPROVED) &&
00094            (niter++ < MAX_BT_ITER));
00095 
00096   // TODO: limit should depend on number of constraints, that is,
00097   // bound transmission should be documented and the cycle should stop
00098   // as soon as no new constraint subgraph has benefited from bound
00099   // transmission. 
00100   //
00101   // BT should be more of a graph spanning procedure that moves from
00102   // one vertex to another when either tightening procedure has given
00103   // some result. This should save some time especially
00104   // w.r.t. applying implied bounds to ALL expressions just because
00105   // one single propagation was found.
00106 
00107   for (int i = 0, j = nVars (); j--; i++) 
00108     if (Var (i) -> Multiplicity () > 0) {
00109 
00110       // final test 
00111       if ((Lb (i) > Ub (i) + COUENNE_EPS * (1 + CoinMin (fabs (Lb (i)), fabs (Ub (i))))) || 
00112           (Ub (i) < - MAX_BOUND) ||
00113           (Lb (i) >   MAX_BOUND)) {
00114 
00115         Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING, "final test: infeasible BT\n");
00116         return false;
00117       }
00118 
00119       // sanity check. Ipopt gives an exception when Lb (i) is above Ub (i)
00120       if (Lb (i) > Ub (i)) {
00121         CouNumber swap = Lb (i);
00122         Lb (i) = Ub (i);
00123         Ub (i) = swap;
00124       }
00125     }
00126 
00127   return true;
00128 }
00129 
00130 
00133 
00134 bool CouenneProblem::boundTightening (t_chg_bounds *chg_bds, 
00135                                       Bonmin::BabInfo * babInfo) const {
00136 
00137   Jnlst()->Printf (J_ITERSUMMARY, J_BOUNDTIGHTENING,
00138                    "Feasibility-based Bound Tightening\n");
00139 
00140   int objInd = Obj (0) -> Body () -> Index ();
00141 
00143 
00144   if ((objInd >= 0) && babInfo && (babInfo -> babPtr ())) {
00145 
00146     CouNumber UB      = babInfo  -> babPtr () -> model (). getObjValue(),
00147               LB      = babInfo  -> babPtr () -> model (). getBestPossibleObjValue (),
00148               primal0 = Ub (objInd), 
00149               dual0   = Lb (objInd);
00150 
00151     // Bonmin assumes minimization. Hence, primal (dual) is an UPPER
00152     // (LOWER) bound.
00153     
00154     if ((UB < COUENNE_INFINITY) && 
00155         (UB < primal0 - COUENNE_EPS)) { // update primal bound (MIP)
00156 
00157       Ub (objInd) = UB;
00158       chg_bds [objInd].setUpper(t_chg_bounds::CHANGED);
00159     }
00160 
00161     if ((LB > - COUENNE_INFINITY) && 
00162         (LB > dual0 + COUENNE_EPS)) { // update dual bound
00163       Lb (objInd) = LB;
00164       chg_bds [objInd].setLower(t_chg_bounds::CHANGED);
00165     }
00166   }
00167 
00168   return btCore (chg_bds);
00169 }
00170 
00171 
00173 int CouenneProblem::redCostBT (const OsiSolverInterface *psi,
00174                                t_chg_bounds *chg_bds) const {
00175   int 
00176     nchanges = 0,
00177     objind   = Obj (0) -> Body () -> Index ();
00178 
00179   assert (objind >= 0);
00180 
00181   CouNumber
00182     UB = getCutOff (), //babInfo -> babPtr () -> model (). getObjValue(), // todo: get cutoff
00183     LB = Lb (objind);  //babInfo -> babPtr () -> model (). getBestPossibleObjValue (); // todo:  w_0^l
00184 
00186 
00187   if ((LB > -COUENNE_INFINITY) && 
00188       (UB <  COUENNE_INFINITY)) {
00189 
00190     const double 
00191       *X  = psi -> getColSolution (),
00192       *L  = psi -> getColLower    (),
00193       *U  = psi -> getColUpper    (),
00194       *RC = psi -> getReducedCost ();
00195 
00196     if (jnlst_ -> ProduceOutput (J_MATRIX, J_BOUNDTIGHTENING)) {
00197       printf ("REDUCED COST BT:\n");
00198       for (int i=0; i < nVars (); i++) 
00199         printf ("%3d %10e [%10e %10e] rc %10e\n", i, X [i], L [i], U [i], RC [i]);
00200       printf ("-----------\n");
00201     }
00202 
00203     int ncols = psi -> getNumCols ();
00204 
00205     for (int i=0; i<ncols; i++)
00206       if ((i != objind) && 
00207           (Var (i) -> Multiplicity () > 0)) {
00208 
00209         CouNumber
00210           x  = X  [i],
00211           l  = L  [i],
00212           u  = U  [i],
00213           rc = RC [i];
00214 
00215         if ((rc < COUENNE_EPS) || (l==u)) // no need to check
00216           continue;
00217 
00218         bool isInt = Var (i) -> isInteger ();
00219 
00220         if (x == l) {
00221           if (LB + (u-l)*rc > UB) {
00222             //printf ("ub [%d]: %g ", i, Ub (i));
00223             Ub (i) = l + (UB-LB) / rc;
00224             if (isInt) 
00225               Ub (i) = floor (Ub (i) + COUENNE_EPS);
00226             //printf ("--> %g\n", Ub (i));
00227             nchanges++;
00228             chg_bds [i].setLower(t_chg_bounds::CHANGED);
00229           }
00230         } else if (x == u) {
00231           if (LB + (u-l) * rc > UB) {
00232             //printf ("lb [%d]: %g ", i, Lb (i));
00233             Lb (i) = u - (UB-LB) / rc;
00234             if (isInt) 
00235               Lb (i) = ceil (Lb (i) - COUENNE_EPS);
00236             //printf ("--> %g\n", Lb (i));
00237             nchanges++;
00238             chg_bds [i].setUpper(t_chg_bounds::CHANGED);
00239           }
00240         }
00241       }
00242 
00243     /*printf ("AFTER reduced cost bt:\n");
00244       for (int i=0; i < nVars (); i++) 
00245         printf ("%3d [%10e %10e]\n", i, Lb (i), Ub (i));
00246         printf ("-----------\n");*/
00247   }
00248 
00249   return nchanges;
00250 }

Generated on Mon May 3 03:05:18 2010 by  doxygen 1.4.7