source: trunk/Bonmin/src/Algorithms/OaGenerators/BonFpForMinlp.cpp @ 1407

Last change on this file since 1407 was 1407, checked in by pbonami, 11 years ago

Change cutoff for iFP

File size: 9.4 KB
Line 
1// (C) Copyright CNRS 2008
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// Authors :
6// P. Bonami, CNRS
7//
8// Date : 02/13/2009
9
10#include "BonFpForMinlp.hpp"
11#include "BonminConfig.h"
12
13#include "OsiClpSolverInterface.hpp"
14
15#include "CbcModel.hpp"
16#include "BonCbcLpStrategy.hpp"
17#ifdef COIN_HAS_CPX
18#include "OsiCpxSolverInterface.hpp"
19#endif
20#include "OsiAuxInfo.hpp"
21#include "BonSolverHelp.hpp"
22
23
24namespace Bonmin
25{
26
27/// Constructor with basic setup
28  MinlpFeasPump::MinlpFeasPump(BabSetupBase & b):
29      OaDecompositionBase(b, true, false)
30  {
31    int ivalue;
32    std::string bonmin="bonmin.";
33    std::string prefix = (b.prefix() == bonmin) ? "" : b.prefix();
34    prefix += "pump_for_minlp.";
35    b.options()->GetEnumValue("milp_solver",ivalue,prefix);
36    if (ivalue <= 0) {//uses cbc
37      //nothing to do?
38    }
39    else if (ivalue == 1) {
40      int nodeS, nStrong, nTrust, mig, mir, probe, cover;
41              b.options()->GetEnumValue("node_comparison",nodeS,prefix);
42              b.options()->GetIntegerValue("number_strong_branch",nStrong, prefix);
43              b.options()->GetIntegerValue("number_before_trust", nTrust, prefix);
44              b.options()->GetIntegerValue("Gomory_cuts", mig, prefix);
45              b.options()->GetIntegerValue("probing_cuts",probe, prefix);
46              b.options()->GetIntegerValue("mir_cuts",mir, prefix);
47              b.options()->GetIntegerValue("cover_cuts",cover,prefix);
48             
49              CbcStrategy * strategy =
50                new CbcOaStrategy(mig, probe, mir, cover, nTrust,
51                    nStrong, nodeS, parameters_.cbcIntegerTolerance_, parameters_.subMilpLogLevel_);
52              setStrategy(*strategy);
53              delete strategy;
54
55            }
56            else if (ivalue == 2) {
57        #ifdef COIN_HAS_CPX
58              OsiCpxSolverInterface * cpxSolver = new OsiCpxSolverInterface;
59              b.nonlinearSolver()->extractLinearRelaxation(*cpxSolver);
60              assignLpInterface(cpxSolver);
61        #else
62              std::cerr << "You have set an option to use CPLEX as the milp\n"
63              << "subsolver in oa decomposition. However, apparently\n"
64              << "CPLEX is not configured to be used in bonmin.\n"
65              << "See the manual for configuring CPLEX\n";
66              throw -1;
67        #endif
68            }
69
70            double oaTime;
71            b.options()->GetNumericValue("time_limit",oaTime,prefix);
72            parameter().localSearchNodeLimit_ = 1000000;
73            parameter().maxLocalSearch_ = 100000;
74            parameter().maxLocalSearchPerNode_ = 10000;
75            parameter().maxLocalSearchTime_ =
76            std::min(b.getDoubleParameter(BabSetupBase::MaxTime), oaTime);
77  }
78  MinlpFeasPump::~MinlpFeasPump()
79  {}
80
81  /// virutal method to decide if local search is performed
82  bool
83  MinlpFeasPump::doLocalSearch(BabInfo * babInfo) const
84  {
85    return (nLocalSearch_<parameters_.maxLocalSearch_ &&
86        parameters_.localSearchNodeLimit_ > 0 &&
87        CoinCpuTime() - timeBegin_ < parameters_.maxLocalSearchTime_ &&
88        numSols_ < parameters_.maxSols_);
89  }
90  /// virtual method which performs the OA algorithm by modifying lp and nlp.
91  double
92  MinlpFeasPump::performOa(OsiCuts &cs,
93      solverManip &lpManip,
94      SubMipSolver * &subMip,
95      BabInfo * babInfo,
96      double & cutoff) const
97  {
98
99    bool interuptOnLimit = false;
100    double lastPeriodicLog= CoinCpuTime();
101
102    const int numcols = nlp_->getNumCols();
103    vector<double> savedColLower(nlp_->getNumCols());
104    CoinCopyN(nlp_->getColLower(), nlp_->getNumCols(), savedColLower());
105    vector<double> savedColUpper(nlp_->getNumCols());
106    CoinCopyN(nlp_->getColUpper(), nlp_->getNumCols(), savedColUpper());
107
108
109    OsiSolverInterface * lp = lpManip.si();
110
111    vector<int> indices;
112    for(int i = 0; i < numcols ; i++) {
113      lp->setObjCoeff(i,0);
114      if(!lp->isInteger(i)) {
115      }
116      else { indices.push_back(i);}
117    }
118
119    // If objective is linear need to add to lp constraint for objective
120    if(lp->getNumCols() == nlp_->getNumCols())
121      nlp_->addObjectiveFunction(*lp, nlp_->getColSolution());
122    lp->setObjCoeff(numcols,0);
123    const double * colsol = NULL;
124    OsiBranchingInformation info(lp, false);
125
126    bool milpOptimal = false;
127    nlp_->resolve();
128    //printf("Time limit is %g", parameters_.maxLocalSearchTime_);
129    if (subMip)//Perform a local search
130    {
131      assert(subMip->solver() == lp);
132      set_fp_objective(*lp, nlp_->getColSolution());
133      lp->initialSolve();
134      lp->setColUpper(numcols, cutoff);
135      subMip->find_good_sol(DBL_MAX, parameters_.subMilpLogLevel_,
136          (parameters_.maxLocalSearchTime_ + timeBegin_ - CoinCpuTime()) /* time limit */,
137          parameters_.localSearchNodeLimit_);
138
139      milpOptimal = subMip -> optimal(); 
140      colsol = subMip->getLastSolution();
141      nLocalSearch_++;
142
143    }
144    int numberPasses = 0;
145
146#ifdef OA_DEBUG
147    bool foundSolution = 0;
148#endif
149    double * nlpSol = NULL;
150
151    while (colsol) {
152      numberPasses++;
153
154      //after a prescribed elapsed time give some information to user
155      double time = CoinCpuTime();
156
157
158      //setup the nlp
159      int numberCutsBefore = cs.sizeRowCuts();
160
161      //Fix the variable which have to be fixed, after having saved the bounds
162      info.solution_ = colsol;
163
164      vector<double> x_bar(indices.size());
165      for(unsigned int i = 0 ; i < indices.size() ; i++){
166         x_bar[i] = colsol[indices[i]];
167      }
168
169      double dist = nlp_->solveFeasibilityProblem(indices.size(), x_bar(), indices(), 1, 0, 2);
170
171      handler_->message(FP_DISTANCE, messages_) 
172      <<dist<<CoinMessageEol;
173
174      if(dist < 1e-06){
175         fixIntegers(*nlp_,info, parameters_.cbcIntegerTolerance_, objects_, nObjects_);
176
177         nlp_->resolve();
178         if (post_nlp_solve(babInfo, cutoff)) {
179           //nlp is solved and feasible
180           // Update the cutoff
181           //cutoff = nlp_->getObjValue() * (1 - 10*parameters_.cbcCutoffIncrement_);
182           cutoff = nlp_->getObjValue() - parameters_.cbcCutoffIncrement_;
183           cutoff = nlp_->getObjValue() - 0.1;
184           numSols_++;
185         }
186         nlpSol = const_cast<double *>(nlp_->getColSolution());
187         nlp_->getOuterApproximation(cs, nlpSol, 1, NULL,
188                                  parameter().global_);
189         nlp_->setColLower(savedColLower());
190         nlp_->setColUpper(savedColUpper());
191         nlp_->resolve();
192      }
193      else {
194         nlpSol = const_cast<double *>(nlp_->getColSolution());
195         nlp_->getOuterApproximation(cs, nlpSol, 1, NULL,
196                                  parameter().global_);
197      }
198
199
200
201      int numberCuts = cs.sizeRowCuts() - numberCutsBefore;
202      assert(numberCuts);
203      installCuts(*lp, cs, numberCuts);
204      numberCutsBefore = cs.sizeRowCuts();
205     
206      //check time
207      if (CoinCpuTime() - timeBegin_ > parameters_.maxLocalSearchTime_)
208        break;
209      //do we perform a new local search ?
210      if (nLocalSearch_ < parameters_.maxLocalSearch_ &&
211          numberPasses < parameters_.maxLocalSearchPerNode_ &&
212          parameters_.localSearchNodeLimit_ > 0 && numSols_ < parameters_.maxSols_) {
213
214        /** do we have a subMip? if not create a new one. */
215        if (subMip == NULL) subMip = new SubMipSolver(lp, parameters_.strategy());
216
217        nLocalSearch_++;
218        set_fp_objective(*lp, nlp_->getColSolution());
219
220        lp->setColUpper(numcols, cutoff);
221
222     
223        subMip->find_good_sol(DBL_MAX, parameters_.subMilpLogLevel_,
224            parameters_.maxLocalSearchTime_ + timeBegin_ - CoinCpuTime(),
225            parameters_.localSearchNodeLimit_);
226        milpOptimal = subMip -> optimal(); 
227        colsol = subMip->getLastSolution();
228      if(colsol)
229        handler_->message(FP_MILP_VAL, messages_) 
230        <<colsol[nlp_->getNumCols()]<<CoinMessageEol;
231         
232      }/** endif localSearch*/
233      else if (subMip!=NULL) {
234        delete subMip;
235        subMip = NULL;
236        colsol = NULL;
237      }
238    }
239    if(colsol || ! milpOptimal)
240      return -DBL_MAX;
241    else
242      return DBL_MAX;
243  }
244
245  /** Register OA options.*/
246  void
247  MinlpFeasPump::registerOptions(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions)
248  {
249    roptions->SetRegisteringCategory("Options for feasibility pump", RegisteredOptions::BonminCategory);
250
251    roptions->AddBoundedIntegerOption("fp_log_level",
252        "specify FP iterations log level.",
253        0,2,1,
254        "Set the level of output of OA decomposition solver : "
255        "0 - none, 1 - normal, 2 - verbose"
256                                     );
257
258    roptions->AddLowerBoundedNumberOption("fp_log_frequency",
259        "display an update on lower and upper bounds in FP every n seconds",
260        0.,1.,100.,
261        "");
262  }
263
264/** Put objective of MIP according to FP scheme. */
265void
266MinlpFeasPump::set_fp_objective(OsiSolverInterface &si, const double * colsol) const{
267  if (objects_) {
268    for (int i = 0 ; i < nObjects_ ; i++) {
269      OsiObject * obj = objects_[i];
270      int colnum = obj->columnNumber();
271      if (colnum >= 0) {//Variable branching object
272        double round = floor(colsol[colnum] + 0.5);
273        double coeff = (colsol[colnum] - round ) < 0;
274        si.setObjCoeff(colnum, 1 - 2 * coeff);
275      }
276      else {
277        throw CoinError("OaDecompositionBase::solverManip",
278                        "setFpObjective",
279                        "Can not use FP on problem with SOS constraints");
280      }
281    }
282  }
283  else {
284    int numcols = nlp_->getNumCols();
285    for (int i = 0; i < numcols ; i++) {
286      if (nlp_->isInteger(i)){
287         double round = floor(colsol[i] + 0.5);
288         double coeff = (colsol[i] - round ) < 0;
289         si.setObjCoeff(i, 1 - 2*coeff);
290      }
291    }
292  }
293  si.initialSolve();
294}
295
296}/* End namespace Bonmin. */
Note: See TracBrowser for help on using the repository browser.