source: trunk/Couenne/src/bound_tightening/obbt_iter.cpp @ 318

Last change on this file since 318 was 318, checked in by pbelotti, 12 years ago

fixed bug that gave 2.982 on minlplib/alan -- wrong OBBT implication on univariate (semi)-aux

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