source: stable/0.2/Couenne/src/bound_tightening/boundTightening.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: 7.1 KB
Line 
1/* $Id: boundTightening.cpp 151 2009-06-16 03:31:06Z pbelotti $
2 *
3 * Name:    boundTightening.cpp
4 * Author:  Pietro Belotti
5 * Purpose: tighten bounds prior to convexification cuts
6 *
7 * (C) Carnegie-Mellon University, 2006.
8 * This file is licensed under the Common Public License (CPL)
9 */
10
11#include "CouenneCutGenerator.hpp"
12#include "CouenneProblem.hpp"
13#include "BonBabInfos.hpp"
14#include "BonCbc.hpp"
15
16// max # bound tightening iterations
17#define MAX_BT_ITER 3
18#define THRES_IMPROVED 0
19
20
21// core of the bound tightening procedure
22
23bool CouenneProblem::btCore (t_chg_bounds *chg_bds) const {
24
25  // Bound propagation and implied bounds ////////////////////
26
27  int   ntightened = 0,
28      nbwtightened = 0,
29      niter = 0;
30
31  bool first = true;
32
33  installCutOff ();
34
35  // check if bt cuts the optimal solution -- now and after bound tightening
36  bool contains_optimum = false;
37
38  if (optimum_ != NULL) {
39    contains_optimum = true;
40    for (int i=0; i < nVars (); i++)
41      if ((optimum_ [i] < Lb (i) * (1 - COUENNE_EPS) - COUENNE_EPS) ||
42          (optimum_ [i] > Ub (i) * (1 + COUENNE_EPS) + COUENNE_EPS)) {
43        /*printf ("won't check BT: %d [%g,%g] (%g) -- %g\n",
44                i, Lb (i), Ub (i), optimum_ [i],
45                CoinMax (- optimum_ [i] + (Lb (i) * (1 - COUENNE_EPS) - COUENNE_EPS),
46                optimum_ [i] - (Ub (i) * (1 + COUENNE_EPS) + COUENNE_EPS)));*/
47        contains_optimum = false;
48        break;
49      }
50  }
51
52  do {
53
54    if (CoinCpuTime () > maxCpuTime_)
55      break;
56
57    // propagate bounds to auxiliary variables
58
59    //    if ((nbwtightened > 0) || (ntightened > 0))
60    ntightened = ((nbwtightened > 0) || first) ? 
61      tightenBounds (chg_bds) : 0;
62
63    // implied bounds. Call also at the beginning, as some common
64    // expression may have non-propagated bounds
65
66    // if last call didn't signal infeasibility
67    nbwtightened = ((ntightened > 0) || ((ntightened==0) && first)) ? impliedBounds (chg_bds) : 0;
68
69    if (first)
70      first = false;
71
72    if ((ntightened < 0) || (nbwtightened < 0)) {
73      Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING, "infeasible BT\n");
74      return false;
75    }
76
77    // continue if EITHER procedures gave (positive) results, as
78    // expression structure is not a tree.
79
80    if (contains_optimum) {
81      for (int i=0; i<nVars (); i++)
82        if ((optimum_[i] < Lb(i) - COUENNE_EPS * (1. + CoinMin (fabs(optimum_ [i]), fabs (Lb(i))))) ||
83            (optimum_[i] > Ub(i) + COUENNE_EPS * (1. + CoinMin (fabs(optimum_ [i]), fabs (Ub(i)))))) {
84          printf ("bound tightening FAIL: %d [%e,%e] (%e) -- %e\n", 
85                  i, Lb (i), Ub (i), optimum_ [i],
86                  CoinMax (- optimum_ [i] + Lb (i),
87                           optimum_ [i] - Ub (i)));
88          contains_optimum = false;
89        }
90    }
91
92  } while (((ntightened > 0) || (nbwtightened > 0)) && 
93           (ntightened + nbwtightened > THRES_IMPROVED) &&
94           (niter++ < MAX_BT_ITER));
95
96  // TODO: limit should depend on number of constraints, that is,
97  // bound transmission should be documented and the cycle should stop
98  // as soon as no new constraint subgraph has benefited from bound
99  // transmission.
100  //
101  // BT should be more of a graph spanning procedure that moves from
102  // one vertex to another when either tightening procedure has given
103  // some result. This should save some time especially
104  // w.r.t. applying implied bounds to ALL expressions just because
105  // one single propagation was found.
106
107  for (int i = 0, j = nVars (); j--; i++) 
108    if (Var (i) -> Multiplicity () > 0) {
109
110      // final test
111      if ((Lb (i) > Ub (i) + COUENNE_EPS * (1 + CoinMin (fabs (Lb (i)), fabs (Ub (i))))) || 
112          (Ub (i) < - MAX_BOUND) ||
113          (Lb (i) >   MAX_BOUND)) {
114
115        Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING, "final test: infeasible BT\n");
116        return false;
117      }
118
119      // sanity check. Ipopt gives an exception when Lb (i) is above Ub (i)
120      if (Lb (i) > Ub (i)) {
121        CouNumber swap = Lb (i);
122        Lb (i) = Ub (i);
123        Ub (i) = swap;
124      }
125    }
126
127  return true;
128}
129
130
131/// procedure to strengthen variable bounds. Return false if problem
132/// turns out to be infeasible with given bounds, true otherwise.
133
134bool CouenneProblem::boundTightening (t_chg_bounds *chg_bds, 
135                                      Bonmin::BabInfo * babInfo) const {
136
137  Jnlst()->Printf (J_ITERSUMMARY, J_BOUNDTIGHTENING,
138                   "Feasibility-based Bound Tightening\n");
139
140  int objInd = Obj (0) -> Body () -> Index ();
141
142  /////////////////////// MIP bound management /////////////////////////////////
143
144  if ((objInd >= 0) && babInfo && (babInfo -> babPtr ())) {
145
146    CouNumber UB      = babInfo  -> babPtr () -> model (). getObjValue(),
147              LB      = babInfo  -> babPtr () -> model (). getBestPossibleObjValue (),
148              primal0 = Ub (objInd), 
149              dual0   = Lb (objInd);
150
151    // Bonmin assumes minimization. Hence, primal (dual) is an UPPER
152    // (LOWER) bound.
153   
154    if ((UB < COUENNE_INFINITY) && 
155        (UB < primal0 - COUENNE_EPS)) { // update primal bound (MIP)
156
157      Ub (objInd) = UB;
158      chg_bds [objInd].setUpper(t_chg_bounds::CHANGED);
159    }
160
161    if ((LB > - COUENNE_INFINITY) && 
162        (LB > dual0 + COUENNE_EPS)) { // update dual bound
163      Lb (objInd) = LB;
164      chg_bds [objInd].setLower(t_chg_bounds::CHANGED);
165    }
166  }
167
168  return btCore (chg_bds);
169}
170
171
172/// reduced cost bound tightening
173int CouenneProblem::redCostBT (const OsiSolverInterface *psi,
174                               t_chg_bounds *chg_bds) const {
175  int 
176    nchanges = 0,
177    objind   = Obj (0) -> Body () -> Index ();
178
179  assert (objind >= 0);
180
181  CouNumber
182    UB = getCutOff (), //babInfo -> babPtr () -> model (). getObjValue(), // todo: get cutoff
183    LB = Lb (objind);  //babInfo -> babPtr () -> model (). getBestPossibleObjValue (); // todo:  w_0^l
184
185  //////////////////////// Reduced cost bound tightening //////////////////////
186
187  if ((LB > -COUENNE_INFINITY) && 
188      (UB <  COUENNE_INFINITY)) {
189
190    const double 
191      *= psi -> getColSolution (),
192      *= psi -> getColLower    (),
193      *= psi -> getColUpper    (),
194      *RC = psi -> getReducedCost ();
195
196    if (jnlst_ -> ProduceOutput (J_MATRIX, J_BOUNDTIGHTENING)) {
197      printf ("REDUCED COST BT:\n");
198      for (int i=0; i < nVars (); i++) 
199        printf ("%3d %10e [%10e %10e] rc %10e\n", i, X [i], L [i], U [i], RC [i]);
200      printf ("-----------\n");
201    }
202
203    int ncols = psi -> getNumCols ();
204
205    for (int i=0; i<ncols; i++)
206      if ((i != objind) && 
207          (Var (i) -> Multiplicity () > 0)) {
208
209        CouNumber
210          x  = X  [i],
211          l  = L  [i],
212          u  = U  [i],
213          rc = RC [i];
214
215        if ((rc < COUENNE_EPS) || (l==u)) // no need to check
216          continue;
217
218        bool isInt = Var (i) -> isInteger ();
219
220        if (x == l) {
221          if (LB + (u-l)*rc > UB) {
222            //printf ("ub [%d]: %g ", i, Ub (i));
223            Ub (i) = l + (UB-LB) / rc;
224            if (isInt) 
225              Ub (i) = floor (Ub (i) + COUENNE_EPS);
226            //printf ("--> %g\n", Ub (i));
227            nchanges++;
228            chg_bds [i].setLower(t_chg_bounds::CHANGED);
229          }
230        } else if (x == u) {
231          if (LB + (u-l) * rc > UB) {
232            //printf ("lb [%d]: %g ", i, Lb (i));
233            Lb (i) = u - (UB-LB) / rc;
234            if (isInt) 
235              Lb (i) = ceil (Lb (i) - COUENNE_EPS);
236            //printf ("--> %g\n", Lb (i));
237            nchanges++;
238            chg_bds [i].setUpper(t_chg_bounds::CHANGED);
239          }
240        }
241      }
242
243    /*printf ("AFTER reduced cost bt:\n");
244      for (int i=0; i < nVars (); i++)
245        printf ("%3d [%10e %10e]\n", i, Lb (i), Ub (i));
246        printf ("-----------\n");*/
247  }
248
249  return nchanges;
250}
Note: See TracBrowser for help on using the repository browser.