00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "SepaHeuristicInnerApproximation.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
00036
00037 namespace Sepa{
00038
00039 HeuristicInnerApproximation::HeuristicInnerApproximation(Bonmin::BonminSetup * setup) :
00040 CbcHeuristic(), setup_(setup), howOften_(100), mip_(NULL),
00041 nbAp_(50), time_limit_(10.) {
00042 Initialize(setup);
00043 }
00044
00045 HeuristicInnerApproximation::HeuristicInnerApproximation(
00046 const HeuristicInnerApproximation ©) :
00047 CbcHeuristic(copy),
00048 setup_(copy.setup_),
00049 howOften_(copy.howOften_),
00050 mip_(new Bonmin::SubMipSolver(*copy.mip_)),
00051 nbAp_(copy.nbAp_),
00052 time_limit_(copy.time_limit_)
00053 {
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 Bonmin::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 Bonmin::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 roptions->AddLowerBoundedIntegerOption("number_inner_approximation_points",
00081 "Set the number of points to use for linear inner approximation of nonlinear functions in heuristic",
00082 1, 20);
00083 roptions->setOptionExtraInfo("number_inner_approximation_points", 63);
00084
00085 roptions->AddLowerBoundedNumberOption("inner_time_limit",
00086 "Time limit for inner approximation",
00087 0, true, 10, "");
00088 roptions->setOptionExtraInfo("number_inner_approximation_points", 63);
00089 }
00090
00091 void
00092 HeuristicInnerApproximation::Initialize(Bonmin::BonminSetup * b) {
00093
00094 delete mip_;
00095 mip_ = new Bonmin::SubMipSolver (*b, "inner_approximation");
00096 b->options()->GetIntegerValue("number_inner_approximation_points",
00097 nbAp_, b->prefix());
00098 b->options()->GetNumericValue("inner_time_limit",
00099 time_limit_, b->prefix());
00100 }
00101
00102 HeuristicInnerApproximation::~HeuristicInnerApproximation() {
00103 delete mip_;
00104 }
00105
00109 int
00110 HeuristicInnerApproximation::solution(double &solutionValue, double *betterSolution)
00111 {
00112 if(model_->getNodeCount() || model_->getCurrentPassNumber() > 1) return 0;
00113 if ((model_->getNodeCount()%howOften_)!=0||model_->getCurrentPassNumber()>1)
00114 return 0;
00115
00116 int returnCode = 0;
00117
00118 Bonmin::OsiTMINLPInterface * nlp = NULL;
00119 if(setup_->getAlgorithm() == Bonmin::B_BB)
00120 nlp = dynamic_cast<Bonmin::OsiTMINLPInterface *>(model_->solver()->clone());
00121 else
00122 nlp = dynamic_cast<Bonmin::OsiTMINLPInterface *>(setup_->nonlinearSolver()->clone());
00123
00124 Bonmin::TMINLP2TNLP* minlp = nlp->problem();
00125
00126
00127
00128
00129 int numberColumns;
00130 int numberRows;
00131 int nnz_jac_g;
00132 int nnz_h_lag;
00133 Ipopt::TNLP::IndexStyleEnum index_style;
00134 minlp->get_nlp_info(numberColumns, numberRows, nnz_jac_g,
00135 nnz_h_lag, index_style);
00136
00137
00138
00139 const double* x_sol = minlp->x_sol();
00140
00141 double* newSolution = new double [numberColumns];
00142 memcpy(newSolution,x_sol,numberColumns*sizeof(double));
00143 double* new_g_sol = new double [numberRows];
00144
00145 bool feasible = true;
00146
00147 #ifdef DEBUG_BON_HEURISTIC
00148 std::cout << "Loading the problem to OSI\n";
00149 #endif
00150 OsiSolverInterface *si = mip_->solver();
00151
00152 bool delete_si = false;
00153 if(si == NULL) {
00154 si = new OsiClpSolverInterface;
00155 mip_->setLpSolver(si);
00156 delete_si = true;
00157 }
00158 CoinMessageHandler * handler = model_->messageHandler()->clone();
00159 si->passInMessageHandler(handler);
00160 si->messageHandler()->setLogLevel(2);
00161 #ifdef DEBUG_BON_HEURISTIC
00162 std::cout << "Loading problem into si\n";
00163 #endif
00164 extractInnerApproximation(*nlp, *si, newSolution, true);
00165 #ifdef DEBUG_BON_HEURISTIC
00166 std::cout << "problem loaded\n";
00167 std::cout << "**** Running optimization ****\n";
00168 #endif
00169 mip_->optimize(DBL_MAX, 2, time_limit_);
00170 #ifdef DEBUG_BON_HEURISTIC
00171 std::cout << "Optimization finished\n";
00172 #endif
00173 if(mip_->getLastSolution()) {
00174 const double* solution = mip_->getLastSolution();
00175 std::copy(solution, solution + numberColumns, newSolution);
00176 }
00177 else
00178 feasible = false;
00179
00180 if(delete_si) {
00181 delete si;
00182 }
00183 delete handler;
00184
00185
00186 #if 0 // Set to 1 if you need to test the feasibility of the returned solution
00187 const double* x_l = minlp->x_l();
00188 const double* x_u = minlp->x_u();
00189 const double* g_l = minlp->g_l();
00190 const double* g_u = minlp->g_u();
00191 double primalTolerance = 1.0e-6;
00192
00193 Bonmin::vector<Ipopt::TNLP::LinearityType> constTypes(numberRows);
00194 minlp->get_constraints_linearity(numberRows, constTypes());
00195 feasible = true;
00196 for (int iColumn=0;iColumn<numberColumns;iColumn++) {
00197 double value=newSolution[iColumn];
00198 if(value - x_l[iColumn] < -1e-8|| value - x_u[iColumn] > 1e-8) {
00199 std::cout<<"Solution found infeasible because: "<<std::endl;
00200 std::cout<<"x_l["<<iColumn<<"]= "<<x_l[iColumn]<<" "
00201 <<"x_sol["<<iColumn<<"]= "<<value<<" "
00202 <<"x_u["<<iColumn<<"]= "<<x_u[iColumn]<<std::endl;
00203 feasible = false;
00204 break;
00205 }
00206 }
00207 minlp->eval_g(numberColumns, newSolution, true,
00208 numberRows, new_g_sol);
00209 for(int iRow=0; iRow<numberRows; iRow++) {
00210 if(new_g_sol[iRow]<g_l[iRow]-primalTolerance ||
00211 new_g_sol[iRow]>g_u[iRow]+primalTolerance) {
00212 std::cout<<"It should be infeasible because: "<<std::endl;
00213 std::cout<<"g_l["<<iRow<<"]= "<<g_l[iRow]<<" "
00214 <<"g_sol["<<iRow<<"]= "<<new_g_sol[iRow]<<" "
00215 <<"g_u["<<iRow<<"]= "<<g_u[iRow]<<std::endl;
00216 std::cout<<"primalTolerance= "<<primalTolerance<<std::endl;
00217 if(constTypes[iRow] == Ipopt::TNLP::NON_LINEAR)
00218 std::cout<<"nonLinear constraint number "<<iRow<<std::endl;
00219 feasible = false;
00220 }
00221 }
00222 std::cout<<"Every thing is feasible"<<std::endl;
00223 #endif
00224
00225 if(feasible) {
00226 double newSolutionValue;
00227 minlp->eval_f(numberColumns, newSolution, true, newSolutionValue);
00228 if(newSolutionValue < solutionValue) {
00229 memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
00230 solutionValue = newSolutionValue;
00231 returnCode = 1;
00232 }
00233 }
00234
00235 delete [] newSolution;
00236 delete [] new_g_sol;
00237
00238 delete nlp;
00239
00240 #ifdef DEBUG_BON_HEURISTIC
00241 std::cout<<"Inner approximation returnCode = "<<returnCode<<std::endl;
00242 #endif
00243 return returnCode;
00244 }
00245
00248 bool
00249 HeuristicInnerApproximation::getMyInnerApproximation(Bonmin::OsiTMINLPInterface &si, OsiCuts &cs, int ind,
00250 const double * x, const double * x2) {
00251
00252 int n, m, nnz_jac_g, nnz_h_lag;
00253 Ipopt::TNLP::IndexStyleEnum index_style;
00254 Bonmin::TMINLP2TNLP * problem = si.problem();
00255 problem->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
00256 double infty = si.getInfinity();
00257
00258
00259 CoinPackedVector cut;
00260 double lb = -infty;
00261 double ub = 0;
00262
00263
00264 double g = 0;
00265 double g2 = 0;
00266 double diff = 0;
00267 double a = 0;
00268 problem->eval_gi(n, x, 1, ind, g);
00269 problem->eval_gi(n, x2, 1, ind, g2);
00270 Bonmin::vector<int> jCol(n);
00271 int nnz;
00272 problem->eval_grad_gi(n, x2, 0, ind, nnz, jCol(), NULL);
00273 Bonmin::vector<double> jValues(nnz);
00274 problem->eval_grad_gi(n, x2, 0, ind, nnz, NULL, jValues());
00275 bool add = false;
00276
00277 for (int i = 0; i < nnz; i++) {
00278 const int &colIdx = jCol[i];
00279 if(index_style == Ipopt::TNLP::FORTRAN_STYLE) jCol[i]--;
00280 diff = x[colIdx] - x2[colIdx];
00281
00282 if (fabs(diff) >= 1e-8) {
00283 a = (g - g2) / diff;
00284 cut.insert(colIdx, a);
00285 ub = (a * x[colIdx] - g) - fabs(a * x[colIdx] - g)*1e-6;
00286
00287 add = true;
00288 } else {
00289 cut.insert(colIdx, jValues[i]);
00290
00291 }
00292 }
00293
00294 if (add) {
00295
00296 OsiRowCut newCut;
00297 newCut.setGloballyValidAsInteger(1);
00298 newCut.setLb(lb);
00299
00300
00301 int binary_id = 0;
00302 const int* ids = problem->get_const_xtra_id();
00303
00304 binary_id = (ids == NULL) ? -1 : ids[ind];
00305 if(binary_id>0) {
00306 cut.insert(binary_id, -ub);
00307 newCut.setUb(0);
00308 }
00309 else
00310 newCut.setUb(ub);
00311
00312
00313 newCut.setRow(cut);
00314 cs.insert(newCut);
00315
00316 return true;
00317 }
00318 return false;
00319 }
00320
00321 void
00322 HeuristicInnerApproximation::extractInnerApproximation(Bonmin::OsiTMINLPInterface & nlp, OsiSolverInterface &si,
00323 const double * x, bool getObj) {
00324 printf("************ Start extracting inner approx");
00325 int n;
00326 int m;
00327 int nnz_jac_g;
00328 int nnz_h_lag;
00329 Ipopt::TNLP::IndexStyleEnum index_style;
00330 Bonmin::TMINLP2TNLP * problem = nlp.problem();
00331
00332 problem->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
00333
00334 Bonmin::vector<int> jRow(nnz_jac_g);
00335 Bonmin::vector<int> jCol(nnz_jac_g);
00336 Bonmin::vector<double> jValues(nnz_jac_g);
00337 problem->eval_jac_g(n, NULL, 0, m, nnz_jac_g, jRow(), jCol(), NULL);
00338 if(index_style == Ipopt::TNLP::FORTRAN_STYLE)
00339 {
00340 for(int i = 0 ; i < nnz_jac_g ; i++){
00341 jRow[i]--;
00342 jCol[i]--;
00343 }
00344 }
00345
00346
00347 problem->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL,
00348 jValues());
00349
00350 Bonmin::vector<double> g(m);
00351 problem->eval_g(n, x, 1, m, g());
00352
00353 Bonmin::vector<int> nonLinear(m);
00354
00355 int numNonLinear = 0;
00356 const double * rowLower = nlp.getRowLower();
00357 const double * rowUpper = nlp.getRowUpper();
00358 const double * colLower = nlp.getColLower();
00359 const double * colUpper = nlp.getColUpper();
00360 assert(m == nlp.getNumRows());
00361 double infty = si.getInfinity();
00362 double nlp_infty = nlp.getInfinity();
00363 Bonmin::vector<Ipopt::TNLP::LinearityType> constTypes(m);
00364 Bonmin::vector<Ipopt::TNLP::LinearityType> varTypes(n);
00365 problem->get_constraints_linearity(m, constTypes());
00366 problem->get_variables_linearity(n, varTypes());
00367 for (int i = 0; i < m; i++) {
00368 if (constTypes[i] == Ipopt::TNLP::NON_LINEAR) {
00369 nonLinear[numNonLinear++] = i;
00370 }
00371 }
00372 Bonmin::vector<double> rowLow(m - numNonLinear);
00373 Bonmin::vector<double> rowUp(m - numNonLinear);
00374 int ind = 0;
00375 for (int i = 0; i < m; i++) {
00376 if (constTypes[i] != Ipopt::TNLP::NON_LINEAR) {
00377 if (rowLower[i] > -nlp_infty) {
00378
00379 rowLow[ind] = (rowLower[i]);
00380 } else
00381 rowLow[ind] = -infty;
00382 if (rowUpper[i] < nlp_infty) {
00383
00384 rowUp[ind] = (rowUpper[i]);
00385 } else
00386 rowUp[ind] = infty;
00387 ind++;
00388 }
00389
00390 }
00391
00392 CoinPackedMatrix mat(true, jRow(), jCol(), jValues(), nnz_jac_g);
00393 mat.setDimensions(m, n);
00394
00395
00396 mat.deleteRows(numNonLinear, nonLinear());
00397
00398 int numcols = nlp.getNumCols();
00399 Bonmin::vector<double> obj(numcols);
00400 for (int i = 0; i < numcols; i++)
00401 obj[i] = 0.;
00402
00403 si.loadProblem(mat, nlp.getColLower(), nlp.getColUpper(),
00404 obj(), rowLow(), rowUp());
00405 const Bonmin::TMINLP::VariableType* variableType = problem->var_types();
00406 for (int i = 0; i < n; i++) {
00407 if ((variableType[i] == Bonmin::TMINLP::BINARY) || (variableType[i] == Bonmin::TMINLP::INTEGER))
00408 si.setInteger(i);
00409 }
00410 if (getObj) {
00411 bool addObjVar = false;
00412 if (problem->hasLinearObjective()) {
00413 double zero;
00414 Bonmin::vector<double> x0(n, 0.);
00415 problem->eval_f(n, x0(), 1, zero);
00416 si.setDblParam(OsiObjOffset, -zero);
00417
00418 problem->eval_grad_f(n, x, 1, obj());
00419 si.setObjective(obj());
00420 } else {
00421 addObjVar = true;
00422 }
00423
00424 if (addObjVar) {
00425 nlp.addObjectiveFunction(si, x);
00426 }
00427 }
00428
00429
00430 int InnerDesc = 1;
00431 if (InnerDesc == 1) {
00432 OsiCuts cs;
00433
00434 double * p = CoinCopyOfArray(colLower, n);
00435 double * pp = CoinCopyOfArray(colLower, n);
00436 double * up = CoinCopyOfArray(colUpper, n);
00437
00438 for (int i = 0; i < n; i++) {
00439 if (p[i] < -1e3){
00440 p[i] = pp[i] = -1e3;
00441 }
00442 if (up[i] > 1e2){
00443 up[i] = 1e2;
00444 }
00445 }
00446
00447 const int& nbAp = nbAp_;
00448 printf("Generating approximation with %i points.\n", nbAp);
00449
00450 std::vector<double> step(n);
00451 int n_lin = 0;
00452
00453 for (int i = 0; i < n; i++) {
00454
00455 if (varTypes[i] == Ipopt::TNLP::LINEAR) {
00456 n_lin ++;
00457 step[i] = 0;
00458 p[i] = pp[i] = up[i] = 0;
00459 }
00460 else {
00461 step[i] = (up[i] - p[i]) / (nbAp);
00462 }
00463 }
00464 printf("Number of linears %i\n", n_lin);
00465
00466 for (int j = 1; j < nbAp; j++) {
00467
00468 for (int i = 0; i < n; i++) {
00469 pp[i] += step[i];
00470 }
00471
00472 for (int i = 0; (i < m ); i++) {
00473 if (constTypes[i] == Ipopt::TNLP::LINEAR) continue;
00474 bool status = getMyInnerApproximation(nlp, cs, i, p, pp);
00475 if(status == false){
00476 printf("Error in generating inner approximation\n");
00477 exit(1);
00478 }
00479 }
00480 std::copy(pp, pp+n, p);
00481
00482 }
00483
00484 for(int i = 0; (i< m); i++) {
00485 if (constTypes[i] == Ipopt::TNLP::LINEAR) continue;
00486 getMyInnerApproximation(nlp, cs, i, p, up);
00487 }
00488
00489 delete [] p;
00490 delete [] pp;
00491 delete [] up;
00492 si.applyCuts(cs);
00493 }
00494 printf("************ Done extracting inner approx ********");
00495 }
00496
00497 }
00498