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 |
---|
29 | const CouNumber SafeCutoff = 1e-4; |
---|
30 | |
---|
31 | // absolute difference |
---|
32 | const CouNumber SafeDelta = 1e-2; |
---|
33 | |
---|
34 | /// initialize auxiliary variables from original variables in the |
---|
35 | /// nonlinear problem |
---|
36 | |
---|
37 | void 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 |
---|
136 | void 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 | |
---|
178 | void 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 |
---|
263 | void 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 |
---|
286 | void 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 |
---|
304 | void 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 |
---|
330 | void 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 | "\ |
---|
395 | If -1, apply at every node (expensive!). \ |
---|
396 | If 0, apply at root node only. \ |
---|
397 | If 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 | "\ |
---|
415 | If -1, apply at every node (expensive!). \ |
---|
416 | If 0, apply at root node only. \ |
---|
417 | If 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 | } |
---|