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

Last change on this file since 788 was 788, checked in by pbelotti, 9 years ago

fplp infeasible -> node infeasible

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