1 /* $Id$
2 *
3 * Name: obbt.cpp
4 * Author: Pietro Belotti
5 * Purpose: Optimality-Based Bound Tightening
6 *
7 * (C) Carnegie-Mellon University, 2006-10.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11 #include "CglCutGenerator.hpp"
12 #include "CouenneCutGenerator.hpp"
13 #include "CouenneProblem.hpp"
14 #include "CouenneProblemElem.hpp"
15 #include "CouenneExprVar.hpp"
16
17 using namespace Ipopt;
18 using namespace Couenne;
19
20 #define OBBT_EPS 1e-3
21
22 // TODO: seems like Clp doesn't like large bounds and crashes on
23 // explicit bounds around 1e200 or so. For now simply use fictitious
24 // bounds around 1e14. Fix.
25
26 int obbt_supplement (const OsiSolverInterface *csi, /// interface to use as a solver
27 int index, /// variable being looked at
28 int sense); /// 1: minimize, -1: maximize
29
30 /// reoptimize and change bound of a variable if needed
obbt_updateBound(OsiSolverInterface * csi,int sense,CouNumber & bound,bool isint)31 static bool obbt_updateBound (OsiSolverInterface *csi, /// interface to use as a solver
32 int sense, /// 1: minimize, -1: maximize
33 CouNumber &bound, /// bound to be updated
34 bool isint) { /// is this variable integer
35
36
37 // NEW: TODO: save min x^\star_i to check at every iteration
38 // (comparison each iteration is less advtgs) when checking if x_i
39 // needs minz for obbt
40
41 //csi -> deleteScaleFactors ();
42 csi -> setDblParam (OsiDualObjectiveLimit, COIN_DBL_MAX);
43 csi -> setDblParam (OsiPrimalObjectiveLimit, (sense==1) ? bound : -bound);
44 //csi -> setObjSense (1); // always minimize, just change the sign of the variable // done in caller
45
46 ////////////////////////////////////////////////////////////////////////
47
48 //csi -> resolve_nobt (); // this is a time-expensive part, be considerate...
49 csi -> resolve (); // this is a time-expensive part, be considerate...
50
51 ////////////////////////////////////////////////////////////////////////
52
53 if (csi -> isProvenOptimal ()) {
54
55 double opt = csi -> getObjValue ();
56
57 if (sense > 0)
58 {if (opt > bound + OBBT_EPS) {bound = (isint ? ceil (opt - COUENNE_EPS) : opt); return true;}}
59 else {if ((opt=-opt) < bound - OBBT_EPS) {bound = (isint ? floor (opt + COUENNE_EPS) : opt); return true;}}
60 }
61
62 return false;
63 }
64
65
66 /// Iteration on one variable
67
obbt_iter(OsiSolverInterface * csi,t_chg_bounds * chg_bds,const CoinWarmStart * warmstart,Bonmin::BabInfo * babInfo,double * objcoe,int sense,int index) const68 int CouenneProblem::obbt_iter (OsiSolverInterface *csi,
69 t_chg_bounds *chg_bds,
70 const CoinWarmStart *warmstart,
71 Bonmin::BabInfo *babInfo,
72 double *objcoe,
73 int sense,
74 int index) const {
75
76 // TODO: do NOT apply OBBT if this is a variable of the form
77 // w2=c*w1, as it suffices to multiply result. More in general, do
78 // not apply if w2 is a unary monotone function of w1. Even more in
79 // general, if w2 is a unary function of w1, apply bound propagation
80 // from w1 to w2 and mark it as exact (depending on whether it is
81 // non-decreasing or non-increasing
82
83 //static int iter = 0;
84
85 // exclude checking known optimal solution if initial bounding box
86 // already excludes it
87
88 CouNumber *knownOptimum = optimum_;
89
90 if (optimum_) {
91
92 for (int i=nVars(); i--; knownOptimum++)
93
94 if (*knownOptimum < Lb (i) ||
95 *knownOptimum > Ub (i)) {
96
97 knownOptimum = NULL;
98 break;
99 }
100
101 if (knownOptimum)
102 knownOptimum -= nVars ();
103 }
104
105 std::set <int> deplist;
106 int deplistsize;
107
108 bool issimple = false;
109
110 exprVar *var = Var (index);
111
112 int
113 objind = Obj (0) -> Body () -> Index (),
114 ncols = csi -> getNumCols (),
115 nImprov = 0;
116
117 if ((var -> Type () == AUX) &&
118 ((deplistsize = var -> Image () -> DepList (deplist, STOP_AT_AUX)) <= 1)) {
119
120 if (!deplistsize) { // funny, the expression is constant...
121
122 CouNumber value = (*(var -> Image ())) ();
123
124 if (csi -> getColLower () [index] < value - COUENNE_EPS) {
125 csi -> setColLower (index, value);
126 chg_bds [index].setLowerBits(t_chg_bounds::CHANGED | t_chg_bounds::EXACT);
127 }
128 else chg_bds [index].setLowerBits(t_chg_bounds::EXACT);
129
130 if (csi -> getColUpper () [index] > value + COUENNE_EPS) {
131 csi -> setColUpper (index, value);
132 chg_bds [index].setUpperBits(t_chg_bounds::CHANGED | t_chg_bounds::EXACT);
133 }
134 else chg_bds [index].setUpperBits(t_chg_bounds::EXACT);
135
136 issimple = true;
137
138 } else { // the expression only depends on one variable, meaning
139 // that bound propagation should be sufficient
140
141 int indInd = *(deplist.begin ());
142
143 // expression *image = var -> Image ();
144 // TODO: write code for monotone functions...
145
146 if // independent variable is exactly bounded in both ways
147 ((((chg_bds [indInd].lower() & t_chg_bounds::EXACT) &&
148 (chg_bds [indInd].upper() & t_chg_bounds::EXACT)) ||
149 // or if this expression is of the form w=cx+d, that is, it
150 // depends on one variable only and it is linear
151 (var -> Image () -> Linearity () <= LINEAR)) &&
152 (var -> sign () == expression::AUX_EQ)) {
153
154 issimple = true;
155
156 CouNumber lb, ub;
157 var -> Image () -> getBounds (lb, ub);
158
159 if (lb < Lb (index)) lb = Lb (index);
160 if (ub > Ub (index)) ub = Ub (index);
161
162 if (var -> isInteger ()) {
163 lb = ceil (lb - COUENNE_EPS);
164 ub = floor (ub + COUENNE_EPS);
165 }
166
167 if (lb > csi -> getColLower () [index] + COUENNE_EPS) {
168 csi -> setColLower (index, lb);
169 Lb (index) = lb;
170 chg_bds [index].setLowerBits(t_chg_bounds::CHANGED | t_chg_bounds::EXACT);
171 } else chg_bds [index].setLowerBits(t_chg_bounds::EXACT);
172
173 if (ub < csi -> getColUpper () [index] - COUENNE_EPS) {
174 csi -> setColUpper (index, ub);
175 Ub (index) = ub;
176 chg_bds [index].setUpperBits(t_chg_bounds::CHANGED | t_chg_bounds::EXACT);
177 } else chg_bds [index].setUpperBits(t_chg_bounds::EXACT);
178 }
179 }
180 }
181
182 // only improve bounds if
183 if (!issimple &&
184 ((Var (index) -> Type () == VAR) || // it is an original variable
185 (Var (index) -> Multiplicity () > 0)) && // or its multiplicity is at least 1
186 (Lb (index) < Ub (index) - COUENNE_EPS) && // in any case, bounds are not equal
187
188 ((index != objind) // this is not the objective
189 // or it is, so we use it for re-solving // TODO: check!
190 || ((sense == 1) && !(chg_bds [index].lower() & t_chg_bounds::EXACT))
191 )) {
192 //((sense==-1) && (psense == MAXIMIZE) && !(chg_bds [index].upper() & t_chg_bounds::EXACT)))) {
193
194 bool isInt = (Var (index) -> isInteger ());
195
196 objcoe [index] = sense;
197
198 csi -> setObjective (objcoe);
199 csi -> setObjSense (1); // minimization
200
201 // TODO: Use something else!
202 #if 0
203 for (int iv=0; iv<csi->getNumCols (); iv++) {
204 if (fabs (csi -> getColLower () [iv]) > 1e7) csi -> setColLower (iv, -1e14);
205 if (fabs (csi -> getColUpper () [iv]) > 1e7) csi -> setColUpper (iv, 1e14);
206 }
207 #endif
208
209 CouNumber &bound =
210 (sense == 1) ?
211 (Lb (index)) :
212 (Ub (index));
213
214 // m{in,ax}imize xi on csi
215
216 /*if (Jnlst()->ProduceOutput(J_MOREVECTOR, J_BOUNDTIGHTENING)) {
217 Jnlst()->Printf(J_MOREVECTOR, J_BOUNDTIGHTENING,
218 "m%simizing x%d [%g,%g] %c= %g\n",
219 (sense==1) ? "in" : "ax", index, Lb (index), Ub (index),
220 (sense==1) ? '>' : '<', bound); fflush (stdout);
221 if (Jnlst()->ProduceOutput(J_MOREMATRIX, J_BOUNDTIGHTENING)) {
222 char fname [20];
223 sprintf (fname, "m%s_w%03d_%03d", (sense == 1) ? "in" : "ax", index, iter);
224 printf ("saving in %s\n", fname);
225 //Jnlst()->Printf(J_MOREVECTOR, J_BOUNDTIGHTENING,"writing %s\n", fname);
226 csi -> writeLp (fname);
227 }
228 }*/
229
230 csi -> setWarmStart (warmstart);
231 //csi -> continuousModel () -> setPerturbation (50);
232
233 /* From ClpSimplex.cpp:
234
235 If you are re-using the same matrix again and again then the
236 setup time to do scaling may be significant. Also you may not
237 want to initialize all values or return all values (especially
238 if infeasible). While an auxiliary model exists it will be
239 faster. If options -1 then model is switched off. Otherwise
240 switched on with following options.
241
242 1 - rhs is constant
243 2 - bounds are constant
244 4 - objective is constant
245 8 - solution in by basis and no djs etc in
246 16 - no duals out (but reduced costs)
247 32 - no output if infeasible
248 */
249
250 //csi -> continuousModel () -> auxiliaryModel (1|8|16|32);
251
252 //Jnlst () -> Printf (J_MATRIX, J_BOUNDTIGHTENING,
253 //"obbt___ index = %d [sen=%d,bd=%g,int=%d]\n",
254 //index, sense, bound, isInt);
255
256 bool has_updated = false;
257
258 if (obbt_updateBound (csi, sense, bound, isInt)) {
259
260 has_updated = true;
261
262 if (knownOptimum) {
263 if (sense == 1) {
264 if (bound > COUENNE_EPS + knownOptimum [index])
265 Jnlst()->Printf(J_STRONGWARNING, J_BOUNDTIGHTENING,
266 "#### OBBT cuts optimum at x%d: lb = %g, opt = %g, new lb = %g\n",
267 index, Lb (index), knownOptimum [index], bound);
268 } else {
269 if (bound < -COUENNE_EPS + knownOptimum [index])
270 Jnlst()->Printf(J_STRONGWARNING, J_BOUNDTIGHTENING,
271 "#### OBBT cuts optimum at x%d: ub = %g, opt = %g, new ub = %g\n",
272 index, Ub (index), knownOptimum [index], bound);
273 }
274 }
275
276 if (sense == 1) {if (bound > Lb (index)) Lb (index) = bound;}
277 else {if (bound < Ub (index)) Ub (index) = bound;}
278
279 // more conservative, only change (and set CHANGED) if improve
280
281 if (sense==1)
282 if (csi -> getColLower () [index] < bound - COUENNE_EPS) {
283 Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING,"l_%d: %g --> %g\n",
284 index, csi -> getColLower () [index], bound);
285 csi -> setColLower (index, bound);
286 chg_bds [index].setLowerBits(t_chg_bounds::CHANGED | t_chg_bounds::EXACT);
287 } else chg_bds [index].setLowerBits(t_chg_bounds::EXACT);
288 else
289 if (csi -> getColUpper () [index] > bound + COUENNE_EPS) {
290 Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING,"u_%d: %g --> %g\n",
291 index, csi -> getColUpper () [index], bound);
292 csi -> setColUpper (index, bound);
293 chg_bds [index].setUpperBits(t_chg_bounds::CHANGED | t_chg_bounds::EXACT);
294 } else chg_bds [index].setUpperBits(t_chg_bounds::EXACT);
295
296 /*
297 if (sense==1) {csi -> setColLower (index, bound); chg_bds [index].lower |= CHANGED | EXACT;}
298 else {csi -> setColUpper (index, bound); chg_bds [index].upper |= CHANGED | EXACT;}
299 */
300
301 #if 0
302 // re-check considering reduced costs (more expensive)
303
304 CouNumber *redcost = NULL;
305
306 // first, compute reduced cost when c = c - e_i, where e_i is
307 // a vector with all zero except a one in position i. This
308 // serves as a base to compute modified reduced cost below.
309
310 for (int j=0; j<ncols; j++)
311 if ((j!=index) && (j!=objind)) {
312
313 // fake a change in the objective function and compute
314 // reduced cost. If resulting vector is all positive
315 // (negative), this solution is also optimal for the
316 // minimization (maximization) of x_j and the corresponding
317 // chg_bds[j].lower (.upper) can be set to EXACT.
318
319 if (!(chg_bds [j].lower & EXACT)) {
320 }
321
322 if (!(chg_bds [j].upper & EXACT)) {
323 }
324 }
325 #endif
326
327 // re-apply bound tightening -- here WITHOUT reduced cost
328 // (first argument =NULL is pointer to solverInterface) as csi
329 // is not our problem
330
331 Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING,
332 " OBBT: x_%d: [%g, %g]\n", index,
333 csi -> getColLower () [index],
334 csi -> getColUpper () [index]);
335
336 // should be faster
337
338 //if (doFBBT_ && !(boundTightening (chg_bds, babInfo))) {
339 if (doFBBT_ && !(btCore (chg_bds))) {
340 Jnlst () -> Printf (J_DETAILED, J_BOUNDTIGHTENING,
341 "node is infeasible after post-OBBT tightening\n");
342 return -1; // tell caller this is infeasible
343 }
344
345 nImprov++;
346 }
347
348 // Check value and bounds of other variables. Do this regardless
349 // of the i-th variable being tightened in this iteration
350
351 const double *sol = csi -> getColSolution ();
352
353 for (int j=0; j<ncols; j++)
354 if ((j!=index) && (j!=objind)) {
355
356 if (sol [j] <= Lb (j) + COUENNE_EPS) {
357
358 // if (!(chg_bds [j].lower() & t_chg_bounds::EXACT) && (!has_updated))
359 // printf ("told you it would be useful: l_i: %g\n", j, Lb (j));
360
361 chg_bds [j].setLowerBits(t_chg_bounds::EXACT);
362 }
363
364 if (sol [j] >= Ub (j) - COUENNE_EPS) {
365
366 // if (!(chg_bds [j].upper() & t_chg_bounds::EXACT) && (!has_updated))
367 // printf ("told you it would be useful: u_i: %g\n", j, Ub (j));
368
369 chg_bds [j].setUpperBits(t_chg_bounds::EXACT);
370 }
371 }
372
373 // check for bounds using dual info
374
375 int result = obbt_supplement (csi, index, sense);
376
377 // if we solved the problem on the objective function's
378 // auxiliary variable (that is, we re-solved the extended
379 // problem), it is worth updating the current point (it will be
380 // used later to generate new cuts).
381
382 // TODO: is it, really? And shouldn't we check the opt sense too?
383 /*
384 if ((objind == index) && (csi -> isProvenOptimal ()) && (sense == 1))
385 update (csi -> getColSolution (), NULL, NULL);
386 */
387 // restore null obj fun
388 objcoe [index] = 0;
389 }
390
391 if (nImprov && jnlst_ -> ProduceOutput (J_ITERSUMMARY, J_BOUNDTIGHTENING)) {
392 Jnlst () -> Printf (J_ITERSUMMARY, J_BOUNDTIGHTENING, "OBBT: tightened ", nImprov);
393 Var (index) -> print ();
394 Jnlst () -> Printf (J_ITERSUMMARY, J_BOUNDTIGHTENING, "\n");
395 }
396
397 return nImprov;
398 }
399