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

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

increase nOrigVars_ in CouenneProblem::addVariables (cleaner)

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