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

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

Remove debug output

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