/home/coin/SVN-release/OS-2.4.0/Couenne/src/heuristics/CouenneFPSolveMILP.cpp

Go to the documentation of this file.
00001 /* $Id: CouenneFPSolveMILP.cpp 720 2011-06-27 13:31:26Z pbelotti $
00002  *
00003  * Name:    CouenneFPSolveMILP.cpp
00004  * Authors: Pietro Belotti
00005  *          Timo Berthold, ZIB Berlin
00006  * Purpose: Solve the MILP within the Feasibility Pump 
00007  * 
00008  * This file is licensed under the Eclipse Public License (EPL)
00009  */
00010 
00011 #include "CbcModel.hpp"
00012 
00013 #include "CouenneConfig.h"
00014 #include "CouenneFeasPump.hpp"
00015 #include "CouenneMINLPInterface.hpp"
00016 #include "CouenneProblem.hpp"
00017 #include "CouenneProblemElem.hpp"
00018 #include "CouenneExprVar.hpp"
00019 
00020 #include "CouenneFPpool.hpp"
00021 
00022 #ifdef COIN_HAS_SCIP
00023 /* general SCIP includes */
00024 #include "scip/scip.h"
00025 #include "scip/cons_linear.h"
00026 #include "scip/scipdefplugins.h"
00027 #endif
00028 
00029 using namespace Couenne;
00030 
00031 #define NUMERICS_THRES 1e19
00032 
00033 #ifdef COIN_HAS_SCIP
00034 void CouenneFeasPump::checkInfinity(SCIP *scip, SCIP_Real val, double infinity){
00035    if( SCIPisInfinity(scip, val) && val < infinity)
00036       printf("Warning: %g will be considered to be Infinity by SCIP\n", val);
00037 }
00038 #endif
00039 
00040 
00042 OsiSolverInterface *createCloneMILP (const CouenneFeasPump *fp, CbcModel *model, bool isMILP);
00043 
00044 
00046 void addDistanceConstraints (const CouenneFeasPump *fp, OsiSolverInterface *lp, double *sol, bool isMILP);
00047 
00048 
00052 CouNumber CouenneFeasPump::solveMILP (CouNumber *nSol0, CouNumber *&iSol, int niter, int* nsuciter) {
00053 
00054   // The problem is of the form
00055   //
00056   // min  f(x)
00057   // s.t. Ax >= b
00058   //      x_i in Z, i in N
00059   //
00060   // where N is the index set of integer variables and f(x) is one of
00061   // the following:
00062   //
00063   // 1) sum {i in Vars} |x_i - x_i^0|
00064   // 2) sum {i in N}    |x_i - x_i^0|
00065   // 3) sum {i in Vars} |P^i (x - x^0)| 
00066   // 4) sum {i in N}    |P^i (x - x^0)|
00067   //
00068   // where is x^0 is the optimal solution of a NLP problem. In the
00069   // last two, the l-1 norm is multiplied by the i-th row of a matrix
00070   // P obtained from the Hessian H of the Lagrangian of the problem. H
00071   // is a positive semidefinite matrix only on the null space of the
00072   // gradient g, i.e.
00073   //
00074   // x' H x >= 0 for all x: g'x = 0
00075   //
00076   // (where x' is x transposed), therefore we add a quadratic term to
00077   // make it positive semidefinite everywhere else. In order to do so,
00078   // the quadratic term (alpha is a scalar)
00079   //
00080   // alpha (g'x)^2
00081   // 
00082   // is zero in the null space of g, and it is strictly positive
00083   // everywhere else. Hence, we need an alpha so that 
00084   //
00085   // P = (H + alpha g g') 
00086   //
00087   // is PSD. In general, we might have a parameter beta in [0,1] such
00088   // that
00089   // 
00090   // P = beta I + (1-beta) (H + alpha g g')
00091   //
00092   // so that we can balance the Hessian and the distance.
00093 
00094   bool firstCall = (milp_ == NULL); // true if this is the first call to
00095                                     // solveMILP; initialization will be
00096                                     // necessary
00097 
00098   if (firstCall) {
00099 
00100     // create MILP
00101 
00102     milp_ = createCloneMILP (this, model_, true);
00103 
00104     // Post-processing LP: persistent if FP_DIST_POST, created on the
00105     // fly if FP_DIST_INT and numerics, not created if FP_DIST_ALL
00106     //
00107     // set up an LP as a copy of the original MILP. Don't do the
00108     // same for FP_DIST_INT as it is only necessary upon numerical
00109     // problems, which might not happen
00110 
00111     if ((compDistInt_ == FP_DIST_POST) && !postlp_)
00112       postlp_ = createCloneMILP (this, model_, false);
00113   }
00114 
00115   int nInitRows = milp_ -> getNumRows ();
00116 
00117   CouNumber * nlpSolExp;
00118 
00119   if (nSol0) {
00120 
00121     nlpSolExp = new CouNumber [problem_ -> nVars ()];
00122 
00123     CoinCopyN (nSol0, problem_ -> nOrigVars (), nlpSolExp);
00124     problem_ -> getAuxs (nlpSolExp);
00125 
00126   } else 
00127     nlpSolExp = CoinCopyOfArray (milp_ -> getColSolution (), 
00128                                  problem_ -> nVars ());
00129 
00130   // create constraints to define l_1 distance objective function
00131   addDistanceConstraints (this, milp_, nlpSolExp, true);
00132 
00133   delete [] nlpSolExp;
00134 
00135   int nFinalRows = milp_ -> getNumRows ();
00136 
00137   // The MILP is complete. We have several ways of solving it, or
00138   // better, to find feasible solutions to it. We have to interface
00139   // with each of them once at the very beginning, and later we loop
00140   // through them in order to find a feasible solution.
00141 
00142   if (firstCall)
00143     init_MILP ();
00144 
00145   if (false) { // should always be false upon commit
00146     static int cntr = 0;
00147     char filename [30];
00148     sprintf (filename, "fp-milp%04d", cntr++);
00149     milp_ -> writeLp (filename);
00150   }
00151 
00152   double obj = findSolution (iSol, niter, nsuciter);
00153 
00154   if ((nSol0 && iSol) &&
00155       (problem_ -> Jnlst () -> ProduceOutput (Ipopt::J_ERROR, J_NLPHEURISTIC))) {
00156 
00157     double dist = 0.;
00158     int nNonint = 0;
00159 
00160     for (int i = 0; i < problem_ -> nVars (); ++i) {
00161 
00162       if (problem_ -> Var (i) -> isInteger () &&
00163           (fabs (iSol [i] - ceil (iSol [i] - .5)) > 1e-4))
00164         ++nNonint;
00165 
00166       dist += 
00167         (iSol [i] - nSol0 [i]) * 
00168         (iSol [i] - nSol0 [i]);
00169     }
00170 
00171     printf ("FP: after MILP, distance %g, %d nonintegers\n", sqrt (dist), nNonint);
00172   }
00173 
00174   //
00175   // POST PROCESSING 
00176   //
00177   // (if we got a solution from MILP, otherwise bail out)
00178   //
00179 
00180   if (iSol &&
00181       (compDistInt_ != FP_DIST_ALL)) {
00182 
00183     // check iSol for numerics (i.e. components whose fabs () is large,
00184     // or >= 1e20) or post-process to obtain a second solution by fixing
00185     // the integer coordinates and solving the resulting LP
00186 
00187       bool numerics = false;
00188 
00189       if (compDistInt_ == FP_DIST_INT) {
00190 
00191         for (std::vector <exprVar *>::iterator i = problem_ -> Variables (). begin (); 
00192              i != problem_ -> Variables (). end (); ++i)
00193 
00194           if ((  (*i) -> Multiplicity () > 0) &&
00195               ! ((*i) -> isInteger    ())     &&
00196               (fabs (iSol [(*i) -> Index ()]) > NUMERICS_THRES)) {
00197 
00198             numerics = true;
00199             break;
00200           }
00201       }
00202 
00203       if (numerics || (compDistInt_ == FP_DIST_POST)) {
00204 
00205         // solve LP where integer variables have been fixed:
00206         //
00207         // a) check if postlp_ exists yet
00208         // 0) save integer bounds
00209         // 1) fix integer variables 
00210         // 2) add variables and inequalities
00211         // 3) solve LP
00212         // 4) if optimal, save solution
00213         // 5) restore IP bounds
00214         // 6) delete variables
00215 
00216         if (!postlp_)
00217           postlp_ = createCloneMILP (this, model_, false);
00218 
00219         int nvars = postlp_ -> getNumCols ();
00220 
00221         // save integer bounds to restore them later
00222         double
00223           *saveLB = CoinCopyOfArray (postlp_ -> getColLower (), nvars),
00224           *saveUB = CoinCopyOfArray (postlp_ -> getColUpper (), nvars),
00225           *newLB  = CoinCopyOfArray (postlp_ -> getColLower (), nvars),
00226           *newUB  = CoinCopyOfArray (postlp_ -> getColUpper (), nvars);
00227 
00228         // fix integer variables
00229 
00230         for (int i = problem_ -> nVars (); i--;)
00231           if (milp_ -> isInteger (i))
00232             newLB [i] = newUB [i] = iSol [i];
00233 
00234         postlp_ -> setColLower (newLB);
00235         postlp_ -> setColUpper (newUB);
00236 
00237         // add inequalities
00238 
00239         int nInitRowsLP  = postlp_ -> getNumRows ();
00240         addDistanceConstraints (this, postlp_, iSol, false);
00241         int nFinalRowsLP = postlp_ -> getNumRows ();
00242 
00243         // Solve the LP, obtain closest point with integer variables fixed
00244 
00245         postlp_ -> initialSolve ();
00246 
00247         // save as new solution
00248 
00249         if (postlp_ -> isProvenOptimal ()) 
00250           CoinCopyN (postlp_ -> getColSolution (), problem_ -> nVars (), iSol);
00251 
00252         postlp_ -> setColLower (saveLB);
00253         postlp_ -> setColUpper (saveUB);
00254 
00255         // delete temp data
00256 
00257         delete [] saveLB;
00258         delete [] saveUB;
00259         delete [] newLB;
00260         delete [] newUB;
00261 
00262         // delete added rows
00263 
00264         int 
00265           nDeleted = nFinalRowsLP - nInitRowsLP,
00266          *deleted  = new int [nDeleted],
00267           nCurRow  = nInitRowsLP;
00268 
00269         for (int i = nDeleted; i--;)
00270           deleted [i] = nCurRow++;
00271 
00272         postlp_ -> deleteRows (nDeleted, deleted);
00273 
00274         delete [] deleted;
00275       }
00276   }
00277 
00278   // delete last rows and add them from scratch (common block below)
00279 
00280   int 
00281     nDeleted = nFinalRows - nInitRows,
00282    *deleted  = new int [nDeleted],
00283     nCurRow  = nInitRows;
00284 
00285   for (int i = nDeleted; i--;)
00286     deleted [i] = nCurRow++;
00287 
00288   milp_ -> deleteRows (nDeleted, deleted);
00289 
00290   delete [] deleted;
00291 
00292   return obj;
00293 }

Generated on Thu Sep 22 03:05:58 2011 by  doxygen 1.4.7