source: trunk/Couenne/src/heuristics/CouenneFPscipSolve.cpp @ 900

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

fixing not-so-integer variable

  • Property svn:keywords set to Id
File size: 15.9 KB
Line 
1/* $Id: CouenneFPscipSolve.cpp 900 2012-08-15 16:23:37Z pbelotti $
2 *
3 * Name:    CouenneFPscipSolve.cpp
4 * Authors: Timo Berthold, ZIB Berlin
5 * Purpose: call SCIP
6 *
7 * This file is licensed under the Eclipse Public License (EPL)
8 */
9
10#include "CouenneProblem.hpp"
11#include "CouenneFeasPump.hpp"
12#include "CouenneExprVar.hpp"
13
14#include "cons_rowcuts.h"
15
16using namespace Couenne;
17
18#ifdef COIN_HAS_SCIP
19
20/* general SCIP includes */
21#include "scip/scip.h"
22#include "scip/cons_linear.h"
23#include "scip/scipdefplugins.h"
24#include "scip/cons_bounddisjunction.h"
25
26SCIP_RETCODE CouenneFeasPump::ScipSolve (double* &sol, int niter, int* nsuciter, CouNumber &obj) {
27
28  static int currentmilpmethod = 0;
29
30  int depth = (model_ -> currentNode ()) ? model_ -> currentNode () -> depth () : 0;
31
32  SCIP* scip;
33
34  SCIP_VAR** vars;
35  const SCIP_Real* lbs;
36  const SCIP_Real* ubs;
37  const SCIP_Real* objs;
38  const char* vartypes;
39  const CoinPackedMatrix * matrix;
40  const CoinBigIndex* rowstarts;
41  const int* rowlengths;
42  const SCIP_Real* coeffs;
43  const SCIP_Real* lhss;
44  const SCIP_Real* rhss;
45  const int* indices; 
46
47  SCIP_Real timelimit;
48
49  SCIP_VAR      **tabuvars;
50  SCIP_Real      *tabubounds;
51  SCIP_BOUNDTYPE *tabuboundtypes;
52
53  double infinity;
54  int nvars;
55  int nconss;
56  int nscipsols;
57
58  bool solveagain;
59
60  // COUENNE_INFINITY , getInfinity()
61
62  // get problem data
63  nvars    = milp_ -> getNumCols ();
64  nconss   = milp_ -> getNumRows ();
65  infinity = milp_ -> getInfinity ();
66
67  // get variable data
68  lbs =      milp_ -> getColLower ();
69  ubs =      milp_ -> getColUpper ();
70  objs =     milp_ -> getObjCoefficients ();
71  vartypes = milp_ -> getColType ();
72
73  // get row data
74  lhss = milp_ -> getRowLower ();
75  rhss = milp_ -> getRowUpper ();
76
77  // get matrix data
78  matrix     = milp_ -> getMatrixByRow();
79  rowstarts  = matrix -> getVectorStarts();
80  rowlengths = matrix -> getVectorLengths();
81  coeffs     = matrix -> getElements();
82  indices    = matrix -> getIndices();
83     
84  // access tabuPool_ in this class
85
86  if (problem_ -> Jnlst () -> ProduceOutput (Ipopt::J_ERROR, J_NLPHEURISTIC)) {
87    SCIPdebugMessage("create SCIP problem instance with %d variables and %d constraints.\n", nvars, nconss);
88  }
89
90  // initialize SCIP
91  SCIP_CALL( SCIPcreate(&scip) );
92  assert(scip != NULL);
93
94   
95  // include default SCIP plugins
96  SCIP_CALL( SCIPincludeDefaultPlugins(scip) );
97
98  // include row cut constraint hanlder
99  if( milpCuttingPlane_ == FP_CUT_INTEGRATED )
100    { 
101      SCIP_CALL( SCIPincludeConshdlrRowcuts(scip, couenneCG_, milp_) );
102    }
103
104  // create problem instance in SCIP
105  SCIP_CALL( SCIPcreateProb(scip, "auxiliary FeasPump MILP", NULL, NULL, NULL, NULL, NULL, NULL, NULL) );
106
107  // allocate local memory for SCIP variable array
108  SCIP_CALL( SCIPallocMemoryArray(scip, &vars, nvars) );
109
110  // Allocating space for tabu stuff
111
112  SCIP_CALL( SCIPallocMemoryArray(scip, &tabuvars      , 2*nvars) ); // fix to nvars + nIntvars
113  SCIP_CALL( SCIPallocMemoryArray(scip, &tabubounds    , 2*nvars) );
114  SCIP_CALL( SCIPallocMemoryArray(scip, &tabuboundtypes, 2*nvars) );
115     
116  // one variable for objective !!!!!!!!!
117
118  // create variables
119  for (int i=0; i<nvars; i++) {
120
121    char varname[SCIP_MAXSTRLEN]; 
122
123    bool neglect = 
124      (i < problem_ -> nVars ()) && 
125      (problem_ -> Var (i) -> Multiplicity () <= 0);
126
127    // check that all data is in valid ranges
128    assert( 0 <= vartypes[i] && vartypes[i] <= 2);
129
130    if (!neglect) {
131
132      checkInfinity(scip, lbs[i], infinity);
133      checkInfinity(scip, ubs[i], infinity);
134    }
135
136    // all variables are named x_i
137    (void) SCIPsnprintf(varname, SCIP_MAXSTRLEN, "x_%d", i);
138    SCIP_CALL( SCIPcreateVar(scip, &vars[i], varname, 
139                             CoinMin (lbs [i], ubs [i]),
140                             CoinMax (lbs [i], ubs [i]),
141                             objs [i], 
142                             vartypes[i] == 0 ? SCIP_VARTYPE_CONTINUOUS : (vartypes[i] == 1 ? SCIP_VARTYPE_BINARY : SCIP_VARTYPE_INTEGER),
143                             TRUE, FALSE, NULL, NULL, NULL, NULL, NULL) );
144
145    // add the variable to SCIP
146    SCIP_CALL( SCIPaddVar(scip, vars[i]) );
147  }
148
149  for (std::set <CouenneFPsolution, compareSol>:: iterator i = tabuPool_ . begin (); 
150       i != tabuPool_ . end (); ++i) {
151
152    const double *x = i -> x ();
153
154    int nEntries = 0;
155
156    SCIP_CONS *tabucons = NULL;
157
158    for (int j = 0; j < i -> n (); ++j) {
159
160      if (problem_ -> Var (j) -> isInteger () && 
161          problem_ -> Var (j) -> Multiplicity () > 0 &&
162          problem_ -> Ub (j) - problem_ -> Lb (j) > .5) {
163
164        // if (fabs (x [j] - floor (x [j] + .5)) >= SCIPfeastol (scip)) {
165        //   printf ("integer var x%d not really integer: %e\n", j, x [j]);
166        // }
167
168        assert (fabs (lbs [j] - problem_ -> Lb (j)) < SCIPfeastol (scip));
169        assert (fabs (ubs [j] - problem_ -> Ub (j)) < SCIPfeastol (scip));
170        assert (fabs (x [j] - floor (x [j] + .5))   < 1e3 * SCIPfeastol (scip));
171
172        assert (nEntries <= 2*nvars - 2);
173
174        double x_rounded = floor (x [j] + .5);
175
176        if        (x [j] >= problem_ -> Ub (j) - COUENNE_EPS) {
177
178          tabuvars       [nEntries]   = vars [j];
179          tabubounds     [nEntries]   = x_rounded - 1.;
180          tabuboundtypes [nEntries++] = SCIP_BOUNDTYPE_UPPER;
181
182        } else if (x [j] <= problem_ -> Lb (j) + COUENNE_EPS) {
183
184          tabuvars       [nEntries]   = vars [j];
185          tabubounds     [nEntries]   = x_rounded + 1.;
186          tabuboundtypes [nEntries++] = SCIP_BOUNDTYPE_LOWER;
187
188        } else {
189
190          tabuvars       [nEntries]   = vars [j];
191          tabubounds     [nEntries]   = x_rounded - 1.;
192          tabuboundtypes [nEntries++] = SCIP_BOUNDTYPE_UPPER;
193
194          tabuvars       [nEntries]   = vars [j];
195          tabubounds     [nEntries]   = x_rounded + 1.;
196          tabuboundtypes [nEntries++] = SCIP_BOUNDTYPE_LOWER;
197        }
198      } 
199    }
200
201    if (nEntries != 0) {
202
203      SCIP_CALL (SCIPcreateConsBounddisjunction (scip, &tabucons, "Tabu Solution", nEntries,
204                                                 tabuvars, tabuboundtypes, tabubounds, 
205                                                 TRUE,  TRUE,  TRUE,  TRUE,  TRUE, 
206                                                 FALSE, FALSE, FALSE, FALSE, FALSE));
207     
208      SCIP_CALL( SCIPaddCons(scip, tabucons) );
209
210      SCIP_CALL (SCIPreleaseCons (scip, &tabucons));
211    }
212  }
213
214  // create constraints
215  for (int i=0; i<nconss; i++) {
216
217    SCIP_CONS* cons;
218       
219    char consname[SCIP_MAXSTRLEN]; 
220    (void) SCIPsnprintf(consname, SCIP_MAXSTRLEN, "row_%d", i);
221
222    // check that all data is in valid ranges
223    checkInfinity(scip, lhss[i], infinity);
224    checkInfinity(scip, rhss[i], infinity);
225
226    // create an empty linear constraint
227    SCIP_CALL( SCIPcreateConsLinear(scip, &cons, consname, 0, NULL, NULL, lhss[i], rhss[i], 
228                                    TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
229
230    // add variables to constraint
231    for(int j=rowstarts[i]; j<rowstarts[i]+rowlengths[i]; j++)       
232      {
233        checkInfinity(scip, coeffs[j], infinity);
234        SCIP_CALL( SCIPaddCoefLinear(scip, cons, vars[indices[j]], coeffs[j]) );
235      }
236
237    // add constraint to SCIP
238    SCIP_CALL( SCIPaddCons(scip, cons) );
239    SCIP_CALL( SCIPreleaseCons(scip, &cons) );       
240  }
241 
242  // determine the method to solve the MILP
243  if (milpMethod_ == 0 && niter == 0)
244    {
245      // initialize currentmilpmethod: at the root node we solve the MILP with SCIP default
246      // deeper in the tree, we will solve the MILP by a FP heuristic
247      if( depth == 0 )
248        currentmilpmethod = 2;
249      else 
250        currentmilpmethod = 4;
251    }
252  else if (milpMethod_ != 0)
253    currentmilpmethod = milpMethod_; // use a fixed method to solve the MILP
254     
255     
256  // MILP solving loop. If the MILP terminates without a solution, it might get resolved with a more expensive atrategy
257  do {
258    solveagain = false;
259       
260    // reset parameters if MILP is solved agian
261    SCIP_CALL( SCIPresetParams(scip) );
262       
263    // deactivate SCIP output
264    if (!(problem_ -> Jnlst () -> ProduceOutput (Ipopt::J_ERROR, J_NLPHEURISTIC))) {
265      SCIP_CALL( SCIPsetIntParam(scip, "display/verblevel", 0) );
266    }
267       
268    // do not abort subproblem on CTRL-C
269    SCIP_CALL( SCIPsetBoolParam(scip, "misc/catchctrlc", FALSE) );
270       
271    // set time limit
272    timelimit = problem_ -> getMaxCpuTime () - CoinCpuTime ();
273
274    if (timelimit < 0) 
275      break;
276
277    SCIP_CALL( SCIPsetRealParam(scip, "limits/time", timelimit) );       
278
279
280    if (problem_ -> Jnlst () -> ProduceOutput (Ipopt::J_ERROR, J_NLPHEURISTIC)) {
281      SCIPinfoMessage(scip, NULL, "using MILP method: %d\n",currentmilpmethod);
282    }
283       
284    // tune SCIP differently, depending on the chosen method to solve the MILP
285    /// -1. Solve the MILP relaxation to proven optimality
286    ///  0. Let Couenne choose
287    ///  1. Partially solve the MILP with emphasis on good solutions
288    ///  2. Solve the MILP relaxation partially, up to a certain node limit
289    ///  3. Apply RENS to 1
290    ///  4. Use Objective FP 2.0 for MILPs
291    switch(currentmilpmethod)
292      {
293      case -1: // solve the MILP completely. SCIP's default setting should be best for this
294        if( milpCuttingPlane_ == FP_CUT_INTEGRATED )
295          { 
296            SCIP_CALL( SCIPsetLongintParam(scip, "constraints/rowcuts/maxcuttingrounds", 0) );
297          }
298        break;
299
300      case 1: // Be aggressive in finding feasible solutions, but lazy about the dual bound.
301        // Set limits on overall nodes and stall nodes (nodes without incumbent improvement).
302        SCIP_CALL( SCIPsetLongintParam(scip, "limits/stallnodes", 1000) );
303        SCIP_CALL( SCIPsetLongintParam(scip, "limits/nodes", 10000) );
304
305        // disable cutting plane separation
306        SCIP_CALL( SCIPsetSeparating(scip, SCIP_PARAMSETTING_OFF, TRUE) );
307       
308        // disable expensive presolving
309        SCIP_CALL( SCIPsetPresolving(scip, SCIP_PARAMSETTING_FAST, TRUE) );
310
311        // use aggressive primal heuristics
312        SCIP_CALL( SCIPsetHeuristics(scip, SCIP_PARAMSETTING_AGGRESSIVE, TRUE) );
313       
314        // use best estimate node selection
315        if( SCIPfindNodesel(scip, "estimate") != NULL )
316          {
317            SCIP_CALL( SCIPsetIntParam(scip, "nodeselection/estimate/stdpriority", INT_MAX/4) ); 
318          }
319       
320        // use inference branching
321        if( SCIPfindBranchrule(scip, "inference") != NULL )
322          {
323            SCIP_CALL( SCIPsetIntParam(scip, "branching/inference/priority", INT_MAX/4) );
324          }
325       
326        // disable conflict analysis
327        SCIP_CALL( SCIPsetBoolParam(scip, "conflict/useprop", FALSE) );
328        SCIP_CALL( SCIPsetBoolParam(scip, "conflict/useinflp", FALSE) );
329        SCIP_CALL( SCIPsetBoolParam(scip, "conflict/useboundlp", FALSE) );
330        SCIP_CALL( SCIPsetBoolParam(scip, "conflict/usesb", FALSE) );
331        SCIP_CALL( SCIPsetBoolParam(scip, "conflict/usepseudo", FALSE) );
332
333        break;
334
335      case 2: // default SCIP with node limits
336        SCIP_CALL( SCIPsetLongintParam(scip, "limits/stallnodes", 500) );
337        SCIP_CALL( SCIPsetLongintParam(scip, "limits/nodes", 5000) );
338        break;
339      case 3: // solve the MILP with RENS. Disable most other features, enable RENS
340        SCIP_CALL( SCIPsetLongintParam(scip, "limits/nodes", 1) );
341
342        // disable cutting plane separation
343        SCIP_CALL( SCIPsetSeparating(scip, SCIP_PARAMSETTING_OFF, TRUE) );
344       
345        // disable expensive presolving
346        SCIP_CALL( SCIPsetPresolving(scip, SCIP_PARAMSETTING_FAST, TRUE) );
347
348        // besides RENS, only use cheap heuristics
349        SCIP_CALL( SCIPsetHeuristics(scip, SCIP_PARAMSETTING_FAST, TRUE) );
350
351        // use inference branching
352        if( SCIPfindBranchrule(scip, "inference") != NULL )
353          {
354            SCIP_CALL( SCIPsetIntParam(scip, "branching/inference/priority", INT_MAX/4) );
355          }
356
357        // ensure that RENS is called
358        if( SCIPfindHeur(scip, "rens") != NULL )
359          {
360            SCIP_CALL( SCIPsetIntParam(scip, "heuristics/rens/freq", 0) );
361            SCIP_CALL( SCIPsetRealParam(scip, "heuristics/rens/minfixingrate", 0.0) );
362          }
363        if( milpCuttingPlane_ == FP_CUT_INTEGRATED )
364          { 
365            SCIP_CALL( SCIPsetLongintParam(scip, "constraints/rowcuts/maxcuttingrounds", 0) );
366          }
367        break;
368
369      case 4: // solve the MILP with Feasibility Pump. Disable most other features, enable stage 3 for feaspump
370        SCIP_CALL( SCIPsetLongintParam(scip, "limits/nodes", 1) );
371
372        // disable cutting plane separation
373        SCIP_CALL( SCIPsetSeparating(scip, SCIP_PARAMSETTING_OFF, TRUE) );
374       
375        // disable expensive presolving
376        SCIP_CALL( SCIPsetPresolving(scip, SCIP_PARAMSETTING_FAST, TRUE) );
377
378        // besides feaspump, only use cheap heuristics
379        SCIP_CALL( SCIPsetHeuristics(scip, SCIP_PARAMSETTING_FAST, TRUE) );
380
381        // use inference branching
382        if( SCIPfindBranchrule(scip, "inference") != NULL )
383          {
384            SCIP_CALL( SCIPsetIntParam(scip, "branching/inference/priority", INT_MAX/4) );
385          }
386
387        // ensure that feasibility pump is called
388        if( SCIPfindHeur(scip, "feaspump") != NULL )
389          {
390            SCIP_CALL( SCIPsetIntParam(scip, "heuristics/feaspump/freq", 0) );
391            SCIP_CALL( SCIPsetIntParam(scip, "heuristics/feaspump/maxsols", -1) );
392
393            // Comment next 8 lines as Workaround for mac/conv/batch
394            // to avoid assertion in stage 3 on nsolutions==0
395
396            if( SCIPversion() <= 2.01 )
397              {
398                SCIP_CALL( SCIPsetBoolParam(scip, "heuristics/feaspump2/stage3", TRUE) );
399              }
400            else
401              {
402                SCIP_CALL( SCIPsetBoolParam(scip, "heuristics/feaspump/stage3", TRUE) );
403              }
404          }
405        if( milpCuttingPlane_ == FP_CUT_INTEGRATED )
406          { 
407            SCIP_CALL( SCIPsetLongintParam(scip, "constraints/rowcuts/maxcuttingrounds", 0) );
408          }
409        break;
410
411      default:
412        printf("invalid MILP method: %d\n", currentmilpmethod);
413        assert(false);
414        break;
415      }
416
417#if 0
418    // writes MILP problem and SCIP settings into a file
419    SCIP_CALL( SCIPwriteOrigProblem(scip, "debug.lp", NULL, FALSE) );
420    SCIP_CALL( SCIPwriteParams(scip, "debug.set", FALSE,TRUE) );
421#endif
422         
423    // solve the MILP
424
425    SCIP_RETCODE retcode = SCIPsolve (scip);
426
427    if (retcode != SCIP_OKAY) {
428      problem_ -> Jnlst () -> Printf (Ipopt::J_ERROR, J_NLPHEURISTIC, "Couenne FP: SCIPsolve did not succeed\n");
429      goto TERMINATION;
430    }
431
432    nscipsols =  SCIPgetNSols(scip);
433     
434    // copy the solution
435    if( nscipsols)
436      {
437        SCIP_SOL** scipsols;
438        SCIP_SOL* bestsol;
439        SCIP_Real cutoffbound;
440
441        int nstoredsols;
442
443        /* get incumbent solution */
444        bestsol = SCIPgetBestSol(scip);
445        assert(bestsol != NULL);
446
447        /* get SCIP solution pool */
448        scipsols = SCIPgetSols(scip);
449        assert(scipsols != NULL);
450
451        if (!sol)
452          sol = new CouNumber [problem_ -> nVars ()];
453
454        // get solution values and objective of incumbent
455        SCIP_CALL( SCIPgetSolVals(scip, bestsol, problem_ -> nVars (), vars, sol) );
456        obj = SCIPgetSolOrigObj(scip, bestsol);
457
458        nstoredsols = 0;
459
460        // if we use an objective feaspump, the objective function might be negative
461        cutoffbound = obj > 0 ? 2*obj : obj / 2.0;
462
463        // insert other SCIP solutions into solution pool
464        // do not store too many or too poor solutions
465        for(int i=1; i<nscipsols && nstoredsols < 10 && 
466              SCIPgetSolOrigObj(scip,scipsols[i]) <= cutoffbound; i++){
467          double* tmpsol;
468
469          tmpsol = new CouNumber [nvars];
470           
471          // get solution values
472          SCIP_CALL( SCIPgetSolVals(scip, scipsols[i], problem_ -> nVars (), vars, tmpsol) );
473          CouenneFPsolution couennesol = CouenneFPsolution (problem_, tmpsol);
474
475          // add solutions to the pool if they are not in the tabu list
476          if (   tabuPool_      . find (couennesol) == tabuPool_      . end () 
477              && pool_ -> Set (). find (couennesol) == pool_ -> Set() . end ()
478              ){
479            pool_ -> Set (). insert (couennesol);
480
481            ++nstoredsols;
482          }
483        }
484
485        ++(*nsuciter);
486
487        // if we succeeded five times in a row, try a cheaper MILP_ solving method next time
488        // TODO: if we want to use time limits, hitting the time limit would be another good reason to switch
489        if( *nsuciter >= 3 && currentmilpmethod < 4 )
490          {
491            ++currentmilpmethod;
492            *nsuciter = 0;
493          }         
494      }
495    //try to use a more aggressive, more expensive way to solve the MILP
496    else if( milpMethod_ == 0 && currentmilpmethod > 1 )
497      {
498        --currentmilpmethod;
499        solveagain = true;
500        *nsuciter = 0;
501
502        // throw away the current solution process
503        SCIP_CALL( SCIPfreeTransform(scip) );
504      }
505    else {
506
507      obj = COIN_DBL_MAX; 
508    }
509
510  } while (solveagain);
511
512  ////////////////////////////////////////////////////////////////
513
514 TERMINATION:
515   
516  // release variables before freeing them
517  for (int i=0; i<nvars; i++) {
518    SCIP_CALL( SCIPreleaseVar(scip, &vars[i]) );
519  }
520
521  // free memory
522  SCIPfreeMemoryArray(scip, &vars);
523
524  // free tabu stuff
525  SCIPfreeMemoryArray(scip, &tabuvars      );
526  SCIPfreeMemoryArray(scip, &tabubounds    );
527  SCIPfreeMemoryArray(scip, &tabuboundtypes);
528
529  SCIP_CALL( SCIPfree(&scip) );
530   
531  BMScheckEmptyMemory();     
532
533  return SCIP_OKAY;
534}
535
536#endif
Note: See TracBrowser for help on using the repository browser.