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 1
00180 if(feasible ) {
00181
00182 std::vector<double> memLow(numberColumns);
00183 std::vector<double> memUpp(numberColumns);
00184 std::copy(minlp->x_l(), minlp->x_l() + numberColumns, memLow.begin());
00185 std::copy(minlp->x_u(), minlp->x_u() + numberColumns, memUpp.begin());
00186
00187 for (int iColumn=0;iColumn<numberColumns;iColumn++) {
00188 if (variableType[iColumn] != Bonmin::TMINLP::CONTINUOUS) {
00189 double value=floor(newSolution[iColumn]+0.5);
00190 minlp->SetVariableUpperBound(iColumn, value);
00191 minlp->SetVariableLowerBound(iColumn, value);
00192 }
00193 }
00194 if(feasible) {
00195 nlp->initialSolve();
00196 if(minlp->optimization_status() != Ipopt::SUCCESS) {
00197 feasible = false;
00198 }
00199 memcpy(newSolution,minlp->x_sol(),numberColumns*sizeof(double));
00200 }
00201
00202
00203 for (int iColumn=0;iColumn<numberColumns;iColumn++) {
00204 if (variableType[iColumn] != Bonmin::TMINLP::CONTINUOUS) {
00205 minlp->SetVariableUpperBound(iColumn, memUpp[iColumn]);
00206 minlp->SetVariableLowerBound(iColumn, memLow[iColumn]);
00207 }
00208 }
00209 }
00210 #endif
00211 #endif
00212
00213 if(feasible) {
00214 double newSolutionValue;
00215 minlp->eval_f(numberColumns, newSolution, true, newSolutionValue);
00216 if(newSolutionValue < solutionValue) {
00217 memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
00218 solutionValue = newSolutionValue;
00219 returnCode = 1;
00220 }
00221 }
00222
00223 delete [] newSolution;
00224 delete [] new_g_sol;
00225
00226 delete nlp;
00227
00228 #ifdef DEBUG_BON_HEURISTIC
00229 std::cout<<"Inner approximation returnCode = "<<returnCode<<std::endl;
00230 #endif
00231 return returnCode;
00232 }
00233
00236 bool
00237 HeuristicInnerApproximation::getMyInnerApproximation(OsiTMINLPInterface &si, OsiCuts &cs, int ind,
00238 const double * x, const double * x2) {
00239
00240 int n, m, nnz_jac_g, nnz_h_lag;
00241 Ipopt::TNLP::IndexStyleEnum index_style;
00242 TMINLP2TNLP * problem = si.problem();
00243 problem->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
00244
00245
00246 CoinPackedVector cut;
00247 double lb = 0;
00248 double ub = 0;
00249
00250 double infty = si.getInfinity();
00251
00252 lb = -infty;
00253
00254 double g = 0;
00255 double g2 = 0;
00256 double diff = 0;
00257 double a = 0;
00258 problem->eval_gi(n, x, 1, ind, g);
00259 problem->eval_gi(n, x2, 1, ind, g2);
00260 vector<int> jCol(n);
00261 int nnz;
00262 problem->eval_grad_gi(n, x2, 0, ind, nnz, jCol(), NULL);
00263 vector<double> jValues(nnz);
00264 problem->eval_grad_gi(n, x2, 0, ind, nnz, NULL, jValues());
00265 bool add = false;
00266 for (int i = 0; i < nnz; i++) {
00267 const int &colIdx = jCol[i];
00268 if(index_style == Ipopt::TNLP::FORTRAN_STYLE) jCol[i]--;
00269 diff = x[colIdx] - x2[colIdx];
00270
00271 if (fabs(diff) >= 1e-8) {
00272 a = (g - g2) / diff;
00273 cut.insert(colIdx, a);
00274 ub = a * x[colIdx] - g;
00275 add = true;
00276 } else
00277 cut.insert(colIdx, jValues[i]);
00278 }
00279
00280 if (add) {
00281
00282 OsiRowCut newCut;
00283 newCut.setGloballyValidAsInteger(1);
00284 newCut.setLb(lb);
00285
00286
00287 int* ids = problem->get_const_xtra_id();
00288 int binary_id = (ids == NULL) ? -1 : ids[ind];
00289 if(binary_id>0) {
00290 cut.insert(binary_id, -ub);
00291 newCut.setUb(0);
00292 }
00293 else
00294 newCut.setUb(ub);
00295
00296
00297 newCut.setRow(cut);
00298 cs.insert(newCut);
00299 return true;
00300 }
00301 return false;
00302 }
00303 void
00304 HeuristicInnerApproximation::extractInnerApproximation(OsiTMINLPInterface & nlp, OsiSolverInterface &si,
00305 const double * x, bool getObj) {
00306 int n;
00307 int m;
00308 int nnz_jac_g;
00309 int nnz_h_lag;
00310 Ipopt::TNLP::IndexStyleEnum index_style;
00311 TMINLP2TNLP * problem = nlp.problem();
00312
00313 problem->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
00314
00315 vector<int> jRow(nnz_jac_g);
00316 vector<int> jCol(nnz_jac_g);
00317 vector<double> jValues(nnz_jac_g);
00318 problem->eval_jac_g(n, NULL, 0, m, nnz_jac_g, jRow(), jCol(), NULL);
00319 if(index_style == Ipopt::TNLP::FORTRAN_STYLE)
00320 {
00321 for(int i = 0 ; i < nnz_jac_g ; i++){
00322 jRow[i]--;
00323 jCol[i]--;
00324 }
00325 }
00326
00327
00328 problem->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL,
00329 jValues());
00330
00331 vector<double> g(m);
00332 problem->eval_g(n, x, 1, m, g());
00333
00334 vector<int> nonLinear(m);
00335
00336 int numNonLinear = 0;
00337 const double * rowLower = nlp.getRowLower();
00338 const double * rowUpper = nlp.getRowUpper();
00339 const double * colLower = nlp.getColLower();
00340 const double * colUpper = nlp.getColUpper();
00341 assert(m == nlp.getNumRows());
00342 double infty = si.getInfinity();
00343 double nlp_infty = nlp.getInfinity();
00344 vector<Ipopt::TNLP::LinearityType> constTypes(m);
00345 problem->get_constraints_linearity(m, constTypes());
00346 for (int i = 0; i < m; i++) {
00347 if (constTypes[i] == Ipopt::TNLP::NON_LINEAR) {
00348 nonLinear[numNonLinear++] = i;
00349 }
00350 }
00351 vector<double> rowLow(m - numNonLinear);
00352 vector<double> rowUp(m - numNonLinear);
00353 int ind = 0;
00354 for (int i = 0; i < m; i++) {
00355 if (constTypes[i] != Ipopt::TNLP::NON_LINEAR) {
00356 if (rowLower[i] > -nlp_infty) {
00357
00358 rowLow[ind] = (rowLower[i]);
00359 } else
00360 rowLow[ind] = -infty;
00361 if (rowUpper[i] < nlp_infty) {
00362
00363 rowUp[ind] = (rowUpper[i]);
00364 } else
00365 rowUp[ind] = infty;
00366 ind++;
00367 }
00368
00369 }
00370
00371 CoinPackedMatrix mat(true, jRow(), jCol(), jValues(), nnz_jac_g);
00372 mat.setDimensions(m, n);
00373
00374
00375 mat.deleteRows(numNonLinear, nonLinear());
00376
00377 int numcols = nlp.getNumCols();
00378 vector<double> obj(numcols);
00379 for (int i = 0; i < numcols; i++)
00380 obj[i] = 0.;
00381
00382 si.loadProblem(mat, nlp.getColLower(), nlp.getColUpper(),
00383 obj(), rowLow(), rowUp());
00384 const Bonmin::TMINLP::VariableType* variableType = problem->var_types();
00385 for (int i = 0; i < n; i++) {
00386 if ((variableType[i] == TMINLP::BINARY) || (variableType[i]
00387 == TMINLP::INTEGER))
00388 si.setInteger(i);
00389 }
00390 if (getObj) {
00391 bool addObjVar = false;
00392 if (problem->hasLinearObjective()) {
00393 double zero;
00394 vector<double> x0(n, 0.);
00395 problem->eval_f(n, x0(), 1, zero);
00396 si.setDblParam(OsiObjOffset, -zero);
00397
00398 problem->eval_grad_f(n, x, 1, obj());
00399 si.setObjective(obj());
00400 } else {
00401 addObjVar = true;
00402 }
00403
00404 if (addObjVar) {
00405 nlp.addObjectiveFunction(si, x);
00406 }
00407 }
00408
00409
00410 int InnerDesc = 1;
00411 if (InnerDesc == 1) {
00412 OsiCuts cs;
00413
00414 double * p = CoinCopyOfArray(colLower, n);
00415 double * pp = CoinCopyOfArray(colLower, n);
00416 double * up = CoinCopyOfArray(colUpper, n);
00417
00418 const int& nbAp = nbAp_;
00419 std::vector<int> nbG(m, 0);
00420
00421 std::vector<double> step(n);
00422
00423 for (int i = 0; i < n; i++) {
00424
00425 if (colUpper[i] > 1e08) {
00426 up[i] = 0;
00427 }
00428
00429 if (colUpper[i] > 1e08 || colLower[i] < -1e08 || (variableType[i]
00430 == TMINLP::BINARY) || (variableType[i] == TMINLP::INTEGER)) {
00431 step[i] = 0;
00432 } else
00433 step[i] = (up[i] - colLower[i]) / (nbAp);
00434
00435 if (colLower[i] < -1e08) {
00436 p[i] = 0;
00437 pp[i] = 0;
00438 }
00439
00440 }
00441 vector<double> g_p(m);
00442 vector<double> g_pp(m);
00443
00444 for (int j = 1; j <= nbAp; j++) {
00445
00446 for (int i = 0; i < n; i++) {
00447 pp[i] += step[i];
00448 }
00449
00450 problem->eval_g(n, p, 1, m, g_p());
00451 problem->eval_g(n, pp, 1, m, g_pp());
00452 double diff = 0;
00453 int varInd = 0;
00454 for (int i = 0; (i < m && constTypes[i] == Ipopt::TNLP::NON_LINEAR); i++) {
00455 if (varInd == n - 1)
00456 varInd = 0;
00457 diff = std::abs(g_p[i] - g_pp[i]);
00458 if (nbG[i] < nbAp - 1) {
00459 getMyInnerApproximation(nlp, cs, i, p, pp);
00460 p[varInd] = pp[varInd];
00461 nbG[i]++;
00462 }
00463 varInd++;
00464 }
00465 }
00466
00467 for(int i = 0; (i< m && constTypes[i] == Ipopt::TNLP::NON_LINEAR); i++) {
00468
00469 getMyInnerApproximation(nlp, cs, i, p, up);
00470 }
00471
00472 delete [] p;
00473 delete [] pp;
00474 delete [] up;
00475 si.applyCuts(cs);
00476 }
00477 }
00478
00479 }
00480