source: stable/0.2/Couenne/src/main/BonCouenneInterface.cpp @ 159

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

created new stable branch 0.2 from trunk (rev. 157)

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