# source:stable/0.2/Couenne/src/problem/problem.cpp@159

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

created new stable branch 0.2 from trunk (rev. 157)

File size: 12.2 KB
Line
1/* \$Id: problem.cpp 154 2009-06-16 18:52:53Z 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.
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_  -> Body ();
182  //int sense = objectives_  -> 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  -> Index ()] =  1; //(sense == MINIMIZE) ?  1 : -1;
195    obj [arglist  -> 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  -> Index ();
237
238        if (index >= 0) obj [index]                      = mulArgList  -> Value ();
239        else            obj [mulArgList  -> Index ()] = mulArgList  -> 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_  -> 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_  -> 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
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
326void CouenneProblem::registerOptions (Ipopt::SmartPtr <Bonmin::RegisteredOptions> roptions) {
327
328  roptions -> SetRegisteringCategory ("Couenne options", Bonmin::RegisteredOptions::CouenneCategory);
329
331    ("art_cutoff",
332     "Artificial cutoff",
333     COIN_DBL_MAX,
334     "Default value is infinity.");
335
337    ("opt_window",
338     "Window around known optimum",
339     COIN_DBL_MAX,
340     "Default value is infinity.");
341
343    ("feas_tolerance",
344     "Tolerance for constraints/auxiliary variables",
345     feas_tolerance_default,
346     "Default value is zero.");
347
349    ("feasibility_bt",
350     "Feasibility-based (cheap) bound tightening",
351     "yes",
352     "no","",
353     "yes","");
354
356    ("redcost_bt",
357     "Reduced cost bound tightening",
358     "yes",
359     "no","",
360     "yes","");
361
365     "no",
366     "no","Use an auxiliary for each bilinear term",
367     "yes","Create one only auxiliary for a quadrati expression");
368
370    ("optimality_bt",
371     "Optimality-based (expensive) bound tightening",
372     "yes",
373     "no","",
374     "yes","");
375
377    ("log_num_obbt_per_level",
378     "Specify the frequency (in terms of nodes) for optimality-based bound tightening.",
379     -1,1,
380     "\
381If -1, apply at every node (expensive!). \
382If 0, apply at root node only. \
383If k>=0, apply with probability 2^(k - level), level being the current depth of the B&B tree.");
384
386    ("aggressive_fbbt",
387     "Aggressive feasibility-based bound tightening (to use with NLP points)",
388     "yes",
389     "no","",
390     "yes","");
391
393    ("log_num_abt_per_level",
394     "Specify the frequency (in terms of nodes) for aggressive bound tightening.",
395     -1,2,
396     "\
397If -1, apply at every node (expensive!). \
398If 0, apply at root node only. \
399If k>=0, apply with probability 2^(k - level), level being the current depth of the B&B tree.");
400
402    ("art_lower",
403     "Artificial lower bound",
404     -COIN_DBL_MAX,
405     "Default value is -COIN_DBL_MAX.");
406
408    ("branching_object",
409     "type of branching object for variable selection",
410     "var_obj",
411     "vt_obj",   "use Violation Transfer from Tawarmalani and Sahinidis",
412     "var_obj",  "use one object for each variable",
413     "expr_obj", "use one object for each nonlinear expression");
414