# source:trunk/Couenne/src/heuristics/CouenneFeasPump.cpp@375

Last change on this file since 375 was 375, checked in by pbelotti, 4 years ago

completed writing Feasibility Pump

• Property svn:keywords set to Id
File size: 5.9 KB
Line
1/* \$Id\$
2 *
3 * Name:    CouenneFeasPump.cpp
4 * Authors: Pietro Belotti, Lehigh University
5 *          Timo Berthold, ZIB Berlin
6 * Purpose: Implement the Feasibility Pump heuristic class
7 * Created: August 5, 2009
8 *
10 */
11
12#include "CbcModel.hpp"
13#include "CoinTime.hpp"
14#include "CoinHelperFunctions.hpp"
15
16#include "CouenneFeasPump.hpp"
17#include "CouenneProblem.hpp"
18#include "CouenneProblemElem.hpp"
19#include "CouenneCutGenerator.hpp"
20#include "CouenneTNLP.hpp"
21
22using namespace Couenne;
23
24// Solve
25int CouenneFeasPump::solution (double &objVal, double *newSolution) {
26
27  if (CoinCpuTime () > problem_ -> getMaxCpuTime ())
28    return 0;
29
30  // This FP works as follows:
31  //
32  // obtain current NLP solution xN or, if none available, current LP solution
33  //
34  // repeat {
35  //
36  //   1) compute MILP solution(s) xI that is H_1-closest to xN
37  //   2) insert them in pool
38  //   3) consider the most promising one in the whole pool
39  //
40  //   if it is MINLP-feasible (hack incumbent callback)
41  //     update cutoff with min obj among MINLP-feasible
42  //     run BT
43  //   else
44  //     apply linearization cuts to MILP
45  //
46  //   select subset of variables (cont or integer), fix them
47  //
48  //   compute NLP solution xN that is H_2-closest to xI
49  //
50  //   if MINLP-feasible
51  //     update cutoff
52  //
53  // } until exit condition satisfied
54  //
55  // fix integer variables
56  // resolve NLP
57
58  CouNumber
59    *nSol = NULL, // solution of the nonlinear problem
60    *iSol = NULL, // solution of the IP problem
61    *best = NULL; // best solution found so far
62
63  int
64    niter  = 0, // current # of iterations
65    retval = 0; // 1 if found a better solution
66
67  /////////////////////////////////////////////////////////////////////////
68  //                      _                   _
69  //                     (_)                 | |
70  //  _ __ ___     __ _   _    _ __          | |    ___     ___    _ __
71  // | '_ ` _ \   / _` | | |  | '_ \         | |   / _ \   / _ \  | '_ `.
72  // | | | | | | | (_| | | |  | | | |        | |  | (_) | | (_) | | |_) |
73  // |_| |_| |_|  \__,_| |_|  |_| |_|        |_|   \___/   \___/  | .__/
74  //                                                              | |
75  /////////////////////////////////////////////////////////////// |_| /////
76
77  milp_ = model_ -> solver () -> clone ();
78
79  expression *origObj = problem_ -> Obj (0) -> Body ();
80
81  // copy bounding box and current solution to the problem (better
82  // linearization cuts)
83  problem_ -> domain () -> push (milp_ -> getNumCols (),
84                                 milp_ -> getColSolution (),
85                                 milp_ -> getColLower (),
86                                 milp_ -> getColUpper ());
87  do {
88
89    // INTEGER PART /////////////////////////////////////////////////////////
90
91    // Solve IP using nSol as the initial point to minimize weighted
92    // l-1 distance from. If nSol==NULL, the MILP is created using the
93    // original milp's LP solution.
94    double z = solveMILP (nSol, iSol);
95
96    if (problem_ -> checkNLP (iSol, z)) {
97
98      // solution is MINLP feasible!
99      // Save and get out of the loop
100
101      if (z < problem_ -> getCutOff ()) {
102
103        retval = 1;
104        objVal = z;
105        best   = iSol;
106        problem_ -> setCutOff (z);
107
108        // if bound tightening clears the whole feasible set, stop
109        if (!(problem_ -> btCore (NULL)))
110          break;
111      }
112    } else {
113
114      // solution non MINLP feasible, it might get cut by
115      // linearization cuts. If so, add a round of cuts and repeat.
116
117      OsiCuts cs;
118      // remaining three arguments at NULL by default
119      couenneCG_ -> genRowCuts (*milp_, cs, 0, NULL);
120
121      if (cs.sizeRowCuts ()) {
122
123        // the (integer, NLP infeasible) solution could be separated
124
125        milp_ -> applyCuts (cs);
126
127        if (milpCuttingPlane_)
128          continue; // found linearization cut, now re-solve MILP (not quite a FP)
129      }
130    }
131
132    // NONLINEAR PART ///////////////////////////////////////////////////////
133
134    // solve the NLP to find a NLP (possibly non-MIP) feasible solution
135    z = solveNLP (iSol, nSol);
136
137    // check if newly found NLP solution is also integer
138    if (problem_ -> checkNLP (nSol, z) &&
139        (z < problem_ -> getCutOff ())) {
140
141      retval = 1;
142      objVal = z;
143      best   = nSol;
144      problem_ -> setCutOff (z);
145
146      // if bound tightening clears the whole feasible set, stop
147      if (!(problem_ -> btCore (NULL)))
148        break;
149    }
150
151  } while ((niter++ < maxIter_) &&
152           (retval == 0));
153
154  // OUT OF THE LOOP ////////////////////////////////////////////////////////
155
156  // If MINLP solution found,
157  //
158  // 1) restore original objective
159  // 2) fix integer variables
160  // 3) resolve NLP
161
162  if (retval > 0) {
163
164    if (!nlp_) // first call (in this call to FP). Create NLP
165      nlp_ = new CouenneTNLP (problem_);
166
167    problem_ -> setObjective (0, origObj);
168
169    fixIntVariables (best);
170    nlp_ -> setInitSol (best);
171
172    ////////////////////////////////////////////////////////////////
173
174    // shamelessly copied from hs071_main.cpp (it's Open Source too!)
175
176    SmartPtr <IpoptApplication> app = IpoptApplicationFactory ();
177
178    ApplicationReturnStatus status = app -> Initialize ();
179    if (status != Solve_Succeeded) printf ("Error in initialization\n");
180
181    status = app -> OptimizeTNLP (nlp_);
182    if (status != Solve_Succeeded) printf ("Error solving problem\n");
183
184    ////////////////////////////////////////////////////////////////
185
186    double z = nlp_ -> solve (nSol);
187
188    problem_ -> domain () -> pop (); // pushed in fixIntVariables
189
190    // check if newly found NLP solution is also integer (unlikely...)
191    if (problem_ -> checkNLP (nSol, z) &&
192        (z < problem_ -> getCutOff ())) {
193
194      problem_ -> setCutOff (z);
195      objVal = z;
196      best   = nSol;
197    }
198  }
199
200  if (retval > 0)
201    CoinCopyN (best, problem_ -> nVars (), newSolution);
202
203  delete [] iSol;
204  delete [] nSol; // best is either iSol or nSol, so don't delete [] it
205
206  // release bounding box
207  problem_ -> domain () -> pop ();
208
209  // milp is deleted at every call since it changes not only in terms
210  // of variable bounds but also in terms of linearization cuts added
211  delete milp_;
212
213  return retval;
214}
Note: See TracBrowser for help on using the repository browser.