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

Go to the documentation of this file.
00001 /* $Id: checkNLP.cpp 298 2010-02-13 14:16:56Z pbelotti $
00002  *
00003  * Name:    checkNLP.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: check NLP feasibility of incumbent integer solution
00006  *
00007  * (C) Carnegie-Mellon University, 2006-10.
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 #include "CoinHelperFunctions.hpp"
00012 #include "CouenneProblem.hpp"
00013 
00014 // check if solution is MINLP feasible
00015 bool CouenneProblem::checkNLP (const double *solution, double &obj, bool recompute) const {
00016 
00017   if (Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM)) {
00018 
00019     printf ("checking solution: [%.12e] ", obj);
00020 
00021     for (int i=0; i<nOrigVars_; i++)
00022       printf ("%.12e ", solution [i]);
00023     printf ("\n");
00024   }
00025 
00026   // pre-check on original variables --- this is done after every LP,
00027   // and should be efficient
00028   for (register int i=0; i < nOrigVars_; i++) {
00029 
00030     CouNumber val = solution [i];
00031 
00032     // check (original and auxiliary) variables' integrality
00033 
00034     exprVar *v = variables_ [i];
00035 
00036     if ((v -> Type ()      == VAR) &&
00037         (v -> Multiplicity () > 0) &&
00038         (v -> isInteger ())        &&
00039         (fabs (val - COUENNE_round (val)) > feas_tolerance_)) {
00040 
00041       Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_PROBLEM,
00042                       "checkNLP: integrality %d violated: %.6f [%g,%g]\n", 
00043                       i, val, domain_.lb (i), domain_.ub (i));
00044 
00045       return false;
00046     }
00047   }
00048 
00049   const int infeasible = 1;
00050   const int wrong_obj  = 2;
00051 
00052   CouNumber *sol = new CouNumber [nVars ()];
00053 
00054   // copy solution, evaluate the corresponding aux, and then replace
00055   // the original variables again for checking
00056   CoinCopyN (solution, nOrigVars_, sol);
00057   getAuxs (sol);
00058   CoinCopyN (solution, nOrigVars_, sol);
00059 
00060   // install NL solution candidate in evaluation structure
00061   domain_.push (nVars (), sol, domain_.lb (), domain_.ub (), false);
00062 
00063 
00064   if (Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM)) {
00065     printf ("checknlp: %d vars -------------------\n", domain_.current () -> Dimension ());
00066     for (int i=0; i<domain_.current () -> Dimension (); i++)
00067       printf ("%4d %.12e [%.12e %.12e]\n", 
00068               i, domain_.x (i), domain_.lb (i), domain_.ub (i));
00069   }
00070 
00071   expression *objBody = Obj (0) -> Body ();
00072 
00073   // BUG: if Ipopt solution violates bounds of original variables and
00074   // objective depends on originals, we may have a "computed object"
00075   // out of bounds
00076 
00077   //CouNumber realobj = (*(objBody -> Image () ? objBody -> Image () : objBody)) ();
00078   CouNumber realobj = obj;
00079 
00080   if (objBody) 
00081     realobj = 
00082       (objBody -> Index () >= 0) ?
00083       sol [objBody -> Index ()] : 
00084       (*(objBody -> Image () ? objBody -> Image () : objBody)) ();
00085 
00086   if (Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM)) {
00087     printf ("%.12e %.12e %.12e ------------------------------\n", 
00088             realobj, sol [objBody -> Index ()], 
00089             (*(objBody -> Image () ? objBody -> Image () : objBody)) ());
00090   }
00091 
00092   bool retval = true;
00093 
00094   try {
00095 
00096     // check if objective corresponds
00097     
00098     if (fabs (realobj - obj) / (1. + fabs (realobj)) > feas_tolerance_) {
00099 
00100       Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_PROBLEM,
00101                       "checkNLP: false objective. %g != %g (diff. %g)\n", 
00102                       realobj, obj, realobj - obj);
00103 
00104       if (!recompute)
00105         throw wrong_obj;
00106     }
00107 
00108     if (recompute)
00109       obj = realobj;
00110 
00111     if (Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM))
00112       printf ("recomputed: %.12e\n", obj);
00113 
00114     for (int i=0; i < nOrigVars_; i++) {
00115 
00116       if (variables_ [i] -> Multiplicity () <= 0) 
00117         continue;
00118 
00119       CouNumber val = domain_.x (i);
00120 
00121       // check bounds
00122 
00123       if ((val > domain_.ub (i) + feas_tolerance_) ||
00124           (val < domain_.lb (i) - feas_tolerance_)) {
00125 
00126         Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_PROBLEM,
00127                         "checkNLP: variable %d out of bounds: %.6f [%g,%g] (diff %g)\n", 
00128                         i, val, domain_.lb (i), domain_.ub (i),
00129                         CoinMax (fabs (val - domain_.lb (i)), 
00130                                  fabs (val - domain_.ub (i))));
00131         throw infeasible;
00132       }
00133 
00134       // check (original and auxiliary) variables' integrality
00135 
00136       if (variables_ [i] -> isInteger () &&
00137           (fabs (val - COUENNE_round (val)) > feas_tolerance_)) {
00138 
00139         Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_PROBLEM,
00140                         "checkNLP: integrality %d violated: %.6f [%g,%g]\n", 
00141                         i, val, domain_.lb (i), domain_.ub (i));
00142 
00143         throw infeasible;
00144       }
00145     }
00146 
00147     // check ALL auxs
00148 
00149     for (int i=0; i < nVars (); i++) {
00150 
00151       if (Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM)) {
00152         if (variables_ [i] -> Type () == AUX) {
00153           printf ("checking aux ");
00154           variables_ [i] -> print (); printf (" := ");
00155           variables_ [i] -> Image () -> print (); 
00156           printf (" --- %.12e = %.12e [%.12e]; {", 
00157                   (*(variables_ [i])) (), 
00158                   (*(variables_ [i] -> Image ())) (),
00159                   (*(variables_ [i])) () -
00160                   (*(variables_ [i] -> Image ())) ());
00161           //for (int j=0; j<nVars (); j++)
00162           //printf ("%.12e ", (*(variables_ [j])) ());
00163           printf ("}\n");
00164         }
00165       }
00166 
00167       CouNumber delta;
00168 
00169       // check if auxiliary has zero infeasibility
00170 
00171       if ((variables_ [i] -> Type () == AUX) && 
00172           (fabs (delta = (*(variables_ [i])) () - 
00173                  (*(variables_ [i] -> Image ())) ()) > feas_tolerance_)) {
00174 
00175         Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_PROBLEM,
00176                         "checkNLP: auxiliarized %d violated (%g)\n", i, delta);
00177 
00178         throw infeasible;
00179       }
00180     }
00181 
00182     // check constraints
00183 
00184     for (int i=0; i < nCons (); i++) {
00185 
00186       CouenneConstraint *c = Con (i);
00187 
00188       CouNumber
00189         body = (*(c -> Body ())) (),
00190         lhs  = (*(c -> Lb   ())) (),
00191         rhs  = (*(c -> Ub   ())) ();
00192 
00193       if (((rhs < COUENNE_INFINITY) &&
00194            (body > rhs + feas_tolerance_ * (1 + CoinMax (fabs (body), fabs (rhs))))) || 
00195           ((lhs > -COUENNE_INFINITY) &&
00196            (body < lhs - feas_tolerance_ * (1 + CoinMax (fabs (body), fabs (lhs)))))) {
00197 
00198         if (Jnlst()->ProduceOutput(Ipopt::J_MOREVECTOR, J_PROBLEM)) {
00199 
00200           printf ("checkNLP: constraint %d violated (lhs=%+e body=%+e rhs=%+e, violation %g): ",
00201                   i, lhs, body, rhs, CoinMax (lhs-body, body-rhs));
00202 
00203           c -> print ();
00204         }
00205 
00206         throw infeasible;
00207       }
00208     }
00209   }
00210 
00211   catch (int exception) {
00212 
00213     switch (exception) {
00214 
00215     case wrong_obj:
00216       retval = false;
00217       break;
00218 
00219     case infeasible:
00220     default:
00221       retval = false;
00222       break;
00223     }
00224   }
00225 
00226   delete [] sol;
00227   domain_.pop ();
00228 
00229   return retval;
00230 }

Generated on Mon May 3 03:05:21 2010 by  doxygen 1.4.7