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

Last change on this file since 620 was 620, checked in by pbelotti, 3 years ago

some cleanup

• Property svn:keywords set to `Author Date Id Revision`
File size: 9.7 KB
Line
1/* \$Id\$
2 *
3 * Name:    CouenneFeasPump.cpp
4 * Authors: Pietro Belotti
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
22#include "CouenneRecordBestSol.hpp"
23
24using namespace Couenne;
25
26/// When the current IP (non-NLP) point is not MINLP feasible, linear
27/// cuts are added and the IP is re-solved not more than this times
28const int numConsCutPasses = 5;
29
30// Solve
31int CouenneFeasPump::solution (double &objVal, double *newSolution) {
32
33  if (problem_ -> nIntVars () <= 0 ||
34      CoinCpuTime () > problem_ -> getMaxCpuTime ())
35    return 0;
36
37#ifdef DEBUG
38  printf ("================= Feasibility Pump =======================\n");
39#endif
40
41  // This FP works as follows:
42  //
43  // obtain current NLP solution xN or, if none available, current LP solution
44  //
45  // repeat {
46  //
47  //   1) compute MILP solution(s) {xI} H_1-closest to xN
48  //   2) insert them in pool
49  //   3) consider the most promising one xI in the whole pool
50  //
51  //   if xI is MINLP-feasible (hack incumbent callback)
52  //     update cutoff with min obj among MINLP-feasible
53  //     run BT
54  //   [else // a little out of the FP scheme
55  //     apply linearization cuts to MILP]
56  //
57  //   select subset of variables (cont or integer), fix them
58  //
59  //   compute NLP solution xN that is H_2-closest to xI
60  //
61  //   if MINLP feasible
62  //     update cutoff with min obj among MINLP-feasible
63  //     run BT
64  //
65  // } until exit condition satisfied
66  //
67  // fix integer variables
68  // resolve NLP
69
70  CouNumber
71    *nSol = NULL, // solution of the nonlinear problem
72    *iSol = NULL, // solution of the IP problem
73    *best = NULL; // best solution found so far
74
75  int
76    niter  = 0, // current # of iterations
77    retval = 0; // 1 if found a better solution
78
79  /////////////////////////////////////////////////////////////////////////
80  //                      _                   _
81  //                     (_)                 | |
82  //  _ __ ___     __ _   _    _ __          | |    ___     ___    _ __
83  // | '_ ` _ \   / _` | | |  | '_ \         | |   / _ \   / _ \  | '_ `.
84  // | | | | | | | (_| | | |  | | | |        | |  | (_) | | (_) | | |_) |
85  // |_| |_| |_|  \__,_| |_|  |_| |_|        |_|   \___/   \___/  | .__/
86  //                                                              | |
87  /////////////////////////////////////////////////////////////// |_| /////
88
89  int
90    objInd = problem_ -> Obj (0) -> Body () -> Index (),
91    nSep = 0;
92
93  do {
94
95    // INTEGER PART /////////////////////////////////////////////////////////
96
97    // Solve IP using nSol as the initial point to minimize weighted
98    // l-1 distance from. If nSol==NULL, the MILP is created using the
99    // original milp's LP solution.
100
101    double z = solveMILP (nSol, iSol);
102
103    bool isChecked = false;
104#ifdef FM_CHECKNLP2
105    isChecked = problem_->checkNLP2(iSol, 0, false, // do not care about obj
106                                    true, // stopAtFirstViol
107                                    true, // checkALL
108                                    problem_->getFeasTol());
109    if(isChecked) {
110      z = problem_->getRecordBestSol()->getModSolVal();
111    }
112#else /* not FM_CHECKNLP2 */
113    isChecked = problem_ -> checkNLP (iSol, z, true);
114#endif  /* not FM_CHECKNLP2 */
115
116    if (isChecked) {
117
118      // solution is MINLP feasible! Save it.
119
120      if (z < problem_ -> getCutOff ()) {
121
122        retval = 1;
123        objVal = z;
124
125#ifdef FM_CHECKNLP2
126#ifdef FM_TRACE_OPTSOL
127        problem_->getRecordBestSol()->update();
128        best = problem_->getRecordBestSol()->getSol();
129        objVal = problem_->getRecordBestSol()->getVal();
130#else /* not FM_TRACE_OPTSOL */
131        best = problem_->getRecordBestSol()->getModSol();
132        objVal = z;
133#endif /* not FM_TRACE_OPTSOL */
134#else /* not FM_CHECKNLP2 */
135#ifdef FM_TRACE_OPTSOL
136        problem_->getRecordBestSol()->update(iSol, problem_->nVars(),
137                                             z, problem_->getFeasTol());
138        best = problem_->getRecordBestSol()->getSol();
139        objVal = problem_->getRecordBestSol()->getVal();
140#else /* not FM_TRACE_OPTSOL */
141        best   = iSol;
142        objVal = z;
143#endif /* not FM_TRACE_OPTSOL */
144#endif /* not FM_CHECKNLP2 */
145
146        problem_ -> setCutOff (objVal);
147
148        t_chg_bounds *chg_bds = NULL;
149
150        if (objInd >= 0) {
151
152          chg_bds = new t_chg_bounds [problem_ -> nVars ()];
153          chg_bds [objInd].setUpper (t_chg_bounds::CHANGED);
154        }
155
156        // if bound tightening clears the whole feasible set, stop
157        bool is_still_feas = problem_ -> btCore (chg_bds);
158
159        if (chg_bds)
160          delete [] chg_bds;
161
162        if (!is_still_feas)
163          break;
164      }
165    } else {
166
167      // solution non MINLP feasible, it might get cut by
168      // linearization cuts. If so, add a round of cuts and repeat.
169
170      OsiCuts cs;
171
172      problem_ -> domain () -> push (milp_);
173      // remaining three arguments at NULL by default
174      couenneCG_ -> genRowCuts (*milp_, cs, 0, NULL);
175      problem_ -> domain () -> pop ();
176
177      if (cs.sizeRowCuts ()) {
178
179        // the (integer, NLP infeasible) solution could be separated
180
181        milp_ -> applyCuts (cs);
182
183        // found linearization cut, now re-solve MILP (not quite a FP)
184        if (milpCuttingPlane_ && (nSep++ < numConsCutPasses))
185          continue;
186      }
187    }
188
189    nSep = 0;
190
191    // NONLINEAR PART ///////////////////////////////////////////////////////
192
193    // fix some variables and solve the NLP to find a NLP (possibly
194    // non-MIP) feasible solution
195
196    z = solveNLP (iSol, nSol);
197
198    // check if newly found NLP solution is also integer
199
200    isChecked = false;
201#ifdef FM_CHECKNLP2
202    isChecked = problem_->checkNLP2(nSol, 0, false, // do not care about obj
203                                    true, // stopAtFirstViol
204                                    true, // checkALL
205                                    problem_->getFeasTol());
206    if(isChecked) {
207      z = problem_->getRecordBestSol()->getModSolVal();
208    }
209#else /* not FM_CHECKNLP2 */
210    isChecked = problem_ -> checkNLP (nSol, z, true);
211#endif  /* not FM_CHECKNLP2 */
212
213    if (isChecked &&
214        (z < problem_ -> getCutOff ())) {
215
216      retval = 1;
217
218#ifdef FM_CHECKNLP2
219#ifdef FM_TRACE_OPTSOL
220      problem_->getRecordBestSol()->update();
221      best = problem_->getRecordBestSol()->getSol();
222      objVal = problem_->getRecordBestSol()->getVal();
223#else /* not FM_TRACE_OPTSOL */
224      best = problem_->getRecordBestSol()->getModSol();
225      objVal = z;
226#endif /* not FM_TRACE_OPTSOL */
227#else /* not FM_CHECKNLP2 */
228#ifdef FM_TRACE_OPTSOL
229      problem_->getRecordBestSol()->update(nSol, problem_->nVars(),
230                                           z, problem_->getFeasTol());
231      best = problem_->getRecordBestSol()->getSol();
232      objVal = problem_->getRecordBestSol()->getVal();
233#else /* not FM_TRACE_OPTSOL */
234      best   = nSol;
235      objVal = z;
236#endif /* not FM_TRACE_OPTSOL */
237#endif /* not FM_CHECKNLP2 */
238
239      problem_ -> setCutOff (objVal);
240
241      t_chg_bounds *chg_bds = NULL;
242
243      if (objInd >= 0) {
244
245        chg_bds = new t_chg_bounds [problem_ -> nVars ()];
246        chg_bds [objInd].setUpper (t_chg_bounds::CHANGED);
247      }
248
249      // if bound tightening clears the whole feasible set, stop
250      bool is_still_feas = problem_ -> btCore (chg_bds);
251
252      if (chg_bds)
253        delete [] chg_bds;
254
255      if (!is_still_feas)
256        break;
257     }
258
259  } while ((niter++ < maxIter_) &&
260           (retval == 0));
261
262  // OUT OF THE LOOP ////////////////////////////////////////////////////////
263
264  // If MINLP solution found,
265  //
266  // 1) restore original objective
267  // 2) fix integer variables
268  // 3) resolve NLP
269
270  if (retval > 0) {
271
272    if (!nlp_) // first call (in this call to FP). Create NLP
273      nlp_ = new CouenneTNLP (problem_);
274
275    //problem_ -> setObjective (0, origObj);
276
277    fixIntVariables (best);
278    nlp_ -> setInitSol (best);
279
280    ////////////////////////////////////////////////////////////////
281
282    // Solve with original objective function
283
284    ApplicationReturnStatus status = app_ -> OptimizeTNLP (nlp_);
285    if (status != Solve_Succeeded)
286      printf ("Feasibility Pump: error solving NLP problem\n");
287
288    ////////////////////////////////////////////////////////////////
289
290    double z = nlp_ -> getSolValue ();
291
292    problem_ -> domain () -> pop (); // pushed in fixIntVariables
293
294    // check if newly found NLP solution is also integer (unlikely...)
295    bool isChecked = false;
296#ifdef FM_CHECKNLP2
297    isChecked = problem_->checkNLP2(nSol, 0, false, // do not care about obj
298                                    true, // stopAtFirstViol
299                                    true, // checkALL
300                                    problem_->getFeasTol());
301    if (isChecked) {
302      z = problem_->getRecordBestSol()->getModSolVal();
303    }
304#else /* not FM_CHECKNLP2 */
305    isChecked = problem_ -> checkNLP (nSol, z, true);
306#endif  /* not FM_CHECKNLP2 */
307
308    if (isChecked &&
309        (z < problem_ -> getCutOff ())) {
310
311#ifdef FM_CHECKNLP2
312#ifdef FM_TRACE_OPTSOL
313      problem_->getRecordBestSol()->update();
314      best = problem_->getRecordBestSol()->getSol();
315      objVal = problem_->getRecordBestSol()->getVal();
316#else /* not FM_TRACE_OPTSOL */
317      best = problem_->getRecordBestSol()->getModSol();
318      objVal = z;
319#endif /* not FM_TRACE_OPTSOL */
320#else /* not FM_CHECKNLP2 */
321#ifdef FM_TRACE_OPTSOL
322      problem_->getRecordBestSol()->update(nSol, problem_->nVars(),
323                                           z, problem_->getFeasTol());
324      best = problem_->getRecordBestSol()->getSol();
325      objVal = problem_->getRecordBestSol()->getVal();
326#else /* not FM_TRACE_OPTSOL */
327      best   = nSol;
328      objVal = z;
329#endif /* not FM_TRACE_OPTSOL */
330#endif /* not FM_CHECKNLP2 */
331
332      problem_ -> setCutOff (objVal);
333    }
334  }
335
336  if (retval > 0)
337    CoinCopyN (best, problem_ -> nVars (), newSolution);
338
339  delete [] iSol;
340  delete [] nSol;
341
342  // release bounding box
343  problem_ -> domain () -> pop (); // pushed in first call to solveMILP
344
345  // deleted at every call from Cbc, since it changes not only in
346  // terms of variable bounds but also in of linearization cuts added
347
348  delete milp_;
349  milp_ = NULL;
350
351  return retval;
352}
Note: See TracBrowser for help on using the repository browser.