source: trunk/Couenne/src/convex/generateCuts.cpp @ 131

Last change on this file since 131 was 131, checked in by pbelotti, 13 years ago

looser checks on bound tightening -- should prevent good upper/lower bound on objective from fathoming a node

File size: 17.0 KB
Line 
1/*
2 * Name:    generateCuts.cpp
3 * Author:  Pietro Belotti
4 * Purpose: the generateCuts() method of the convexification class CouenneCutGenerator
5 *
6 * (C) Carnegie-Mellon University, 2006-08.
7 * This file is licensed under the Common Public License (CPL)
8 */
9
10#include "BonAuxInfos.hpp"
11#include "CglCutGenerator.hpp"
12
13#include "CouenneCutGenerator.hpp"
14#include "CouenneProblem.hpp"
15#include "CouenneSolverInterface.hpp"
16
17
18// checks bad cuts against known optimum
19bool isOptimumCut (const CouNumber *opt, OsiCuts &cs, CouenneProblem *p);
20
21
22// set and lift bound for auxiliary variable associated with objective
23// function
24void fictitiousBound (OsiCuts &cs,
25                      CouenneProblem *p, 
26                      bool action) {     // true before convexifying, false afterwards
27
28  // fictitious bound for initial unbounded lp relaxations
29  const CouNumber large_bound =  9.999e12;
30  const CouNumber large_tol = (large_bound / 1e6);
31
32  // set trivial dual bound to objective function, if there is none
33
34  int ind_obj = p -> Obj (0) -> Body () -> Index ();
35
36  if (ind_obj < 0) return;
37
38  // we have a single variable objective function
39
40  //int sense = -1; //(p -> Obj (0) -> Sense () == MINIMIZE) ? -1 : 1;
41
42  if (action)
43    //if (sense<0)
44      {if (p -> Lb (ind_obj) < - large_bound) p -> Lb (ind_obj) = - large_bound;}
45  //else         {if (p -> Ub (ind_obj) >   large_bound) p -> Ub (ind_obj) =   large_bound;}
46  else
47    //if (sense>0) {if (fabs (p->Ub(ind_obj)-large_bound)<large_tol) p->Ub(ind_obj)=COUENNE_INFINITY;}
48    //else         
49      {if (fabs (p->Lb(ind_obj)+large_bound)<large_tol) p->Lb(ind_obj) =-COUENNE_INFINITY;}
50}
51
52
53// translate changed bound sparse array into a dense one
54void sparse2dense (int ncols, t_chg_bounds *chg_bds, int *&changed, int &nchanged) {
55
56  // convert sparse chg_bds in something handier
57
58  changed  = (int *) realloc (changed, ncols * sizeof (int));
59  nchanged = 0;
60
61  for (register int i=ncols, j=0; i--; j++, chg_bds++)
62    if (chg_bds -> lower() != t_chg_bounds::UNCHANGED ||
63        chg_bds -> upper() != t_chg_bounds::UNCHANGED ) {
64      *changed++ = j;
65      nchanged++;
66    }
67
68  changed -= nchanged;
69  //changed = (int *) realloc (changed, nchanged * sizeof (int));
70}
71
72
73/// get new bounds from parents' bounds + branching rules
74void updateBranchInfo (const OsiSolverInterface &si, CouenneProblem *p, 
75                       t_chg_bounds *chg, const CglTreeInfo &info);
76
77/// a convexifier cut generator
78
79void CouenneCutGenerator::generateCuts (const OsiSolverInterface &si,
80                                        OsiCuts &cs, 
81                                        const CglTreeInfo info) const {
82  const int infeasible = 1;
83
84  int nInitCuts = cs.sizeRowCuts ();
85
86  CouNumber
87    *&realOpt = problem_ -> bestSol (),
88    *saveOptimum = realOpt;
89
90  if (!firstcall_ && realOpt) { 
91
92    // have a debug optimal solution. Check if current bounds
93    // contain it, otherwise pretend it does not exist
94
95    CouNumber *opt = realOpt;
96
97    const CouNumber
98      *sol = si.getColSolution (),
99      *lb  = si.getColLower (),
100      *ub  = si.getColUpper ();
101
102    for (int i=problem_ -> nVars (); i--; opt++, lb++, ub++)
103      if ((*opt < *lb - COUENNE_EPS) || 
104          (*opt > *ub + COUENNE_EPS)) {
105        jnlst_ -> Printf (J_VECTOR, J_CONVEXIFYING, 
106                          "out of bounds, ignore x%d = %g [%g,%g] opt = %g\n", 
107                          problem_ -> nVars () - i - 1, *sol, *lb, *ub, *opt);
108
109        // optimal point is not in current bounding box,
110        // pretend realOpt is NULL until we return from this procedure
111        realOpt = NULL;
112        break;
113      }
114  }
115
116  /*static int count = 0;
117  char fname [20];
118  sprintf (fname, "relax_%d", count++);
119  si.writeLp (fname);
120  printf ("writing %s\n", fname);*/
121
122  jnlst_ -> Printf (J_DETAILED, J_CONVEXIFYING,
123                    "generateCuts: level = %d, pass = %d, intree = %d\n",
124                    info.level, info.pass, info.inTree);
125
126  Bonmin::BabInfo * babInfo = dynamic_cast <Bonmin::BabInfo *> (si.getAuxiliaryInfo ());
127
128  if (babInfo)
129    babInfo -> setFeasibleNode ();
130
131  double now   = CoinCpuTime ();
132  int    ncols = problem_ -> nVars ();
133
134  // This vector contains variables whose bounds have changed due to
135  // branching, reduced cost fixing, or bound tightening below. To be
136  // used with malloc/realloc/free
137
138  t_chg_bounds *chg_bds = new t_chg_bounds [ncols];
139
140  /*for (int i=0; i < ncols; i++)
141    if (problem_ -> Var (i) -> Multiplicity () <= 0) {
142      chg_bds [i].setLower (t_chg_bounds::UNCHANGED);
143      chg_bds [i].setUpper (t_chg_bounds::UNCHANGED);
144      }*/
145
146  problem_ -> installCutOff (); // install upper bound
147
148  if (firstcall_) {
149
150    // First convexification //////////////////////////////////////
151
152    // OsiSolverInterface is empty yet, no information can be obtained
153    // on variables or bounds -- and none is needed since our
154    // constructor populated *problem_ with variables and bounds. We
155    // only need to update the auxiliary variables and bounds with
156    // their current value.
157
158    for (int i=0; i < ncols; i++) 
159      if (problem_ -> Var (i) -> Multiplicity () > 0) {
160        chg_bds [i].setLower (t_chg_bounds::CHANGED);
161        chg_bds [i].setUpper (t_chg_bounds::CHANGED);
162      }
163
164    // start with FBBT, should take advantage of cutoff found by NLP
165    // run AFTER initial FBBT...
166    if (problem_ -> doFBBT () &&
167        (! (problem_ -> boundTightening (chg_bds, babInfo))))
168      printf ("Couenne: WARNING, first convexification is infeasible\n");
169
170    // For each auxiliary variable replacing the original (nonlinear)
171    // constraints, check if corresponding bounds are violated, and
172    // add cut to cs
173
174    int nnlc = problem_ -> nCons ();
175
176    for (int i=0; i<nnlc; i++) {
177
178      // for each constraint
179      CouenneConstraint *con = problem_ -> Con (i);
180
181      // (which has an aux as its body)
182      int index = con -> Body () -> Index ();
183
184      if ((index >= 0) && 
185          ((con -> Body () -> Type () == AUX) ||
186           (con -> Body () -> Type () == VAR))) {
187
188        // get the auxiliary that is at the lhs
189        exprVar *conaux = problem_ -> Var (index);
190
191        if (conaux &&
192            (conaux -> Type () == AUX) &&
193            (conaux -> Image ()) && 
194            (conaux -> Image () -> Linearity () <= LINEAR)) {
195
196          // reduce density of problem by adding w >= l rather than
197          // ax + b >= l for any linear auxiliary defined as w := ax+b
198
199          double 
200            lb = (*(con -> Lb ())) (), 
201            ub = (*(con -> Ub ())) ();
202
203          OsiColCut newBound;
204          if (lb > -COUENNE_INFINITY) newBound.setLbs (1, &index, &lb);
205          if (ub <  COUENNE_INFINITY) newBound.setUbs (1, &index, &ub);
206
207          cs.insert (newBound);
208
209          // the auxiliary w of constraint w <= b is associated with a
210          // linear expression w = ax: add constraint ax <= b
211          /*conaux -> Image () -> generateCuts (conaux, si, cs, this, chg_bds,
212                                              conaux -> Index (),
213                                              (*(con -> Lb ())) (),
214                                              (*(con -> Ub ())) ());*/
215
216          // take it from the list of the variables to be linearized
217          //
218          // DO NOT decrease multiplicity. Even if it is a linear
219          // term, its bounds can still be used in implied bounds
220          //
221          // Are we sure? That will happen only if its multiplicity is
222          // nonzero, for otherwise this aux is only used here, and is
223          // useless elsewhere
224          //
225          //conaux -> decreaseMult (); // !!!
226        }
227
228        // also, add constraint w <= b
229
230        // not now, do it later
231
232//      // if there exists violation, add constraint
233//      CouNumber l = con -> Lb () -> Value (),
234//                u = con -> Ub () -> Value ();
235
236//      // tighten bounds in Couenne's problem representation
237//      problem_ -> Lb (index) = CoinMax (l, problem_ -> Lb (index));
238//      problem_ -> Ub (index) = CoinMin (u, problem_ -> Ub (index));
239
240      } else { // body is more than just a variable, but it should be
241               // linear. If so, generate equivalent linear cut
242
243        assert (false); // TODO
244      }
245    }
246
247    if (jnlst_ -> ProduceOutput (J_ITERSUMMARY, J_CONVEXIFYING)) {
248      if (cs.sizeRowCuts ()) {
249        jnlst_ -> Printf (J_ITERSUMMARY, J_CONVEXIFYING,"Couenne: %d constraint row cuts\n",
250                          cs.sizeRowCuts ());
251        for (int i=0; i<cs.sizeRowCuts (); i++) 
252          cs.rowCutPtr (i) -> print ();
253      }
254      if (cs.sizeColCuts ()) {
255        jnlst_ -> Printf (J_ITERSUMMARY, J_CONVEXIFYING,"Couenne: %d constraint col cuts\n",
256                          cs.sizeColCuts ());
257        for (int i=0; i<cs.sizeColCuts (); i++) 
258          cs.colCutPtr (i) -> print ();
259      }
260    }
261  } else {
262
263    // use new optimum as lower bound for variable associated w/objective
264    int indobj = problem_ -> Obj (0) -> Body () -> Index ();
265
266    assert (indobj >= 0);
267
268    // transmit solution from OsiSolverInterface to problem
269    problem_ -> domain () -> push
270      (problem_ -> nVars (),
271       si. getColSolution (), 
272       si. getColLower    (),
273       si. getColUpper    ());
274
275    if (indobj >= 0) {
276
277      // Use current value of objvalue's x as a lower bound for bound
278      // tightening
279      double lp_bound = problem_ -> domain () -> x (indobj);
280
281      //if (problem_ -> Obj (0) -> Sense () == MINIMIZE)
282      {if (lp_bound > problem_ -> Lb (indobj)) problem_ -> Lb (indobj) = lp_bound;}
283           //else {if (lp_bound < problem_ -> Ub (indobj)) problem_ -> Ub (indobj) = lp_bound;}
284    }
285
286    updateBranchInfo (si, problem_, chg_bds, info); // info.depth >= 0 || info.pass >= 0
287  }
288
289  // restore constraint bounds before tightening and cut generation
290  for (int i = problem_ -> nCons (); i--;) {
291
292    // for each constraint
293    CouenneConstraint *con = problem_ -> Con (i);
294
295    // (which has an aux as its body)
296    int index = con -> Body () -> Index ();
297
298    if ((index >= 0) && 
299        ((con -> Body () -> Type () == AUX) ||
300         (con -> Body () -> Type () == VAR))) {
301
302      // if there exists violation, add constraint
303      CouNumber l = con -> Lb () -> Value (),   
304        u = con -> Ub () -> Value ();
305
306      // tighten bounds in Couenne's problem representation
307      problem_ -> Lb (index) = CoinMax (l, problem_ -> Lb (index));
308      problem_ -> Ub (index) = CoinMin (u, problem_ -> Ub (index));
309    }
310  }
311
312  problem_ -> installCutOff (); // install upper bound
313
314  fictitiousBound (cs, problem_, false); // install finite lower bound, if currently -inf
315
316  int *changed = NULL, nchanged;
317
318  // Bound tightening ///////////////////////////////////////////
319
320  // do bound tightening only at first pass of cutting plane in a node
321  // of BB tree (info.pass == 0) or if first call (creation of RLT,
322  // info.pass == -1)
323
324  try {
325
326    // Bound tightening ////////////////////////////////////
327
328    /*printf ("== BT ================\n");
329    for (int i = 0; i < problem_ -> nVars (); i++)
330      if (problem_ -> Var (i) -> Multiplicity () > 0)
331        printf ("%4d %+20.8g [%+20.8g,%+20.8g]\n", i,
332                problem_ -> X  (i), problem_ -> Lb (i), problem_ -> Ub (i));
333                printf("=============================\n");*/
334
335    // Reduced Cost BT -- to be done first to use rcost correctly
336    if (!firstcall_  &&                         // have a linearization already
337        problem_ -> doRCBT () &&                // authorized to do reduced cost tightening
338        problem_ -> redCostBT (&si, chg_bds) && // some variables were tightened with reduced cost
339        !(problem_ -> btCore (chg_bds)))        // in this case, do another round of FBBT
340      throw infeasible;
341
342    // FBBT
343    if (problem_ -> doFBBT () && 
344        //(info.pass <= 0) && // do it in subsequent rounds too
345        (! (problem_ -> boundTightening (chg_bds, babInfo))))
346      throw infeasible;
347
348    // OBBT
349    if (!firstcall_ && // no obbt if first call (there is no LP to work with)
350        problem_ -> obbt (this, si, cs, info, babInfo, chg_bds) < 0)
351      throw infeasible;
352
353    // Bound tightening done /////////////////////////////
354
355    if ((problem_ -> doFBBT () ||
356         problem_ -> doOBBT () ||
357         problem_ -> doABT  ()) &&
358        (jnlst_ -> ProduceOutput (J_VECTOR, J_CONVEXIFYING))) {
359
360      jnlst_ -> Printf(J_VECTOR, J_CONVEXIFYING,"== after bt =============\n");
361      for (int i = 0; i < problem_ -> nVars (); i++)
362        if (problem_ -> Var (i) -> Multiplicity () > 0)
363          jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"%4d %+20.8g [%+20.8g,%+20.8g]\n", i,
364                         problem_ -> X  (i), problem_ -> Lb (i), problem_ -> Ub (i));
365      jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"=============================\n");
366    }
367
368    // Generate convexification cuts //////////////////////////////
369
370    sparse2dense (ncols, chg_bds, changed, nchanged);
371
372    double *nlpSol;
373
374    //--------------------------------------------
375
376    if (babInfo && ((nlpSol = const_cast <double *> (babInfo -> nlpSolution ())))) {
377
378      // Aggressive Bound Tightening ////////////////////////////////
379
380      int logAbtLev = problem_ -> logAbtLev ();
381
382      if (problem_ -> doABT () &&           // flag is checked, AND
383          ((logAbtLev != 0) ||                // (parameter is nonzero OR
384           (info.level == 0)) &&              //  we are at root node), AND
385          (info.pass == 0) &&               // at first round of cuts, AND
386          ((logAbtLev < 0) ||                 // (logAbtLev = -1, OR
387           (info.level <= logAbtLev) ||       //  depth is lower than COU_OBBT_CUTOFF_LEVEL, OR
388           (CoinDrand48 () <                  //  probability inversely proportional to the level)
389            pow (2., (double) logAbtLev - (info.level + 1))))) {
390
391        jnlst_ -> Printf(J_ITERSUMMARY, J_CONVEXIFYING,"  performing ABT\n");
392        if (! (problem_ -> aggressiveBT (nlp_, chg_bds, babInfo)))
393          throw infeasible;
394
395        sparse2dense (ncols, chg_bds, changed, nchanged);
396      }
397
398      // obtain solution just found by nlp solver
399
400      // Auxiliaries should be correct. solution should be the one found
401      // at the node even if not as good as the best known.
402
403      // save violation flag and disregard it while adding cut at NLP
404      // point (which are not violated by the current, NLP, solution)
405      bool save_av = addviolated_;
406      addviolated_ = false;
407
408      // save values
409      problem_ -> domain () -> push
410        (problem_ -> nVars (), 
411         problem_ -> domain () -> x  (), 
412         problem_ -> domain () -> lb (), 
413         problem_ -> domain () -> ub (), false);
414
415      // fill originals with nlp values
416      CoinCopyN (nlpSol, problem_ -> nOrigVars (), problem_ -> domain () -> x ());
417      problem_ -> initAuxs ();
418
419      if (jnlst_ -> ProduceOutput (J_VECTOR, J_CONVEXIFYING)) {
420        jnlst_ -> Printf(J_VECTOR, J_CONVEXIFYING,"== genrowcuts on NLP =============\n");
421        for (int i = 0; i < problem_ -> nVars (); i++)
422          if (problem_ -> Var (i) -> Multiplicity () > 0)
423            jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"%4d %+20.8g [%+20.8g,%+20.8g]\n", i,
424                           problem_ -> X  (i),
425                           problem_ -> Lb (i),
426                           problem_ -> Ub (i));
427        jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"=============================\n");
428      }
429
430      problem_ -> domain () -> current () -> isNlp () = true;
431      genRowCuts (si, cs, nchanged, changed, chg_bds);  // add cuts
432
433      problem_ -> domain () -> pop (); // restore point
434
435      addviolated_ = save_av;     // restore previous value
436
437      //    if (!firstcall_) // keep solution if called from extractLinearRelaxation()
438      babInfo -> setHasNlpSolution (false); // reset it after use //AW HERE
439    } else {
440
441      if (jnlst_ -> ProduceOutput (J_VECTOR, J_CONVEXIFYING)) {
442        jnlst_ -> Printf(J_VECTOR, J_CONVEXIFYING,"== genrowcuts on LP =============\n");
443        for (int i = 0; i < problem_ -> nVars (); i++)
444          if (problem_ -> Var (i) -> Multiplicity () > 0)
445            jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"%4d %+20.8g [%+20.8g,%+20.8g]\n", i,
446                           problem_ -> X  (i),
447                           problem_ -> Lb (i),
448                           problem_ -> Ub (i));
449        jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"=============================\n");
450      }
451
452      genRowCuts (si, cs, nchanged, changed, chg_bds);
453    }
454
455    // change tightened bounds through OsiCuts
456    if (nchanged)
457      genColCuts (si, cs, nchanged, changed);
458
459    if (firstcall_ && (cs.sizeRowCuts () >= 1))
460      jnlst_->Printf(J_ITERSUMMARY, J_CONVEXIFYING,
461                     "Couenne: %d initial row cuts\n", cs.sizeRowCuts ());
462
463    if (realOpt && // this is a good time to check if we have messed up with the optimal solution
464        isOptimumCut (realOpt, cs, problem_))
465      jnlst_->Printf(J_ITERSUMMARY, J_CONVEXIFYING,
466                     "\n\n CUT OPTIMUM\n\n");
467  }
468
469  catch (int exception) {
470
471    if ((exception == infeasible) && (!firstcall_)) {
472
473      jnlst_ -> Printf (J_ITERSUMMARY, J_CONVEXIFYING,
474                        "Couenne: Infeasible node\n");
475
476      OsiColCut *infeascut = new OsiColCut;
477
478      if (infeascut) {
479        int i=0;
480        double upper = -1., lower = +1.;
481        infeascut -> setLbs (1, &i, &lower);
482        infeascut -> setUbs (1, &i, &upper);
483        cs.insert (infeascut);
484        delete infeascut;
485      }
486    }
487
488    if (babInfo) // set infeasibility to true in order to skip NLP heuristic
489      babInfo -> setInfeasibleNode ();
490  }
491
492  delete [] chg_bds;
493
494  if (changed) 
495    free (changed);
496
497  if (firstcall_) {
498
499    jnlst_ -> Printf (J_SUMMARY, J_CONVEXIFYING, 
500                      "Couenne: %d cuts (%d row, %d col) for linearization\n", 
501                      cs.sizeRowCuts () + cs.sizeColCuts (),
502                      cs.sizeRowCuts (),  cs.sizeColCuts ());
503
504    fictitiousBound (cs, problem_, true);
505    firstcall_  = false;
506    ntotalcuts_ = nrootcuts_ = cs.sizeRowCuts ();
507  } else { 
508    problem_ -> domain () -> pop ();
509
510    ntotalcuts_ += (cs.sizeRowCuts () - nInitCuts);
511
512    if (saveOptimum)
513      realOpt = saveOptimum; // restore debug optimum
514  }
515
516  septime_ += CoinCpuTime () - now;
517
518  if (jnlst_ -> ProduceOutput (J_ITERSUMMARY, J_CONVEXIFYING)) {
519
520    if (cs.sizeColCuts ()) {
521      jnlst_ -> Printf (J_ITERSUMMARY, J_CONVEXIFYING,"Couenne col cuts:\n");
522      for (int i=0; i<cs.sizeColCuts (); i++) 
523        cs.colCutPtr (i) -> print ();
524    }
525  }
526
527  if (!(info.inTree)) 
528    rootTime_ = CoinCpuTime ();
529}
Note: See TracBrowser for help on using the repository browser.