source: trunk/Couenne/src/bound_tightening/boundTightening.cpp @ 112

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

isolate restoreUnused to when there is at least one unused -- takes care of reformulated stockcycle problem (thanks to C. D'Ambrosio for the bug report). Check linear terms before adding a full exprGroup -- need a constructor pilot. Quicker solution feasibility check (with initial integrality check).

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