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

Last change on this file since 267 was 267, checked in by stefan, 11 years ago

some fixes to options documentation; integrate some docu from manual

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