source: stable/0.3/Couenne/src/problem/problem.cpp @ 580

Last change on this file since 580 was 580, checked in by pbelotti, 9 years ago

merging change in cutoff saving from trunk

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