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

Last change on this file since 961 was 961, checked in by pbelotti, 7 years ago

compiling to use SCIP without error

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