source: stable/0.1/Couenne/src/main/BonCouenneInterface.cpp @ 77

Last change on this file since 77 was 77, checked in by stefan, 13 years ago

merge changeset 70 from trunk into stable: compile Couenne library without ASL

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