00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "BonCutStrengthener.hpp"
00010 #include "IpBlas.hpp"
00011
00012 namespace Bonmin
00013 {
00014 using namespace Ipopt;
00015
00016 CutStrengthener::CutStrengthener(SmartPtr<TNLPSolver> tnlp_solver,
00017 SmartPtr<OptionsList> options)
00018 :
00019 tnlp_solver_(tnlp_solver)
00020 {
00021 options->GetIntegerValue("oa_log_level", oa_log_level_, "bonmin.");
00022 options->GetEnumValue("cut_strengthening_type", cut_strengthening_type_,
00023 "bonmin.");
00024 options->GetEnumValue("disjunctive_cut_type", disjunctive_cut_type_,
00025 "bonmin.");
00026
00027 tnlp_solver_->options()->clear();
00028 if (!tnlp_solver_->Initialize("strength.opt")) {
00029 throw CoinError("CutStrengthener","CutStrengthener","Error during initialization of tnlp_solver_");
00030 }
00031 tnlp_solver_->options()->SetStringValue("hessian_approximation","limited-memory");
00032 tnlp_solver_->options()->SetStringValue("mu_strategy", "adaptive");
00033 }
00034
00035 CutStrengthener::~CutStrengthener()
00036 {
00037 }
00038
00039 bool CutStrengthener::HandleOneCut(bool is_tight, TMINLP* tminlp,
00040 TMINLP2TNLP* problem,
00041 const double* minlp_lb,
00042 const double* minlp_ub,
00043 const int gindex, CoinPackedVector& cut,
00044 double& cut_lb, double& cut_ub,
00045 int n, const double* x,
00046 double infty)
00047 {
00048 bool retval = true;
00049 const int cut_nele = cut.getNumElements();
00050 const int* cut_indices = cut.getIndices();
00051 const TMINLP::VariableType* vartypes = problem->var_types();
00052 const double* cut_elements = cut.getElements();
00053 switch (disjunctive_cut_type_) {
00054 case DC_None:
00055 if (!is_tight) {
00056 retval = StrengthenCut(tminlp, gindex, cut, n, x, minlp_lb, minlp_ub,
00057 cut_lb, cut_ub);
00058 }
00059 break;
00060 case DC_MostFractional: {
00061
00062 int imostfra = -1;
00063 double viol = 1e-6;
00064 for (int i=0; i<cut_nele; i++) {
00065 const int& idx = cut_indices[i];
00066 if (idx < n && (vartypes[idx] == TMINLP::BINARY ||
00067 vartypes[idx] == TMINLP::INTEGER)) {
00068 const double& xi = x[idx];
00069 const double this_viol = CoinMin(xi-floor(xi),ceil(xi)-xi);
00070 if (this_viol > viol) {
00071 imostfra = i;
00072 viol = this_viol;
00073 }
00074 }
00075 }
00076 if (imostfra == -1) {
00077
00078 if (!is_tight) {
00079 retval = StrengthenCut(tminlp, gindex, cut, n, x, minlp_lb, minlp_ub,
00080 cut_lb, cut_ub);
00081 }
00082 }
00083 else {
00084
00085 const int& idx = cut_indices[imostfra];
00086 const double& xi = x[idx];
00087 if (oa_log_level_>=2) {
00088 printf("Doing disjunction for constr %d on x[%d] = %e\n", gindex, idx, xi);
00089 }
00090 const double down_xi = floor(xi);
00091 double* changed_bnds = new double[n];
00092
00093 CoinCopyN(minlp_ub, n, changed_bnds);
00094 changed_bnds[idx] = down_xi;
00095 double cut_lb_down = cut_lb;
00096 double cut_ub_down = cut_ub;
00097 retval = StrengthenCut(tminlp, gindex, cut, n, x, minlp_lb,
00098 changed_bnds, cut_lb_down, cut_ub_down);
00099 double cut_lb_up = cut_lb;
00100 double cut_ub_up = cut_ub;
00101 if (retval) {
00102 CoinCopyN(minlp_lb, n, changed_bnds);
00103 changed_bnds[idx] = down_xi + 1.;
00104 retval = StrengthenCut(tminlp, gindex, cut, n, x, changed_bnds,
00105 minlp_ub, cut_lb_up, cut_ub_up);
00106 }
00107 delete [] changed_bnds;
00108 if (retval) {
00109
00110 double coeff_xi = cut_elements[imostfra];
00111 const double old_coeff = coeff_xi;
00112 if (cut_lb <= -infty) {
00113
00114 coeff_xi += cut_ub_down - cut_ub_up;
00115 cut_ub = cut_ub_up + (down_xi+1.)*(cut_ub_down - cut_ub_up);
00116 }
00117 else {
00118
00119 coeff_xi += cut_lb_down - cut_lb_up;
00120 cut_lb = cut_lb_up + (down_xi+1.)*(cut_lb_down - cut_lb_up);
00121 }
00122 cut.setElement(imostfra, coeff_xi);
00123 printf("old coeff = %e new = %e\n", old_coeff, coeff_xi);
00124 }
00125 }
00126 }
00127 break;
00128 default:
00129 std::cerr << "Invalid case for disjunctive_cut_type_ in CutStrengthener HandleOneCut\n";
00130 exit(-2);
00131 break;
00132 }
00133
00134 return retval;
00135 }
00136
00137 bool CutStrengthener::ComputeCuts(OsiCuts &cs,
00138 TMINLP* tminlp,
00139 TMINLP2TNLP* problem,
00140 const int gindex, CoinPackedVector& cut,
00141 double& cut_lb, double& cut_ub,
00142 const double g_val, const double g_lb,
00143 const double g_ub,
00144 int n, const double* x,
00145 double infty)
00146 {
00147
00148
00149 bool is_tight = false;
00150 if (gindex==-1) {
00151 is_tight = true;
00152 }
00153 else {
00154 const Number tight_tol = 1e-8;
00155 if (cut_lb <= -infty && g_ub - g_val <= tight_tol) {
00156 is_tight = true;
00157 }
00158 else if (cut_ub >= infty && g_val - g_lb <= tight_tol) {
00159 is_tight = true;
00160 }
00161 }
00162 if (cut_strengthening_type_ == CS_StrengthenedGlobal ||
00163 cut_strengthening_type_ == CS_StrengthenedGlobal_StrengthenedLocal) {
00164 const double orig_lb = cut_lb;
00165 const double orig_ub = cut_ub;
00166 bool retval = HandleOneCut(is_tight, tminlp, problem,
00167 problem->orig_x_l(),
00168 problem->orig_x_u(), gindex, cut,
00169 cut_lb, cut_ub, n, x, infty);
00170 if (!retval) {
00171 if (oa_log_level_ >= 1) {
00172 printf(" Error during strengthening of global cut for constraint %d\n", gindex);
00173 }
00174 }
00175 else {
00176 if (oa_log_level_ >=2 && (fabs(orig_lb-cut_lb)>1e-4 ||
00177 fabs(orig_ub-cut_ub)>1e-4)) {
00178 if (orig_ub < infty) {
00179 printf(" Strengthening ub of global cut for constraint %d from %e to %e\n", gindex, orig_ub, cut_ub);
00180 }
00181 else {
00182 printf(" Strengthening lb of global cut for constraint %d from %e to %e\n", gindex, orig_lb, cut_lb);
00183 }
00184 }
00185 }
00186 }
00187 if (cut_strengthening_type_ == CS_UnstrengthenedGlobal_StrengthenedLocal ||
00188 cut_strengthening_type_ == CS_StrengthenedGlobal_StrengthenedLocal) {
00189 Number lb2 = cut_lb;
00190 Number ub2 = cut_ub;
00191 CoinPackedVector cut2(cut);
00192 bool retval = HandleOneCut(is_tight, tminlp, problem, problem->x_l(),
00193 problem->x_u(), gindex, cut2,
00194 lb2, ub2, n, x, infty);
00195 if (!retval) {
00196 if (oa_log_level_ >= 1) {
00197 printf(" Error during strengthening of local cut for constraint %d\n", gindex);
00198 }
00199 }
00200 else {
00201 const Number localCutTol = 1e-4;
00202 if (fabs(lb2-cut_lb) >= localCutTol || fabs(cut_ub-ub2) >= localCutTol) {
00203 if (ub2 < infty) {
00204 printf(" Strengthening ub of local cut for constraint %d from %e to %e\n", gindex, cut_ub, ub2);
00205 }
00206 else {
00207 printf(" Strengthening ub of local cut for constraint %d from %e to %e\n", gindex, cut_lb, lb2);
00208 }
00209
00210 OsiRowCut newCut2;
00211 newCut2.setEffectiveness(99.99e99);
00212 newCut2.setLb(lb2);
00213 newCut2.setUb(ub2);
00214 newCut2.setRow(cut2);
00215 cs.insert(newCut2);
00216 }
00217 }
00218 }
00219 return true;
00220 }
00221
00222 bool CutStrengthener::StrengthenCut(SmartPtr<TMINLP> tminlp,
00223 int constr_index,
00224 const CoinPackedVector& row,
00225 int n,
00226 const double* x,
00227 const double* x_l,
00228 const double* x_u,
00229 double& lb,
00230 double& ub)
00231 {
00232
00233 Index nele_grad_gi;
00234 Index* jCol = new Index[n+1];
00235 bool new_x = true;
00236 if (constr_index == -1) {
00237
00238
00239 double* x_rand = new double[n];
00240 for (int i=0; i<n; i++) {
00241 const double radius = CoinMin(1., x_u[i]-x_l[i]);
00242 const double p = CoinMax(CoinMin(x[i]-0.5*radius, x_u[i]-radius),
00243 x_l[i]);
00244 x_rand[i] = p + radius*CoinDrand48();
00245 }
00246 Number* grad_f = new Number[n];
00247 bool retval = tminlp->eval_grad_f(n, x_rand, new_x, grad_f);
00248 delete [] x_rand;
00249 if (!retval) {
00250 delete [] grad_f;
00251 delete [] jCol;
00252 return false;
00253 }
00254 nele_grad_gi = 0;
00255 for (int i=0; i<n; i++) {
00256 if (grad_f[i] != 0.) {
00257 jCol[nele_grad_gi++] = i;
00258 }
00259 }
00260 delete [] grad_f;
00261 jCol[nele_grad_gi++] = n;
00262 }
00263 else {
00264 if (!tminlp->eval_grad_gi(n, x, new_x, constr_index, nele_grad_gi,
00265 jCol, NULL)) {
00266 delete [] jCol;
00267 return false;
00268 }
00269 }
00270
00271 bool lower_bound;
00272 if (lb <= -COIN_DBL_MAX) {
00273 assert(ub < COIN_DBL_MAX);
00274 lower_bound = false;
00275 }
00276 else {
00277 assert(ub >= COIN_DBL_MAX);
00278 lower_bound = true;
00279 }
00280 SmartPtr<StrengtheningTNLP> stnlp =
00281 new StrengtheningTNLP(tminlp, row, lower_bound, n, x, x_l, x_u,
00282 constr_index, nele_grad_gi, jCol);
00283
00284 delete [] jCol;
00285
00286 TNLPSolver::ReturnStatus status =
00287 tnlp_solver_->OptimizeTNLP(GetRawPtr(stnlp));
00288
00289 if (status == TNLPSolver::solvedOptimal ||
00290 status == TNLPSolver::solvedOptimalTol) {
00291 const Number tiny_move = 0e-8;
00292 const Number final_bound = stnlp->StrengthenedBound();
00293 if (lower_bound) {
00294 lb = final_bound - tiny_move;
00295 }
00296 else {
00297 ub = final_bound + tiny_move;
00298 }
00299 }
00300 else {
00301 return false;
00302 }
00303
00304 return true;
00305 }
00306
00307 CutStrengthener::StrengtheningTNLP::
00308 StrengtheningTNLP(SmartPtr<TMINLP> tminlp,
00309 const CoinPackedVector& cut,
00310 bool lower_bound,
00311 Index n,
00312 const Number* starting_point,
00313 const double* x_l_orig,
00314 const double* x_u_orig,
00315 Index constr_index,
00316 Index nvar_constr ,
00317 const Index* jCol)
00318 :
00319 tminlp_(tminlp),
00320 n_orig_(n),
00321 constr_index_(constr_index),
00322 nvar_constr_(nvar_constr),
00323 lower_bound_(lower_bound),
00324 have_final_bound_(false),
00325 grad_f_(NULL)
00326 {
00327 starting_point_ = new Number[n_orig_];
00328 x_full_ = new Number[n_orig_];
00329 IpBlasDcopy(n_orig_, starting_point, 1, starting_point_, 1);
00330 IpBlasDcopy(n_orig_, starting_point, 1, x_full_, 1);
00331
00332 obj_grad_ = new Number[nvar_constr_];
00333 x_l_ = new Number[nvar_constr_];
00334 x_u_ = new Number[nvar_constr_];
00335 const Number zero = 0.;
00336 IpBlasDcopy(nvar_constr_, &zero, 0, obj_grad_, 1);
00337
00338 const int cut_nele = cut.getNumElements();
00339 const int* cut_indices = cut.getIndices();
00340 const double* cut_elements = cut.getElements();
00341
00342 for (int i = 0; i<cut_nele; i++) {
00343 const int& idx = cut_indices[i];
00344
00345 Index jidx = -1;
00346 for (int j=0; j<nvar_constr_; j++) {
00347 if (idx == jCol[j]) {
00348 jidx = j;
00349 break;
00350 }
00351 }
00352 if (jidx < 0) {
00353 printf("There is an index (%d) in the cut that does not appear in the constraint.\n",idx);
00354 exit(-99);
00355 }
00356
00357 if (lower_bound) {
00358 obj_grad_[jidx] = cut_elements[i];
00359 }
00360 else {
00361 obj_grad_[jidx] = -cut_elements[i];
00362 }
00363
00364 }
00365
00366 var_indices_ = new Index[nvar_constr_];
00367 for (int i=0; i<nvar_constr_; i++) {
00368 const Index& j = jCol[i];
00369 var_indices_[i] = j;
00370 if (j < n) {
00371 x_l_[i] = x_l_orig[j];
00372 x_u_[i] = x_u_orig[j];
00373 }
00374 else {
00375 x_l_[i] = -1e100;
00376 x_u_[i] = 1e100;
00377 }
00378 }
00379
00380 if (constr_index_ == -1) {
00381 grad_f_ = new Number[n_orig_];
00382 }
00383 }
00384
00385 CutStrengthener::StrengtheningTNLP::~StrengtheningTNLP()
00386 {
00387 delete [] obj_grad_;
00388 delete [] x_l_;
00389 delete [] x_u_;
00390 delete [] var_indices_;
00391 delete [] starting_point_;
00392 delete [] x_full_;
00393 delete [] grad_f_;
00394 }
00395
00396 bool CutStrengthener::StrengtheningTNLP::
00397 get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
00398 Index& nnz_h_lag, IndexStyleEnum& index_style)
00399 {
00400 n = nvar_constr_;
00401 m = 1;
00402 nnz_jac_g = nvar_constr_;
00403 nnz_h_lag = 0;
00404 index_style = C_STYLE;
00405
00406 Index n_orig;
00407 Index nnz_jac_g_orig;
00408 Index nnz_h_lag_orig;
00409 TNLP::IndexStyleEnum index_style_orig;
00410 if(!tminlp_->get_nlp_info(n_orig, m_orig_, nnz_jac_g_orig, nnz_h_lag_orig,
00411 index_style_orig)) {
00412 return false;
00413 }
00414 if (n_orig_ != n_orig) {
00415 std::cerr << "Number of variables inconsistent in StrengtheningTNLP::get_nlp_info\n";
00416 return false;
00417 }
00418
00419 return true;
00420 }
00421
00422 bool CutStrengthener::StrengtheningTNLP::
00423 get_bounds_info(Index n, Number* x_l, Number* x_u,
00424 Index m, Number* g_l, Number* g_u)
00425 {
00426 if (constr_index_ == -1) {
00427 g_l[0] = -1e100;
00428 g_u[0] = 0.;
00429 }
00430 else {
00431 Number* x_l_orig = new Number[n_orig_];
00432 Number* x_u_orig = new Number[n_orig_];
00433 Number* g_l_orig = new Number[m_orig_];
00434 Number* g_u_orig = new Number[m_orig_];
00435
00436 if (!tminlp_->get_bounds_info(n_orig_, x_l_orig, x_u_orig,
00437 m_orig_, g_l_orig, g_u_orig)) {
00438 delete [] x_l_orig;
00439 delete [] x_u_orig;
00440 delete [] g_l_orig;
00441 delete [] g_u_orig;
00442 return false;
00443 }
00444
00445 g_l[0] = g_l_orig[constr_index_];
00446 g_u[0] = g_u_orig[constr_index_];
00447
00448 delete [] x_l_orig;
00449 delete [] x_u_orig;
00450 delete [] g_l_orig;
00451 delete [] g_u_orig;
00452 }
00453
00454 for (Index i=0; i<nvar_constr_; i++) {
00455 x_l[i] = x_l_[i];
00456 x_u[i] = x_u_[i];
00457 }
00458
00459 return true;
00460 }
00461
00462 bool CutStrengthener::StrengtheningTNLP::
00463 get_starting_point(Index n, bool init_x, Number* x,
00464 bool init_z, Number* z_L, Number* z_U,
00465 Index m, bool init_lambda,
00466 Number* lambda)
00467 {
00468 assert(!init_z && !init_lambda);
00469 assert(n = nvar_constr_);
00470 if (init_x) {
00471 if (constr_index_ == -1) {
00472 for (Index i=0; i<n-1; i++) {
00473 x[i] = starting_point_[var_indices_[i]];
00474 }
00475 x[n-1] = 0.;
00476 }
00477 else {
00478 for (Index i=0; i<n; i++) {
00479 x[i] = starting_point_[var_indices_[i]];
00480 }
00481 }
00482 }
00483
00484 return true;
00485 }
00486
00487 bool CutStrengthener::StrengtheningTNLP::
00488 eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
00489 {
00490 obj_value = 0.;
00491 for (Index i=0; i<n; i++) {
00492 obj_value += obj_grad_[i]*x[i];
00493 }
00494 return true;
00495 }
00496
00497 bool CutStrengthener::StrengtheningTNLP::
00498 eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
00499 {
00500 IpBlasDcopy(n, obj_grad_, 1, grad_f, 1);
00501 return true;
00502 }
00503
00504 bool CutStrengthener::StrengtheningTNLP::
00505 eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
00506 {
00507 update_x_full(x);
00508 bool retval;
00509 if (constr_index_ == -1) {
00510 retval = tminlp_->eval_f(n_orig_, x_full_, new_x, g[0]);
00511 g[0] -= x[n-1];
00512 }
00513 else {
00514 retval = tminlp_->eval_gi(n_orig_, x_full_, new_x, constr_index_, g[0]);
00515 }
00516 return retval;
00517 }
00518
00519 bool CutStrengthener::StrengtheningTNLP::
00520 eval_jac_g(Index n, const Number* x, bool new_x,
00521 Index m, Index nele_jac, Index* iRow, Index *jCol,
00522 Number* values)
00523 {
00524 bool retval = true;
00525 if (iRow) {
00526 DBG_ASSERT(jCol && !values);
00527 for (Index i=0; i<nele_jac; i++) {
00528 iRow[i] = 0;
00529 jCol[i] = i;
00530 }
00531 }
00532 else {
00533 DBG_ASSERT(!iRow && values);
00534 update_x_full(x);
00535 if (constr_index_ == -1) {
00536 retval = tminlp_->eval_grad_f(n_orig_, x_full_, new_x, grad_f_);
00537 if (retval) {
00538 for (Index i=0; i<n-1; i++) {
00539 values[i] = grad_f_[var_indices_[i]];
00540 }
00541 values[n-1] = -1.;
00542 }
00543 }
00544 else {
00545 retval = tminlp_->eval_grad_gi(n_orig_, x_full_, new_x, constr_index_,
00546 nele_jac, NULL, values);
00547 }
00548 }
00549 return retval;
00550 }
00551
00552 bool CutStrengthener::StrengtheningTNLP::
00553 eval_h(Index n, const Number* x, bool new_x,
00554 Number obj_factor, Index m, const Number* lambda,
00555 bool new_lambda, Index nele_hess,
00556 Index* iRow, Index* jCol, Number* values)
00557 {
00558 std::cerr << "At the moment, StrengtheningTNLP::eval_h is not yet implemented\n";
00559 return false;
00560 }
00561
00562 void CutStrengthener::StrengtheningTNLP::
00563 update_x_full(const Number *x)
00564 {
00565 if (constr_index_ == -1) {
00566 for (Index i=0; i<nvar_constr_-1; i++) {
00567 x_full_[var_indices_[i]] = x[i];
00568 }
00569 }
00570 else {
00571 for (Index i=0; i<nvar_constr_; i++) {
00572 x_full_[var_indices_[i]] = x[i];
00573 }
00574 }
00575 }
00576
00577 void CutStrengthener::StrengtheningTNLP::
00578 finalize_solution(SolverReturn status,
00579 Index n, const Number* x, const Number* z_L, const Number* z_U,
00580 Index m, const Number* g, const Number* lambda,
00581 Number obj_value,
00582 const IpoptData* ip_data,
00583 IpoptCalculatedQuantities* ip_cq)
00584 {
00585 if (status == SUCCESS || status == STOP_AT_ACCEPTABLE_POINT) {
00586 have_final_bound_ = true;
00587 strengthened_bound_ = obj_value;
00588 }
00589 else {
00590 have_final_bound_ = false;
00591 }
00592 }
00593
00594 Number CutStrengthener::StrengtheningTNLP::StrengthenedBound() const
00595 {
00596 DBG_ASSERT(have_final_bound_);
00597 if (lower_bound_) {
00598 return strengthened_bound_;
00599 }
00600 else {
00601 return -strengthened_bound_;
00602 }
00603 }
00604
00605
00606 }