source: trunk/Couenne/src/heuristics/CouenneFeasPumpConstructors.cpp @ 898

Last change on this file since 898 was 898, checked in by pbelotti, 8 years ago

fixing two ninja bugs affecting iSol assignment and problem passing (a push with NULL bounds sets them to COUENNE_INFINITY)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.5 KB
Line 
1/* $Id: CouenneFeasPumpConstructors.cpp 898 2012-08-15 12:32:55Z pbelotti $
2 *
3 * Name:    CouenneFeasPump.cpp
4 * Authors: Pietro Belotti
5 *          Timo Berthold, ZIB Berlin
6 * Purpose: Constructors and service methods of the Feasibility Pump class
7 *
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11#include <string>
12
13#include "CouenneConfig.h"
14#include "CouenneFeasPump.hpp"
15#include "CouenneFPpool.hpp"
16#include "CouenneMINLPInterface.hpp"
17#include "CouenneObject.hpp"
18#include "CouenneProblemElem.hpp"
19#include "CouenneProblem.hpp"
20#include "CouenneExprClone.hpp"
21#include "CouenneExprSub.hpp"
22#include "CouenneExprPow.hpp"
23#include "CouenneExprSum.hpp"
24#include "CouenneTNLP.hpp"
25#include "CouenneSparseMatrix.hpp"
26
27using namespace Ipopt;
28using namespace Couenne;
29
30// common code for initializing ipopt application
31void CouenneFeasPump::initIpoptApp () {
32
33  // Although app_ is only used in CouenneFPSolveNLP, we need to have
34  // an object lasting the program's lifetime as otherwise it appears
35  // to delete the nlp pointer at deletion.
36
37  if (!app_)
38    app_ = IpoptApplicationFactory ();
39
40  ApplicationReturnStatus status = app_ -> Initialize ();
41
42  app_ -> Options () -> SetIntegerValue ("max_iter", 1000);
43  app_ -> Options () -> SetIntegerValue // 0 for none, 4 for summary, 5 for iteration output
44    ("print_level", (problem_ -> Jnlst () -> ProduceOutput (J_ITERSUMMARY,  J_NLPHEURISTIC) ? 4 : 
45                     problem_ -> Jnlst () -> ProduceOutput (J_MOREDETAILED, J_NLPHEURISTIC) ? 5 : 0));
46
47  app_ -> Options () -> SetStringValue ("fixed_variable_treatment", "make_parameter");
48
49  // Suppress iteration output from nonlinear solver
50  app_ -> Options () -> SetStringValue ("sb", "yes", false, true);
51
52  if (status != Solve_Succeeded)
53    printf ("FP: Error in initialization\n");
54}
55
56
57// Constructor //////////////////////////////////////////////////
58CouenneFeasPump::CouenneFeasPump (CouenneProblem *couenne,
59                                  CouenneCutGenerator *cg,
60                                  Ipopt::SmartPtr<Ipopt::OptionsList> options):
61  CbcHeuristic         (),
62
63  problem_             (couenne),
64  couenneCG_           (cg),
65  nlp_                 (NULL),
66  app_                 (NULL),
67  milp_                (NULL),
68  postlp_              (NULL),
69  pool_                (NULL),
70
71  numberSolvePerLevel_ (5), // if options are not valid, don't overuse FP
72
73  multDistNLP_         (1.), // settings for classical FP
74  multHessNLP_         (0.),
75  multObjFNLP_         (0.),
76
77  multDistMILP_        (1.),
78  multHessMILP_        (0.),
79  multObjFMILP_        (0.),
80
81  compDistInt_         (FP_DIST_INT),
82  milpCuttingPlane_    (FP_CUT_NONE),
83  nSepRounds_          (0),
84  maxIter_             (COIN_INT_MAX),
85  useSCIP_             (false),
86  milpMethod_          (0),
87  tabuMgt_             (FP_TABU_NONE) {
88
89  if (IsValid (options)) {
90
91    std::string s;
92
93    int compareTerm;
94
95    options -> GetIntegerValue ("feas_pump_iter",       maxIter_,             "couenne.");
96    options -> GetIntegerValue ("feas_pump_level",      numberSolvePerLevel_, "couenne.");
97    options -> GetIntegerValue ("feas_pump_milpmethod", milpMethod_,          "couenne.");
98
99    options -> GetNumericValue ("feas_pump_mult_dist_nlp",  multDistNLP_,     "couenne.");
100    options -> GetNumericValue ("feas_pump_mult_hess_nlp",  multHessNLP_,     "couenne.");
101    options -> GetNumericValue ("feas_pump_mult_objf_nlp",  multObjFNLP_,     "couenne.");
102
103    options -> GetNumericValue ("feas_pump_mult_dist_milp", multDistMILP_,    "couenne.");
104    options -> GetNumericValue ("feas_pump_mult_hess_milp", multHessMILP_,    "couenne.");
105    options -> GetNumericValue ("feas_pump_mult_objf_milp", multObjFMILP_,    "couenne.");
106
107    options -> GetStringValue  ("feas_pump_convcuts", s, "couenne."); 
108
109    milpCuttingPlane_ = 
110      (s == "none")       ? FP_CUT_NONE       :
111      (s == "integrated") ? FP_CUT_INTEGRATED : 
112      (s == "postcut")    ? FP_CUT_POST       : FP_CUT_EXTERNAL;
113
114    options -> GetIntegerValue ("feas_pump_nseprounds", nSepRounds_, "couenne."); 
115
116    options -> GetStringValue  ("feas_pump_vardist",  s, "couenne."); 
117
118    compDistInt_ = 
119      (s == "integer") ? FP_DIST_INT : 
120      (s == "all")     ? FP_DIST_ALL : FP_DIST_POST;
121
122    options -> GetIntegerValue ("feas_pump_milpmethod", milpMethod_, "couenne."); 
123    options -> GetIntegerValue ("feas_pump_poolcomp",   compareTerm, "couenne."); 
124
125    pool_ = new CouenneFPpool (problem_, (enum what_to_compare) compareTerm);
126
127    options -> GetStringValue  ("feas_pump_tabumgt", s, "couenne.");
128
129    tabuMgt_ = 
130      (s == "pool")    ? FP_TABU_POOL    :
131      (s == "perturb") ? FP_TABU_PERTURB : 
132      (s == "cut")     ? FP_TABU_CUT     : FP_TABU_NONE;
133
134    options -> GetStringValue  ("feas_pump_usescip", s, "couenne."); 
135
136#ifdef COIN_HAS_SCIP
137    useSCIP_ = (s == "yes");
138    if (milpMethod_ < 0)
139      milpMethod_ = 0;
140#else
141    if (s == "yes") 
142      problem_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "Warning: you have set feas_pump_usescip to true, but SCIP is not installed.\n");
143#endif
144
145  } else
146    //pool_ = new CouenneFPpool (SUM_NINF);
147    pool_ = new CouenneFPpool (problem_, INTEGER_VARS);
148
149  setHeuristicName ("Couenne Feasibility Pump");
150
151  initIpoptApp ();
152}
153
154 
155// Copy constructor /////////////////////////////////////////////
156CouenneFeasPump::CouenneFeasPump (const CouenneFeasPump &other) {
157
158  operator= (other);
159}
160
161
162// Clone ////////////////////////////////////////////////////////
163CbcHeuristic *CouenneFeasPump::clone () const
164{return new CouenneFeasPump (*this);}
165
166
167// Assignment operator //////////////////////////////////////////
168CouenneFeasPump &CouenneFeasPump::operator= (const CouenneFeasPump & rhs) {
169
170  if (this != &rhs) {
171
172    CbcHeuristic::operator= (rhs);
173
174    problem_             = rhs. problem_;
175    couenneCG_           = rhs. couenneCG_;
176    nlp_                 = rhs. nlp_;
177    app_                 = NULL;
178    milp_                = rhs. milp_   ? rhs. milp_   -> clone () : NULL;
179    postlp_              = rhs. postlp_ ? rhs. postlp_ -> clone () : NULL;
180    pool_                = NULL;
181
182    numberSolvePerLevel_ = rhs. numberSolvePerLevel_;
183
184    multDistNLP_         = rhs. multDistNLP_;
185    multHessNLP_         = rhs. multHessNLP_;
186    multObjFNLP_         = rhs. multObjFNLP_;
187                                       
188    multDistMILP_        = rhs. multDistMILP_;
189    multHessMILP_        = rhs. multHessMILP_;
190    multObjFMILP_        = rhs. multObjFMILP_;
191
192    compDistInt_         = rhs. compDistInt_;
193    milpCuttingPlane_    = rhs. milpCuttingPlane_;
194    nSepRounds_          = rhs. nSepRounds_;
195    maxIter_             = rhs. maxIter_;
196    useSCIP_             = rhs. useSCIP_;
197    milpMethod_          = rhs. milpMethod_;
198    tabuMgt_             = rhs. tabuMgt_;
199
200    if (rhs. pool_)
201      pool_ = new CouenneFPpool (*(rhs. pool_));
202
203    for (std::set <CouenneFPsolution, compareSol>::const_iterator i = rhs.tabuPool_.begin (); 
204         i != rhs.tabuPool_.end ();
205         ++i)
206      tabuPool_. insert (CouenneFPsolution (*i));
207
208    initIpoptApp ();
209  }
210
211  return *this;
212}
213
214
215// Destructor ///////////////////////////////////////////////////
216CouenneFeasPump::~CouenneFeasPump () {
217
218  if (pool_)   delete pool_;
219  if (app_)    delete app_;
220  if (milp_)   delete milp_;
221  if (postlp_) delete postlp_;
222
223  //if (nlp_) delete nlp_; // already deleted by "delete app_;"
224}
225
226
227/// Set new expression as the NLP objective function using
228/// argument as point to minimize distance from. Return new
229/// objective function.
230expression *CouenneFeasPump::updateNLPObj (const double *iSol) {
231
232  expression **list = NULL; 
233
234  int nTerms = 0;
235
236  const double *iS = iSol;
237
238  if ((multHessNLP_ == 0.) || 
239      (nlp_ -> optHessian () == NULL)) {
240
241    list = new expression * [1 + problem_ -> nVars ()];
242
243    // here the objective function is ||x-x^0||_2^2
244
245    // create the argument list (x_i - x_i^0)^2 for all i's
246    for (int i=0; i<problem_ -> nVars (); ++i, ++iS) {
247
248      if (problem_ -> Var (i) -> Multiplicity () <= 0)
249        continue;
250
251      if (compDistInt_ == FP_DIST_INT && 
252          !(problem_ -> Var (i) -> isInteger ()))
253        continue;
254
255      expression *base;
256
257      if      (*iS == 0.) base =              new exprClone (problem_ -> Var (i));
258      else if (*iS <  0.) base = new exprSum (new exprClone (problem_ -> Var (i)), new exprConst (-*iS));
259      else                base = new exprSub (new exprClone (problem_ -> Var (i)), new exprConst  (*iS));
260
261      list [nTerms++] = new exprPow (base, new exprConst (2.));
262    }
263
264  } else {
265
266    // possibly a much larger set of operands
267
268    list = new expression * [problem_ -> nVars () * 
269                             problem_ -> nVars ()];
270
271    // here the objective function is
272    //
273    // ||P(x-x^0)||_2^2 = (x-x^0)' P'P (x-x^0)
274    //
275    // with P'P positive semidefinite stored in CouenneTNLP::optHessian_
276
277    // P is a convex combination, with weights multDistMILP_ and
278    // multHessMILP_, of the distance and the Hessian respectively
279
280    bool *diag = NULL;
281
282    if (multDistNLP_ > 0.) { // only use this if distance is used
283
284      diag = new bool [problem_ -> nVars ()];
285      CoinFillN (diag, problem_ -> nVars (), false);
286    }
287
288    int    *row = nlp_ -> optHessian () -> row ();
289    int    *col = nlp_ -> optHessian () -> col ();
290    double *val = nlp_ -> optHessian () -> val ();
291
292    int     num = nlp_ -> optHessian () -> num ();
293
294    // Add Hessian part -- only lower triangular part
295    for (int i=0; i<num; ++i, ++row, ++col, ++val) {
296
297      if ((problem_ -> Var (*row) -> Multiplicity () <= 0) ||
298          (problem_ -> Var (*col) -> Multiplicity () <= 0))
299        continue;
300
301      // check if necessary given options
302
303      if (compDistInt_ == FP_DIST_INT && 
304          !(problem_ -> Var (*row) -> isInteger () &&
305            problem_ -> Var (*col) -> isInteger ()))
306        continue;
307
308      // second, only do subdiagonal elements
309
310      if (*col < *row) { // that is, lower triangular
311
312        if (2. * *val == 1.) // check if this would have trivial coefficient when doubled (off-diagonal element)
313
314          list [nTerms++] = new exprMul (new exprSub (new exprClone (problem_ -> Var (*row)), new exprConst (iSol [*row])),
315                                         new exprSub (new exprClone (problem_ -> Var (*col)), new exprConst (iSol [*col])));
316
317        else if (fabs (*val) > COUENNE_EPS) { // we don't need extreme precision...
318
319          expression **mlist = new expression * [3];
320
321          mlist [0] = new exprConst (2. * *val);  // twice elements off diagonal
322          mlist [1] = new exprSub (new exprClone (problem_ -> Var (*row)), new exprConst (iSol [*row]));
323          mlist [2] = new exprSub (new exprClone (problem_ -> Var (*col)), new exprConst (iSol [*col]));
324
325          list [nTerms++] = new exprMul (mlist, 3);
326        }
327
328      } else if (*col == *row) { // or diagonal elements
329
330        if (multDistNLP_ > 0.)
331          diag [*col] = true;
332
333        if (*val + multDistNLP_ == 1.)
334
335          list [nTerms++] = new exprPow (new exprSub (new exprClone (problem_ -> Var (*row)), 
336                                                      new exprConst (iSol [*row])), 
337                                         new exprConst (2.));
338
339        else if (fabs (*val + multDistNLP_) > COUENNE_EPS)
340
341          list [nTerms++] = new exprMul (new exprConst (*val + multDistNLP_),
342                                         new exprPow (new exprSub (new exprClone (problem_ -> Var (*row)), 
343                                                                   new exprConst (iSol [*row])), 
344                                                      new exprConst (2.)));
345      }
346    }
347
348    // third, add missing diagonal elements
349
350    if (multDistNLP_ > 0.) {
351
352      // create the argument list (x_i - x_i^0)^2 for all i's
353      for (int i=0; i<problem_ -> nVars (); ++i, ++iS) {
354
355        if (problem_ -> Var (i) -> Multiplicity () <= 0)
356          continue;
357         
358        if ((compDistInt_ == FP_DIST_INT && 
359             !(problem_ -> Var (i) -> isInteger ())) ||
360            diag [i])
361          continue;
362         
363        expression *base;
364
365        if      (*iS == 0.) base =              new exprClone (problem_ -> Var (i));
366        else if (*iS <  0.) base = new exprSum (new exprClone (problem_ -> Var (i)), new exprConst (-*iS));
367        else                base = new exprSub (new exprClone (problem_ -> Var (i)), new exprConst  (*iS));
368
369        list [nTerms++] = new exprPow (base, new exprConst (2.));
370      }
371
372      delete [] diag;
373    }
374  }
375
376  if (multObjFNLP_ != 0.) 
377    list [nTerms++] = new exprMul (new exprConst (multObjFNLP_),
378                                   new exprClone (problem_ -> Obj (0) -> Body ()));
379
380  // resize list
381
382  expression **tmp = list;
383  list = CoinCopyOfArray (tmp, nTerms);
384  delete [] tmp;
385
386  expression *retexpr = new exprSum (list, nTerms);
387
388  //printf ("new objective: "); retexpr -> print (); printf ("\n");
389
390  return retexpr;
391}
392
393
394/// Reads a (possibly fractional) solution and fixes the integer
395/// components in the nonlinear problem for later re-solve
396bool CouenneFeasPump::fixIntVariables (const double *sol) {
397
398  assert (sol);
399
400  t_chg_bounds *chg_bds = new t_chg_bounds [problem_ -> nVars ()];
401
402  for (int i = problem_ -> nVars (); i--;)
403
404    if ((problem_ -> Var (i) -> isInteger ()) &&
405        (problem_ -> Var (i) -> Multiplicity () > 0)) {
406
407      double 
408        value = sol [i],
409        rUp   = ceil  (value - COUENNE_EPS),
410        rDn   = floor (value + COUENNE_EPS);
411
412      // If numerics or sol[i] fractional, set to closest
413
414      value = 
415        (rUp < rDn + 0.5)           ? rUp : 
416        (rUp - value < value - rDn) ? rUp : rDn;
417
418#define INT_NLP_BRACKET 1e-6
419
420      problem_ -> Lb (i) = value - INT_NLP_BRACKET;
421      problem_ -> Ub (i) = value + INT_NLP_BRACKET;
422
423      chg_bds [i].setLower (t_chg_bounds::CHANGED); 
424      chg_bds [i].setUpper (t_chg_bounds::CHANGED); 
425    }
426
427  // Now, to restrict the bounding box even more (and hopefully make
428  // it easier) apply BT
429
430  bool retval = problem_ -> btCore (chg_bds); // maybe fixing makes the nlp infeasible
431
432  delete [] chg_bds;
433
434  return retval;
435}
436
437
438/// initialize options
439void CouenneFeasPump::registerOptions (Ipopt::SmartPtr <Bonmin::RegisteredOptions> roptions) {
440
441  roptions -> AddStringOption2
442    ("feas_pump_heuristic",
443     "Apply the nonconvex Feasibility Pump",
444     "no",
445     "no","",
446     "yes","",
447     "An implementation of the Feasibility Pump for nonconvex MINLPs");
448
449  roptions -> AddLowerBoundedIntegerOption
450    ("feas_pump_level",
451     "Specify the logarithm of the number of feasibility pumps to perform" 
452     " on average for each level of given depth of the tree.",
453     -1,
454     3, "Solve as many nlp's at the nodes for each level of the tree. "
455     "Nodes are randomly selected. If for a "
456     "given level there are less nodes than this number nlp are solved for every nodes. "
457     "For example, if parameter is 8 NLPs are solved for all node until level 8, " 
458     "then for half the node at level 9, 1/4 at level 10.... "
459     "Set to -1 to perform at all nodes.");
460
461  roptions -> AddLowerBoundedIntegerOption
462    ("feas_pump_iter",
463     "Number of iterations in the main Feasibility Pump loop",
464     -1,
465     10, "-1 means no limit");
466
467  // six options
468
469  char option [40];
470  char help   [250];
471
472  std::string terms [] = {"dist", "hess", "objf"};
473  std::string types [] = {"nlp",  "milp"};
474
475  for   (int j=0; j<3; j++) 
476    for (int i=0; i<2; i++) {
477
478      sprintf (option, "feas_pump_mult_%s_%s",                                              terms [j].c_str (), types [i].c_str ());
479      sprintf (help,       "Weight of the %s in the distance function of the %s problem", 
480               !(strcmp ("dist", terms [j].c_str ())) ? "distance" : 
481               !(strcmp ("hess", terms [j].c_str ())) ? "Hessian"  : "original objective function",             types [i].c_str ());
482
483      roptions -> AddBoundedNumberOption
484        (option, help,
485         0., false,
486         1., false,
487         0., "0: no weight, 1: full weight");
488    }
489
490  roptions -> AddStringOption3
491    ("feas_pump_vardist",
492     "Distance computed on integer-only or on both types of variables, in different flavors.",
493     "integer",
494     "integer",         "Only compute the distance based on integer coordinates (use post-processing if numerical errors occur)",
495     "all",             "Compute the distance using continuous and integer variables",
496     "int-postprocess", "Use a post-processing fixed-IP LP to determine a closest-point solution");
497
498  roptions -> AddStringOption4
499    ("feas_pump_convcuts",
500     "Separate MILP-feasible, MINLP-infeasible solution during or after MILP solver.",
501     "none",
502     "integrated", "Done within the MILP solver in a branch-and-cut fashion",
503     "external",   "Done after the MILP solver, in a Benders-like fashion",
504     "postcut",    "Do one round of cuts and proceed with NLP",
505     "none",       "Just proceed to the NLP");
506
507  roptions -> AddBoundedIntegerOption
508    ("feas_pump_nseprounds",
509     "Number of rounds that separate convexification cuts. Must be at least 1",
510     1, 1e5, 4,
511     "");
512
513  roptions -> AddStringOption4
514    ("feas_pump_tabumgt",
515     "Retrieval of MILP solutions when the one returned is unsatisfactory",
516     "pool",
517     "pool",       "Use a solution pool and replace unsatisfactory solution with Euclidean-closest in pool",
518     "perturb",    "Randomly perturb unsatisfactory solution",
519     "cut",        "Separate convexification cuts",
520     "none",       "Bail out of feasibility pump");
521
522  roptions -> AddStringOption2
523    ("feas_pump_usescip",
524     "Should SCIP be used to solve the MILPs?",
525#ifdef COIN_HAS_SCIP
526     "yes", // we want it by default if SCIP is available
527#else
528     "no",  // otherwise switch it off and warn the user if turned on
529#endif
530     "no",  "Use Cbc's branch-and-cut to solve the MILP",
531     "yes", "Use SCIP's branch-and-cut or heuristics (see feas_pump_milpmethod option) to solve the MILP",
532     "");
533
534  roptions -> AddBoundedIntegerOption
535    ("feas_pump_milpmethod",
536     "How should the integral solution be constructed?",
537     -1, 4, -1, 
538       "0: automatic, 1: aggressive heuristics, large node limit, 2: default, node limit, 3: RENS, 4: Objective Feasibility Pump,  -1: solve MILP completely");
539
540  roptions -> AddBoundedIntegerOption
541    ("feas_pump_poolcomp",
542     "Priority field to compare solutions in FP pool",
543     0, 2, 0, 
544       "0: total number of infeasible objects (integer and nonlinear), 1: maximum infeasibility (integer or nonlinear), 2: objective value.");
545}
Note: See TracBrowser for help on using the repository browser.