source: trunk/Couenne/src/bound_tightening/FixPointGenCuts.cpp @ 945

Last change on this file since 945 was 945, checked in by stefan, 8 years ago

generateCuts in Cgl >= 0.58 will not be const anymore

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.4 KB
Line 
1/* $Id: FixPointGenCuts.cpp 945 2013-04-06 20:25:21Z stefan $
2 *
3 * Name:    FixPointGenCuts.cpp
4 * Author:  Pietro Belotti
5 * Purpose: Fix point bound tightener -- separator
6 *
7 * (C) Pietro Belotti, 2010.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11#include "CoinHelperFunctions.hpp"
12#include "OsiClpSolverInterface.hpp"
13#include "OsiCuts.hpp"
14#include "CoinTime.hpp"
15
16#include "CouenneFixPoint.hpp"
17
18#include "CouenneProblem.hpp"
19#include "CouennePrecisions.hpp"
20#include "CouenneExprVar.hpp"
21#include "CouenneInfeasCut.hpp"
22
23using namespace Ipopt;
24using namespace Couenne;
25
26/// Cut Generator for FBBT fixpoint
27
28void CouenneFixPoint::generateCuts (const OsiSolverInterface &si,
29                                    OsiCuts &cs,
30                                    const CglTreeInfo treeInfo)
31#if CGL_VERSION_MAJOR == 0 && CGL_VERSION_MINOR <= 57
32  const
33#endif
34  {
35
36  /// Only run this if the latest FBBT terminated on the iteration
37  /// limit, as this suggest that the FPLP might be of some help.
38  /// Termination before iteration limit reached implies that a
39  /// relaxation (on which the FPLP is based) won't generate better
40  /// bounds.
41  ///
42  /// However, we do run the first time as otherwise it would be
43  /// nixed for the whole branch-and-bound.
44
45  if (firstCall_) 
46    firstCall_ = false;
47  else
48    if (!(problem_ -> fbbtReachedIterLimit ()))
49      return;
50
51  if (isWiped (cs))
52    return;
53
54  if (treeInfo.inTree && 
55      treeInfo.level > 0 &&
56      treeInfo.pass > 1)
57    return;
58
59  double startTime = CoinCpuTime ();
60
61  perfIndicator_. setOldBounds (si. getColLower (), si. getColUpper ());
62
63  int nInitTightened = nTightened_;
64
65  if (treeInfo.inTree && 
66      treeInfo.level <= 0) {
67
68    problem_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "Fixed Point FBBT: "); 
69    fflush (stdout);
70  }
71
72  ++nRuns_;
73
74  double now = CoinCpuTime ();
75
76  problem_ -> domain () -> push (&si, &cs);
77
78  /// An LP relaxation of a MINLP problem is available in the first
79  /// parameter passed. Let us suppose that this LP relaxation is of
80  /// the form
81  ///
82  /// LP = {x in R^n: Ax <= b}
83  ///
84  /// for suitable nxm matrix A, rhs vector b, and variable vector
85  /// x. Our purpose is that of creating a much larger LP that will
86  /// help us find the interval [l,u] corresponding to the fixpoint of
87  /// an FBBT algorithm. To this purpose, consider a single constraint
88  /// of the above system:
89  ///
90  /// sum {i=1..n} a_ji x_i <= b_j
91  ///
92  /// According to two schools of thought (Leo's and mine), this
93  /// single constraint can give rise to a number of FBBT
94  /// constraints. The two schools of thoughts differ in the meaning
95  /// of b: in mine, it is constant. In Leo's, it is a variable.
96
97  OsiSolverInterface *fplp = NULL;
98
99  if (true) { // placeholder for later selection of LP solver among
100              // those available
101
102    fplp = new OsiClpSolverInterface;
103  }
104
105  /// We need to perform the following steps:
106  ///
107  /// define variables xL and xU
108  /// define variables gL and gU for constraints (downward variables)
109  ///
110  /// add objective function sum_j (u_j - l_j)
111  ///
112  /// for each constraint a^j x <= b_j in Ax <= b:
113  ///   for each variable x_i contained:
114  ///     depending on sign of a_ji, add constraint on x_i^L or x_i^U
115  ///     (*) add constraints on g_j as well
116  ///
117  /// solve LP
118  ///
119  /// if new bounds are better than si's old bounds
120  ///   add OsiColCuts
121
122  /// Get the original problem's coefficient matrix and rhs vector, A and b
123
124  const CoinPackedMatrix *A = si. getMatrixByRow ();
125
126  const int
127       n     = si.  getNumCols  (),
128       m     = si.  getNumRows  (),
129       nCuts = cs.sizeRowCuts   (),
130    *ind     = A -> getIndices  ();
131
132  const double
133    *lb  = si.  getColLower (),
134    *ub  = si.  getColUpper (),
135    *rlb = si.  getRowLower (),
136    *rub = si.  getRowUpper (),
137    *coe = A -> getElements ();
138
139  if (problem_ -> Jnlst () -> ProduceOutput (J_ERROR, J_BOUNDTIGHTENING))
140    for (int i=0; i<n; i++) 
141      printf ("----------- x_%d in [%g,%g]\n", i, lb [i], ub [i]);
142
143  // turn off logging
144  fplp -> messageHandler () -> setLogLevel (0);
145
146  // add lvars and uvars to the new problem
147  for (int i=0; i<n; i++) {
148    bool isActive = problem_ -> Var (i) -> Multiplicity () > 0;
149    fplp -> addCol (0, NULL, NULL, lb [i], ub [i], isActive ? -1. : 0.); // xL_i
150  }
151
152  for (int i=0; i<n; i++) {
153    bool isActive = problem_ -> Var (i) -> Multiplicity () > 0;
154    fplp -> addCol (0, NULL, NULL, lb [i], ub [i], isActive ? +1. : 0.); // xU_i
155  }
156
157  if (extendedModel_) {
158
159    for (int j=0; j<m; j++) fplp -> addCol (0, NULL, NULL, rlb [j],       COIN_DBL_MAX, 0.); // bL_j
160    for (int j=0; j<m; j++) fplp -> addCol (0, NULL, NULL, -COIN_DBL_MAX, rub [j],      0.); // bU_j
161  }
162
163  // Scan each row of the matrix
164 
165  for (int j=0; j<m; j++) { // for each row
166
167    int nEl = A -> getVectorSize (j); // # elements in each row
168
169    if (!nEl)
170      continue;
171
172    if (problem_ -> Jnlst () -> ProduceOutput (J_ERROR, J_BOUNDTIGHTENING)) {
173
174      printf ("row %4d, %4d elements: ", j, nEl);
175
176      for (int i=0; i<nEl; i++) {
177        printf ("%+g x%d ", coe [i], ind [i]);
178        fflush (stdout);
179      }
180
181      printf ("\n");
182    }
183
184    // create cuts for the xL and xU elements //////////////////////
185
186    if (extendedModel_ || rlb [j] > -COUENNE_INFINITY) 
187      for (int i=0; i<nEl; i++) 
188        createRow (-1, ind [i], n, fplp, ind, coe, rlb [j], nEl, extendedModel_, j, m + nCuts); // downward constraints -- on x_i
189
190    if (extendedModel_ || rub [j] <  COUENNE_INFINITY) 
191      for (int i=0; i<nEl; i++) 
192        createRow (+1, ind [i], n, fplp, ind, coe, rub [j], nEl, extendedModel_, j, m + nCuts); // downward constraints -- on x_i
193
194    // create (at most 2) cuts for the bL and bU elements //////////////////////
195
196    if (extendedModel_) {
197      createRow (-1, 2*n     + j, n, fplp, ind, coe, rlb [j], nEl, extendedModel_, j, m + nCuts); // upward constraints -- on bL_i
198      createRow (+1, 2*n + m + j, n, fplp, ind, coe, rub [j], nEl, extendedModel_, j, m + nCuts); // upward constraints -- on bU_i
199    }
200
201    ind += nEl;
202    coe += nEl;
203  }
204
205  // similarly, scan previous cuts in cs //////////////////////////////////////
206
207  for (int j = 0, jj = nCuts; jj--; j++) {
208
209    // create cuts for the xL and xU elements //////////////////////
210
211    OsiRowCut *cut = cs.rowCutPtr (j);
212
213    const CoinPackedVector &row = cut -> row ();
214
215    const int 
216      nEl = row.getNumElements (),
217      *ind = row.getIndices ();
218
219    const double *coe = row.getElements ();
220
221    if (extendedModel_) {
222      fplp -> addCol (0, NULL, NULL, cut -> lb (),       COIN_DBL_MAX, 0.); // bL_j
223      fplp -> addCol (0, NULL, NULL, -COIN_DBL_MAX, cut -> ub (),      0.); // bU_j
224    }
225
226    if (extendedModel_ || cut -> lb () > -COUENNE_INFINITY) 
227      for (int i=0; i<nEl; i++) 
228        createRow (-1, ind [i], n, fplp, ind, coe, cut -> lb (), nEl, extendedModel_, m + j, m + nCuts); // downward constraints -- on x_i
229
230    if (extendedModel_ || cut -> ub () <  COUENNE_INFINITY) 
231      for (int i=0; i<nEl; i++) 
232        createRow (+1, ind [i], n, fplp, ind, coe, cut -> ub (), nEl, extendedModel_, m + j, m + nCuts); // downward constraints -- on x_i
233
234    // create (at most 2) cuts for the bL and bU elements
235
236    if (extendedModel_) {
237      createRow (-1, 2*n             + j, n, fplp, ind, coe, cut -> lb (), nEl, extendedModel_, m + j, m + nCuts); // upward constraints -- on bL_i
238      createRow (+1, 2*n + m + nCuts + j, n, fplp, ind, coe, cut -> ub (), nEl, extendedModel_, m + j, m + nCuts); // upward constraints -- on bU_i
239    }
240
241    ind += nEl;
242    coe += nEl;
243  }
244
245  // finally, add consistency cuts, bL <= bU
246
247  if (extendedModel_)
248
249    for (int j=0; j<m; j++) { // for each row
250
251      int    ind [2] = {2*n + j, 2*n + m + j};
252      double coe [2] = {1.,      -1.};
253
254      CoinPackedVector row (2, ind, coe);
255      fplp -> addRow (row, -COIN_DBL_MAX, 0.);
256    }
257
258  /// Now we have an fbbt-fixpoint LP problem. Solve it to get
259  /// (possibly) better bounds
260
261  fplp -> setObjSense (-1.); // we want to maximize
262
263  //printf ("(writing lp) ");
264  //fplp -> writeLp ("fplp");
265
266  fplp -> initialSolve ();
267
268  if (fplp -> isProvenOptimal ()) {
269
270    // if problem not solved to optimality, bounds are useless
271
272    const double 
273      *newLB = fplp -> getColSolution (),
274      *newUB = newLB + n,
275      *oldLB = si. getColLower (),
276      *oldUB = si. getColUpper ();
277
278    // check old and new bounds
279
280    int 
281      *indLB = new int [n],
282      *indUB = new int [n],
283      ntightenedL = 0,
284      ntightenedU = 0;
285
286    double 
287      *valLB = new double [n],
288      *valUB = new double [n];
289
290    for (int i=0; i<n; i++) {
291
292      if (problem_ -> Jnlst () -> ProduceOutput (J_ERROR, J_BOUNDTIGHTENING))
293        printf ("x%d: [%g,%g] --> [%g,%g]\n", i, 
294                oldLB [i], oldUB [i], 
295                newLB [i], newUB [i]);
296
297      if (newLB [i] > oldLB [i] + COUENNE_EPS) {
298
299        indLB [ntightenedL]   = i;
300        valLB [ntightenedL++] = newLB [i];
301
302        ++nTightened_;
303      }
304
305      if (newUB [i] < oldUB [i] - COUENNE_EPS) {
306
307        indUB [ntightenedU]   = i;
308        valUB [ntightenedU++] = newUB [i];
309
310        ++nTightened_;
311      }
312    }
313
314    if (ntightenedL || ntightenedU) {
315
316      OsiColCut newBound;
317
318      newBound.setLbs (ntightenedL, indLB, valLB);
319      newBound.setUbs (ntightenedU, indUB, valUB);
320
321      cs.insert (newBound);
322    }
323
324    delete [] indLB;
325    delete [] indUB;
326    delete [] valLB;
327    delete [] valUB;
328
329    CPUtime_ += CoinCpuTime () - now;
330
331    if (treeInfo.inTree && 
332        treeInfo.level <= 0) {
333
334      problem_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "%d bounds tightened (%g seconds)\n", 
335                                      nTightened_ - nInitTightened, CoinCpuTime () - now); 
336    }
337
338    perfIndicator_. update (problem_ -> Lb (), problem_ -> Ub (), treeInfo.level);
339    perfIndicator_. addToTimer (CoinCpuTime () - startTime);
340
341  } else {
342    if (treeInfo.inTree && 
343        treeInfo.level <= 0)
344      problem_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, " FPLP infeasible or unbounded.\n");
345
346    WipeMakeInfeas (cs);
347  }
348
349  delete fplp;
350
351  problem_ -> domain () -> pop ();
352}
353
354
355// single cut creation. Parameters:
356//
357//  1) sign:     tells us whether this is a <= or a >= (part of a) constraint.
358//  2) indexVar: index of variable we want to do upward or downward bt
359//  3) nVars:    number of variables in the original problems (original +
360//               auxiliaries). Used to understand if we are adding an
361//               up or a down constraint
362//  4) p:        solver interface to which we are adding constraints
363//  5) indices:  vector containing indices of the linearization constraint (the    i's)
364//  6) coe:                        coeffs                                       a_ji's
365//  7) rhs:      right-hand side of constraint
366//  8) nEl:      number of elements of this linearization cut
367//  9) extMod:   extendedModel_
368// 10) indCon:   index of constraint being treated (and corresponding bL, bU)
369// 11) nCon:     number of constraints
370
371void CouenneFixPoint::createRow (int sign,
372                                 int indexVar,
373                                 int nVars,
374                                 OsiSolverInterface *p,
375                                 const int *indices,
376                                 const double *coe,
377                                 const double rhs,
378                                 const int nEl,
379                                 bool extMod,
380                                 int indCon,
381                                 int nCon) const {
382
383  ///////////////////////////////////////////////////////////////////////////////////////////////////////
384  ///
385  /// Define
386  ///
387  /// I+ the subset of {1..n} such that a_ji > 0 and i != indexVar;
388  /// I- the subset of {1..n} such that a_ji < 0 and i != indexVar.
389  ///
390  /// Consider
391  ///
392  /// sum {i=1..n} a_ji x_i = b_j,      (=)
393  ///
394  /// equivalent to the two constraints
395  ///
396  /// sum {i=1..n} a_ji x_i >= b_j.     (>) -- sign will be -1 (rlb)
397  /// sum {i=1..n} a_ji x_i <= b_j      (<) -- sign will be +1 (rub)
398  ///
399  /// Hence we should consider both (<) and (>) when we have an
400  /// equality constraint. The resulting set of upward FBBT is as
401  /// follows:
402  ///
403  /// sum {i in I+} a_ji xU_i + sum {i in I-} a_ji xL_i >= bU_j  (only if (<) enforced, i.e., sign ==  1)
404  /// sum {i in I+} a_ji xL_i + sum {i in I-} a_ji xU_i <= bL_j  (only if (>) enforced, i.e., sign == -1)
405  ///
406  /// together with the constraints defining the initial bounding
407  /// interval of the auxiliary variable (already included):
408  ///
409  /// bU_j <= bU0_j (<)
410  /// bL_j >= bL0_j (>)
411  ///
412  /// The set of downward FBBT constraints is instead:
413  ///
414  /// xL_i >= (bL_j - sum {k in I+} a_jk xU_k - sum {k in I-} a_jk xL_k) / a_ji   (if a_ji > 0 and (>))
415  /// xU_i <= (bU_j - sum {k in I+} a_jk xL_k - sum {k in I-} a_jk xU_k) / a_ji   (if a_ji > 0 and (<))
416  ///
417  /// xU_i <= (bL_j - sum {k in I+} a_jk xU_k - sum {k in I-} a_jk xL_k) / a_ji   (if a_ji < 0 and (>))
418  /// xL_i >= (bU_j - sum {k in I+} a_jk xL_k - sum {k in I-} a_jk xU_k) / a_ji   (if a_ji < 0 and (<))
419  ///
420  /// probably better rewritten, to avoid (further) numerical issues, as
421  ///
422  ///   a_ji xL_i >=   bL_j - sum {k in I+} a_jk xU_k - sum {k in I-} a_jk xL_k   (if a_ji > 0 and (>))
423  ///   a_ji xU_i <=   bU_j - sum {k in I+} a_jk xL_k - sum {k in I-} a_jk xU_k   (if a_ji > 0 and (<))
424  ///
425  /// - a_ji xU_i <= - bL_j + sum {k in I+} a_jk xU_k + sum {k in I-} a_jk xL_k   (if a_ji < 0 and (>))
426  /// - a_ji xL_i >= - bU_j + sum {k in I+} a_jk xL_k + sum {k in I-} a_jk xU_k   (if a_ji < 0 and (<))
427  ///
428  /// or even better, to keep the old coefficients (but not the indices!), like this:
429  ///
430  /// a_ji xL_i + sum {k in I+} a_jk xU_k + sum {k in I-} a_jk xL_k - bL_j >= 0  (if a_ji > 0 and (>))
431  /// a_ji xU_i + sum {k in I+} a_jk xL_k + sum {k in I-} a_jk xU_k - bU_j <= 0  (if a_ji > 0 and (<))
432  ///
433  /// a_ji xU_i + sum {k in I+} a_jk xU_k + sum {k in I-} a_jk xL_k - bL_j >= 0  (if a_ji < 0 and (>))
434  /// a_ji xL_i + sum {k in I+} a_jk xL_k + sum {k in I-} a_jk xU_k - bU_j <= 0  (if a_ji < 0 and (<))
435  ///
436  /// and finally we need initial lower/upper bounds:
437  ///
438  /// xL_i >= xL0_i
439  /// xU_i <= xU0_i
440  ///
441  /// and some consistency constraints
442  ///
443  /// bL_i <= bU_i
444  ///
445  /// (these and the two bound constraints above are already added in
446  /// the main function above).
447  ///
448  /////////////////////////////////////////////////////////////////////////////////////////////////////
449  ///
450  /// According to my school of thought, instead, there is no
451  /// upward/downward directions to simulate. Hence, considering again
452  ///
453  /// sum {i=1..n} a_ji x_i >= b_j      (>) -- sign will be -1 (rlb)
454  /// sum {i=1..n} a_ji x_i <= b_j      (<) -- sign will be +1 (rub)
455  ///
456  /// we'll have similar constraints, where bL and bU are replaced by
457  /// the original rhs.
458  ///
459  /// xL_i >= (b_j - sum {k in I+} a_jk xU_k - sum {k in I-} a_jk xL_k) / a_ji   (if a_ji > 0 and (>))
460  /// xU_i <= (b_j - sum {k in I+} a_jk xL_k - sum {k in I-} a_jk xU_k) / a_ji   (if a_ji > 0 and (<))
461  ///                                                                               
462  /// xU_i <= (b_j - sum {k in I+} a_jk xU_k - sum {k in I-} a_jk xL_k) / a_ji   (if a_ji < 0 and (>))
463  /// xL_i >= (b_j - sum {k in I+} a_jk xL_k - sum {k in I-} a_jk xU_k) / a_ji   (if a_ji < 0 and (<))
464  ///
465  /// once again rewritten as
466  ///
467  ///   a_ji xL_i >=   b_j - sum {k in I+} a_jk xU_k - sum {k in I-} a_jk xL_k   (if a_ji > 0 and (>))
468  ///   a_ji xU_i <=   b_j - sum {k in I+} a_jk xL_k - sum {k in I-} a_jk xU_k   (if a_ji > 0 and (<))
469  ///                                                                               
470  /// - a_ji xU_i <= - b_j + sum {k in I+} a_jk xU_k + sum {k in I-} a_jk xL_k   (if a_ji < 0 and (>))
471  /// - a_ji xL_i >= - b_j + sum {k in I+} a_jk xL_k + sum {k in I-} a_jk xU_k   (if a_ji < 0 and (<))
472  ///
473  /// or even better:
474  ///
475  /// a_ji xL_i + sum {k in I+} a_jk xU_k + sum {k in I-} a_jk xL_k >= b_j       (if a_ji > 0 and (>))
476  /// a_ji xU_i + sum {k in I+} a_jk xL_k + sum {k in I-} a_jk xU_k <= b_j       (if a_ji > 0 and (<))
477  ///                                                                               
478  /// a_ji xU_i + sum {k in I+} a_jk xU_k + sum {k in I-} a_jk xL_k >= b_j       (if a_ji < 0 and (>))
479  /// a_ji xL_i + sum {k in I+} a_jk xL_k + sum {k in I-} a_jk xU_k <= b_j       (if a_ji < 0 and (<))
480  ///
481  /// No other cuts are needed.
482  ///
483  /////////////////////////////////////////////////////////////////////////////////////////////////////
484
485
486  if (problem_ -> Jnlst () -> ProduceOutput (J_ERROR, J_BOUNDTIGHTENING)) {
487
488    printf ("creating constraint from: ");
489
490    for (int i=0; i<nEl; i++)
491      printf ("%+g x%d ", coe [i], indices [i]);
492
493    printf ("%c= %g for variable index %d: ", sign > 0 ? '<' : '>', rhs, indexVar);
494  }
495
496  int nTerms = nEl;
497
498  if (extMod) 
499    nTerms++; // always add one element when using extended model
500
501  int    *iInd = new int    [nTerms];
502  double *elem = new double [nTerms];
503
504  // coefficients are really easy
505
506  CoinCopyN (coe, nEl, elem);
507
508  if (extMod) {
509    elem [nEl] = -1.; // extended model, coefficient for bL or bU
510    iInd [nEl] = 2*nVars + indCon + ((sign > 0) ? nCon : 0);
511  }
512
513  // indices are not so easy...
514
515  for (int k=0; k<nEl; k++) {
516
517    int curInd = indices [k];
518
519    iInd [k] = curInd; // Begin with xL_i, same index as x_i in the
520                       // original model. Should add n depending on a
521                       // few things...
522
523    if (curInd == indexVar) { // x_k is x_i itself
524      if (((sign > 0) && (coe [k] > 0.)) || 
525          ((sign < 0) && (coe [k] < 0.)))
526
527      iInd [k] += nVars;
528
529    } else if (((coe [k] > 0.) && (sign < 0)) ||
530               ((coe [k] < 0.) && (sign > 0)))
531
532      iInd [k] += nVars;
533  }
534
535  CoinPackedVector vec (nTerms, iInd, elem);
536
537  double 
538    lb = sign > 0 ? -COIN_DBL_MAX : extMod ? 0. : rhs,
539    ub = sign < 0 ? +COIN_DBL_MAX : extMod ? 0. : rhs;
540
541  p -> addRow (vec, lb, ub);
542
543  // Update time spent doing this
544
545  if (problem_ -> Jnlst () -> ProduceOutput (J_ERROR, J_BOUNDTIGHTENING)) {
546
547    for (int i=0; i<nEl; i++)
548      printf ("%+g x%d ", elem [i], iInd [i]);
549
550    printf ("in [%g,%g]\n", lb, ub);
551  }
552
553  // OsiRowCut *cut = new OsiRowCut (lb, ub, nTerms, nTerms, iInd, elem);
554  // cut -> print ();
555  // delete cut;
556
557  delete [] iInd;
558  delete [] elem;
559}
Note: See TracBrowser for help on using the repository browser.