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. |
---|
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 |
---|
29 | const CouNumber SafeCutoff = 1e-4; |
---|
30 | |
---|
31 | /// initialize auxiliary variables from original variables in the |
---|
32 | /// nonlinear problem |
---|
33 | |
---|
34 | void 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 |
---|
133 | void 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 | |
---|
175 | void 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 |
---|
260 | void 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 |
---|
282 | void 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 |
---|
300 | void 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 |
---|
326 | void 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", |
---|
351 | "yes", |
---|
352 | "no","", |
---|
353 | "yes",""); |
---|
354 | |
---|
355 | roptions -> AddStringOption2 |
---|
356 | ("redcost_bt", |
---|
357 | "Reduced cost bound tightening", |
---|
358 | "yes", |
---|
359 | "no","", |
---|
360 | "yes",""); |
---|
361 | |
---|
362 | roptions -> AddStringOption2 |
---|
363 | ("use_quadratic", |
---|
364 | "Use quadratic expressions and related exprQuad class", |
---|
365 | "no", |
---|
366 | "no","Use an auxiliary for each bilinear term", |
---|
367 | "yes","Create one only auxiliary for a quadrati expression"); |
---|
368 | |
---|
369 | roptions -> AddStringOption2 |
---|
370 | ("optimality_bt", |
---|
371 | "Optimality-based (expensive) bound tightening", |
---|
372 | "yes", |
---|
373 | "no","", |
---|
374 | "yes",""); |
---|
375 | |
---|
376 | roptions -> AddLowerBoundedIntegerOption |
---|
377 | ("log_num_obbt_per_level", |
---|
378 | "Specify the frequency (in terms of nodes) for optimality-based bound tightening.", |
---|
379 | -1,1, |
---|
380 | "\ |
---|
381 | If -1, apply at every node (expensive!). \ |
---|
382 | If 0, apply at root node only. \ |
---|
383 | If k>=0, apply with probability 2^(k - level), level being the current depth of the B&B tree."); |
---|
384 | |
---|
385 | roptions -> AddStringOption2 |
---|
386 | ("aggressive_fbbt", |
---|
387 | "Aggressive feasibility-based bound tightening (to use with NLP points)", |
---|
388 | "yes", |
---|
389 | "no","", |
---|
390 | "yes",""); |
---|
391 | |
---|
392 | roptions -> AddLowerBoundedIntegerOption |
---|
393 | ("log_num_abt_per_level", |
---|
394 | "Specify the frequency (in terms of nodes) for aggressive bound tightening.", |
---|
395 | -1,2, |
---|
396 | "\ |
---|
397 | If -1, apply at every node (expensive!). \ |
---|
398 | If 0, apply at root node only. \ |
---|
399 | If k>=0, apply with probability 2^(k - level), level being the current depth of the B&B tree."); |
---|
400 | |
---|
401 | roptions -> AddNumberOption |
---|
402 | ("art_lower", |
---|
403 | "Artificial lower bound", |
---|
404 | -COIN_DBL_MAX, |
---|
405 | "Default value is -COIN_DBL_MAX."); |
---|
406 | |
---|
407 | roptions -> AddStringOption3 |
---|
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 | |
---|
415 | roptions -> AddStringOption2 |
---|
416 | ("delete_redundant", |
---|
417 | "Eliminate redundant variables, which appear in the problem as x_k = x_h", |
---|
418 | "yes", |
---|
419 | "no","Keep redundant variables, making the problem a bit larger", |
---|
420 | "yes","Eliminate redundant (the problem will be equivalent, only smaller)"); |
---|
421 | } |
---|