00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <vector>
00012
00013 #include "BonRegisteredOptions.hpp"
00014
00015 #include "CoinHelperFunctions.hpp"
00016 #include "CoinTime.hpp"
00017
00018 #include "BonBabSetupBase.hpp"
00019
00020 #include "CouenneTypes.hpp"
00021
00022 #include "CouenneExpression.hpp"
00023 #include "CouenneExprConst.hpp"
00024 #include "CouenneExprGroup.hpp"
00025 #include "CouenneExprClone.hpp"
00026 #include "CouenneExprAux.hpp"
00027 #include "CouenneProblem.hpp"
00028 #include "CouenneGlobalCutOff.hpp"
00029 #include "CouenneProblemElem.hpp"
00030 #include "CouenneLQelems.hpp"
00031
00032 using namespace Couenne;
00033
00034 #define MAX_FBBT_ITER 3
00035
00037 void CouenneProblem::setBase (Bonmin::BabSetupBase *base) {
00038 bonBase_ = base;
00039 jnlst_ = base -> journalist ();
00040 }
00041
00042
00043 const CouNumber SafeCutoff = 1e-4;
00044
00045
00046 const CouNumber SafeDelta = 1e-2;
00047
00050
00051 void CouenneProblem::initAuxs () const {
00052
00053 domain_.current () -> resize (nVars ());
00054
00055
00056
00057
00058 int nvars = nVars ();
00059
00060 for (int i=0; i < nvars; i++) {
00061
00062 int indvar = variables_ [i] -> Index ();
00063
00064 if (((variables_ [i] -> Type () == AUX) &&
00065 (indvar >= nOrigVars_)) ||
00066 (variables_ [i] -> Multiplicity () == 0))
00067
00068 Lb (indvar) = - (Ub (indvar) = COIN_DBL_MAX);
00069 }
00070
00071
00072
00073
00074
00075 for (std::vector <CouenneConstraint *>::const_iterator con = constraints_.begin ();
00076 con != constraints_.end (); ++con) {
00077
00078 CouNumber
00079 lb = (*((*con) -> Lb ())) (),
00080 ub = (*((*con) -> Ub ())) ();
00081
00082 int index = (*con) -> Body () -> Index ();
00083
00084 assert (index >= 0);
00085
00086 if ((Lb (index) = CoinMax (Lb (index), lb)) <= -COUENNE_INFINITY) Lb (index) = -COIN_DBL_MAX;
00087 if ((Ub (index) = CoinMin (Ub (index), ub)) >= COUENNE_INFINITY) Ub (index) = COIN_DBL_MAX;
00088 }
00089
00090
00091
00092
00093 Jnlst () -> Printf (Ipopt::J_MOREMATRIX, J_PROBLEM, "InitAux -- assigning bounds\n");
00094
00095 for (int j=0, i=nVars (); i--; j++) {
00096
00097 int ord = numbering_ [j];
00098
00099
00100 if (variables_ [ord] -> Multiplicity () == 0) {
00101 Lb (ord) = - (Ub (ord) = COIN_DBL_MAX);
00102 X (ord) = 0.;
00103 continue;
00104 }
00105
00106 exprVar *var = variables_ [ord];
00107
00108
00109 if (var -> Type () == AUX) {
00110
00111 Jnlst () -> Printf (Ipopt::J_MOREMATRIX, J_PROBLEM,
00112 "w_%04d [%10g,%10g] ", ord, Lb (ord), Ub (ord));
00113
00114 CouNumber l, u;
00115
00116 var -> Image () -> getBounds (l, u);
00117
00118
00119
00120
00121
00122 Jnlst () -> Printf (Ipopt::J_MOREMATRIX, J_PROBLEM,
00123 " ( --> w_%04d [%10g,%10g] ) vs [%10g %10g]",
00124 ord, l, u, Lb (ord), Ub (ord));
00125
00126
00127 if ((var -> sign () != expression::AUX_LEQ) && ((Lb (ord) = CoinMax (Lb (ord), l)) <= -COUENNE_INFINITY)) Lb (ord) = -COIN_DBL_MAX;
00128 if ((var -> sign () != expression::AUX_GEQ) && ((Ub (ord) = CoinMin (Ub (ord), u)) >= COUENNE_INFINITY)) Ub (ord) = COIN_DBL_MAX;
00129
00130
00131
00132
00133 bool integer = var -> isDefinedInteger ();
00134
00135 if (integer) {
00136 if (var -> sign () != expression::AUX_GEQ) Lb (ord) = ceil (Lb (ord) - COUENNE_EPS);
00137 if (var -> sign () != expression::AUX_LEQ) Ub (ord) = floor (Ub (ord) + COUENNE_EPS);
00138 }
00139
00140 X (ord) = CoinMax (Lb (ord), CoinMin (Ub (ord), (*(var -> Image ())) ()));
00141
00142 Jnlst () -> Printf (Ipopt::J_MOREMATRIX, J_PROBLEM,
00143 " --> [%10g,%10g] (%g)\n", Lb (ord), Ub (ord), X (ord));
00144 }
00145 }
00146
00147 restoreUnusedOriginals ();
00148 }
00149
00150
00153 void CouenneProblem::getAuxs (CouNumber * x) const {
00154
00155
00156 domain_.push (nVars (), x, domain_.lb (), domain_.ub (), false);
00157
00158
00159
00160 if (ndefined_ > 0)
00161 for (int i = 0; i < nVars (); ++i) {
00162 int ii = numbering_ [i];
00163 if (ii >= nOrigVars_ - ndefined_ &&
00164 ii < nOrigVars_)
00165 X (ii) = (*(commonexprs_ [ii - nOrigVars_ + ndefined_])) ();
00166 }
00167
00168
00169
00170
00171
00172
00173 for (int j=0, i=nVars (); i--; j++) {
00174
00175 int index = numbering_ [j];
00176 exprVar *var = variables_ [index];
00177
00178 if (var -> Multiplicity () > 0) {
00179
00180 CouNumber l, u;
00181
00182 if (var -> Type () == AUX)
00183 var -> Image () -> getBounds (l,u);
00184 else {
00185 l = Lb (index);
00186 u = Ub (index);
00187 }
00188
00189 if (var -> Type () == AUX) {
00190
00191 CouNumber &x = X (index);
00192
00193 bool isInt = var -> isDefinedInteger ();
00194
00195
00196
00197
00198
00199
00200 if ((var -> sign () == expression::AUX_EQ) &&
00201 ((index >= nOrigVars_) ||
00202 (index < nOrigVars_ - ndefined_)))
00203 x = (*(var -> Image ())) ();
00204
00205 x =
00206 CoinMax ((var -> sign () != expression::AUX_LEQ) ? (isInt ? ceil (l - COUENNE_EPS) : l) : -COIN_DBL_MAX,
00207 CoinMin ((var -> sign () != expression::AUX_GEQ) ? (isInt ? floor (u + COUENNE_EPS) : u) : COIN_DBL_MAX, x));
00208
00209
00210
00211 if (isInt) {
00212 if (var -> sign () == expression::AUX_GEQ) x = ceil (x - COUENNE_EPS);
00213 if (var -> sign () == expression::AUX_LEQ) x = floor (x + COUENNE_EPS);
00214 }
00215
00216
00217 }
00218 } else X (index) = 0.;
00219 }
00220
00221 restoreUnusedOriginals ();
00222
00223 domain_.pop ();
00224 }
00225
00226
00229
00230 void CouenneProblem::fillObjCoeff (double *&obj) {
00231
00232
00233
00234
00235
00236 expression *body = objectives_ [0] -> Body ();
00237
00238
00239 switch (body -> code ()) {
00240
00241 case COU_EXPRVAR:
00242 obj [body -> Index ()] = 1;
00243 break;
00244
00245 case COU_EXPRSUB: {
00246
00247 expression **arglist = body -> ArgList ();
00248
00249 obj [arglist [0] -> Index ()] = 1;
00250 obj [arglist [1] -> Index ()] = -1;
00251
00252 } break;
00253
00254 case COU_EXPRGROUP: {
00255
00256 exprGroup *eg = dynamic_cast <exprGroup *> (body -> isaCopy () ?
00257 body -> Copy () :
00258 body);
00259
00260 const exprGroup::lincoeff &lcoe = eg -> lcoeff ();
00261
00262
00263
00264
00265 for (int n = (int) lcoe.size (), i=0; n--; i++)
00266
00267 obj [lcoe [i]. first -> Index ()] = lcoe [i]. second;
00268
00269
00270
00271
00272 }
00273
00274 case COU_EXPRSUM: {
00275
00276 expression **arglist = body -> ArgList ();
00277
00278 for (int i = body -> nArgs (); i--;)
00279 switch ((arglist [i]) -> code ()) {
00280
00281 case COU_EXPRCONST:
00282 break;
00283
00284 case COU_EXPRVAR:
00285 obj [arglist [i] -> Index ()] = 1;
00286 break;
00287
00288 case COU_EXPRMUL: {
00289
00290 expression **mulArgList = arglist [i] -> ArgList ();
00291 int index = mulArgList [0] -> Index ();
00292
00293 if (index >= 0) obj [index] = mulArgList [1] -> Value ();
00294 else obj [mulArgList [1] -> Index ()] = mulArgList [0] -> Value ();
00295 } break;
00296
00297 default:
00298 Jnlst()->Printf(Ipopt::J_ERROR, J_PROBLEM,
00299 "Couenne: invalid element of sum\nAborting\n");
00300 exit (-1);
00301 }
00302 } break;
00303
00304 case COU_EXPRCONST: break;
00305
00306 default:
00307 Jnlst()->Printf(Ipopt::J_WARNING, J_PROBLEM,
00308 "Couenne: warning, objective function not recognized\n");
00309 break;
00310 }
00311 }
00312
00313
00315 void CouenneProblem::setCutOff (CouNumber cutoff, const double *s) const {
00316
00317 int indobj = objectives_ [0] -> Body () -> Index ();
00318
00319
00320
00321 if ((indobj >= 0) && (cutoff < pcutoff_ -> getCutOff () - COUENNE_EPS)) {
00322
00323
00324 Jnlst () -> Printf (Ipopt::J_WARNING, J_PROBLEM, "Couenne: new MINLP solution, value %.10e\n", cutoff);
00325
00326
00327 if (Var (indobj) -> isInteger ())
00328 pcutoff_ -> setCutOff (this, floor (cutoff + COUENNE_EPS), s);
00329 else pcutoff_ -> setCutOff (this, cutoff, s);
00330 }
00331 }
00332
00334 void CouenneProblem::resetCutOff (CouNumber value) const {
00335
00336 int indobj = objectives_ [0] -> Body () -> Index ();
00337
00338 if ((indobj >= 0)) {
00339
00340 if (Var (indobj) -> isInteger ())
00341 pcutoff_ -> setCutOff (this, floor (value + COUENNE_EPS), NULL);
00342 else pcutoff_ -> setCutOff (this, value, NULL);
00343 }
00344 }
00345
00346
00349 void CouenneProblem::installCutOff () const {
00350
00351
00352 double cutoff = pcutoff_ -> getCutOff();
00353
00354 if (cutoff > COUENNE_INFINITY)
00355 return;
00356
00357 int indobj = objectives_ [0] -> Body () -> Index ();
00358
00359 assert (indobj >= 0);
00360
00361
00362
00363
00364
00365 cutoff = (Var (indobj) -> isInteger ()) ?
00366 floor (cutoff + COUENNE_EPS) :
00367 (cutoff + CoinMin (SafeDelta, SafeCutoff * (1. + fabs (cutoff))));
00368
00369 if (cutoff < Ub (indobj))
00370 Ub (indobj) = cutoff;
00371 }
00372
00373
00374
00375 void CouenneProblem::realign () {
00376
00377
00378 for (std::vector <exprVar *>::iterator i = variables_.begin ();
00379 i != variables_.end (); ++i) {
00380
00381 (*i) -> linkDomain (&domain_);
00382 (*i) -> realign (this);
00383 if ((*i) -> Type () == AUX)
00384 (*i) -> Image () -> realign (this);
00385 }
00386
00387
00388 for (std::vector <CouenneObjective *>::iterator i = objectives_.begin ();
00389 i != objectives_.end (); ++i)
00390 (*i) -> Body () -> realign (this);
00391
00392
00393
00394 for (std::vector <CouenneConstraint *>::iterator i = constraints_.begin ();
00395 i != constraints_.end (); ++i)
00396 (*i) -> Body () -> realign (this);
00397 }
00398
00399
00401 void CouenneProblem::registerOptions (Ipopt::SmartPtr <Bonmin::RegisteredOptions> roptions) {
00402
00403 roptions -> AddNumberOption
00404 ("art_cutoff",
00405 "Artificial cutoff",
00406 COIN_DBL_MAX,
00407 "Default value is infinity.");
00408
00409 roptions -> AddNumberOption
00410 ("opt_window",
00411 "Window around known optimum",
00412 COIN_DBL_MAX,
00413 "Default value is infinity.");
00414
00415 roptions -> AddStringOption2
00416 ("use_semiaux",
00417 "Use semiauxiliaries, i.e. auxiliaries defined as w >= f(x) rather than w := f(x))",
00418 "yes",
00419 "no", "Only use auxiliaries assigned with \"=\" ",
00420 "yes", "Use auxiliaries defined by w <= f(x), w >= f(x), and w = f(x)"
00421 );
00422
00423 roptions -> AddStringOption2
00424 ("use_auxcons",
00425 "Use constraints-defined auxiliaries, i.e. auxiliaries w = f(x) defined by original constraints f(x) - w = 0",
00426 "yes",
00427 "no", "",
00428 "yes", ""
00429 );
00430
00431 roptions -> AddStringOption2
00432 ("redcost_bt",
00433 "Reduced cost bound tightening",
00434 "yes",
00435 "no","",
00436 "yes","",
00437 "This bound reduction technique uses the reduced costs of the LP in order to infer better variable bounds.");
00438
00439 roptions -> AddStringOption2
00440 ("use_quadratic",
00441 "Use quadratic expressions and related exprQuad class",
00442 "no",
00443 "no","Use an auxiliary for each bilinear term",
00444 "yes","Create only one auxiliary for a quadratic expression",
00445 "If enabled, then quadratic forms are not reformulated and therefore decomposed as a sum of auxiliary variables, each associated with a bilinear term, but rather taken as a whole expression. "
00446 "Envelopes for these expressions are generated through alpha-convexification."
00447 );
00448
00449 roptions -> AddStringOption2
00450 ("optimality_bt",
00451 "Optimality-based (expensive) bound tightening (OBBT)",
00452 "yes",
00453 "no","",
00454 "yes","",
00455 "This is another bound reduction technique aiming at reducing the solution set by looking at the initial LP relaxation. "
00456 "This technique is computationally expensive, and should be used only when necessary."
00457 );
00458
00459 roptions -> AddLowerBoundedIntegerOption
00460 ("log_num_obbt_per_level",
00461 "Specify the frequency (in terms of nodes) for optimality-based bound tightening.",
00462 -1,1,
00463 "\
00464 If -1, apply at every node (expensive!). \
00465 If 0, apply at root node only. \
00466 If k>=0, apply with probability 2^(k - level), level being the current depth of the B&B tree.");
00467
00468 roptions -> AddLowerBoundedIntegerOption
00469 ("max_fbbt_iter",
00470 "Number of FBBT iterations before stopping even with tightened bounds.",
00471 -1, MAX_FBBT_ITER,
00472 "Set to -1 to impose no upper limit");
00473
00474 roptions -> AddStringOption2
00475 ("aggressive_fbbt",
00476 "Aggressive feasibility-based bound tightening (to use with NLP points)",
00477 "yes",
00478 "no","",
00479 "yes","",
00480 "Aggressive FBBT is a version of probing that also allows to reduce the solution set, although it is not as quick as FBBT. "
00481 "It can be applied up to a certain depth of the B&B tree -- see ``log_num_abt_per_level''. "
00482 "In general, this option is useful but can be switched off if a problem is too large and seems not to benefit from it."
00483 );
00484
00485 roptions -> AddLowerBoundedIntegerOption
00486 ("log_num_abt_per_level",
00487 "Specify the frequency (in terms of nodes) for aggressive bound tightening.",
00488 -1,2,
00489 "\
00490 If -1, apply at every node (expensive!). \
00491 If 0, apply at root node only. \
00492 If k>=0, apply with probability 2^(k - level), level being the current depth of the B&B tree.");
00493
00494 roptions -> AddNumberOption
00495 ("art_lower",
00496 "Artificial lower bound",
00497 -COIN_DBL_MAX,
00498 "Default value is -COIN_DBL_MAX.");
00499
00500 roptions -> AddStringOption3
00501 ("branching_object",
00502 "type of branching object for variable selection",
00503 "var_obj",
00504 "vt_obj", "use Violation Transfer from Tawarmalani and Sahinidis",
00505 "var_obj", "use one object for each variable",
00506 "expr_obj", "use one object for each nonlinear expression");
00507
00508 roptions -> AddStringOption2
00509 ("delete_redundant",
00510 "Eliminate redundant variables, which appear in the problem as x_k = x_h",
00511 "yes",
00512 "no","Keep redundant variables, making the problem a bit larger",
00513 "yes","Eliminate redundant variables (the problem will be equivalent, only smaller)");
00514
00515 roptions -> AddStringOption4
00516 ("quadrilinear_decomp",
00517 "type of decomposition for quadrilinear terms (see work by Cafieri, Lee, Liberti)",
00518 "rAI",
00519 "rAI", "Recursive decomposition in bilinear terms (as in Ryoo and Sahinidis): x5 = ((x1 x2) x3) x4)",
00520 "tri+bi", "Trilinear and bilinear term: x5 = (x1 (x2 x3 x4))",
00521 "bi+tri", "Bilinear, THEN trilinear term: x5 = ((x1 x2) x3 x4))",
00522 "hier-bi", "Hierarchical decomposition: x5 = ((x1 x2) (x3 x4))");
00523 }