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

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

giving Couenne its own name...

File size: 10.6 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, 0.);
124          setColUpper (i, 0.);
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          p -> checkNLP (solution, obj, true) && // true for recomputing obj
231          (obj < p -> getCutOff ())) {
232
233        // tell caller there is an initial solution to be fed to the initHeuristic
234        have_nlp_solution_ = true;
235
236        // set cutoff to take advantage of bound tightening
237        p -> setCutOff (obj);
238
239        OsiAuxInfo * auxInfo = si.getAuxiliaryInfo ();
240        BabInfo * babInfo = dynamic_cast <BabInfo *> (auxInfo);
241
242        if (babInfo) {
243          babInfo -> setNlpSolution (solution, getNumCols (), obj);
244          babInfo -> setHasNlpSolution (true);
245        }
246      }
247    }
248  }
249
250  if (!is_feasible) // nothing else to do, problem infeasible
251    return;
252
253  int 
254    numcols     = getNumCols (), // # original               variables
255    numcolsconv = p -> nVars (); // # original + # auxiliary variables
256
257  const double
258    *lb = getColLower (),
259    *ub = getColUpper ();
260
261   // add original and auxiliary variables to the new problem
262   for (int i=0; i<numcols; i++) 
263     if (p -> Var (i) -> Multiplicity () > 0) si.addCol (0, NULL,NULL, lb [i],       ub [i],      0);
264     else                                     si.addCol (0, NULL,NULL, 0.,           0.,          0);
265   for (int i=numcols; i<numcolsconv; i++)    si.addCol (0, NULL,NULL, -COIN_DBL_MAX,COIN_DBL_MAX,0);
266
267   // get initial relaxation
268   OsiCuts cs;
269   couenneCg.generateCuts (si, cs);
270
271   // store all (original + auxiliary) bounds in the relaxation
272   CouNumber * colLower = new CouNumber [numcolsconv];
273   CouNumber * colUpper = new CouNumber [numcolsconv];
274
275   CouNumber *ll = p -> Lb ();
276   CouNumber *uu = p -> Ub ();
277
278   // overwrite original bounds, could have improved within generateCuts
279   for (int i = numcolsconv; i--;) 
280     if (p -> Var (i) -> Multiplicity () > 0) {
281       colLower [i] = (ll [i] > - COUENNE_INFINITY) ? ll [i] : -COIN_DBL_MAX;
282       colUpper [i] = (uu [i] <   COUENNE_INFINITY) ? uu [i] :  COIN_DBL_MAX;
283     } else {
284       colLower [i] = 0.;
285       colUpper [i] = 0.;
286     }
287
288   int numrowsconv = cs.sizeRowCuts ();
289
290   // create matrix and other stuff
291   CoinBigIndex * start = new CoinBigIndex [numrowsconv + 1];
292
293   int    * length   = new int    [numrowsconv];
294   double * rowLower = new double [numrowsconv];
295   double * rowUpper = new double [numrowsconv];
296
297   start[0] = 0;
298   int nnz = 0;
299   /* fill the four arrays. */
300   for(int i = 0 ; i < numrowsconv ; i++)
301   {
302     OsiRowCut * cut = cs.rowCutPtr (i);
303
304     const CoinPackedVector &v = cut->row();
305     start[i+1] = start[i] + v.getNumElements();
306     nnz += v.getNumElements();
307     length[i] = v.getNumElements();
308
309     rowLower[i] = cut->lb();
310     rowUpper[i] = cut->ub();
311   }
312
313   assert (nnz == start [numrowsconv]);
314   /* Now fill the elements arrays. */
315   int * ind = new int[start[numrowsconv]];
316   double * elem = new double[start[numrowsconv]];
317   for(int i = 0 ; i < numrowsconv ; i++) {
318
319     OsiRowCut * cut = cs.rowCutPtr (i);
320
321     const CoinPackedVector &v = cut->row();
322
323     if(v.getNumElements() != length[i])
324       std::cout<<"Empty row"<<std::endl;
325     //     cut->print();
326     CoinCopyN (v.getIndices(),  length[i], ind  + start[i]);
327     CoinCopyN (v.getElements(), length[i], elem + start[i]);
328   }
329
330   // Ok everything done now create interface
331   CoinPackedMatrix A;
332   A.assignMatrix(false, numcolsconv, numrowsconv,
333                  start[numrowsconv], elem, ind,
334                  start, length);
335   if(A.getNumCols() != numcolsconv || A.getNumRows() != numrowsconv){
336     std::cout<<"Error in row number"<<std::endl;
337   }
338   assert(A.getNumElements() == nnz);
339   // Objective function
340   double * obj = new double[numcolsconv];
341   CoinFillN(obj,numcolsconv,0.);
342
343   // some instances have no (or null) objective function, check it here
344   if (p -> nObjs () > 0)
345     p -> fillObjCoeff (obj);
346
347   // Finally, load interface si with the initial LP relaxation
348   si.loadProblem (A, colLower, colUpper, obj, rowLower, rowUpper);
349
350   delete [] rowLower; 
351   delete [] rowUpper;
352   delete [] colLower;
353   delete [] colUpper;
354   delete [] obj;
355
356   for (int i=0; i<numcolsconv; i++)
357     if ((p -> Var (i) -> Multiplicity () > 0) &&
358         (p -> Var (i) -> isDefinedInteger ()))
359       si.setInteger (i);
360 
361   //si.writeMpsNative("toto",NULL,NULL,1);
362   //si.writeLp ("toto");
363   app_ -> enableWarmStart();
364
365   // restored check. With "TOO FEW DEGREES OF FREEDOM" exception, x_sol() is null
366   if (problem () -> x_sol ()) {
367     setColSolution (problem () -> x_sol     ());
368     setRowPrice    (problem () -> duals_sol ());
369   }
370}
371
372
373/** To set some application specific defaults. */
374void CouenneInterface::setAppDefaultOptions(Ipopt::SmartPtr<Ipopt::OptionsList> Options){
375  Options->SetStringValue("bonmin.algorithm", "B-Couenne", true, true);
376  Options->SetIntegerValue("bonmin.filmint_ecp_cuts", 1, true, true);
377}
378} /** End Bonmin namespace. */
Note: See TracBrowser for help on using the repository browser.