00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "BonHeuristicInnerApproximation.hpp"
00014 #include "CoinHelperFunctions.hpp"
00015 #include "CbcModel.hpp"
00016 #include "BonSubMipSolver.hpp"
00017 #include "BonCbcLpStrategy.hpp"
00018
00019 #ifdef COIN_HAS_CPX
00020 #include "OsiCpxSolverInterface.hpp"
00021 #endif
00022
00023 #include "OsiClpSolverInterface.hpp"
00024
00025 #include "OsiAuxInfo.hpp"
00026
00027 #include "CoinTime.hpp"
00028
00029 #include <fstream>
00030
00031 #include <iomanip>
00032
00033 #include "CoinHelperFunctions.hpp"
00034
00035 #define DEBUG_BON_HEURISTIC
00036
00037 using namespace std;
00038
00039 namespace Bonmin {
00040
00041 HeuristicInnerApproximation::HeuristicInnerApproximation(BonminSetup * setup) :
00042 CbcHeuristic(), setup_(setup), howOften_(100), mip_(NULL),
00043 nbAp_(50) {
00044 Initialize(setup);
00045 }
00046
00047 HeuristicInnerApproximation::HeuristicInnerApproximation(
00048 const HeuristicInnerApproximation ©) :
00049 CbcHeuristic(copy),
00050 setup_(copy.setup_),
00051 howOften_(copy.howOften_),
00052 mip_(new SubMipSolver(*copy.mip_)),
00053 nbAp_(copy.nbAp_) {
00054 }
00055
00056 HeuristicInnerApproximation &
00057 HeuristicInnerApproximation::operator=(const HeuristicInnerApproximation& rhs) {
00058 if (this != &rhs) {
00059 CbcHeuristic::operator=(rhs);
00060 setup_ = rhs.setup_;
00061 howOften_ = rhs.howOften_;
00062 nbAp_ = rhs.nbAp_;
00063 delete mip_;
00064 if (rhs.mip_)
00065 mip_ = new SubMipSolver(*rhs.mip_);
00066 }
00067 return *this;
00068 }
00069
00070 void HeuristicInnerApproximation::registerOptions(Ipopt::SmartPtr<
00071 Bonmin::RegisteredOptions> roptions) {
00072 roptions->SetRegisteringCategory("Initial Approximations descriptions",
00073 RegisteredOptions::UndocumentedCategory);
00074 roptions->AddStringOption2("heuristic_inner_approximation",
00075 "if yes runs the InnerApproximation heuristic", "yes", "no",
00076 "don't run it", "yes", "runs the heuristic", "");
00077
00078 roptions->setOptionExtraInfo("heuristic_inner_approximation", 63);
00079 }
00080
00081 void
00082 HeuristicInnerApproximation::Initialize(BonminSetup * b) {
00083
00084 delete mip_;
00085 mip_ = new SubMipSolver (*b, "inner_approximation");
00086 b->options()->GetIntegerValue("number_approximations_initial_outer",
00087 nbAp_, b->prefix());
00088 }
00089
00090 HeuristicInnerApproximation::~HeuristicInnerApproximation() {
00091 delete mip_;
00092 }
00093
00097 int
00098 HeuristicInnerApproximation::solution(double &solutionValue, double *betterSolution)
00099 {
00100 if(model_->getNodeCount() || model_->getCurrentPassNumber() > 1) return 0;
00101 if ((model_->getNodeCount()%howOften_)!=0||model_->getCurrentPassNumber()>1)
00102 return 0;
00103
00104 int returnCode = 0;
00105
00106 OsiTMINLPInterface * nlp = NULL;
00107 if(setup_->getAlgorithm() == B_BB)
00108 nlp = dynamic_cast<OsiTMINLPInterface *>(model_->solver()->clone());
00109 else
00110 nlp = dynamic_cast<OsiTMINLPInterface *>(setup_->nonlinearSolver()->clone());
00111
00112 TMINLP2TNLP* minlp = nlp->problem();
00113
00114 double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
00115
00116 int numberColumns;
00117 int numberRows;
00118 int nnz_jac_g;
00119 int nnz_h_lag;
00120 Ipopt::TNLP::IndexStyleEnum index_style;
00121 minlp->get_nlp_info(numberColumns, numberRows, nnz_jac_g,
00122 nnz_h_lag, index_style);
00123
00124 const Bonmin::TMINLP::VariableType* variableType = minlp->var_types();
00125
00126 const double* x_sol = minlp->x_sol();
00127
00128 double* newSolution = new double [numberColumns];
00129 memcpy(newSolution,x_sol,numberColumns*sizeof(double));
00130 double* new_g_sol = new double [numberRows];
00131
00132 bool feasible = true;
00133
00134 #ifdef DEBUG_BON_HEURISTIC
00135 cout << "Loading the problem to OSI\n";
00136 #endif
00137 OsiSolverInterface *si = mip_->solver();
00138
00139 bool delete_si = false;
00140 if(si == NULL) {
00141 si = new OsiClpSolverInterface;
00142 mip_->setLpSolver(si);
00143 delete_si = true;
00144 }
00145 CoinMessageHandler * handler = model_->messageHandler()->clone();
00146 si->passInMessageHandler(handler);
00147 si->messageHandler()->setLogLevel(2);
00148 #ifdef DEBUG_BON_HEURISTIC
00149 cout << "Loading problem into si\n";
00150 #endif
00151 extractInnerApproximation(*nlp, *si, newSolution, true);
00152 #ifdef DEBUG_BON_HEURISTIC
00153 cout << "problem loaded\n";
00154 cout << "**** Running optimization ****\n";
00155 #endif
00156 mip_->optimize(DBL_MAX, 2, 180);
00157 #ifdef DEBUG_BON_HEURISTIC
00158 cout << "Optimization finished\n";
00159 #endif
00160 if(mip_->getLastSolution()) {
00161 const double* solution = mip_->getLastSolution();
00162 for (size_t iLCol=0;iLCol<numberColumns;iLCol++) {
00163 newSolution[iLCol] = solution[iLCol];
00164 }
00165 }
00166 else
00167 feasible = false;
00168
00169 if(delete_si) {
00170 delete si;
00171 }
00172 delete handler;
00173
00174 const double* x_l = minlp->x_l();
00175 const double* x_u = minlp->x_u();
00176 const double* g_l = minlp->g_l();
00177 const double* g_u = minlp->g_u();
00178 double primalTolerance = 1.0e-6;
00179 #if 0 // Set to 1 if you need to test the feasibility of the returned solution
00180 vector<Ipopt::TNLP::LinearityType> constTypes(numberRows);
00181 minlp->get_constraints_linearity(numberRows, constTypes());
00182 feasible = true;
00183 for (int iColumn=0;iColumn<numberColumns;iColumn++) {
00184 double value=newSolution[iColumn];
00185 if(value - x_l[iColumn] < -1e-8|| value - x_u[iColumn] > 1e-8) {
00186 #ifdef DEBUG_BON_HEURISTIC
00187 cout<<"It should be infeasible because: "<<endl;
00188 cout<<"x_l["<<iColumn<<"]= "<<x_l[iColumn]<<" "
00189 <<"x_sol["<<iColumn<<"]= "<<value<<" "
00190 <<"x_u["<<iColumn<<"]= "<<x_u[iColumn]<<endl;
00191 #endif
00192 feasible = false;
00193 break;
00194 }
00195 if (variableType[iColumn] != Bonmin::TMINLP::CONTINUOUS) {
00196 if (fabs(floor(value+0.5)-value)>integerTolerance) {
00197 #ifdef DEBUG_BON_HEURISTIC
00198 cout<<"It should be infeasible because: "<<endl;
00199 cout<<"x["<<iColumn<<"]= "<<value<<" Should be integer!"<<"\\Integer tolerance = "<<integerTolerance<<"\n";
00200 #endif
00201 feasible = false;
00202 break;
00203 }
00204 }
00205 }
00206 minlp->eval_g(numberColumns, newSolution, true,
00207 numberRows, new_g_sol);
00208 for(int iRow=0; iRow<numberRows; iRow++) {
00209 if(new_g_sol[iRow]<g_l[iRow]-primalTolerance ||
00210 new_g_sol[iRow]>g_u[iRow]+primalTolerance) {
00211 if(minlp->optimization_status() != SUCCESS) {
00212 feasible = false;
00213 break;
00214 } else {
00215 #ifdef DEBUG_BON_HEURISTIC
00216 cout<<"It should be infeasible because: "<<endl;
00217 cout<<"g_l["<<iRow<<"]= "<<g_l[iRow]<<" "
00218 <<"g_sol["<<iRow<<"]= "<<new_g_sol[iRow]<<" "
00219 <<"g_u["<<iRow<<"]= "<<g_u[iRow]<<endl;
00220 cout<<"primalTolerance= "<<primalTolerance<<endl;
00221 if(constTypes[iRow] == TNLP::NON_LINEAR)
00222 cout<<"nonLinear constraint number "<<iRow<<endl;
00223 #endif
00224 feasible = false;
00225 }
00226 }
00227 }
00228 #ifdef DEBUG_BON_HEURISTIC
00229 cout<<"Every thing is feasible"<<endl;
00230 #endif
00231
00232
00233
00234 #if 0
00235 if(feasible ) {
00236
00237 std::vector<double> memLow(numberColumns);
00238 std::vector<double> memUpp(numberColumns);
00239 std::copy(minlp->x_l(), minlp->x_l() + numberColumns, memLow.begin());
00240 std::copy(minlp->x_u(), minlp->x_u() + numberColumns, memUpp.begin());
00241
00242 for (int iColumn=0;iColumn<numberColumns;iColumn++) {
00243 double value=floor(newSolution[iColumn]+0.5);
00244 minlp->SetVariableUpperBound(iColumn, value);
00245 minlp->SetVariableLowerBound(iColumn, value);
00246 }
00247 if(feasible) {
00248 nlp->initialSolve();
00249 if(minlp->optimization_status() != SUCCESS) {
00250 feasible = false;
00251 }
00252 memcpy(newSolution,minlp->x_sol(),numberColumns*sizeof(double));
00253 }
00254
00255
00256 for (int iColumn=0;iColumn<numberColumns;iColumn++) {
00257 if (variableType[iColumn] != Bonmin::TMINLP::CONTINUOUS) {
00258 minlp->SetVariableUpperBound(iColumn, memUpp[iColumn]);
00259 minlp->SetVariableLowerBound(iColumn, memLow[iColumn]);
00260 }
00261 }
00262 }
00263 #endif
00264 #endif
00265
00266 if(feasible) {
00267 double newSolutionValue;
00268 minlp->eval_f(numberColumns, newSolution, true, newSolutionValue);
00269 if(newSolutionValue < solutionValue) {
00270 memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
00271 solutionValue = newSolutionValue;
00272 returnCode = 1;
00273 }
00274 }
00275
00276 delete [] newSolution;
00277 delete [] new_g_sol;
00278
00279 delete nlp;
00280
00281 #ifdef DEBUG_BON_HEURISTIC
00282 std::cout<<"Inner approximation returnCode = "<<returnCode<<std::endl;
00283 #endif
00284 return returnCode;
00285 }
00286
00289 bool
00290 HeuristicInnerApproximation::getMyInnerApproximation(OsiTMINLPInterface &si, OsiCuts &cs, int ind,
00291 const double * x, const double * x2) {
00292
00293 int n, m, nnz_jac_g, nnz_h_lag;
00294 Ipopt::TNLP::IndexStyleEnum index_style;
00295 TMINLP2TNLP * problem = si.problem();
00296 problem->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
00297
00298
00299 CoinPackedVector cut;
00300 double lb = 0;
00301 double ub = 0;
00302
00303 double infty = si.getInfinity();
00304
00305 lb = -infty;
00306
00307 double g = 0;
00308 double g2 = 0;
00309 double diff = 0;
00310 double a = 0;
00311 problem->eval_gi(n, x, 1, ind, g);
00312 problem->eval_gi(n, x2, 1, ind, g2);
00313 vector<int> jCol(n);
00314 int nnz;
00315 problem->eval_grad_gi(n, x2, 0, ind, nnz, jCol(), NULL);
00316 vector<double> jValues(nnz);
00317 problem->eval_grad_gi(n, x2, 0, ind, nnz, NULL, jValues());
00318
00319 bool add = false;
00320 for (int i = 0; i < nnz; i++) {
00321 const int &colIdx = jCol[i];
00322 if(index_style == Ipopt::TNLP::FORTRAN_STYLE) jCol[i]--;
00323 diff = x[colIdx] - x2[colIdx];
00324
00325 if (fabs(diff) >= 1e-8) {
00326 a = (g - g2) / diff;
00327 cut.insert(colIdx, a);
00328 ub = a * x[colIdx] - g;
00329 add = true;
00330 } else
00331 cut.insert(colIdx, jValues[i]);
00332 }
00333
00334 if (add) {
00335
00336 OsiRowCut newCut;
00337 newCut.setGloballyValidAsInteger(1);
00338 newCut.setLb(lb);
00339 newCut.setUb(ub);
00340 newCut.setRow(cut);
00341 cs.insert(newCut);
00342 return true;
00343 }
00344 return false;
00345 }
00346 void
00347 HeuristicInnerApproximation::extractInnerApproximation(OsiTMINLPInterface & nlp, OsiSolverInterface &si,
00348 const double * x, bool getObj) {
00349 int n;
00350 int m;
00351 int nnz_jac_g;
00352 int nnz_h_lag;
00353 Ipopt::TNLP::IndexStyleEnum index_style;
00354 TMINLP2TNLP * problem = nlp.problem();
00355
00356 problem->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
00357
00358 vector<int> jRow(nnz_jac_g);
00359 vector<int> jCol(nnz_jac_g);
00360 vector<double> jValues(nnz_jac_g);
00361 problem->eval_jac_g(n, NULL, 0, m, nnz_jac_g, jRow(), jCol(), NULL);
00362 if(index_style == Ipopt::TNLP::FORTRAN_STYLE)
00363 {
00364 for(int i = 0 ; i < nnz_jac_g ; i++){
00365 jRow[i]--;
00366 jCol[i]--;
00367 }
00368 }
00369
00370
00371 problem->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL,
00372 jValues());
00373
00374 vector<double> g(m);
00375 problem->eval_g(n, x, 1, m, g());
00376
00377 vector<int> nonLinear(m);
00378
00379 int numNonLinear = 0;
00380 const double * rowLower = nlp.getRowLower();
00381 const double * rowUpper = nlp.getRowUpper();
00382 const double * colLower = nlp.getColLower();
00383 const double * colUpper = nlp.getColUpper();
00384 assert(m == nlp.getNumRows());
00385 double infty = si.getInfinity();
00386 double nlp_infty = nlp.getInfinity();
00387 vector<Ipopt::TNLP::LinearityType> constTypes(m);
00388 problem->get_constraints_linearity(m, constTypes());
00389 for (int i = 0; i < m; i++) {
00390 if (constTypes[i] == Ipopt::TNLP::NON_LINEAR) {
00391 nonLinear[numNonLinear++] = i;
00392 }
00393 }
00394 vector<double> rowLow(m - numNonLinear);
00395 vector<double> rowUp(m - numNonLinear);
00396 int ind = 0;
00397 for (int i = 0; i < m; i++) {
00398 if (constTypes[i] != Ipopt::TNLP::NON_LINEAR) {
00399 if (rowLower[i] > -nlp_infty) {
00400
00401 rowLow[ind] = (rowLower[i]);
00402 } else
00403 rowLow[ind] = -infty;
00404 if (rowUpper[i] < nlp_infty) {
00405
00406 rowUp[ind] = (rowUpper[i]);
00407 } else
00408 rowUp[ind] = infty;
00409 ind++;
00410 }
00411
00412 }
00413
00414 CoinPackedMatrix mat(true, jRow(), jCol(), jValues(), nnz_jac_g);
00415 mat.setDimensions(m, n);
00416
00417
00418 mat.deleteRows(numNonLinear, nonLinear());
00419
00420 int numcols = nlp.getNumCols();
00421 vector<double> obj(numcols);
00422 for (int i = 0; i < numcols; i++)
00423 obj[i] = 0.;
00424
00425 si.loadProblem(mat, nlp.getColLower(), nlp.getColUpper(),
00426 obj(), rowLow(), rowUp());
00427 const Bonmin::TMINLP::VariableType* variableType = problem->var_types();
00428 for (int i = 0; i < n; i++) {
00429 if ((variableType[i] == TMINLP::BINARY) || (variableType[i]
00430 == TMINLP::INTEGER))
00431 si.setInteger(i);
00432 }
00433 if (getObj) {
00434 bool addObjVar = false;
00435 if (problem->hasLinearObjective()) {
00436 double zero;
00437 vector<double> x0(n, 0.);
00438 problem->eval_f(n, x0(), 1, zero);
00439 si.setDblParam(OsiObjOffset, -zero);
00440
00441 problem->eval_grad_f(n, x, 1, obj());
00442 si.setObjective(obj());
00443 } else {
00444 addObjVar = true;
00445 }
00446
00447 if (addObjVar) {
00448 nlp.addObjectiveFunction(si, x);
00449 }
00450 }
00451
00452
00453 int InnerDesc = 1;
00454 if (InnerDesc == 1) {
00455 OsiCuts cs;
00456
00457 double * p = CoinCopyOfArray(colLower, n);
00458 double * pp = CoinCopyOfArray(colLower, n);
00459 double * up = CoinCopyOfArray(colUpper, n);
00460
00461 const int& nbAp = nbAp_;
00462 std::vector<int> nbG(m, 0);
00463
00464 std::vector<double> step(n);
00465
00466 for (int i = 0; i < n; i++) {
00467
00468 if (colUpper[i] > 1e08) {
00469 up[i] = 0;
00470 }
00471
00472 if (colUpper[i] > 1e08 || colLower[i] < -1e08 || (variableType[i]
00473 == TMINLP::BINARY) || (variableType[i] == TMINLP::INTEGER)) {
00474 step[i] = 0;
00475 } else
00476 step[i] = (up[i] - colLower[i]) / (nbAp);
00477
00478 if (colLower[i] < -1e08) {
00479 p[i] = 0;
00480 pp[i] = 0;
00481 }
00482
00483 }
00484 vector<double> g_p(m);
00485 vector<double> g_pp(m);
00486
00487 for (int j = 1; j <= nbAp; j++) {
00488
00489 for (int i = 0; i < n; i++) {
00490 pp[i] += step[i];
00491 }
00492
00493 problem->eval_g(n, p, 1, m, g_p());
00494 problem->eval_g(n, pp, 1, m, g_pp());
00495 double diff = 0;
00496 int varInd = 0;
00497 for (int i = 0; (i < m && constTypes[i] == Ipopt::TNLP::NON_LINEAR); i++) {
00498 if (varInd == n - 1)
00499 varInd = 0;
00500 diff = std::abs(g_p[i] - g_pp[i]);
00501 if (nbG[i] < nbAp - 1) {
00502 getMyInnerApproximation(nlp, cs, i, p, pp);
00503 p[varInd] = pp[varInd];
00504 nbG[i]++;
00505 }
00506 varInd++;
00507 }
00508 }
00509
00510 for(int i = 0; (i< m && constTypes[i] == Ipopt::TNLP::NON_LINEAR); i++) {
00511
00512 getMyInnerApproximation(nlp, cs, i, p, up);
00513 }
00514
00515 delete [] p;
00516 delete [] pp;
00517 delete [] up;
00518 si.applyCuts(cs);
00519 }
00520 }
00521
00522 }
00523