/home/coin/SVN-release/OS-2.4.1/Couenne/src/problem/problem.cpp

Go to the documentation of this file.
00001 /* $Id: problem.cpp 752 2011-08-08 03:45:07Z pbelotti $
00002  *
00003  * Name:    problem.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: methods of the class CouenneProblem
00006  *
00007  * (C) Carnegie-Mellon University, 2006-10.
00008  * This file is licensed under the Eclipse Public License (EPL)
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 // tricky... smaller values cut the optimum in OS unitTest
00043 const CouNumber SafeCutoff = 1e-4; 
00044 
00045 // absolute difference
00046 const CouNumber SafeDelta = 1e-2; 
00047 
00050 
00051 void CouenneProblem::initAuxs () const {
00052 
00053   domain_.current () -> resize (nVars ());
00054 
00055   // initially, auxiliary variables are unbounded, their bounds only
00056   // depending on their function
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) &&     // this is an auxiliary
00065          (indvar >= nOrigVars_))            ||     // and not an original, originally
00066         (variables_ [i] -> Multiplicity () == 0))  // or a useless one
00067       //int index = variables_ [i] -> Index ();
00068       Lb (indvar) = - (Ub (indvar) = COIN_DBL_MAX);
00069   }
00070 
00071   // first initialize with values from constraints
00072 
00073   //Jnlst()->Printf(Ipopt::J_VECTOR, J_PROBLEM, "Initial bounds for aux (initAuxs):\n");
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   // only one loop is sufficient here, since auxiliary variable are
00091   // defined in such a way that w_i does NOT depend on w_j if i<j.
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     // ignore these variables!
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     // and handle only those with nonzero multiplicity
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       /*printf ("printing bounds: [%g %g]\n", Lb (ord), Ub (ord));
00119       var -> Lb () -> print (); printf ("\n");
00120       var -> Ub () -> print (); printf ("\n");*/
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       // set bounds 
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       //if ((lb_ [ord] = (*(aux -> Lb ())) ()) <= -COUENNE_INFINITY) lb_ [ord] = -DBL_MAX;
00131       //if ((ub_ [ord] = (*(aux -> Ub ())) ()) >=  COUENNE_INFINITY) ub_ [ord] =  DBL_MAX;
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 (); // no argument as the local x vector will be used
00148 }
00149 
00150 
00153 void CouenneProblem::getAuxs (CouNumber * x) const {
00154 
00155   // set point at x, don't copy
00156   domain_.push (nVars (), x, domain_.lb (), domain_.ub (), false);
00157 
00158   // if there are common expressions, set them (they usually don't get
00159   // set by the AMPL interface as they are not really variables)
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   // set auxiliary w to f(x). This procedure is exact even though the
00169   // auxiliary variables have an incomplete image, i.e. they have been
00170   // decomposed previously, since they are updated with increasing
00171   // index.
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         /*printf ("checking "); var -> print ();
00196         printf (" = %g = %g. Sign: %d, int: %d, [%g,%g]", 
00197                 X (index), (*(var -> Image ())) (),
00198                 var -> sign (), isInt, l, u);*/
00199 
00200         if ((var -> sign () == expression::AUX_EQ) &&
00201             ((index >= nOrigVars_) ||
00202              (index  < nOrigVars_ - ndefined_)))
00203           x = (*(var -> Image ())) ();  // addresses of x[] and X() are equal
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         // heuristic feasibility heuristic: if semiaux, round value to nearest integer if variable is integer
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         //printf (" -> %g\n", X (index));
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   // linearized objective can be an exprAux, an exprSub, an exprGroup,
00233   // or an exprSum. In the last two cases, the components are
00234   // variables or constants
00235 
00236   expression *body = objectives_ [0] -> Body ();
00237   //int sense = objectives_ [0] -> Sense ();
00238 
00239   switch (body -> code ()) {
00240 
00241   case COU_EXPRVAR:   //
00242     obj [body -> Index ()] = 1; //(sense == MINIMIZE) ? 1 : -1;
00243     break;
00244 
00245   case COU_EXPRSUB: { // 
00246 
00247     expression **arglist = body -> ArgList ();
00248 
00249     obj [arglist [0] -> Index ()] =  1; //(sense == MINIMIZE) ?  1 : -1;
00250     obj [arglist [1] -> Index ()] = -1; //(sense == MINIMIZE) ? -1 :  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     //    if (sense == MINIMIZE) while (*index >= 0) obj [*index++] =  *coeff++;
00263     //    else                   while (*index >= 0) obj [*index++] = -*coeff++;      
00264 
00265     for (int n = (int) lcoe.size (), i=0; n--; i++)
00266       //exprGroup::lincoeff::iterator el = lcoe.begin (); el != lcoe.end (); ++el)
00267       obj [lcoe [i]. first -> Index ()] = lcoe [i]. second;
00268     //(sense == MINIMIZE) ? 
00269     //(lcoe [i]. second) : 
00270     //-(lcoe [i]. second);
00271 
00272   } // no break, as exprGroup is derived from exprSum
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; //(sense == MINIMIZE) ? 1 : -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; // a constant objective
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   // AW: Should we use the value of the objective variable computed by 
00320   //     Couenne here?
00321   if ((indobj >= 0) && (cutoff < pcutoff_ -> getCutOff () - COUENNE_EPS)) {
00322 
00323     //if (fabs (cutoff - pcutoff_ -> getCutOff ()) > (1 + fabs (cutoff)) * 2 * SafeCutoff) // avoid too many printouts
00324     Jnlst () -> Printf (Ipopt::J_WARNING, J_PROBLEM, "Couenne: new MINLP solution, value %.10e\n", cutoff);
00325                         //pcutoff_ -> getCutOff ());
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   // all problem are assumed to be minimization
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   //Jnlst () -> Printf (Ipopt::J_WARNING, J_PROBLEM,
00362   //"installing cutoff %.10e vs current ub %.10e\n",
00363   //cutoff, Ub (indobj));
00364 
00365   cutoff = (Var (indobj) -> isInteger ()) ?
00366     floor (cutoff + COUENNE_EPS) :
00367     (cutoff + CoinMin (SafeDelta, SafeCutoff * (1. + fabs (cutoff))));  // tolerance needed to retain feasibility
00368 
00369   if (cutoff < Ub (indobj))
00370     Ub (indobj) = cutoff;
00371 }
00372 
00373 
00374 // clear all spurious variables pointers not referring to the variables_ vector
00375 void CouenneProblem::realign () {
00376 
00377   // link variables to problem's domain
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   // link variables to problem's domain
00388   for (std::vector <CouenneObjective *>::iterator i = objectives_.begin ();
00389        i != objectives_.end (); ++i) 
00390     (*i) -> Body () -> realign (this);
00391 
00392 
00393   // link variables to problem's domain
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 }

Generated on Thu Nov 10 03:05:46 2011 by  doxygen 1.4.7