source: trunk/Couenne/src/main/BonCouenneInterface.cpp @ 115

Last change on this file since 115 was 115, checked in by pbelotti, 11 years ago

fix right value of objective computed in checkNLP -- numerics can make a bound interval infeasible. Some cleanup

File size: 10.9 KB
Line 
1// (C) Copyright International Business Machines Corporation (IBM) 2006, 2007
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// Authors :
6// Pietro Belotti, Carnegie Mellon University
7// Pierre Bonami, International Business Machines Corporation
8//
9// Date : 12/19/2006
10
11
12#include "BonCouenneInterface.hpp"
13#include "CoinHelperFunctions.hpp"
14#include "CouenneProblem.hpp"
15
16namespace Bonmin {
17
18/** Default constructor. */
19CouenneInterface::CouenneInterface():
20  AmplInterface(),
21  have_nlp_solution_ (false)
22{}
23
24/** Copy constructor. */
25CouenneInterface::CouenneInterface(const CouenneInterface &other):
26  AmplInterface(other),
27  have_nlp_solution_ (false)
28{}
29
30/** virtual copy constructor. */
31CouenneInterface * CouenneInterface::clone(bool CopyData){
32  return new CouenneInterface(*this);
33}
34
35/** Destructor. */
36CouenneInterface::~CouenneInterface(){
37}
38
39#ifdef COIN_HAS_ASL   
40void 
41CouenneInterface::readAmplNlFile(char **& argv, Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,
42                                 Ipopt::SmartPtr<Ipopt::OptionsList> options,
43                                 Ipopt::SmartPtr<Ipopt::Journalist> journalist){
44  //  if (!IsValid (app_))
45  //createApplication (roptions, options, journalist, "couenne.");
46  AmplInterface::readAmplNlFile(argv, roptions, options, journalist);
47}
48#endif
49
50/** \name Overloaded methods to build outer approximations */
51  //@{
52  /** \brief Extract a linear relaxation of the MINLP.
53   *
54   * Solve the continuous relaxation and takes first-order
55   * outer-approximation constraints at the optimum.  Then put
56   * everything in an OsiSolverInterface.
57   *
58   * The OsiSolverInterface si is empty and has to be populated with
59   * the initial linear relaxation.
60   */
61
62void
63CouenneInterface::extractLinearRelaxation
64(OsiSolverInterface &si, CouenneCutGenerator & couenneCg, bool getObj, bool solveNlp) {
65
66  CouenneProblem *p = couenneCg.Problem ();
67  bool is_feasible = true;
68
69  if (solveNlp) {
70
71    int nvars = p -> nVars();
72
73    if (p -> doFBBT ()) {
74
75      // include the rhs of auxiliary-based constraints into the FBBT
76      // (should be useful with Vielma's problems, for example)
77
78      for (int i=0; i < p -> nCons (); i++) {
79
80        // for each constraint
81        CouenneConstraint *con = p -> Con (i);
82
83        // (which has an aux as its body)
84        int index = con -> Body () -> Index ();
85
86        if ((index >= 0) && (con -> Body () -> Type () == AUX)) {
87
88          // if there exists violation, add constraint
89          CouNumber
90            l = con -> Lb () -> Value (),       
91            u = con -> Ub () -> Value ();
92
93          // tighten bounds in Couenne's problem representation
94          p -> Lb (index) = CoinMax (l, p -> Lb (index));
95          p -> Ub (index) = CoinMin (u, p -> Ub (index));
96        }
97      }
98
99      t_chg_bounds *chg_bds = new t_chg_bounds [nvars];
100
101      for (int i=0; i<nvars; i++) {
102        chg_bds [i].setLower(t_chg_bounds::CHANGED);
103        chg_bds [i].setUpper(t_chg_bounds::CHANGED);
104      }
105
106      if (!(p -> boundTightening (chg_bds, NULL))) {
107        is_feasible = false;
108        printf ("warning, tightened NLP is infeasible\n");
109      }
110
111      delete [] chg_bds;
112
113      const double 
114        *nlb = getColLower (),
115        *nub = getColUpper ();
116
117      for (int i=0; i < p -> nOrigVars (); i++) 
118        if (p -> Var (i) -> Multiplicity () > 0) {
119          if (nlb [i] < p -> Lb (i) - COUENNE_EPS) setColLower (i, p -> Lb (i));
120          if (nub [i] > p -> Ub (i) + COUENNE_EPS) setColUpper (i, p -> Ub (i));
121        } else { 
122          // if not enabled, fix them in the NLP solver
123          setColLower (i, -COIN_DBL_MAX);
124          setColUpper (i,  COIN_DBL_MAX);
125        }
126    }
127
128    if (is_feasible) 
129      initialSolve ();
130    else {
131      OsiAuxInfo * auxInfo = si.getAuxiliaryInfo ();
132      BabInfo * babInfo = dynamic_cast <BabInfo *> (auxInfo);
133
134      if (babInfo) 
135        babInfo -> setInfeasibleNode ();
136    }
137
138    if (is_feasible && isProvenOptimal ()) {
139
140      CouNumber obj             = getObjValue    ();
141      const CouNumber *solution = getColSolution ();
142
143      if (getNumIntegers () > 0) {
144
145        int
146          norig = p -> nOrigVars (),
147          nvars = p -> nVars ();
148
149        bool fractional = false;
150
151        // if problem is integer, check if any integral variable is
152        // fractional. If so, round them and re-optimize
153
154        for (int i=0; i<norig; i++)
155          if (-> Var (i) -> isInteger () &&
156              (p -> Var (i) -> Multiplicity () > 0) &&
157              (!::isInteger (solution [i]))) {
158            fractional = true;
159            break;
160          }
161
162        if (fractional) { // try again if solution found by Ipopt is fractional
163
164          double 
165            *lbSave = new double [norig],
166            *ubSave = new double [norig],
167
168            *lbCur  = new double [nvars],
169            *ubCur  = new double [nvars],
170
171            *Y      = new double [nvars];
172
173          CoinCopyN (getColLower (), norig, lbSave);
174          CoinCopyN (getColUpper (), norig, ubSave);
175
176          CoinFillN (Y,     nvars, 0.);
177          CoinFillN (lbCur, nvars, -COUENNE_INFINITY);
178          CoinFillN (ubCur, nvars,  COUENNE_INFINITY);
179
180          CoinCopyN (getColLower (), norig, lbCur);
181          CoinCopyN (getColUpper (), norig, ubCur);
182
183          if (p -> getIntegerCandidate (solution, Y, lbCur, ubCur) >= 0) {
184
185            for (int i = getNumCols (); i--;) {
186
187              if (lbCur [i] > ubCur [i]) {
188                double swap = lbCur [i];
189                lbCur [i] = ubCur [i];
190                ubCur [i] = swap;
191              }
192
193              if      (Y [i] < lbCur [i]) Y [i] = lbCur [i];
194              else if (Y [i] > ubCur [i]) Y [i] = ubCur [i];
195            }
196
197            for (int i=0; i<norig; i++)
198              if ((p -> Var (i) -> Multiplicity () > 0) &&
199                  p  -> Var (i) -> isInteger ()) {
200                setColLower (i, lbCur [i]);
201                setColUpper (i, ubCur [i]);
202              }
203
204            setColSolution (Y); // use initial solution given
205
206            resolve (); // solve with integer variables fixed
207
208            obj      = getObjValue ();
209            solution = getColSolution ();
210
211            // restore previous bounds on integer variables
212            for (int i=0; i<norig; i++)
213              if ((p -> Var (i) -> Multiplicity () > 0) &&
214                  p  -> Var (i) -> isInteger ()) {
215                setColLower (i, lbSave [i]);
216                setColUpper (i, ubSave [i]);
217              }
218          }
219
220          delete [] Y;
221          delete [] lbSave;
222          delete [] ubSave;
223          delete [] lbCur;
224          delete [] ubCur;
225        } 
226      }
227
228      // re-check optimality in case resolve () was called
229      if (isProvenOptimal () && 
230          (obj < p -> getCutOff ()) && // first check (before re-computing)
231          p -> checkNLP (solution, obj, true) && // true for recomputing obj
232          (obj < p -> getCutOff ())) { // second check (real object might be different)
233
234        // tell caller there is an initial solution to be fed to the initHeuristic
235        have_nlp_solution_ = true;
236
237        // set cutoff to take advantage of bound tightening
238        //printf ("new cutoff %g from BonCouenneInterface\n", obj);
239        p -> setCutOff (obj);
240
241        OsiAuxInfo * auxInfo = si.getAuxiliaryInfo ();
242        BabInfo * babInfo = dynamic_cast <BabInfo *> (auxInfo);
243
244        if (babInfo) {
245          babInfo -> setNlpSolution (solution, getNumCols (), obj);
246          babInfo -> setHasNlpSolution (true);
247        }
248      }
249    }
250  }
251
252  if (!is_feasible) // nothing else to do, problem infeasible
253    return;
254
255  int 
256    numcols     = getNumCols (), // # original               variables
257    numcolsconv = p -> nVars (); // # original + # auxiliary variables
258
259  const double
260    *lb = getColLower (),
261    *ub = getColUpper ();
262
263   // add original and auxiliary variables to the new problem
264   for (int i=0; i<numcols; i++) 
265     if (p -> Var (i) -> Multiplicity () > 0) si.addCol (0, NULL,NULL, lb [i],       ub [i],      0);
266     else                                     si.addCol (0, NULL,NULL, -COIN_DBL_MAX,COIN_DBL_MAX,0);
267   for (int i=numcols; i<numcolsconv; i++)    si.addCol (0, NULL,NULL, -COIN_DBL_MAX,COIN_DBL_MAX,0);
268
269   // get initial relaxation
270   OsiCuts cs;
271   couenneCg.generateCuts (si, cs);
272
273   // store all (original + auxiliary) bounds in the relaxation
274   CouNumber * colLower = new CouNumber [numcolsconv];
275   CouNumber * colUpper = new CouNumber [numcolsconv];
276
277   CouNumber *ll = p -> Lb ();
278   CouNumber *uu = p -> Ub ();
279
280   // overwrite original bounds, could have improved within generateCuts
281   for (int i = numcolsconv; i--;) 
282     if (p -> Var (i) -> Multiplicity () > 0) {
283       colLower [i] = (ll [i] > - COUENNE_INFINITY) ? ll [i] : -COIN_DBL_MAX;
284       colUpper [i] = (uu [i] <   COUENNE_INFINITY) ? uu [i] :  COIN_DBL_MAX;
285     } else {
286       colLower [i] = -COIN_DBL_MAX;
287       colUpper [i] =  COIN_DBL_MAX;
288     }
289
290   int numrowsconv = cs.sizeRowCuts ();
291
292   // create matrix and other stuff
293   CoinBigIndex * start = new CoinBigIndex [numrowsconv + 1];
294
295   int    * length   = new int    [numrowsconv];
296   double * rowLower = new double [numrowsconv];
297   double * rowUpper = new double [numrowsconv];
298
299   start[0] = 0;
300   int nnz = 0;
301   /* fill the four arrays. */
302   for(int i = 0 ; i < numrowsconv ; i++)
303   {
304     OsiRowCut * cut = cs.rowCutPtr (i);
305
306     const CoinPackedVector &v = cut->row();
307     start[i+1] = start[i] + v.getNumElements();
308     nnz += v.getNumElements();
309     length[i] = v.getNumElements();
310
311     rowLower[i] = cut->lb();
312     rowUpper[i] = cut->ub();
313   }
314
315   assert (nnz == start [numrowsconv]);
316   /* Now fill the elements arrays. */
317   int * ind = new int[start[numrowsconv]];
318   double * elem = new double[start[numrowsconv]];
319   for(int i = 0 ; i < numrowsconv ; i++) {
320
321     OsiRowCut * cut = cs.rowCutPtr (i);
322
323     const CoinPackedVector &v = cut->row();
324
325     if(v.getNumElements() != length[i])
326       std::cout<<"Empty row"<<std::endl;
327     //     cut->print();
328     CoinCopyN (v.getIndices(),  length[i], ind  + start[i]);
329     CoinCopyN (v.getElements(), length[i], elem + start[i]);
330   }
331
332   // Ok everything done now create interface
333   CoinPackedMatrix A;
334   A.assignMatrix(false, numcolsconv, numrowsconv,
335                  start[numrowsconv], elem, ind,
336                  start, length);
337   if(A.getNumCols() != numcolsconv || A.getNumRows() != numrowsconv){
338     std::cout<<"Error in row number"<<std::endl;
339   }
340   assert(A.getNumElements() == nnz);
341   // Objective function
342   double * obj = new double[numcolsconv];
343   CoinFillN(obj,numcolsconv,0.);
344
345   // some instances have no (or null) objective function, check it here
346   if (p -> nObjs () > 0)
347     p -> fillObjCoeff (obj);
348
349   // Finally, load interface si with the initial LP relaxation
350   si.loadProblem (A, colLower, colUpper, obj, rowLower, rowUpper);
351
352   delete [] rowLower; 
353   delete [] rowUpper;
354   delete [] colLower;
355   delete [] colUpper;
356   delete [] obj;
357
358   for (int i=0; i<numcolsconv; i++)
359     if ((p -> Var (i) -> Multiplicity () > 0) &&
360         (p -> Var (i) -> isDefinedInteger ()))
361       si.setInteger (i);
362 
363   //si.writeMpsNative("toto",NULL,NULL,1);
364   //si.writeLp ("toto");
365   app_ -> enableWarmStart();
366
367   // restored check. With "TOO FEW DEGREES OF FREEDOM" exception, x_sol() is null
368   if (problem () -> x_sol ()) {
369     setColSolution (problem () -> x_sol     ());
370     setRowPrice    (problem () -> duals_sol ());
371   }
372}
373
374
375/** To set some application specific defaults. */
376void CouenneInterface::setAppDefaultOptions(Ipopt::SmartPtr<Ipopt::OptionsList> Options){
377  Options->SetStringValue("bonmin.algorithm", "B-Couenne", true, true);
378  Options->SetIntegerValue("bonmin.filmint_ecp_cuts", 1, true, true);
379}
380} /** End Bonmin namespace. */
Note: See TracBrowser for help on using the repository browser.