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

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

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

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