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

Last change on this file since 1536 was 1536, checked in by pbonami, 10 years ago

Port Claudia's code to trunk

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