1 /* $Id$
2 *
3 * Name: boundTightening.cpp
4 * Author: Pietro Belotti
5 * Purpose: tighten bounds prior to convexification cuts
6 *
7 * (C) Carnegie-Mellon University, 2006-10.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11 #include "CouenneCutGenerator.hpp"
12 #include "CouenneProblem.hpp"
13 #include "CouenneExprVar.hpp"
14 #include "CouenneProblemElem.hpp"
15 #include "CouenneBTPerfIndicator.hpp"
16 #include "BonBabInfos.hpp"
17 #include "BonCbc.hpp"
18
19 // max # bound tightening iterations
20 #define THRES_IMPROVED 0
21
22 using namespace Couenne;
23
24 // core of the bound tightening procedure
25
btCore(t_chg_bounds * chg_bds) const26 bool CouenneProblem::btCore (t_chg_bounds *chg_bds) const {
27
28 fbbtReachedIterLimit_ = false;
29
30 bool null_chg_bds = (NULL == chg_bds);
31
32 if (!chg_bds) {
33
34 chg_bds = new t_chg_bounds [nVars ()];
35
36 for (int i=0; i < nVars (); i++)
37
38 if (Var (i) -> Multiplicity () > 0) {
39
40 chg_bds [i].setLower (t_chg_bounds::CHANGED);
41 chg_bds [i].setUpper (t_chg_bounds::CHANGED);
42 }
43 }
44
45 // Bound propagation and implied bounds ////////////////////
46
47 int ntightened = 0,
48 nbwtightened = 0,
49 niter = 0;
50
51 bool first = true;
52
53 installCutOff ();
54
55 // check if bt cuts the optimal solution -- now and after bound tightening
56 bool contains_optimum = false;
57
58 if (optimum_ != NULL) {
59 contains_optimum = true;
60 for (int i=0; i < nVars (); i++)
61 if ((optimum_ [i] < Lb (i) * (1 - COUENNE_EPS) - COUENNE_EPS) ||
62 (optimum_ [i] > Ub (i) * (1 + COUENNE_EPS) + COUENNE_EPS)) {
63 /*printf ("won't check BT: %d [%g,%g] (%g) -- %g\n",
64 i, Lb (i), Ub (i), optimum_ [i],
65 CoinMax (- optimum_ [i] + (Lb (i) * (1 - COUENNE_EPS) - COUENNE_EPS),
66 optimum_ [i] - (Ub (i) * (1 + COUENNE_EPS) + COUENNE_EPS)));*/
67 contains_optimum = false;
68 break;
69 }
70 }
71
72 if (max_fbbt_iter_)
73
74 do {
75
76 if (CoinCpuTime () > maxCpuTime_)
77 break;
78
79 // propagate bounds to auxiliary variables
80
81 // if ((nbwtightened > 0) || (ntightened > 0))
82 ntightened = ((nbwtightened > 0) || first) ?
83 tightenBounds (chg_bds) : 0;
84
85 // implied bounds. Call also at the beginning, as some common
86 // expression may have non-propagated bounds
87
88 // if last call didn't signal infeasibility
89 nbwtightened = ((ntightened > 0) || ((ntightened==0) && first)) ? impliedBounds (chg_bds) : 0;
90
91 if (first)
92 first = false;
93
94 if ((ntightened < 0) || (nbwtightened < 0)) {
95 Jnlst () -> Printf (Ipopt::J_ITERSUMMARY, J_BOUNDTIGHTENING, "infeasible BT\n");
96 if (null_chg_bds)
97 delete [] chg_bds;
98 return false;
99 }
100
101 // continue if EITHER procedures gave (positive) results, as
102 // expression structure is not a tree.
103
104 if (contains_optimum) {
105 for (int i=0; i<nVars (); i++)
106 if ((optimum_[i] < Lb(i) - COUENNE_EPS * (1. + CoinMin (fabs(optimum_ [i]), fabs (Lb(i))))) ||
107 (optimum_[i] > Ub(i) + COUENNE_EPS * (1. + CoinMin (fabs(optimum_ [i]), fabs (Ub(i)))))) {
108 printf ("bound tightening CUTS optimum: x%d [%e,%e] val = %e, violation = %e\n",
109 i, Lb (i), Ub (i), optimum_ [i],
110 CoinMax (- optimum_ [i] + Lb (i),
111 optimum_ [i] - Ub (i)));
112 contains_optimum = false;
113 }
114 }
115
116 // double width = 0.;
117
118 // int
119 // ninf = 0,
120 // ndinf = 0;
121
122 // for (int i=nOrigVars (); i--;)
123
124 // if ((Lb (i) < -COUENNE_INFINITY/1e3) &&
125 // (Ub (i) > COUENNE_INFINITY/1e3)) ndinf++;
126 // else if ((Lb (i) < -COUENNE_INFINITY/1e3) ||
127 // (Ub (i) > COUENNE_INFINITY/1e3)) ninf++;
128 // else width += Ub (i) - Lb (i);
129
130 // printf ("pass %5d: %5d fwd, %5d bwd, %5d inf, %5d double inf, width: %g\n",
131 // niter, ntightened, nbwtightened, ninf, ndinf, width);
132
133 } while (((ntightened > 0) || (nbwtightened > 0)) &&
134 (ntightened + nbwtightened > THRES_IMPROVED) &&
135 ((max_fbbt_iter_ < 0) || (niter++ < max_fbbt_iter_)));
136
137 if (null_chg_bds)
138 delete [] chg_bds;
139
140 fbbtReachedIterLimit_ = ((max_fbbt_iter_ > 0) && (niter >= max_fbbt_iter_));
141
142 // TODO: limit should depend on number of constraints, that is,
143 // bound transmission should be documented and the cycle should stop
144 // as soon as no new constraint subgraph has benefited from bound
145 // transmission.
146 //
147 // BT should be more of a graph spanning procedure that moves from
148 // one vertex to another when either tightening procedure has given
149 // some result. This should save some time especially
150 // w.r.t. applying implied bounds to ALL expressions just because
151 // one single propagation was found.
152
153 for (int i = 0, j = nVars (); j--; i++)
154 if (Var (i) -> Multiplicity () > 0) {
155
156 // final test
157 if ((Lb (i) > Ub (i) + COUENNE_BOUND_PREC * (1 + CoinMin (fabs (Lb (i)), fabs (Ub (i))))) ||
158 (Ub (i) < - MAX_BOUND) ||
159 (Lb (i) > MAX_BOUND)) {
160
161 Jnlst () -> Printf (Ipopt::J_ITERSUMMARY, J_BOUNDTIGHTENING, "final test: infeasible BT\n");
162 return false;
163 }
164
165 // sanity check. Ipopt gives an exception when Lb (i) is above Ub (i)
166 if (Lb (i) > Ub (i)) {
167 CouNumber swap = Lb (i);
168 Lb (i) = Ub (i);
169 Ub (i) = swap;
170 }
171 }
172
173 return true;
174 }
175
176
177 /// procedure to strengthen variable bounds. Return false if problem
178 /// turns out to be infeasible with given bounds, true otherwise.
179
boundTightening(t_chg_bounds * chg_bds,const CglTreeInfo info,Bonmin::BabInfo * babInfo) const180 bool CouenneProblem::boundTightening (t_chg_bounds *chg_bds,
181 const CglTreeInfo info,
182 Bonmin::BabInfo * babInfo) const {
183
184 double startTime = CoinCpuTime ();
185
186 FBBTperfIndicator_ -> setOldBounds (Lb (), Ub ());
187
188 //
189 // #define SMALL_BOUND 1e4
190 // for (int i=nOrigVars (); i--;) {
191 // if (Lb (i) < -SMALL_BOUND) Lb (i) = -SMALL_BOUND;
192 // if (Ub (i) > SMALL_BOUND) Ub (i) = SMALL_BOUND;
193 // }
194
195 Jnlst () -> Printf (Ipopt::J_ITERSUMMARY, J_BOUNDTIGHTENING,
196 "Feasibility-based Bound Tightening\n");
197
198 int objInd = Obj (0) -> Body () -> Index ();
199
200 /////////////////////// MIP bound management /////////////////////////////////
201
202 if ((objInd >= 0) && babInfo && (babInfo -> babPtr ())) {
203
204 CouNumber UB = babInfo -> babPtr () -> model (). getObjValue(),
205 LB = babInfo -> babPtr () -> model (). getBestPossibleObjValue (),
206 primal0 = Ub (objInd),
207 dual0 = Lb (objInd);
208
209 // Bonmin assumes minimization. Hence, primal (dual) is an UPPER
210 // (LOWER) bound.
211
212 if ((UB < COUENNE_INFINITY) &&
213 (UB < primal0 - COUENNE_EPS)) { // update primal bound (MIP)
214
215 Ub (objInd) = UB;
216 chg_bds [objInd].setUpper(t_chg_bounds::CHANGED);
217 }
218
219 if ((LB > - COUENNE_INFINITY) &&
220 (LB > dual0 + COUENNE_EPS)) { // update dual bound
221 Lb (objInd) = LB;
222 chg_bds [objInd].setLower(t_chg_bounds::CHANGED);
223 }
224 }
225
226 bool retval = btCore (chg_bds);
227
228 FBBTperfIndicator_ -> update (Lb (), Ub (), info.level);
229 FBBTperfIndicator_ -> addToTimer (CoinCpuTime () - startTime);
230
231 return retval;
232
233 //printf ("Total cpu time = %e\n", CoinCpuTime () - startTime);
234 //exit (-1);
235 //return retval;
236 }
237
238
239 /// reduced cost bound tightening
redCostBT(const OsiSolverInterface * psi,t_chg_bounds * chg_bds) const240 int CouenneProblem::redCostBT (const OsiSolverInterface *psi,
241 t_chg_bounds *chg_bds) const {
242
243 int
244 nchanges = 0,
245 objind = Obj (0) -> Body () -> Index ();
246
247 if (objind < 0)
248 return 0;
249
250 CouNumber
251 UB = getCutOff (), //babInfo -> babPtr () -> model (). getObjValue(),
252 LB = Lb (objind); //babInfo -> babPtr () -> model (). getBestPossibleObjValue ();
253
254 //////////////////////// Reduced cost bound tightening //////////////////////
255
256 if ((LB > -COUENNE_INFINITY) &&
257 (UB < COUENNE_INFINITY)) {
258
259 const double
260 *X = psi -> getColSolution (),
261 *L = psi -> getColLower (),
262 *U = psi -> getColUpper (),
263 *RC = psi -> getReducedCost ();
264
265 if (jnlst_ -> ProduceOutput (Ipopt::J_MATRIX, J_BOUNDTIGHTENING)) {
266 printf ("REDUCED COST BT (LB=%g, UB=%g):\n", LB, UB);
267 for (int i=0, j=0; i < nVars (); i++) {
268 if (Var (i) -> Multiplicity () <= 0)
269 continue;
270 printf ("%3d %7e [%7e %7e] c %7e ", i, X [i], L [i], U [i], RC [i]);
271 if (!(++j % 3))
272 printf ("\n");
273 }
274 printf ("-----------\n");
275 }
276
277 int ncols = psi -> getNumCols ();
278
279 for (int i=0; i<ncols; i++)
280 if ((i != objind) &&
281 (Var (i) -> Multiplicity () > 0)) {
282
283 CouNumber
284 x = X [i],
285 l = L [i],
286 u = U [i],
287 rc = RC [i];
288
289 #define CLOSE_TO_BOUNDS 1.e-15
290
291 if ((fabs (rc) < CLOSE_TO_BOUNDS) ||
292 (fabs (l-u) < CLOSE_TO_BOUNDS)) // no need to check
293 continue;
294
295 bool isInt = Var (i) -> isInteger ();
296
297 if ((rc >= 0.) &&
298 (fabs (x-l) <= CLOSE_TO_BOUNDS)) {
299
300 if (LB + (u-l)*rc > UB) {
301
302 CouNumber newUb = l + (UB-LB) / rc; // which is surely < u
303 newUb = !isInt ? newUb : floor (newUb + COUENNE_EPS);
304
305 if (newUb < Ub (i)) {
306
307 Ub (i) = newUb;
308
309 nchanges++;
310 chg_bds [i].setLower (t_chg_bounds::CHANGED);
311 }
312 }
313
314 } else if ((rc <= 0.) &&
315 (fabs (x-u) <= CLOSE_TO_BOUNDS)) {
316
317 if (LB - (u-l) * rc > UB) {
318
319 CouNumber newLb = u + (UB-LB) / rc; // recall rc < 0 here
320
321 newLb = !isInt ? newLb : ceil (newLb - COUENNE_EPS);
322
323 if (newLb > Lb (i)) {
324
325 Lb (i) = newLb;
326
327 nchanges++;
328 chg_bds [i].setUpper (t_chg_bounds::CHANGED);
329 }
330 }
331 }
332 }
333
334 if (jnlst_ -> ProduceOutput (Ipopt::J_MATRIX, J_BOUNDTIGHTENING)) {
335 printf ("AFTER reduced cost bt:\n");
336 for (int i=0, j=0; i < nVars (); ++i) {
337 if (Var (i) -> Multiplicity () <= 0)
338 continue;
339 printf ("%3d [%7e %7e] ", i, Lb (i), Ub (i));
340 if (!(++j % 4))
341 printf ("\n");
342 }
343 printf ("-----------\n");
344 }
345 }
346
347 return nchanges;
348 }
349