source: trunk/Couenne/src/problem/problem.cpp @ 125

Last change on this file since 125 was 125, checked in by pbelotti, 13 years ago

less tight cutoff -- avoids cutting some optima, for instance in OS unitTest (thanks to Kipp for pointing that out)

File size: 11.7 KB
Line 
1/*
2 * Name:    problem.cpp
3 * Author:  Pietro Belotti
4 * Purpose: methods of the class CouenneProblem
5 *
6 * (C) Carnegie-Mellon University, 2006-08.
7 * This file is licensed under the Common Public License (CPL)
8 */
9
10#include <vector>
11
12#include "CoinHelperFunctions.hpp"
13#include "CoinTime.hpp"
14
15#include "CouenneTypes.hpp"
16
17#include "expression.hpp"
18#include "exprConst.hpp"
19#include "exprGroup.hpp"
20#include "exprClone.hpp"
21#include "exprAux.hpp"
22#include "lqelems.hpp"
23#include "CouenneProblem.hpp"
24#include "CouenneProblemElem.hpp"
25#include "lqelems.hpp"
26
27// tricky... smaller values cut the optimum in OS unitTest
28const CouNumber SafeCutoff = 1e-4; 
29
30/// initialize auxiliary variables from original variables in the
31/// nonlinear problem
32
33void CouenneProblem::initAuxs () const {
34
35  domain_.current () -> resize (nVars ());
36
37  // initially, auxiliary variables are unbounded, their bounds only
38  // depending on their function
39
40  int nvars = nVars ();
41
42  for (int i=0; i < nvars; i++) {
43
44    int indvar = variables_ [i] -> Index ();
45
46    if ((variables_ [i] -> Type () == AUX) &&                   // this is an auxiliary
47        (indvar >= nOrigVars_) || // and not an original, originally
48        (variables_ [i] -> Multiplicity () == 0))               // or a useless one
49      //int index = variables_ [i] -> Index ();
50      Lb (indvar) = - (Ub (indvar) = COIN_DBL_MAX);
51  }
52
53  // first initialize with values from constraints
54
55  //Jnlst()->Printf(Ipopt::J_VECTOR, J_PROBLEM, "Initial bounds for aux (initAuxs):\n");
56
57  for (std::vector <CouenneConstraint *>::const_iterator con = constraints_.begin ();
58       con != constraints_.end (); ++con) {
59
60    CouNumber
61      lb = (*((*con) -> Lb ())) (),
62      ub = (*((*con) -> Ub ())) ();
63
64    int index = (*con) -> Body () -> Index ();
65
66    assert (index >= 0);
67
68    if ((Lb (index) = CoinMax (Lb (index), lb)) <= -COUENNE_INFINITY) Lb (index) = -COIN_DBL_MAX;
69    if ((Ub (index) = CoinMin (Ub (index), ub)) >=  COUENNE_INFINITY) Ub (index) =  COIN_DBL_MAX;
70  }
71
72  // only one loop is sufficient here, since auxiliary variable are
73  // defined in such a way that w_i does NOT depend on w_j if i<j.
74
75  Jnlst () -> Printf (Ipopt::J_MOREMATRIX, J_PROBLEM, "InitAux -- assigning bounds\n");
76
77  for (int j=0, i=nVars (); i--; j++) {
78
79    int ord = numbering_ [j];
80
81    // ignore these variables!
82    if (variables_ [ord] -> Multiplicity () == 0) {
83      Lb (ord) = - (Ub (ord) = COIN_DBL_MAX);
84      X (ord) = 0.;
85      continue;
86    }
87
88    // and handle only those with nonzero multiplicity
89    if (variables_ [ord] -> Type () == AUX) {
90
91      Jnlst () -> Printf (Ipopt::J_MOREMATRIX, J_PROBLEM, 
92                          "w_%04d [%10g,%10g] ", ord, Lb (ord), Ub (ord));
93
94      CouNumber l, u;
95
96      variables_ [ord] -> Image () -> getBounds (l, u);
97
98      /*printf ("printing bounds: [%g %g]\n", Lb (ord), Ub (ord));
99      variables_ [ord] -> Lb () -> print (); printf ("\n");
100      variables_ [ord] -> Ub () -> print (); printf ("\n");*/
101
102      Jnlst () -> Printf (Ipopt::J_MOREMATRIX, J_PROBLEM, 
103                          " [ --> w_%04d [%10g,%10g] ] vs [%10g %10g]", 
104                          ord, l, u, Lb (ord), Ub (ord));
105
106      // set bounds
107      if ((Lb (ord) = CoinMax (Lb (ord), l)) <= -COUENNE_INFINITY) Lb (ord) = -COIN_DBL_MAX;
108      if ((Ub (ord) = CoinMin (Ub (ord), u)) >=  COUENNE_INFINITY) Ub (ord) =  COIN_DBL_MAX;
109      //if ((lb_ [ord] = (*(aux -> Lb ())) ()) <= -COUENNE_INFINITY) lb_ [ord] = -DBL_MAX;
110      //if ((ub_ [ord] = (*(aux -> Ub ())) ()) >=  COUENNE_INFINITY) ub_ [ord] =  DBL_MAX;
111
112      Jnlst () -> Printf (Ipopt::J_MOREMATRIX, J_PROBLEM, 
113                          " --> [%10g,%10g]\n", Lb (ord), Ub (ord));
114
115      bool integer = variables_ [ord] -> isInteger ();
116
117      if (integer) {
118        Lb (ord) = ceil  (Lb (ord) - COUENNE_EPS);
119        Ub (ord) = floor (Ub (ord) + COUENNE_EPS);
120      }
121
122      X (ord) = CoinMax (Lb (ord), CoinMin (Ub (ord), (*(variables_ [ord] -> Image ())) ()));
123    }
124  }
125
126  restoreUnusedOriginals (); // no argument as the local x vector will be used
127}
128
129
130/// get auxiliary variables from original variables in the nonlinear
131/// problem
132void CouenneProblem::getAuxs (CouNumber * x) const {
133
134  // set point at x, don't copy
135  domain_.push (nVars (), x, domain_.lb (), domain_.ub (), false);
136
137  // set auxiliary w to f(x). This procedure is exact even though the
138  // auxiliary variables have an incomplete image, i.e. they have been
139  // decomposed previously, since they are updated with increasing
140  // index.
141
142  for (int j=0, i=nVars (); i--; j++) {
143
144    int index = numbering_ [j];
145    exprVar *var = variables_ [index];
146
147    if (var -> Multiplicity () > 0) {
148
149      CouNumber l, u;
150
151      if (var -> Type () == AUX)
152        var -> Image () -> getBounds (l,u);
153      else {
154        l = Lb (index);
155        u = Ub (index);
156      }
157
158      if (var -> Type () == AUX) {
159        X (index) =  // addresses of x[] and X() are equal
160          CoinMax (l, CoinMin (u, (*(var -> Image ())) ())); 
161      }
162    } else X (index) = 0.;
163  }
164
165  restoreUnusedOriginals ();
166
167  domain_.pop ();
168}
169
170
171/// fill obj vector with coefficient of the (linearized) obj function
172/// (depends on sense of optimization -- invert if sense()==MAXIMIZE)
173
174void CouenneProblem::fillObjCoeff (double *&obj) {
175
176  // linearized objective can be an exprAux, an exprSub, an exprGroup,
177  // or an exprSum. In the last two cases, the components are
178  // variables or constants
179
180  expression *body = objectives_ [0] -> Body ();
181  //int sense = objectives_ [0] -> Sense ();
182
183  switch (body -> code ()) {
184
185  case COU_EXPRVAR:   //
186    obj [body -> Index ()] = 1; //(sense == MINIMIZE) ? 1 : -1;
187    break;
188
189  case COU_EXPRSUB: { //
190
191    expression **arglist = body -> ArgList ();
192
193    obj [arglist [0] -> Index ()] =  1; //(sense == MINIMIZE) ?  1 : -1;
194    obj [arglist [1] -> Index ()] = -1; //(sense == MINIMIZE) ? -1 :  1;
195
196  } break;
197
198  case COU_EXPRGROUP: { //
199
200    exprGroup *eg    = dynamic_cast <exprGroup *> (body -> isaCopy () ? 
201                                                   body -> Copy () :
202                                                   body);
203
204    const exprGroup::lincoeff &lcoe = eg -> lcoeff ();
205
206    //    if (sense == MINIMIZE) while (*index >= 0) obj [*index++] =  *coeff++;
207    //    else                   while (*index >= 0) obj [*index++] = -*coeff++;     
208
209    for (int n = lcoe.size (), i=0; n--; i++)
210      //exprGroup::lincoeff::iterator el = lcoe.begin (); el != lcoe.end (); ++el)
211      obj [lcoe [i]. first -> Index ()] = lcoe [i]. second;
212    //(sense == MINIMIZE) ?
213    //(lcoe [i]. second) :
214    //-(lcoe [i]. second);
215
216  } // no break, as exprGroup is derived from exprSum
217
218  case COU_EXPRSUM: { //
219
220    expression **arglist = body -> ArgList ();
221
222    for (int i = body -> nArgs (); i--;)
223      switch ((arglist [i]) -> code ()) {
224
225      case COU_EXPRCONST: 
226        break;
227
228      case COU_EXPRVAR: 
229        obj [arglist [i] -> Index ()] = 1; //(sense == MINIMIZE) ? 1 : -1;
230        break;
231
232      case COU_EXPRMUL: {
233
234        expression **mulArgList = arglist [i] -> ArgList ();
235        int index = mulArgList [0] -> Index ();
236
237        if (index >= 0) obj [index]                      = mulArgList [1] -> Value ();
238        else            obj [mulArgList [1] -> Index ()] = mulArgList [0] -> Value ();
239      } break;
240
241      default: 
242        Jnlst()->Printf(Ipopt::J_ERROR, J_PROBLEM,
243                        "Couenne: invalid element of sum\nAborting\n");
244        exit (-1);
245      }
246  } break;
247
248  case COU_EXPRCONST: break; // a constant objective
249
250  default:
251    Jnlst()->Printf(Ipopt::J_WARNING, J_PROBLEM,
252                    "Couenne: warning, objective function not recognized\n");
253    break;
254  }
255}
256
257
258/// set cutoff from NLP solution
259void CouenneProblem::setCutOff (CouNumber cutoff) const {
260
261  int indobj = objectives_ [0] -> Body () -> Index ();
262
263  // AW: Should we use the value of the objective variable computed by
264  //     Couenne here?
265  if ((indobj >= 0) && (cutoff < pcutoff_ -> getCutOff () - COUENNE_EPS)) {
266
267    Jnlst () -> Printf (Ipopt::J_DETAILED, J_PROBLEM,
268                        "Setting new cutoff %.10e for optimization variable index %d val = %.10e\n",
269                        cutoff, indobj,
270                        pcutoff_ -> getCutOff ());
271
272    if (Var (indobj) -> isInteger ())
273      pcutoff_    -> setCutOff (floor (cutoff + COUENNE_EPS));
274    else pcutoff_ -> setCutOff (cutoff + SafeCutoff * (1. + fabs(cutoff)));
275  }
276} // tolerance needed to retain feasibility
277
278
279/// Tell problem that auxiliary related to obj has a cutoff, to be
280/// used in bound tightening
281void CouenneProblem::installCutOff () const {
282
283  int indobj = objectives_ [0] -> Body () -> Index ();
284
285  if (indobj >= 0) {
286
287    // all problem are assumed to be minimization
288    double cutoff = pcutoff_ -> getCutOff();
289
290    if (cutoff < Ub (indobj))
291      Ub (indobj) = cutoff;
292
293  } else jnlst_ -> Printf (J_SUMMARY, J_PROBLEM, 
294                           "Warning, could not install cutoff - negative objective index\n");
295}
296
297
298// clear all spurious variables pointers not referring to the variables_ vector
299void CouenneProblem::realign () {
300
301  // link variables to problem's domain
302  for (std::vector <exprVar *>::iterator i = variables_.begin ();
303       i != variables_.end (); ++i) {
304
305    (*i) -> linkDomain (&domain_);
306    (*i) -> realign (this);
307    if ((*i) -> Type () == AUX)
308      (*i) -> Image () -> realign (this);
309  }
310
311  // link variables to problem's domain
312  for (std::vector <CouenneObjective *>::iterator i = objectives_.begin ();
313       i != objectives_.end (); ++i) 
314    (*i) -> Body () -> realign (this);
315
316
317  // link variables to problem's domain
318  for (std::vector <CouenneConstraint *>::iterator i = constraints_.begin ();
319       i != constraints_.end (); ++i)
320    (*i) -> Body () -> realign (this);
321}
322
323
324/// Add list of options to be read from file
325void CouenneProblem::registerOptions (Ipopt::SmartPtr <Bonmin::RegisteredOptions> roptions) {
326
327  roptions -> SetRegisteringCategory ("Couenne options", Bonmin::RegisteredOptions::CouenneCategory);
328
329  roptions -> AddNumberOption
330    ("art_cutoff",
331     "Artificial cutoff",
332     COIN_DBL_MAX,
333     "Default value is infinity.");
334
335  roptions -> AddNumberOption
336    ("opt_window",
337     "Window around known optimum",
338     COIN_DBL_MAX,
339     "Default value is infinity.");
340
341  roptions -> AddNumberOption
342    ("feas_tolerance",
343     "Tolerance for constraints/auxiliary variables",
344     feas_tolerance_default,
345     "Default value is zero.");
346
347  roptions -> AddStringOption2
348    ("feasibility_bt",
349     "Feasibility-based (cheap) bound tightening",
350     "yes",
351     "no","",
352     "yes","");
353
354  roptions -> AddStringOption2
355    ("use_quadratic",
356     "Use quadratic expressions and related exprQuad class",
357     "no",
358     "no","Use an auxiliary for each bilinear term",
359     "yes","Create one only auxiliary for a quadrati expression");
360
361  roptions -> AddStringOption2
362    ("optimality_bt",
363     "Optimality-based (expensive) bound tightening",
364     "yes",
365     "no","",
366     "yes","");
367
368  roptions -> AddLowerBoundedIntegerOption
369    ("log_num_obbt_per_level",
370     "Specify the frequency (in terms of nodes) for optimality-based bound tightening.",
371     -1,1,
372     "\
373If -1, apply at every node (expensive!). \
374If 0, apply at root node only. \
375If k>=0, apply with probability 2^(k - level), level being the current depth of the B&B tree.");
376
377  roptions -> AddStringOption2
378    ("aggressive_fbbt",
379     "Aggressive feasibility-based bound tightening (to use with NLP points)",
380     "yes",
381     "no","",
382     "yes","");
383
384  roptions -> AddLowerBoundedIntegerOption
385    ("log_num_abt_per_level",
386     "Specify the frequency (in terms of nodes) for aggressive bound tightening.",
387     -1,2,
388     "\
389If -1, apply at every node (expensive!). \
390If 0, apply at root node only. \
391If k>=0, apply with probability 2^(k - level), level being the current depth of the B&B tree.");
392
393  roptions -> AddNumberOption
394    ("art_lower",
395     "Artificial lower bound",
396     -COIN_DBL_MAX,
397     "Default value is -COIN_DBL_MAX.");
398
399  roptions -> AddStringOption3
400    ("branching_object",
401     "type of branching object for variable selection",
402     "var_obj",
403     "vt_obj",   "use Violation Transfer from Tawarmalani and Sahinidis",
404     "var_obj",  "use one object for each variable",
405     "expr_obj", "use one object for each nonlinear expression");
406}
Note: See TracBrowser for help on using the repository browser.