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

Last change on this file since 694 was 694, checked in by stefan, 10 years ago

Bonmin headers do not do 'using namespace Ipopt' anymore - which is good, but Couenne needs to adapt to this

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