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

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

Add MSVisualStudio dir

File size: 10.5 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      subMip->optimize(DBL_MAX, parameters_.subMilpLogLevel_,
137          (parameters_.maxLocalSearchTime_ + timeBegin_ - CoinCpuTime()) /* time limit */,
138          parameters_.localSearchNodeLimit_);
139
140      milpOptimal = subMip -> optimal(); 
141      colsol = subMip->getLastSolution();
142      nLocalSearch_++;
143      if(milpOptimal)
144        handler_->message(SOLVED_LOCAL_SEARCH, messages_)
145        <<subMip->nodeCount()<<subMip->iterationCount()<<CoinMessageEol;
146      else
147        handler_->message(LOCAL_SEARCH_ABORT, messages_)
148        <<subMip->nodeCount()<<subMip->iterationCount()<<CoinMessageEol;
149    }
150    int numberPasses = 0;
151
152#ifdef OA_DEBUG
153    bool foundSolution = 0;
154#endif
155    double * nlpSol = NULL;
156
157    while (colsol) {
158      numberPasses++;
159
160      //after a prescribed elapsed time give some information to user
161      double time = CoinCpuTime();
162
163
164      //setup the nlp
165      int numberCutsBefore = cs.sizeRowCuts();
166
167      //Fix the variable which have to be fixed, after having saved the bounds
168      info.solution_ = colsol;
169
170      vector<double> x_bar(indices.size());
171      for(unsigned int i = 0 ; i < indices.size() ; i++){
172         x_bar[i] = colsol[indices[i]];
173      }
174
175      double dist = nlp_->solveFeasibilityProblem(indices.size(), x_bar(), indices(), 1, 0, 2);
176
177      handler_->message(FP_DISTANCE, messages_) 
178      <<dist<<CoinMessageEol;
179
180      if(dist < 1e-05){
181         fixIntegers(*nlp_,info, parameters_.cbcIntegerTolerance_, objects_, nObjects_);
182
183         nlp_->resolve();
184         bool restart = false;
185         if (post_nlp_solve(babInfo, cutoff)) {
186           restart = true;
187           //nlp is solved and feasible
188           // Update the cutoff
189           printf("Cutoff increment %g\n", parameters_.cbcCutoffIncrement_);
190           cutoff = nlp_->getObjValue() - 
191                    parameters_.cbcCutoffIncrement_;
192           cutoff = nlp_->getObjValue() - 0.1;
193           printf("New cutoff value %g\n", cutoff);
194           numSols_++;
195         }
196         else{
197           //nlp_->setColLower(savedColLower());
198           //nlp_->setColUpper(savedColUpper());
199           //dist = nlp_->solveFeasibilityProblem(indices.size(), x_bar(), indices(), 1, 0, 2);
200         }
201         nlpSol = const_cast<double *>(nlp_->getColSolution());
202         nlp_->getOuterApproximation(cs, nlpSol, 1, NULL,
203                                  parameter().global_);
204         //if(restart){
205           nlp_->setColLower(savedColLower());
206           nlp_->setColUpper(savedColUpper());
207         if(restart)
208           nlp_->resolve();
209         //}
210      }
211      else {
212         nlpSol = const_cast<double *>(nlp_->getColSolution());
213         nlp_->getOuterApproximation(cs, nlpSol, 1, NULL,
214                                  parameter().global_);
215      }
216
217
218
219      int numberCuts = cs.sizeRowCuts() - numberCutsBefore;
220      assert(numberCuts);
221      installCuts(*lp, cs, numberCuts);
222      numberCutsBefore = cs.sizeRowCuts();
223     
224      //check time
225      if (CoinCpuTime() - timeBegin_ > parameters_.maxLocalSearchTime_)
226        break;
227      //do we perform a new local search ?
228      if (nLocalSearch_ < parameters_.maxLocalSearch_ &&
229          numberPasses < parameters_.maxLocalSearchPerNode_ &&
230          parameters_.localSearchNodeLimit_ > 0 && numSols_ < parameters_.maxSols_) {
231
232        /** do we have a subMip? if not create a new one. */
233        if (subMip == NULL) subMip = new SubMipSolver(lp, parameters_.strategy());
234
235        nLocalSearch_++;
236        set_fp_objective(*lp, nlp_->getColSolution());
237
238        lp->setColUpper(numcols, cutoff);
239
240     
241        //subMip->find_good_sol(DBL_MAX, parameters_.subMilpLogLevel_,
242        subMip->optimize(DBL_MAX, parameters_.subMilpLogLevel_,
243                         parameters_.maxLocalSearchTime_ + timeBegin_ - CoinCpuTime(),
244                         parameters_.localSearchNodeLimit_);
245        milpOptimal = subMip -> optimal(); 
246        colsol = subMip->getLastSolution();
247      if(milpOptimal)
248        handler_->message(SOLVED_LOCAL_SEARCH, messages_)<<subMip->nodeCount()<<subMip->iterationCount()<<CoinMessageEol;
249      else
250        handler_->message(LOCAL_SEARCH_ABORT, messages_)<<subMip->nodeCount()<<subMip->iterationCount()<<CoinMessageEol;
251      if(colsol)
252        handler_->message(FP_MILP_VAL, messages_) 
253        <<colsol[nlp_->getNumCols()]<<CoinMessageEol;
254         
255      }/** endif localSearch*/
256      else if (subMip!=NULL) {
257        delete subMip;
258        subMip = NULL;
259        colsol = NULL;
260      }
261    }
262    if(colsol || ! milpOptimal)
263      return -DBL_MAX;
264    else
265      return DBL_MAX;
266  }
267
268  /** Register OA options.*/
269  void
270  MinlpFeasPump::registerOptions(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions)
271  {
272    roptions->SetRegisteringCategory("Options for feasibility pump", RegisteredOptions::BonminCategory);
273
274    roptions->AddBoundedIntegerOption("fp_log_level",
275        "specify FP iterations log level.",
276        0,2,1,
277        "Set the level of output of OA decomposition solver : "
278        "0 - none, 1 - normal, 2 - verbose"
279                                     );
280
281    roptions->AddLowerBoundedNumberOption("fp_log_frequency",
282        "display an update on lower and upper bounds in FP every n seconds",
283        0.,1.,100.,
284        "");
285  }
286
287/** Put objective of MIP according to FP scheme. */
288void
289MinlpFeasPump::set_fp_objective(OsiSolverInterface &si, const double * colsol) const{
290  if (objects_) {
291    for (int i = 0 ; i < nObjects_ ; i++) {
292      OsiObject * obj = objects_[i];
293      int colnum = obj->columnNumber();
294      if (colnum >= 0) {//Variable branching object
295        double round = floor(colsol[colnum] + 0.5);
296        double coeff = (colsol[colnum] - round ) < 0;
297        si.setObjCoeff(colnum, 1 - 2 * coeff);
298      }
299      else {
300        throw CoinError("OaDecompositionBase::solverManip",
301                        "setFpObjective",
302                        "Can not use FP on problem with SOS constraints");
303      }
304    }
305  }
306  else {
307    int numcols = nlp_->getNumCols();
308    for (int i = 0; i < numcols ; i++) {
309      if (nlp_->isInteger(i)){
310         double round = floor(colsol[i] + 0.5);
311         double coeff = (colsol[i] - round ) < 0;
312         si.setObjCoeff(i, 1 - 2*coeff);
313      }
314    }
315  }
316  si.initialSolve();
317}
318
319}/* End namespace Bonmin. */
Note: See TracBrowser for help on using the repository browser.