/home/coin/SVN-release/OS-2.4.0/Bonmin/src/CbcBonmin/Heuristics/BonMilpRounding.cpp

Go to the documentation of this file.
00001 // Copyright (C) 2007, International Business Machines Corporation and others. 
00002 // All Rights Reserved.
00003 // This code is published under the Common Public License.
00004 //
00005 // Authors :
00006 // Pierre Bonami CNRS
00007 //
00008 // Date : May, 26 2010
00009 
00010 #include "BonMilpRounding.hpp"
00011 #include "CoinHelperFunctions.hpp"
00012 #include "CbcModel.hpp"
00013 #include "BonSubMipSolver.hpp"
00014 
00015 #include "CoinTime.hpp"
00016 
00017 #include <fstream>
00018 
00019 #include <iomanip>
00020 
00021 #include "CoinHelperFunctions.hpp"
00022 #include "OsiClpSolverInterface.hpp"
00023 
00024 //#define DEBUG_BON_HEURISTIC_DIVE_MIP
00025 
00026 using namespace std;
00027 
00028 namespace Bonmin
00029 {
00030   MilpRounding::MilpRounding(BonminSetup * setup)
00031     :
00032     CbcHeuristic(),
00033     setup_(setup),
00034     howOften_(20),
00035     mip_(NULL)
00036   {
00037     Initialize(setup);
00038   }
00039 
00040   void
00041   MilpRounding::Initialize(BonminSetup * b){
00042     delete mip_;
00043     mip_ = new SubMipSolver (*b, b->prefix());
00044    
00045   }
00046 
00047   MilpRounding::MilpRounding(const MilpRounding &copy)
00048     :
00049     CbcHeuristic(copy),
00050     setup_(copy.setup_),
00051     howOften_(copy.howOften_),
00052     mip_(new SubMipSolver(*copy.mip_))
00053   {
00054   }
00055 
00056   MilpRounding &
00057   MilpRounding::operator=(const MilpRounding & rhs)
00058   {
00059     if(this != &rhs) {
00060       CbcHeuristic::operator=(rhs);
00061       setup_ = rhs.setup_;
00062       howOften_ = rhs.howOften_;
00063       delete mip_;
00064       if(rhs.mip_)
00065         mip_ = new SubMipSolver(*rhs.mip_);
00066     }
00067     return *this;
00068   }
00069 
00070   MilpRounding::~MilpRounding(){
00071     delete mip_;
00072   }
00073 
00074   struct MatComp{
00075     const int * iRow;
00076     const int * jCol;
00078     bool operator()(int i,int j){
00079       return (jCol[i] < jCol[j]) || (jCol[i] == jCol[j] && iRow[i] < iRow[j]);
00080     }
00081   };
00082 
00083 
00084   int
00085   MilpRounding::solution(double &solutionValue, double *betterSolution)
00086   {
00087     printf("Entering heuristic\n");
00088     if(model_->getCurrentPassNumber() > 1) return 0;
00089     if ((model_->getNodeCount()%howOften_)!=0||model_->getCurrentPassNumber()>1)
00090       return 0;
00091  
00092     int returnCode = 0; // 0 means it didn't find a feasible solution
00093 
00094     OsiTMINLPInterface * nlp = NULL;
00095     if(setup_->getAlgorithm() == B_BB)
00096       nlp = dynamic_cast<OsiTMINLPInterface *>(model_->solver()->clone());
00097     else
00098       nlp = dynamic_cast<OsiTMINLPInterface *>(setup_->nonlinearSolver()->clone());
00099 
00100     TMINLP2TNLP* minlp = nlp->problem();
00101  
00102     // set tolerances
00103     double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
00104     double primalTolerance = 1.0e-6;
00105 
00106     int n;
00107     int m;
00108     int nnz_jac_g;
00109     int nnz_h_lag;
00110     Ipopt::TNLP::IndexStyleEnum index_style;
00111     minlp->get_nlp_info(n, m, nnz_jac_g,
00112                         nnz_h_lag, index_style);
00113 
00114     const Bonmin::TMINLP::VariableType* variableType = minlp->var_types();
00115     const double* x_sol = minlp->x_sol();
00116     const double* g_l = minlp->g_l();
00117     const double* g_u = minlp->g_u();
00118 
00119     const double * colsol = model_->solver()->getColSolution();
00120 
00121 
00122     // Get information about the linear and nonlinear part of the instance
00123     TMINLP* tminlp = nlp->model();
00124     vector<Ipopt::TNLP::LinearityType> c_lin(m);
00125     tminlp->get_constraints_linearity(m, c_lin());
00126     vector<int> c_idx(m);
00127     int n_lin = 0;
00128     for (int i=0;i<m;i++) {
00129       if (c_lin[i]==Ipopt::TNLP::LINEAR)
00130         c_idx[i] = n_lin++;
00131       else
00132         c_idx[i] = -1;
00133     }
00134 
00135 
00136     // Get the structure of the jacobian
00137     vector<int> indexRow(nnz_jac_g);
00138     vector<int> indexCol(nnz_jac_g);
00139     minlp->eval_jac_g(n, x_sol, false,
00140                       m, nnz_jac_g,
00141                       indexRow(), indexCol(), 0);
00142 
00143     // get the jacobian values 
00144     vector<double> jac_g(nnz_jac_g);
00145     minlp->eval_jac_g(n, x_sol, false,
00146                       m, nnz_jac_g,
00147                       NULL, NULL, jac_g());
00148 
00149     // Sort the matrix to column ordered
00150     vector<int> sortedIndex(nnz_jac_g);
00151     CoinIotaN(sortedIndex(), nnz_jac_g, 0);
00152     MatComp c;
00153     c.iRow = indexRow();
00154     c.jCol = indexCol();
00155     std::sort(sortedIndex.begin(), sortedIndex.end(), c);
00156 
00157     vector<int> row (nnz_jac_g);
00158     vector<double> value (nnz_jac_g);
00159     vector<int> columnStart(n,0.); 
00160     vector<int> columnLength(n,0.);
00161     int indexCorrection = (index_style == Ipopt::TNLP::C_STYLE) ? 0 : 1;
00162     int iniCol = -1;
00163     int nnz = 0;
00164     for(int i=0; i<nnz_jac_g; i++) {
00165       int thisIndexCol = indexCol[sortedIndex[i]]-indexCorrection;
00166       int thisIndexRow = c_idx[indexRow[sortedIndex[i]] - indexCorrection];
00167       if(thisIndexCol != iniCol) {
00168         iniCol = thisIndexCol;
00169         columnStart[thisIndexCol] = nnz;
00170         columnLength[thisIndexCol] = 0;
00171       }
00172       if(thisIndexRow == -1) continue;
00173       columnLength[thisIndexCol]++;
00174       row[nnz] = thisIndexRow;
00175       value[nnz] = jac_g[i];
00176       nnz++;
00177     }
00178 
00179     // Build the row lower and upper bounds
00180     vector<double> newRowLower(n_lin);
00181     vector<double> newRowUpper(n_lin);
00182     for(int i = 0 ; i < m ; i++){
00183       if(c_idx[i] == -1) continue;
00184       newRowLower[c_idx[i]] = g_l[i];
00185       newRowUpper[c_idx[i]] = g_u[i];
00186     }
00187 
00188     // Get solution array for heuristic solution
00189     vector<double> newSolution(n);
00190     std::copy(x_sol, x_sol + n, newSolution.begin());
00191 
00192     // Define the constraint matrix for MILP
00193     CoinPackedMatrix matrix(true,n_lin,n, nnz, value(), row(), columnStart(), columnLength());
00194 
00195       // create objective function and columns lower and upper bounds for MILP
00196       // and create columns for matrix in MILP
00197       double alpha = 0;
00198       double beta = 1;
00199       vector<double> objective(n);
00200       vector<int> idxIntegers;
00201       idxIntegers.reserve(n);
00202       for(int i = 0 ; i < n ; i++){
00203          if(variableType[i] != Bonmin::TMINLP::CONTINUOUS){
00204             idxIntegers.push_back(i);
00205             objective[i] = beta*(1 - 2*colsol[i]);
00206          }
00207       }
00208 
00209 #if 0
00210       // Get dual multipliers and build gradient of the lagrangean
00211       const double * duals = nlp->getRowPrice() + 2 *n;
00212       vector<double> grad(n, 0); 
00213       vector<int> indices(n, 0);
00214       tminlp->eval_grad_f(n, x_sol, false, grad());
00215       for(int i = 0 ; i < m ; i++){
00216         if(c_lin[i] == Ipopt::TNLP::LINEAR) continue;
00217         int nnz;
00218         tminlp->eval_grad_gi(n, x_sol, false, i, nnz, indices(), NULL);  
00219         tminlp->eval_grad_gi(n, x_sol, false, i, nnz, NULL, grad());
00220         for(int k = 0 ; k < nnz ; k++){
00221           objective[indices[k]] += alpha *duals[i] * grad[k];
00222         } 
00223       }
00224       for(int i = 0 ; i < n ; i++){
00225          if(variableType[i] != Bonmin::TMINLP::CONTINUOUS)
00226          objective[i] += alpha * grad[i];
00227          //if(fabs(objective[i]) < 1e-4) objective[i] = 0;
00228          else objective[i] = 0;
00229       }
00230       std::copy(objective.begin(), objective.end(), std::ostream_iterator<double>(std::cout, " "));
00231       std::cout<<std::endl;
00232 #endif
00233 
00234       // load the problem to OSI
00235       OsiSolverInterface *si = mip_->solver();
00236       assert(si != NULL);
00237       CoinMessageHandler * handler = model_->messageHandler()->clone();
00238       si->passInMessageHandler(handler);
00239       si->messageHandler()->setLogLevel(1);
00240 
00241       si->loadProblem(matrix, model_->solver()->getColLower(), model_->solver()->getColUpper(), objective(), 
00242                       newRowLower(), newRowUpper());
00243       si->setInteger(idxIntegers(), static_cast<int>(idxIntegers.size()));
00244       si->applyCuts(noGoods);
00245       printf("Done creating mip, start solving\n"); 
00246 
00247       bool hasFractionnal = true;
00248       while(hasFractionnal){
00249         mip_->optimize(DBL_MAX, 0, 60);
00250         printf("Done solving\n"); 
00251         hasFractionnal = false;
00252 #if 0
00253         bool feasible = false;
00254         if(mip_->getLastSolution()) {
00255         const double* solution = mip_->getLastSolution();
00256           std::copy(solution, solution + n, newSolution.begin());
00257         feasible = true;
00258   
00259         }
00260 
00261     if(feasible) {
00262       // fix the integer variables and solve the NLP
00263       // also add no good cut
00264       CoinPackedVector v;
00265       double lb = 1;
00266       for (int iColumn=0;iColumn<n;iColumn++) {
00267         if (variableType[iColumn] != Bonmin::TMINLP::CONTINUOUS) {
00268           double value=newSolution[iColumn];
00269           if (fabs(floor(value+0.5)-value)>integerTolerance) {
00270 #ifdef DEBUG_BON_HEURISTIC_DIVE_MIP
00271             cout<<"It should be infeasible because: "<<endl;
00272             cout<<"variable "<<iColumn<<" is not integer"<<endl;
00273 #endif
00274             feasible = false;
00275             break;
00276           }
00277           else {
00278             value=floor(newSolution[iColumn]+0.5);
00279             if(value > 0.5){
00280               v.insert(iColumn, -1);
00281               lb -= value;
00282             }
00283             minlp->SetVariableUpperBound(iColumn, value);
00284             minlp->SetVariableLowerBound(iColumn, value);
00285           }
00286         }
00287       }
00288       }
00289 #endif
00290       }
00291       bool feasible = false;
00292       if(mip_->getLastSolution()) {
00293         const double* solution = mip_->getLastSolution();
00294         std::copy(solution, solution + n, newSolution.begin());
00295         feasible = true;
00296 
00297         delete handler;
00298       }
00299       
00300 
00301     if(feasible) {
00302       // fix the integer variables and solve the NLP
00303       // also add no good cut
00304       printf("We found a solution");
00305       CoinPackedVector v;
00306       double lb = 1;
00307       for (int iColumn=0;iColumn<n;iColumn++) {
00308         if (variableType[iColumn] != Bonmin::TMINLP::CONTINUOUS) {
00309           double value=newSolution[iColumn];
00310           if (fabs(floor(value+0.5)-value)>integerTolerance) {
00311 #ifdef DEBUG_BON_HEURISTIC_DIVE_MIP
00312             cout<<"It should be infeasible because: "<<endl;
00313             cout<<"variable "<<iColumn<<" is not integer"<<endl;
00314 #endif
00315             feasible = false;
00316             break;
00317           }
00318           else {
00319             value=floor(newSolution[iColumn]+0.5);
00320             if(value > 0.5){
00321               v.insert(iColumn, -1);
00322               lb -= value;
00323             }
00324             minlp->SetVariableUpperBound(iColumn, value);
00325             minlp->SetVariableLowerBound(iColumn, value);
00326           }
00327         }
00328       }
00329       OsiRowCut c;
00330       c.setRow(v);
00331       c.setLb(lb);
00332       c.setUb(DBL_MAX);
00333       noGoods.insert(c);
00334       if(feasible) {
00335         nlp->initialSolve();
00336         if(minlp->optimization_status() != Ipopt::SUCCESS) {
00337           printf("Infeasible for NLP");
00338           feasible = false;
00339         }
00340         std::copy(x_sol,x_sol+n, newSolution.begin());
00341       }
00342     }
00343     else {printf("No solution found\n");}
00344     if(feasible) {
00345       double newSolutionValue;
00346       minlp->eval_f(n, newSolution(), true, newSolutionValue); 
00347       if(newSolutionValue < solutionValue) {
00348         std::copy(newSolution.begin(), newSolution.end(), betterSolution);
00349         solutionValue = newSolutionValue;
00350         returnCode = 1;
00351       }
00352     }
00353 
00354     delete nlp;
00355 
00356 #ifdef DEBUG_BON_HEURISTIC_DIVE_MIP
00357     std::cout<<"DiveMIP returnCode = "<<returnCode<<std::endl;
00358 #endif
00359 
00360     return returnCode;
00361   }
00362   void
00363   MilpRounding::registerOptions(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions){
00364     roptions->SetRegisteringCategory("Undocumented Heuristics", RegisteredOptions::UndocumentedCategory);
00365    roptions->AddStringOption2(
00366      "MILP_rounding_heuristic",
00367      "if yes runs the heuristic",
00368      "no",
00369      "no", "don't run it",
00370      "yes", "runs the heuristic",
00371      "");
00372   }
00373 }

Generated on Thu Sep 22 03:05:54 2011 by  doxygen 1.4.7