Scippy

SCIP

Solving Constraint Integer Programs

cons_abspower.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2020 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file cons_abspower.c
17  * @ingroup DEFPLUGINS_CONS
18  * @brief Constraint handler for absolute power constraints \f$\textrm{lhs} \leq \textrm{sign}(x+a) |x+a|^n + c z \leq \textrm{rhs}\f$
19  * @author Stefan Vigerske
20  */
21 
22 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
23 
24 #define SCIP_PRIVATE_ROWPREP
25 
26 #include <ctype.h>
27 #include "nlpi/pub_expr.h"
28 #include "nlpi/type_expr.h"
29 #include "nlpi/type_nlpi.h"
30 #include "scip/cons_abspower.h"
31 #include "scip/cons_indicator.h"
32 #include "scip/cons_linear.h"
33 #include "scip/cons_nonlinear.h"
34 #include "scip/cons_quadratic.h"
35 #include "scip/cons_varbound.h"
36 #include "scip/debug.h"
37 #include "scip/heur_subnlp.h"
38 #include "scip/heur_trysol.h"
39 #include "scip/intervalarith.h"
40 #include "scip/pub_cons.h"
41 #include "scip/pub_event.h"
42 #include "scip/pub_heur.h"
43 #include "scip/pub_lp.h"
44 #include "scip/pub_message.h"
45 #include "scip/pub_misc.h"
46 #include "scip/pub_nlp.h"
47 #include "scip/pub_sol.h"
48 #include "scip/pub_tree.h"
49 #include "scip/pub_var.h"
50 #include "scip/scip_branch.h"
51 #include "scip/scip_conflict.h"
52 #include "scip/scip_cons.h"
53 #include "scip/scip_copy.h"
54 #include "scip/scip_cut.h"
55 #include "scip/scip_event.h"
56 #include "scip/scip_general.h"
57 #include "scip/scip_heur.h"
58 #include "scip/scip_lp.h"
59 #include "scip/scip_mem.h"
60 #include "scip/scip_message.h"
61 #include "scip/scip_nlp.h"
62 #include "scip/scip_numerics.h"
63 #include "scip/scip_param.h"
64 #include "scip/scip_prob.h"
65 #include "scip/scip_probing.h"
66 #include "scip/scip_sepa.h"
67 #include "scip/scip_sol.h"
68 #include "scip/scip_tree.h"
69 #include "scip/scip_var.h"
70 #include <string.h>
71 
72 /* constraint handler properties */
73 #define CONSHDLR_NAME "abspower"
74 #define CONSHDLR_DESC "constraint handler for absolute power constraints lhs <= sign(x+offset)abs(x+offset)^n + c*z <= rhs"
75 #define CONSHDLR_SEPAPRIORITY 0 /**< priority of the constraint handler for separation */
76 #define CONSHDLR_ENFOPRIORITY -30 /**< priority of the constraint handler for constraint enforcing */
77 #define CONSHDLR_CHECKPRIORITY -3500000 /**< priority of the constraint handler for checking feasibility */
78 #define CONSHDLR_SEPAFREQ 1 /**< frequency for separating cuts; zero means to separate only in the root node */
79 #define CONSHDLR_PROPFREQ 1 /**< frequency for propagating domains; zero means only preprocessing propagation */
80 #define CONSHDLR_EAGERFREQ 100 /**< frequency for using all instead of only the useful constraints in separation,
81  * propagation and enforcement, -1 for no eager evaluations, 0 for first only */
82 #define CONSHDLR_MAXPREROUNDS -1 /**< maximal number of presolving rounds the constraint handler participates in (-1: no limit) */
83 #define CONSHDLR_DELAYSEPA FALSE /**< should separation method be delayed, if other separators found cuts? */
84 #define CONSHDLR_DELAYPROP FALSE /**< should propagation method be delayed, if other propagators found reductions? */
85 #define CONSHDLR_NEEDSCONS TRUE /**< should the constraint handler be skipped, if no constraints are available? */
86 
87 #define CONSHDLR_PRESOLTIMING SCIP_PRESOLTIMING_FAST | SCIP_PRESOLTIMING_MEDIUM
88 #define CONSHDLR_PROP_TIMING SCIP_PROPTIMING_ALWAYS /**< when should the constraint handlers propagation routines be called? */
89 
90 #define QUADCONSUPGD_PRIORITY 50000 /**< priority of the constraint handler for upgrading of quadratic constraints */
91 #define NONLINCONSUPGD_PRIORITY 50000 /**< priority of the constraint handler for upgrading of nonlinear constraints and reformulating expression graph nodes */
92 
93 /*
94  * Local defines
95  */
96 
97 #define PROPVARTOL SCIPepsilon(scip) /**< tolerance to add to variable bounds in domain propagation */
98 #define PROPSIDETOL SCIPepsilon(scip) /**< tolerance to add to constraint sides in domain propagation */
99 #define INITLPMAXVARVAL 1000.0 /**< maximal absolute value of variable for still generating a linearization cut at that point in initlp */
101 /** power function type to be used by a constraint instead of the general pow */
102 #define DECL_MYPOW(x) SCIP_Real x (SCIP_Real base, SCIP_Real exponent)
104 /** sign of a value (-1 or +1)
105  *
106  * 0.0 has sign +1
107  */
108 #define SIGN(x) ((x) >= 0.0 ? 1.0 : -1.0)
110 
111 /*
112  * Data structures
113  */
114 
115 #define ROOTS_KNOWN 10 /**< up to which (integer) exponents precomputed roots have been stored */
117 /** The positive root of the polynomial (n-1) y^n + n y^(n-1) - 1 is needed in separation.
118  * Here we store these roots for small integer values of n.
119  */
120 static
122  -1.0, /* no root for n=0 */
123  -1.0, /* no root for n=1 */
124  0.41421356237309504880, /* root for n=2 (-1+sqrt(2)) */
125  0.5, /* root for n=3 */
126  0.56042566045031785945, /* root for n=4 */
127  0.60582958618826802099, /* root for n=5 */
128  0.64146546982884663257, /* root for n=6 */
129  0.67033204760309682774, /* root for n=7 */
130  0.69428385661425826738, /* root for n=8 */
131  0.71453772716733489700, /* root for n=9 */
132  0.73192937842370733350 /* root for n=10 */
133 };
134 
135 /** constraint data for absolute power constraints */
136 struct SCIP_ConsData
137 {
138  SCIP_VAR* x; /**< variable x in sign(x+offset)|x+offset|^n term */
139  SCIP_VAR* z; /**< linear variable in constraint */
140  SCIP_Real exponent; /**< exponent n of |x+offset| */
141  SCIP_Real xoffset; /**< offset in x+offset */
142  SCIP_Real zcoef; /**< coefficient of linear variable z */
143  SCIP_Real lhs; /**< left hand side of constraint */
144  SCIP_Real rhs; /**< right hand side of constraint */
145 
146  SCIP_Real root; /**< root of polynomial */
147  DECL_MYPOW ((*power)); /**< function for computing power*/
148 
149  SCIP_Real lhsviol; /**< current violation of left hand side */
150  SCIP_Real rhsviol; /**< current violation of right hand side */
151 
152  int xeventfilterpos; /**< position of x var event in SCIP event filter */
153  int zeventfilterpos; /**< position of z var event in SCIP event filter */
154  unsigned int propvarbounds:1; /**< have variable bounds been propagated? */
155 
156  SCIP_NLROW* nlrow; /**< nonlinear row representation of constraint */
157 };
158 
159 /** constraint handler data */
160 struct SCIP_ConshdlrData
161 {
162  SCIP_Real cutmaxrange; /**< maximal coef range (maximal abs coef / minimal abs coef) of a cut in order to be added to LP */
163  SCIP_Bool projectrefpoint; /**< whether to project the reference point when linearizing a absolute power constraint in a convex region */
164  int preferzerobranch; /**< how much we prefer to branch on 0.0 first */
165  SCIP_Bool branchminconverror; /**< whether to compute branching point such that the convexification error is minimized after branching on 0.0 */
166  SCIP_Bool addvarboundcons; /**< should variable bound constraints be added? */
167  SCIP_Bool linfeasshift; /**< try linear feasibility shift heuristic in CONSCHECK */
168  SCIP_Bool dualpresolve; /**< should dual presolve be applied? */
169  SCIP_Bool sepainboundsonly; /**< should tangents only be generated in variable bounds during separation? */
170  SCIP_Real sepanlpmincont; /**< minimal required fraction of continuous variables in problem to use solution of NLP relaxation in root for separation */
171  SCIP_Bool enfocutsremovable; /**< are cuts added during enforcement removable from the LP in the same node? */
172 
173  SCIP_HEUR* subnlpheur; /**< a pointer to the subnlp heuristic */
174  SCIP_HEUR* trysolheur; /**< a pointer to the trysol heuristic */
175  SCIP_EVENTHDLR* eventhdlr; /**< our handler for bound change events on variable x */
176  SCIP_CONSHDLR* conshdlrindicator; /**< a pointer to the indicator constraint handler */
177  int newsoleventfilterpos;/**< filter position of new solution event handler, if catched */
178  SCIP_Bool comparedpairwise; /**< did we compare absolute power constraints pairwise in this run? */
179  SCIP_Bool sepanlp; /**< where linearization of the NLP relaxation solution added? */
180  SCIP_NODE* lastenfonode; /**< the node for which enforcement was called the last time (and some constraint was violated) */
181  int nenforounds; /**< counter on number of enforcement rounds for the current node */
182  unsigned int nsecantcuts; /**< number of secant cuts created so far */
183  unsigned int ncuts; /**< number of linearization cuts created so far */
184 };
185 
186 /*
187  * Propagation rules
188  */
189 
190 enum Proprule
191 {
192  PROPRULE_1, /**< left hand side and bounds on z -> lower bound on x */
193  PROPRULE_2, /**< left hand side and upper bound on x -> bound on z */
194  PROPRULE_3, /**< right hand side and bounds on z -> upper bound on x */
195  PROPRULE_4, /**< right hand side and lower bound on x -> bound on z */
196  PROPRULE_INVALID /**< propagation was applied without a specific propagation rule */
197 };
198 typedef enum Proprule PROPRULE;
200 /*
201  * Local methods
202  */
203 
204 /** computes bounds on x in a absolute power constraints for given bounds on z */
205 static
206 void computeBoundsX(
207  SCIP* scip, /**< SCIP data structure */
208  SCIP_CONS* cons, /**< constraint */
209  SCIP_INTERVAL zbnds, /**< bounds on x that are to be propagated */
210  SCIP_INTERVAL* xbnds /**< buffer to store corresponding bounds on z */
211  );
212 
213 /** power function for square, that should be faster than using pow(x, 2.0) */
214 static
216 {
217  assert(exponent == 2.0);
218  return base*base;
219 }
220 
221 /** process variable event */
222 static
223 SCIP_DECL_EVENTEXEC(processVarEvent)
224 {
225  SCIP_CONS* cons;
226 
227  assert(scip != NULL);
228  assert(event != NULL);
230 
231  cons = (SCIP_CONS*) eventdata;
232  assert(cons != NULL);
233 
234  assert(SCIPconsGetData(cons) != NULL);
235  assert(SCIPeventGetVar(event) == SCIPconsGetData(cons)->x || SCIPeventGetVar(event) == SCIPconsGetData(cons)->z);
236 
238 
239  return SCIP_OKAY;
240 } /*lint !e715*/
241 
242 /** catch variable bound tightening events */
243 static
245  SCIP* scip, /**< SCIP data structure */
246  SCIP_EVENTHDLR* eventhdlr, /**< event handler for variables */
247  SCIP_CONS* cons /**< constraint for which to catch bound change events */
248  )
249 {
250  SCIP_CONSDATA* consdata;
251  SCIP_EVENTTYPE eventtype;
252 
253  assert(scip != NULL);
254  assert(cons != NULL);
255  assert(eventhdlr != NULL);
256 
257  consdata = SCIPconsGetData(cons);
258  assert(consdata != NULL);
259 
260  /* if z is multiaggregated, then bound changes on x could not be propagated, so we do not need to catch them */
261  if( SCIPvarGetStatus(consdata->z) != SCIP_VARSTATUS_MULTAGGR )
262  {
263  eventtype = SCIP_EVENTTYPE_DISABLED;
264  if( !SCIPisInfinity(scip, -consdata->lhs) )
265  eventtype |= SCIP_EVENTTYPE_UBTIGHTENED;
266  if( !SCIPisInfinity(scip, consdata->rhs) )
267  eventtype |= SCIP_EVENTTYPE_LBTIGHTENED;
268 
269  SCIP_CALL( SCIPcatchVarEvent(scip, consdata->x, eventtype, eventhdlr, (SCIP_EVENTDATA*)cons, &consdata->xeventfilterpos) );
270 
271  SCIP_CALL( SCIPmarkConsPropagate(scip, cons) );
272  }
273 
274  /* if x is multiaggregated, then bound changes on z could not be propagated, so we do not need to catch them */
275  if( SCIPvarGetStatus(consdata->x) != SCIP_VARSTATUS_MULTAGGR )
276  {
277  eventtype = SCIP_EVENTTYPE_DISABLED;
278  if( consdata->zcoef > 0.0 )
279  {
280  if( !SCIPisInfinity(scip, -consdata->lhs) )
281  eventtype |= SCIP_EVENTTYPE_UBTIGHTENED;
282  if( !SCIPisInfinity(scip, consdata->rhs) )
283  eventtype |= SCIP_EVENTTYPE_LBTIGHTENED;
284  }
285  else
286  {
287  if( !SCIPisInfinity(scip, -consdata->lhs) )
288  eventtype |= SCIP_EVENTTYPE_LBTIGHTENED;
289  if( !SCIPisInfinity(scip, consdata->rhs) )
290  eventtype |= SCIP_EVENTTYPE_UBTIGHTENED;
291  }
292 
293  SCIP_CALL( SCIPcatchVarEvent(scip, consdata->z, eventtype, eventhdlr, (SCIP_EVENTDATA*)cons, &consdata->zeventfilterpos) );
294 
295  SCIP_CALL( SCIPmarkConsPropagate(scip, cons) );
296  }
297 
298  return SCIP_OKAY;
299 }
300 
301 /** drop variable bound tightening events */
302 static
304  SCIP* scip, /**< SCIP data structure */
305  SCIP_EVENTHDLR* eventhdlr, /**< event handler for variables */
306  SCIP_CONS* cons /**< constraint for which to drop bound change events */
307  )
308 {
309  SCIP_CONSDATA* consdata;
310  SCIP_EVENTTYPE eventtype;
311 
312  assert(scip != NULL);
313  assert(cons != NULL);
314  assert(eventhdlr != NULL);
315 
316  consdata = SCIPconsGetData(cons);
317  assert(consdata != NULL);
318 
319  if( SCIPvarGetStatus(consdata->z) != SCIP_VARSTATUS_MULTAGGR )
320  {
321  eventtype = SCIP_EVENTTYPE_DISABLED;
322  if( !SCIPisInfinity(scip, -consdata->lhs) )
323  eventtype |= SCIP_EVENTTYPE_UBTIGHTENED;
324  if( !SCIPisInfinity(scip, consdata->rhs) )
325  eventtype |= SCIP_EVENTTYPE_LBTIGHTENED;
326 
327  SCIP_CALL( SCIPdropVarEvent(scip, consdata->x, eventtype, eventhdlr, (SCIP_EVENTDATA*)cons, consdata->xeventfilterpos) );
328  consdata->xeventfilterpos = -1;
329  }
330 
331  if( SCIPvarGetStatus(consdata->x) != SCIP_VARSTATUS_MULTAGGR )
332  {
333  eventtype = SCIP_EVENTTYPE_DISABLED;
334  if( consdata->zcoef > 0.0 )
335  {
336  if( !SCIPisInfinity(scip, -consdata->lhs) )
337  eventtype |= SCIP_EVENTTYPE_UBTIGHTENED;
338  if( !SCIPisInfinity(scip, consdata->rhs) )
339  eventtype |= SCIP_EVENTTYPE_LBTIGHTENED;
340  }
341  else
342  {
343  if( !SCIPisInfinity(scip, -consdata->lhs) )
344  eventtype |= SCIP_EVENTTYPE_LBTIGHTENED;
345  if( !SCIPisInfinity(scip, consdata->rhs) )
346  eventtype |= SCIP_EVENTTYPE_UBTIGHTENED;
347  }
348 
349  SCIP_CALL( SCIPdropVarEvent(scip, consdata->z, eventtype, eventhdlr, (SCIP_EVENTDATA*)cons, consdata->zeventfilterpos) );
350  consdata->zeventfilterpos = -1;
351  }
352 
353  assert(consdata->xeventfilterpos == -1);
354  assert(consdata->zeventfilterpos == -1);
355 
356  return SCIP_OKAY;
357 }
358 
359 /** get key of hash element */
360 static
361 SCIP_DECL_HASHGETKEY(presolveFindDuplicatesGetKey)
362 {
363  return elem;
364 } /*lint !e715*/
365 
366 /** checks if two constraints have the same x variable, the same exponent, and either the same offset or the same linear variable and are both equality constraint */
367 static
368 SCIP_DECL_HASHKEYEQ(presolveFindDuplicatesKeyEQ)
369 {
370  SCIP_CONSDATA* consdata1;
371  SCIP_CONSDATA* consdata2;
372 
373  consdata1 = SCIPconsGetData((SCIP_CONS*)key1);
374  consdata2 = SCIPconsGetData((SCIP_CONS*)key2);
375  assert(consdata1 != NULL);
376  assert(consdata2 != NULL);
377 
378  if( consdata1->x != consdata2->x )
379  return FALSE;
380 
381  if( consdata1->exponent != consdata2->exponent ) /*lint !e777*/
382  return FALSE;
383 
384  if( consdata1->xoffset != consdata2->xoffset && consdata1->z != consdata2->z ) /*lint !e777*/
385  return FALSE;
386 
387  return TRUE;
388 } /*lint !e715*/
389 
390 /** get value of hash element when comparing on x */
391 static
392 SCIP_DECL_HASHKEYVAL(presolveFindDuplicatesKeyVal)
393 {
394  SCIP_CONSDATA* consdata;
395 
396  consdata = SCIPconsGetData((SCIP_CONS*)key);
397  assert(consdata != NULL);
398 
399  return SCIPhashTwo(SCIPvarGetIndex(consdata->x),
400  SCIPrealHashCode(consdata->exponent));
401 } /*lint !e715*/
402 
403 /** checks if two constraints have the same z variable and the same exponent */
404 static
405 SCIP_DECL_HASHKEYEQ(presolveFindDuplicatesKeyEQ2)
406 {
407  SCIP_CONSDATA* consdata1;
408  SCIP_CONSDATA* consdata2;
409 
410  consdata1 = SCIPconsGetData((SCIP_CONS*)key1);
411  consdata2 = SCIPconsGetData((SCIP_CONS*)key2);
412  assert(consdata1 != NULL);
413  assert(consdata2 != NULL);
414 
415  if( consdata1->z != consdata2->z )
416  return FALSE;
417 
418  if( consdata1->exponent != consdata2->exponent ) /*lint !e777*/
419  return FALSE;
420 
421  return TRUE;
422 } /*lint !e715*/
423 
424 /** get value of hash element when comparing on z */
425 static
426 SCIP_DECL_HASHKEYVAL(presolveFindDuplicatesKeyVal2)
427 {
428  SCIP_CONSDATA* consdata;
429 
430  consdata = SCIPconsGetData((SCIP_CONS*)key);
431  assert(consdata != NULL);
432 
433  return SCIPhashTwo(SCIPvarGetIndex(consdata->z),
434  SCIPrealHashCode(consdata->exponent));
435 } /*lint !e715*/
436 
437 /** upgrades a signpower constraint to a linear constraint if a second signpower constraint with same nonlinear term is available */
438 static
440  SCIP* scip, /**< SCIP data structure */
441  SCIP_CONS* cons1, /**< constraint to upgrade to a linear constraint */
442  SCIP_CONS* cons2, /**< constraint which defines a relation for x|x|^{n-1} */
443  SCIP_Bool* infeas, /**< buffer where to indicate if infeasibility has been detected */
444  int* nupgdconss, /**< buffer where to add number of upgraded conss */
445  int* ndelconss, /**< buffer where to add number of deleted conss */
446  int* naggrvars /**< buffer where to add number of aggregated variables */
447  )
448 {
449  SCIP_CONSDATA* consdata1;
450  SCIP_CONSDATA* consdata2;
451  SCIP_CONS* lincons;
452  SCIP_Real lhs;
453  SCIP_Real rhs;
454  SCIP_VAR* vars[2];
455  SCIP_Real coefs[2];
456 
457  assert(scip != NULL);
458  assert(cons1 != NULL);
459  assert(cons2 != NULL);
460  assert(infeas != NULL);
461  assert(nupgdconss != NULL);
462  assert(ndelconss != NULL);
463  assert(naggrvars != NULL);
464 
465  consdata1 = SCIPconsGetData(cons1);
466  consdata2 = SCIPconsGetData(cons2);
467  assert(consdata1 != NULL);
468  assert(consdata2 != NULL);
469 
470  assert(SCIPisEQ(scip, consdata2->lhs, consdata2->rhs));
471  assert(!SCIPisInfinity(scip, consdata2->lhs));
472  assert(consdata1->x == consdata2->x);
473  assert(consdata1->exponent == consdata2->exponent); /*lint !e777*/
474  assert(consdata1->xoffset == consdata2->xoffset); /*lint !e777*/
475 
476  lhs = consdata1->lhs;
477  if( !SCIPisInfinity(scip, -lhs) )
478  lhs -= consdata2->lhs;
479  rhs = consdata1->rhs;
480  if( !SCIPisInfinity(scip, rhs) )
481  rhs -= consdata2->lhs;
482 
483  vars[0] = consdata1->z;
484  vars[1] = consdata2->z;
485 
486  coefs[0] = consdata1->zcoef;
487  coefs[1] = -consdata2->zcoef;
488 
489  if( SCIPisEQ(scip, lhs, rhs) )
490  {
491  SCIP_Bool redundant;
492  SCIP_Bool aggregated;
493 
494  /* try aggregation */
495  SCIP_CALL( SCIPaggregateVars(scip, consdata1->z, consdata2->z, consdata1->zcoef, -consdata2->zcoef, rhs, infeas, &redundant, &aggregated) );
496 
497  /* if infeasibility has been detected, stop here */
498  if( *infeas )
499  return SCIP_OKAY;
500 
501  if( redundant )
502  {
503  /* if redundant is TRUE, then either the aggregation has been done, or it was redundant */
504  if( aggregated )
505  ++*naggrvars;
506 
507  ++*ndelconss;
508 
509  SCIP_CALL( SCIPdelCons(scip, cons1) );
510  return SCIP_OKAY;
511  }
512  }
513 
514  /* if aggregation did not succeed, then either because some variable is multi-aggregated or due to numerics or because lhs != rhs
515  * we then add a linear constraint instead
516  */
517  vars[0] = consdata1->z;
518  vars[1] = consdata2->z;
519  coefs[0] = consdata1->zcoef;
520  coefs[1] = -consdata2->zcoef;
521 
522  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, SCIPconsGetName(cons1), 2, vars, coefs, lhs, rhs,
526  SCIPconsIsStickingAtNode(cons1)) );
527  SCIP_CALL( SCIPaddCons(scip, lincons) );
528  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
529 
530  SCIP_CALL( SCIPdelCons(scip, cons1) );
531  ++*nupgdconss;
532 
533  return SCIP_OKAY;
534 }
535 
536 /** solves a system of two absolute power equations
537  * Given: (x+xoffset1)|x+xoffset1|^{exponent-1} + zcoef1 * z == rhs1
538  * and (x+xoffset2)|x+xoffset2|^{exponent-1} + zcoef2 * z == rhs2
539  * with xoffset1 != xoffset2 and zcoef1 * rhs2 == zcoef2 * rhs1 and exponent == 2,
540  * finds values for x and z that satisfy these equations, or reports infeasibility if no solution exists.
541  *
542  * Multiplying the second equation by -zcoef1/zcoef2 and adding it to the first one gives
543  * (x+xoffset1)|x+xoffset1| - zcoef1/zcoef2 (x+offset2)|x+offset2| == 0
544  *
545  * If zcoef1 == zcoef2, then there exists, due to monotonicity of x|x|, no x such that
546  * (x+xoffset1)|x+xoffset1| == (x+xoffset2)|x+xoffset2|.
547  *
548  * In general, for zcoef1 / zcoef2 > 0.0, we get
549  * x = (xoffset2 - xoffset1) / (sqrt(zcoef2 / zcoef1) - 1.0) - xoffset1,
550  * and for zcoef1 / zcoef2 < 0.0, we get
551  * x = (xoffset2 - xoffset1) / (-sqrt(-zcoef2 / zcoef1) - 1.0) - xoffset1.
552  *
553  * This then yields z = (rhs1 - (x+xoffset1)|x+xoffset1|) / zcoef1.
554  */
555 static
557  SCIP* scip, /**< SCIP data structure */
558  SCIP_Bool* infeas, /**< buffer to indicate if the system of equations has no solution */
559  SCIP_Real* xval, /**< buffer to store value of x in the solution, if any */
560  SCIP_Real* zval, /**< buffer to store value of z in the solution, if any */
561  SCIP_Real exponent, /**< exponent in absolute power equations */
562  SCIP_Real xoffset1, /**< offset for x in first absolute power equation */
563  SCIP_Real zcoef1, /**< coefficient of z in first absolute power equation */
564  SCIP_Real rhs1, /**< right-hand-side in first absolute power equation */
565  SCIP_Real xoffset2, /**< offset for x in second absolute power equation */
566  SCIP_Real zcoef2, /**< coefficient of z in second absolute power equation */
567  SCIP_Real rhs2 /**< right-hand-side in second absolute power equation */
568  )
569 {
570  assert(scip != NULL);
571  assert(infeas != NULL);
572  assert(xval != NULL);
573  assert(zval != NULL);
574  assert(exponent == 2.0);
575  assert(!SCIPisEQ(scip, xoffset1, xoffset2));
576  assert(SCIPisEQ(scip, zcoef1 * rhs2, zcoef2 * rhs1));
577  assert(zcoef1 != 0.0);
578  assert(zcoef2 != 0.0);
579 
580  if( xoffset2 < xoffset1 )
581  {
582  presolveFindDuplicatesSolveEquations(scip, infeas, xval, zval, exponent, xoffset2, zcoef2, rhs2, xoffset1, zcoef1, rhs1);
583  return;
584  }
585 
586  if( SCIPisEQ(scip, zcoef1, zcoef2) )
587  {
588  *infeas = TRUE;
589  return;
590  }
591 
592  *infeas = FALSE;
593 
594  if( SCIPisEQ(scip, zcoef1, -zcoef2) )
595  {
596  *xval = - (xoffset1 + xoffset2) / 2.0;
597  }
598  else if( zcoef2 * zcoef1 > 0.0 )
599  {
600  *xval = (xoffset2 - xoffset1) / (sqrt(zcoef2 / zcoef1) - 1.0) - xoffset1;
601  }
602  else
603  {
604  assert(zcoef2 * zcoef1 < 0.0);
605  *xval = (xoffset2 - xoffset1) / (-sqrt(-zcoef2 / zcoef1) - 1.0) - xoffset1;
606  }
607 
608  *zval = rhs1 - (*xval + xoffset1) * REALABS(*xval + xoffset1);
609  *zval /= zcoef1;
610 
611  assert(SCIPisFeasEQ(scip, (*xval + xoffset1) * REALABS(*xval + xoffset1) + zcoef1 * *zval, rhs1));
612  assert(SCIPisFeasEQ(scip, (*xval + xoffset2) * REALABS(*xval + xoffset2) + zcoef2 * *zval, rhs2));
613 }
614 
615 /** finds and removes duplicates in a set of absolute power constraints */
616 static
618  SCIP* scip, /**< SCIP data structure */
619  SCIP_CONSHDLR* conshdlr, /**< constraint handler for absolute power constraints */
620  SCIP_CONS** conss, /**< constraints */
621  int nconss, /**< number of constraints */
622  int* nupgdconss, /**< pointer where to add number of upgraded constraints */
623  int* ndelconss, /**< pointer where to add number of deleted constraints */
624  int* naddconss, /**< pointer where to add number of added constraints */
625  int* nfixedvars, /**< pointer where to add number of fixed variables */
626  int* naggrvars, /**< pointer where to add number of aggregated variables */
627  SCIP_Bool* success, /**< pointer to store whether a duplicate was found (and removed) */
628  SCIP_Bool* infeas /**< pointer to store whether infeasibility was detected */
629  )
630 {
631  SCIP_MULTIHASH* multihash;
632  SCIP_MULTIHASHLIST* multihashlist;
633  SCIP_CONSHDLRDATA* conshdlrdata;
634  int c;
635 
636  assert(scip != NULL);
637  assert(conshdlr != NULL);
638  assert(conss != NULL || nconss == 0);
639  assert(nupgdconss != NULL);
640  assert(ndelconss != NULL);
641  assert(naddconss != NULL);
642  assert(nfixedvars != NULL);
643  assert(naggrvars != NULL);
644  assert(success != NULL);
645  assert(infeas != NULL);
646 
647  *success = FALSE;
648  *infeas = FALSE;
649 
650  if( nconss <= 1 )
651  return SCIP_OKAY;
652 
653  conshdlrdata = SCIPconshdlrGetData(conshdlr);
654  assert(conshdlrdata != NULL);
655 
656  /* check all constraints in the given set for duplicates, dominance, or possible simplifications w.r.t. the x variable */
657 
659  presolveFindDuplicatesGetKey, presolveFindDuplicatesKeyEQ, presolveFindDuplicatesKeyVal, (void*)scip) );
660 
661  for( c = 0; c < nconss && !*infeas; ++c )
662  {
663  SCIP_CONS* cons0;
664  SCIP_CONS* cons1;
665 
666  cons0 = conss[c]; /*lint !e613*/
667 
668  assert(!SCIPconsIsModifiable(cons0)); /* absolute power constraints aren't modifiable */
669  assert(!SCIPconsIsLocal(cons0)); /* shouldn't have local constraints in presolve */
670  assert(SCIPconsIsActive(cons0)); /* shouldn't get inactive constraints here */
671 
672  multihashlist = NULL;
673 
674  do
675  {
676  SCIP_CONSDATA* consdata0;
677  SCIP_CONSDATA* consdata1;
678 
679  /* get constraint from current hash table with same x variable as cons0 and same exponent */
680  cons1 = (SCIP_CONS*)(SCIPmultihashRetrieveNext(multihash, &multihashlist, (void*)cons0));
681  if( cons1 == NULL )
682  {
683  /* processed all constraints like cons0 from hash table, so insert cons0 and go to conss[c+1] */
684  SCIP_CALL( SCIPmultihashInsert(multihash, (void*) cons0) );
685  break;
686  }
687 
688  assert(cons0 != cons1);
689 
690  consdata0 = SCIPconsGetData(cons0);
691  consdata1 = SCIPconsGetData(cons1);
692  assert(consdata0 != NULL);
693  assert(consdata1 != NULL);
694 
695  SCIPdebugPrintCons(scip, cons0, NULL);
696  SCIPdebugPrintCons(scip, cons1, NULL);
697 
698  assert(consdata0->x == consdata1->x);
699  assert(consdata0->exponent == consdata1->exponent); /*lint !e777*/
700 
701  if( SCIPisEQ(scip, consdata0->xoffset, consdata1->xoffset) )
702  {
703  /* we have two constraints with the same (x+offset)|x+offset|^n term */
704 
705  /* if both constraints have the same functions; strengthen sides of cons1 and throw cons0 away */
706  if( consdata0->z == consdata1->z && SCIPisEQ(scip, consdata0->zcoef, consdata1->zcoef) )
707  {
708  /* check if side strenghtening would result in inconsistency */
709  if( SCIPisGT(scip, consdata0->lhs, consdata1->rhs) || SCIPisGT(scip, consdata1->lhs, consdata0->rhs) )
710  {
711  SCIPdebugMsg(scip, "<%s> and <%s> are contradictory; declare infeasibility\n", SCIPconsGetName(cons0), SCIPconsGetName(cons1));
712  *infeas = TRUE;
713  break;
714  }
715 
716  SCIPdebugMsg(scip, "<%s> and <%s> are equivalent; dropping the first\n", SCIPconsGetName(cons0), SCIPconsGetName(cons1));
717 
718  /* if a side of cons1 gets finite via merging with cons0, then this changes locks and events */
719  if( (SCIPisInfinity(scip, -consdata1->lhs) && !SCIPisInfinity(scip, -consdata0->lhs)) ||
720  ( SCIPisInfinity(scip, consdata1->rhs) && !SCIPisInfinity(scip, consdata0->rhs)) )
721  {
722  SCIP_CALL( dropVarEvents(scip, conshdlrdata->eventhdlr, cons1) );
723  SCIP_CALL( SCIPunlockVarCons(scip, consdata1->x, cons1, !SCIPisInfinity(scip, -consdata1->lhs), !SCIPisInfinity(scip, consdata1->rhs)) );
724  if( consdata1->zcoef > 0.0 )
725  SCIP_CALL( SCIPunlockVarCons(scip, consdata1->z, cons1, !SCIPisInfinity(scip, -consdata1->lhs), !SCIPisInfinity(scip, consdata1->rhs)) );
726  else
727  SCIP_CALL( SCIPunlockVarCons(scip, consdata1->z, cons1, !SCIPisInfinity(scip, consdata1->rhs), !SCIPisInfinity(scip, -consdata1->lhs)) );
728 
729  consdata1->lhs = MAX(consdata0->lhs, consdata1->lhs);
730  consdata1->rhs = MIN(consdata0->rhs, consdata1->rhs);
731 
732  SCIP_CALL( catchVarEvents(scip, conshdlrdata->eventhdlr, cons1) );
733  SCIP_CALL( SCIPlockVarCons(scip, consdata1->x, cons1, !SCIPisInfinity(scip, -consdata1->lhs), !SCIPisInfinity(scip, consdata1->rhs)) );
734  if( consdata1->zcoef > 0.0 )
735  SCIP_CALL( SCIPlockVarCons(scip, consdata1->z, cons1, !SCIPisInfinity(scip, -consdata1->lhs), !SCIPisInfinity(scip, consdata1->rhs)) );
736  else
737  SCIP_CALL( SCIPlockVarCons(scip, consdata1->z, cons1, !SCIPisInfinity(scip, consdata1->rhs), !SCIPisInfinity(scip, -consdata1->lhs)) );
738  }
739  else
740  {
741  consdata1->lhs = MAX(consdata0->lhs, consdata1->lhs);
742  consdata1->rhs = MIN(consdata0->rhs, consdata1->rhs);
743  }
744 
745  SCIP_CALL( SCIPdelCons(scip, cons0) );
746  ++*ndelconss;
747  *success = TRUE;
748 
749  break;
750  }
751 
752  /* if cons1 defines a linear expression for sign(x+offset)|x+offset|^n, use it to replace cons0 by a linear constraint */
753  if( SCIPisEQ(scip, consdata1->lhs, consdata1->rhs) )
754  {
755  SCIPdebugMsg(scip, "substitute <%s> in <%s> to make linear constraint\n", SCIPconsGetName(cons1), SCIPconsGetName(cons0));
756  SCIP_CALL( presolveFindDuplicatesUpgradeCons(scip, cons0, cons1, infeas, nupgdconss, ndelconss, naggrvars) );
757 
758  *success = TRUE;
759  break;
760  }
761 
762  /* if cons0 defines a linear expression for sign(x+offset)|x+offset|^n, use it to replace cons1 by a linear constraint */
763  if( SCIPisEQ(scip, consdata0->lhs, consdata0->rhs) )
764  {
765  SCIPdebugMsg(scip, "substitute <%s> in <%s> to make linear constraint\n", SCIPconsGetName(cons0), SCIPconsGetName(cons1));
766  SCIP_CALL( presolveFindDuplicatesUpgradeCons(scip, cons1, cons0, infeas, nupgdconss, ndelconss, naggrvars) );
767 
768  SCIP_CALL( SCIPmultihashRemove(multihash, cons1) );
769  *success = TRUE;
770 
771  if( *infeas )
772  break;
773  }
774  else
775  {
776  /* introduce a new equality constraint for sign(x+offset)|x+offset|^n and use it to replace cons0 and cons1 */
777  /* @todo maybe we could be more clever by looking which constraint sides are finite */
778  SCIP_VAR* auxvar;
779  SCIP_CONS* auxcons;
780  char name[SCIP_MAXSTRLEN];
781  SCIP_VAR* vars[2];
782  SCIP_Real coefs[2];
783 
784  SCIPdebugMsg(scip, "introduce new auxvar for signpower(%s+%g, %g) to make <%s> and <%s> linear constraint\n", SCIPvarGetName(consdata0->x), consdata0->exponent, consdata0->xoffset, SCIPconsGetName(cons0), SCIPconsGetName(cons1));
785 
786  /* create auxiliary variable to represent sign(x+offset)|x+offset|^n */
787  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "auxvar_abspower%s_%g_%g", SCIPvarGetName(consdata0->x), consdata0->exponent, consdata0->xoffset);
788  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS,
789  TRUE, TRUE, NULL, NULL, NULL, NULL, NULL) );
790  SCIP_CALL( SCIPaddVar(scip, auxvar) );
791 
792  /* create auxiliary constraint auxvar = sign(x+offset)|x+offset|^n
793  * as we introduced a new variable, the constraint that "defines" the value for this variable need to be enforced, that is, is not redundent
794  */
795  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "auxcons_abspower%s_%g_%g", SCIPvarGetName(consdata0->x), consdata0->exponent, consdata0->xoffset);
796  SCIP_CALL( SCIPcreateConsAbspower(scip, &auxcons, name, consdata0->x, auxvar, consdata0->exponent, consdata0->xoffset, -1.0, 0.0, 0.0,
797  SCIPconsIsInitial(cons0) || SCIPconsIsInitial(cons1),
798  SCIPconsIsSeparated(cons0) || SCIPconsIsSeparated(cons1),
799  TRUE,
800  TRUE,
802  FALSE,
803  FALSE,
804  SCIPconsIsDynamic(cons0) || SCIPconsIsDynamic(cons1),
805  SCIPconsIsRemovable(cons0) || SCIPconsIsRemovable(cons1),
807  ) );
808  SCIP_CALL( SCIPaddCons(scip, auxcons) );
809  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
810  ++*naddconss;
811 
812 #ifdef WITH_DEBUG_SOLUTION
813  if( SCIPdebugIsMainscip(scip) )
814  {
815  SCIP_Real xval;
816 
817  SCIP_CALL( SCIPdebugGetSolVal(scip, consdata0->x, &xval) );
818  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, SIGN(xval + consdata0->xoffset) * pow(REALABS(xval + consdata0->xoffset), consdata0->exponent)) );
819  }
820 #endif
821 
822  /* create linear constraint equivalent for cons0 */
823  vars[0] = auxvar;
824  vars[1] = consdata0->z;
825  coefs[0] = 1.0;
826  coefs[1] = consdata0->zcoef;
827  SCIP_CALL( SCIPcreateConsLinear(scip, &auxcons, SCIPconsGetName(cons0), 2, vars, coefs, consdata0->lhs, consdata0->rhs,
831  SCIPconsIsStickingAtNode(cons0)) );
832  SCIP_CALL( SCIPaddCons(scip, auxcons) );
833  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
834  ++*nupgdconss;
835 
836  /* create linear constraint equivalent for cons1 */
837  vars[1] = consdata1->z;
838  coefs[1] = consdata1->zcoef;
839  SCIP_CALL( SCIPcreateConsLinear(scip, &auxcons, SCIPconsGetName(cons1), 2, vars, coefs, consdata1->lhs, consdata1->rhs,
843  SCIPconsIsStickingAtNode(cons1)) );
844  SCIP_CALL( SCIPaddCons(scip, auxcons) );
845  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
846  ++*nupgdconss;
847 
848  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
849 
850  SCIP_CALL( SCIPdelCons(scip, cons0) );
851  SCIP_CALL( SCIPdelCons(scip, cons1) );
852  SCIP_CALL( SCIPmultihashRemove(multihash, cons1) );
853  *success = TRUE;
854 
855  break;
856  }
857  }
858  else if( consdata0->z == consdata1->z &&
859  consdata0->exponent == 2.0 &&
860  !SCIPisZero(scip, consdata0->zcoef) &&
861  !SCIPisZero(scip, consdata1->zcoef) &&
862  SCIPisEQ(scip, consdata0->lhs, consdata0->rhs) &&
863  SCIPisEQ(scip, consdata1->lhs, consdata1->rhs) &&
864  SCIPisEQ(scip, consdata0->lhs * consdata1->zcoef, consdata1->lhs * consdata0->zcoef) )
865  {
866  /* If we have two equality constraints with the same variables and the same exponent and compatible constants,
867  * then this system of equations should have either no or a single solution.
868  * Thus, we can report cutoff or fix the variables to this solution, and forget about the constraints.
869  * @todo think about inequalities, differing exponents, and exponents != 2
870  */
871 
872  SCIP_Real xval;
873  SCIP_Real zval;
874 
875  assert(consdata0->x == consdata1->x);
876  assert(consdata0->exponent == consdata1->exponent); /*lint !e777*/
877  assert(!SCIPisEQ(scip, consdata0->xoffset, consdata1->xoffset));
878 
879  presolveFindDuplicatesSolveEquations(scip, infeas, &xval, &zval,
880  consdata0->exponent,
881  consdata0->xoffset, consdata0->zcoef, consdata0->lhs,
882  consdata1->xoffset, consdata1->zcoef, consdata1->lhs);
883 
884  if( *infeas )
885  {
886  SCIPdebugMsg(scip, "infeasibility detected while solving the equations, no solution exists\n");
887  SCIPdebugPrintCons(scip, cons0, NULL);
888  SCIPdebugPrintCons(scip, cons1, NULL);
889  break;
890  }
891 
892  SCIPdebugMsg(scip, "fixing variables <%s>[%g, %g] to %g and <%s>[%g, %g] to %g due to equations\n",
893  SCIPvarGetName(consdata0->x), SCIPvarGetLbLocal(consdata0->x), SCIPvarGetUbLocal(consdata0->x), xval,
894  SCIPvarGetName(consdata0->z), SCIPvarGetLbLocal(consdata0->z), SCIPvarGetUbLocal(consdata0->z), zval);
895  SCIPdebugPrintCons(scip, cons0, NULL);
896  SCIPdebugPrintCons(scip, cons1, NULL);
897 
899  {
900  SCIP_Bool fixed;
901 
902  SCIP_CALL( SCIPfixVar(scip, consdata0->x, xval, infeas, &fixed) );
903  ++*ndelconss;
904 
905  if( fixed )
906  ++*nfixedvars;
907 
908  if( *infeas )
909  {
910  SCIPdebugMsg(scip, "infeasibility detected after fixing <%s>\n", SCIPvarGetName(consdata0->x));
911  break;
912  }
913  }
914  else
915  {
916  SCIP_CONS* lincons;
917  SCIP_Real one;
918 
919  one = 1.0;
920  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, SCIPconsGetName(cons0), 1, &consdata0->x, &one, xval, xval,
922  SCIP_CALL( SCIPaddCons(scip, lincons) );
923  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
924  ++*nupgdconss;
925  }
926 
928  {
929  SCIP_Bool fixed;
930 
931  SCIP_CALL( SCIPfixVar(scip, consdata0->z, zval, infeas, &fixed) );
932  ++*ndelconss;
933 
934  if( fixed )
935  ++*nfixedvars;
936 
937  if( *infeas )
938  {
939  SCIPdebugMsg(scip, "infeasibility detected after fixing <%s>\n", SCIPvarGetName(consdata0->z));
940  break;
941  }
942  }
943  else
944  {
945  SCIP_CONS* lincons;
946  SCIP_Real one;
947 
948  one = 1.0;
949  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, SCIPconsGetName(cons1), 1, &consdata0->z, &one, zval, zval,
951  SCIP_CALL( SCIPaddCons(scip, lincons) );
952  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
953  ++*nupgdconss;
954  }
955 
956  SCIP_CALL( SCIPdelCons(scip, cons0) );
957  SCIP_CALL( SCIPdelCons(scip, cons1) );
958  SCIP_CALL( SCIPmultihashRemove(multihash, cons1) );
959  *success = TRUE;
960 
961  break;
962  }
963 
964  if( multihashlist == NULL )
965  {
966  /* processed all constraints like cons0 from hash table, but cons0 could not be removed, so insert cons0 into hashmap and go to conss[c+1] */
967  SCIP_CALL( SCIPmultihashInsert(multihash, (void*) cons0) );
968  break;
969  }
970  }
971  while( TRUE ); /*lint !e506*/
972  }
973 
974  /* free hash table */
975  SCIPmultihashFree(&multihash);
976 
977  if( *infeas )
978  return SCIP_OKAY;
979 
980  /* check all constraints in the given set for duplicates, dominance, or possible simplifications w.r.t. the z variable */
981 
983  presolveFindDuplicatesGetKey, presolveFindDuplicatesKeyEQ2, presolveFindDuplicatesKeyVal2, (void*) scip) );
984 
985  for( c = 0; c < nconss && !*infeas; ++c )
986  {
987  SCIP_CONS* cons0;
988  SCIP_CONS* cons1;
989  SCIP_CONSDATA* consdata0;
990 
991  cons0 = conss[c]; /*lint !e613*/
992 
993  assert(!SCIPconsIsModifiable(cons0)); /* absolute power constraints aren't modifiable */
994  assert(!SCIPconsIsLocal(cons0)); /* shouldn't have local constraints in presolve */
995 
996  /* do not consider constraints that we have deleted in the above loop */
997  if( SCIPconsIsDeleted(cons0) )
998  continue;
999  assert(SCIPconsIsActive(cons0)); /* shouldn't get inactive constraints here */
1000 
1001  consdata0 = SCIPconsGetData(cons0);
1002  assert(consdata0 != NULL);
1003 
1004  /* consider only equality constraints so far
1005  * @todo do also something with inequalities
1006  */
1007  if( !SCIPisEQ(scip, consdata0->lhs, consdata0->rhs) )
1008  continue;
1009 
1010  multihashlist = NULL;
1011 
1012  do
1013  {
1014  SCIP_CONSDATA* consdata1;
1015 
1016  /* get constraint from current hash table with same z variable as cons0 and same exponent */
1017  cons1 = (SCIP_CONS*)(SCIPmultihashRetrieveNext(multihash, &multihashlist, (void*)cons0));
1018  if( cons1 == NULL )
1019  {
1020  /* processed all constraints like cons0 from hash table, so insert cons0 and go to conss[c+1] */
1021  SCIP_CALL( SCIPmultihashInsert(multihash, (void*) cons0) );
1022  break;
1023  }
1024 
1025  assert(cons0 != cons1);
1026  assert(!SCIPconsIsDeleted(cons1));
1027 
1028  consdata1 = SCIPconsGetData(cons1);
1029  assert(consdata1 != NULL);
1030 
1031  SCIPdebugPrintCons(scip, cons0, NULL);
1032  SCIPdebugPrintCons(scip, cons1, NULL);
1033 
1034  assert(consdata0->z == consdata1->z);
1035  assert(consdata0->exponent == consdata1->exponent); /*lint !e777*/
1036  assert(SCIPisEQ(scip, consdata1->lhs, consdata1->rhs));
1037  assert(!SCIPisZero(scip, consdata1->zcoef));
1038 
1039  if( SCIPisEQ(scip, consdata0->lhs*consdata1->zcoef, consdata1->lhs*consdata0->zcoef) )
1040  {
1041  /* have two absolute power equations with same z and compatible constants
1042  * we can then reduce this to one absolute power and one linear equation
1043  * -> x0 + xoffset0 = signpower(zcoef0/zcoef1, 1/exponent) (x1 + xoffset1)
1044  * -> keep cons1
1045  * the latter can be realized as an aggregation (if x0 and x1 are not multiaggregated) or linear constraint
1046  */
1047  SCIP_Bool redundant;
1048  SCIP_Bool aggregated;
1049  SCIP_Real coef;
1050  SCIP_Real rhs;
1051 
1052  SCIPdebugMsg(scip, "<%s> and <%s> can be reformulated to one abspower and one aggregation\n", SCIPconsGetName(cons0), SCIPconsGetName(cons1));
1053  SCIPdebugPrintCons(scip, cons0, NULL);
1054  SCIPdebugPrintCons(scip, cons1, NULL);
1055 
1056  if( consdata0->exponent == 2.0 )
1057  coef = SIGN(consdata0->zcoef / consdata1->zcoef) * sqrt(REALABS(consdata0->zcoef / consdata1->zcoef));
1058  else
1059  coef = SIGN(consdata0->zcoef / consdata1->zcoef) * pow(REALABS(consdata0->zcoef / consdata1->zcoef), 1.0/consdata0->exponent);
1060  rhs = coef * consdata1->xoffset - consdata0->xoffset;
1061 
1062  /* try aggregation */
1063  SCIP_CALL( SCIPaggregateVars(scip, consdata0->x, consdata1->x, 1.0, -coef, rhs, infeas, &redundant, &aggregated) );
1064  if( *infeas )
1065  {
1066  /* if infeasibility has been detected, stop here */
1067  break;
1068  }
1069  else if( redundant )
1070  {
1071  /* if redundant is TRUE, then either the aggregation has been done, or it was redundant */
1072  if( aggregated )
1073  ++*naggrvars;
1074 
1075  ++*ndelconss;
1076  }
1077  else
1078  {
1079  /* if aggregation did not succeed, then either because some variable is multi-aggregated or due to numerics
1080  * we then add a linear constraint instead
1081  */
1082  SCIP_CONS* auxcons;
1083  SCIP_VAR* vars[2];
1084  SCIP_Real coefs[2];
1085 
1086  vars[0] = consdata0->x;
1087  vars[1] = consdata1->x;
1088  coefs[0] = 1.0;
1089  coefs[1] = -coef;
1090 
1091  /* create linear constraint equivalent for cons0 */
1092  SCIP_CALL( SCIPcreateConsLinear(scip, &auxcons, SCIPconsGetName(cons0), 2, vars, coefs, rhs, rhs,
1096  SCIPconsIsStickingAtNode(cons0)) );
1097  SCIP_CALL( SCIPaddCons(scip, auxcons) );
1098  SCIPdebugPrintCons(scip, auxcons, NULL);
1099  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
1100 
1101  ++*nupgdconss;
1102  }
1103  SCIP_CALL( SCIPdelCons(scip, cons0) );
1104 
1105  *success = TRUE;
1106  break;
1107  }
1108 
1109  if( multihashlist == NULL )
1110  {
1111  /* processed all constraints like cons0 from hash table, but cons0 could not be removed, so insert cons0 into hashmap and go to conss[c+1] */
1112  SCIP_CALL( SCIPmultihashInsert(multihash, (void*) cons0) );
1113  break;
1114  }
1115  }
1116  while( TRUE ); /*lint !e506*/
1117  }
1118 
1119  /* free hash table */
1120  SCIPmultihashFree(&multihash);
1121 
1122  return SCIP_OKAY;
1123 }
1124 
1125 /** fix variables not appearing in any other constraint
1126  *
1127  * @todo generalize to inequalities
1128  * @todo generalize to support discrete variables
1129  * @todo generalize to arbitrary exponents also if z is in objective
1130  */
1131 static
1133  SCIP* scip, /**< SCIP data structure */
1134  SCIP_CONS* cons, /**< constraint */
1135  SCIP_Bool* cutoff, /**< buffer to indicate whether a cutoff was detected */
1136  int* ndelconss, /**< buffer to increase with the number of deleted constraint */
1137  int* nfixedvars /**< buffer to increase with the number of fixed variables */
1138  )
1139 {
1140  SCIP_CONSDATA* consdata;
1141  SCIP_Bool lhsexists;
1142  SCIP_Bool rhsexists;
1143 
1144  assert(scip != NULL);
1145  assert(cons != NULL);
1146  assert(cutoff != NULL);
1147  assert(nfixedvars != NULL);
1148  assert(ndelconss != NULL);
1149 
1150  /* only process checked constraints (for which the locks are increased);
1151  * otherwise we would have to check for variables with nlocks == 0, and these are already processed by the
1152  * dualfix presolver
1153  */
1154  if( !SCIPconsIsChecked(cons) )
1155  return SCIP_OKAY;
1156 
1157  consdata = SCIPconsGetData(cons);
1158  assert(consdata != NULL);
1159 
1160  /* skip dual presolve if multiaggregated variables are present for now (bounds are not updated, difficult to fix) */
1161  if( SCIPvarGetStatus(consdata->x) == SCIP_VARSTATUS_MULTAGGR )
1162  return SCIP_OKAY;
1163  if( SCIPvarGetStatus(consdata->z) == SCIP_VARSTATUS_MULTAGGR )
1164  return SCIP_OKAY;
1165 
1166  /* skip dual presolve if discrete variables are present for now (more difficult to compute fixing value) */
1167  if( SCIPvarGetType(consdata->x) <= SCIP_VARTYPE_INTEGER )
1168  return SCIP_OKAY;
1169  if( SCIPvarGetType(consdata->z) <= SCIP_VARTYPE_INTEGER )
1170  return SCIP_OKAY;
1171 
1172  /* we assume that domain propagation has been run and fixed variables were removed if possible */
1173  assert(!SCIPconsIsMarkedPropagate(cons));
1174  assert(consdata->zcoef != 0.0);
1175 
1176  lhsexists = !SCIPisInfinity(scip, -consdata->lhs);
1177  rhsexists = !SCIPisInfinity(scip, consdata->rhs);
1178 
1179  if( SCIPvarGetNLocksDownType(consdata->x, SCIP_LOCKTYPE_MODEL) == (lhsexists ? 1 : 0) &&
1180  SCIPvarGetNLocksUpType(consdata->x, SCIP_LOCKTYPE_MODEL) == (rhsexists ? 1 : 0) &&
1181  (consdata->zcoef > 0.0 ? SCIPvarGetNLocksDownType(consdata->z, SCIP_LOCKTYPE_MODEL) :
1182  SCIPvarGetNLocksUpType(consdata->z, SCIP_LOCKTYPE_MODEL)) == (lhsexists ? 1 : 0) &&
1183  (consdata->zcoef > 0.0 ? SCIPvarGetNLocksUpType(consdata->z, SCIP_LOCKTYPE_MODEL) :
1184  SCIPvarGetNLocksDownType(consdata->z, SCIP_LOCKTYPE_MODEL)) == (rhsexists ? 1 : 0) )
1185  {
1186  /* x and z are only locked by cons, so we can fix them to an optimal solution of
1187  * min xobj * x + zobj * z
1188  * s.t. lhs <= sign(x+offset)*abs(x+offset)^exponent + zcoef * z <= rhs
1189  * xlb <= x <= xub
1190  * zlb <= z <= zub
1191  */
1192  if( SCIPisEQ(scip, consdata->lhs, consdata->rhs) )
1193  {
1194  /* much simpler case where we can substitute z:
1195  * min xobj * x + zobj/zcoef * (rhs - sign(x+offset)*abs(x+offset)^exponent)
1196  * s.t. xlb <= x <= xub
1197  */
1198  SCIP_Real xfix;
1199  SCIP_Real xlb;
1200  SCIP_Real xub;
1201  SCIP_Real zfix;
1202  SCIP_INTERVAL xbnds;
1203  SCIP_INTERVAL zbnds;
1204  SCIP_Bool fixed;
1205 
1206  /* Since domain propagation has been applied, we would like to assume that for any valid value for x,
1207  * also the corresponding z value is valid. However, domain propagation only applies sufficiently
1208  * strong bound tightenings, so we better recompute the bounds on x.
1209  */
1210  SCIPintervalSetBounds(&zbnds, SCIPvarGetLbGlobal(consdata->z), SCIPvarGetUbGlobal(consdata->z));
1211  computeBoundsX(scip, cons, zbnds, &xbnds);
1212  xlb = MAX(SCIPvarGetLbGlobal(consdata->x), xbnds.inf); /*lint !e666*/
1213  xub = MIN(SCIPvarGetUbGlobal(consdata->x), xbnds.sup); /*lint !e666*/
1214 
1215  /* with our own "local" boundtightening, xlb might end slightly above xub,
1216  * which can result in xfix being outside bounds below, see also #2202
1217  */
1218  assert(SCIPisFeasLE(scip, xlb, xub));
1219  if( xub < xlb )
1220  xlb = xub = (xlb + xub)/2.0;
1221 
1222  if( SCIPisZero(scip, SCIPvarGetObj(consdata->z)) )
1223  {
1224  /* even simpler case where objective is linear in x */
1225  if( SCIPisZero(scip, SCIPvarGetObj(consdata->x)) )
1226  {
1227  /* simplest case where objective is zero:
1228  * if zero is within bounds, fix to zero, otherwise
1229  * fix x to middle of bounds for numerical stability. */
1230  if(SCIPisLT(scip, xlb, 0.0) && SCIPisGT(scip, xub, 0.0))
1231  xfix = 0.0;
1232  else
1233  xfix = 0.5 * (xlb + xub);
1234  }
1235  else
1236  {
1237  /* fix x to best bound */
1238  xfix = (SCIPvarGetObj(consdata->x) >= 0.0) ? xlb : xub;
1239  }
1240  }
1241  else if( consdata->exponent == 2.0 )
1242  {
1243  /* consider cases x <= -offset and x >= -offset separately */
1244  SCIP_Real a;
1245  SCIP_Real b;
1246  SCIP_Real c;
1247  SCIP_Real cand;
1248  SCIP_Real xfixobjval;
1249 
1250  xfix = SCIP_INVALID;
1251  xfixobjval = SCIP_INVALID;
1252 
1253  if( SCIPisLT(scip, xlb, -consdata->xoffset) )
1254  {
1255  /* For x <= -offset, the objective is equivalent to
1256  * zobj/zcoef * x^2 + (xobj + 2 offset zobj/zcoef) * x + offset^2 * zobj/zcoef + other constant
1257  * <-> a * x^2 + b * x + c
1258  *
1259  * critical values for x are xlb, MIN(xub,-offset), and -b/(2*a)
1260  */
1261  a = SCIPvarGetObj(consdata->z) / consdata->zcoef;
1262  b = SCIPvarGetObj(consdata->x) + 2 * consdata->xoffset * SCIPvarGetObj(consdata->z) / consdata->zcoef;
1263  c = consdata->xoffset * consdata->xoffset * SCIPvarGetObj(consdata->z) / consdata->zcoef;
1264 
1265  if( a < 0.0 && SCIPisInfinity(scip, -xlb) )
1266  {
1267  /* if a < 0.0, then a*x^2 is unbounded for x -> -infinity, thus fix x to -infinity */
1268  xfix = -SCIPinfinity(scip);
1269  xfixobjval = -SCIPinfinity(scip);
1270  }
1271  else
1272  {
1273  /* initialize with value for x=xlb */
1274  xfix = xlb;
1275  xfixobjval = a * xlb * xlb + b * xlb + c;
1276 
1277  /* compare with value for x=MIN(-offset,xub) */
1278  cand = MIN(-consdata->xoffset, xub);
1279  if( xfixobjval > a * cand * cand + b * cand + c )
1280  {
1281  xfix = cand;
1282  xfixobjval = a * cand * cand + b * cand + c;
1283  }
1284 
1285  /* compare with value for x=-b/(2*a), if within bounds */
1286  cand = -b/(2.0*a);
1287  if( cand > xlb && cand < -consdata->xoffset && cand < xub && xfixobjval > -b*b/(4.0*a) + c )
1288  {
1289  xfix = cand;
1290  xfixobjval = -b*b/(4.0*a) + c;
1291  }
1292  }
1293  }
1294 
1295  if( SCIPisGT(scip, xub, -consdata->xoffset) )
1296  {
1297  /* For x >= -offset, the objective is equivalent to
1298  * -zobj/zcoef * x^2 + (xobj - 2 offset zobj/zcoef) * x - offset^2 * zobj/zcoef + constants
1299  * <-> a * x^2 + b * x + c
1300  *
1301  * critical values for x are xub, MAX(xlb,-offset), and -b/(2*a)
1302  */
1303  a = -SCIPvarGetObj(consdata->z) / consdata->zcoef;
1304  b = SCIPvarGetObj(consdata->x) - 2 * consdata->xoffset * SCIPvarGetObj(consdata->z) / consdata->zcoef;
1305  c = -consdata->xoffset * consdata->xoffset * SCIPvarGetObj(consdata->z) / consdata->zcoef;
1306 
1307  if( a < 0.0 && SCIPisInfinity(scip, xub) )
1308  {
1309  /* if a < 0.0, then a*x^2 is unbounded for x -> infinity, thus fix x to infinity */
1310  xfix = SCIPinfinity(scip);
1311  /* not needed: xfixobjval = SCIPinfinity(scip); */
1312  }
1313  else
1314  {
1315  if( xfix == SCIP_INVALID ) /*lint !e777*/
1316  {
1317  /* initialize with value for x=xub */
1318  xfix = xub;
1319  xfixobjval = a * xub * xub + b * xub + c;
1320  }
1321  else
1322  {
1323  /* compare with value for x=xub */
1324  cand = xub;
1325  if( xfixobjval > a * cand * cand + b * cand + c )
1326  {
1327  xfix = cand;
1328  xfixobjval = a * cand * cand + b * cand + c;
1329  }
1330  }
1331 
1332  /* compare with value for x=MAX(xlb,-offset) */
1333  cand = MAX(xlb, -consdata->xoffset);
1334  if( xfixobjval > a * cand * cand + b * cand + c )
1335  {
1336  xfix = cand;
1337  xfixobjval = a * cand * cand + b * cand + c;
1338  }
1339 
1340  /* compare with value for x=-b/(2*a), if within bounds */
1341  cand = -b/(2.0*a);
1342  if( cand > xlb && cand > -consdata->xoffset && cand < xub && xfixobjval > -b*b/(4.0*a) + c )
1343  {
1344  xfix = cand;
1345  /* not needed: xfixobjval = -b*b/(4.0*a) + c; */
1346  }
1347  }
1348  }
1349  assert(xfix != SCIP_INVALID); /*lint !e777*/
1350  assert(SCIPisInfinity(scip, -xlb) || SCIPisLE(scip, xlb, xfix));
1351  assert(SCIPisInfinity(scip, xub) || SCIPisGE(scip, xub, xfix));
1352  }
1353  else
1354  {
1355  /* skip dual presolve for exponents != 2 and z in objective for now */
1356  return SCIP_OKAY;
1357  }
1358 
1359  /* compute fixing value for z */
1360  if( SCIPisInfinity(scip, xfix) )
1361  {
1362  if( consdata->zcoef > 0.0 )
1363  {
1364  assert(SCIPisInfinity(scip, -SCIPvarGetLbGlobal(consdata->z)));
1365  zfix = -SCIPinfinity(scip);
1366  }
1367  else
1368  {
1369  assert(SCIPisInfinity(scip, SCIPvarGetUbGlobal(consdata->z)));
1370  zfix = SCIPinfinity(scip);
1371  }
1372  }
1373  else if( SCIPisInfinity(scip, -xfix) )
1374  {
1375  if( consdata->zcoef > 0.0 )
1376  {
1377  assert(SCIPisInfinity(scip, SCIPvarGetUbGlobal(consdata->z)));
1378  zfix = SCIPinfinity(scip);
1379  }
1380  else
1381  {
1382  assert(SCIPisInfinity(scip, -SCIPvarGetLbGlobal(consdata->z)));
1383  zfix = -SCIPinfinity(scip);
1384  }
1385  }
1386  else
1387  {
1388  SCIP_Real zlb;
1389  SCIP_Real zub;
1390 
1391  zlb = SCIPvarGetLbGlobal(consdata->z);
1392  zub = SCIPvarGetUbGlobal(consdata->z);
1393  zfix = consdata->rhs - SIGN(xfix + consdata->xoffset) * consdata->power(ABS(xfix + consdata->xoffset), consdata->exponent);
1394  zfix /= consdata->zcoef;
1395 
1396  /* project zfix into box, it should be at least very close */
1397  assert(SCIPisFeasLE(scip, zlb, zfix));
1398  assert(SCIPisFeasGE(scip, zub, zfix));
1399  zfix = MAX(zlb, MIN(zub, zfix));
1400  }
1401 
1402  /* fix variables according to x=xfix */
1403  SCIPdebugMsg(scip, "dual presolve fixes x=<%s>[%g,%g] to %g and z=<%s>[%g,%g] to %g in cons <%s>\n",
1404  SCIPvarGetName(consdata->x), xlb, xub, xfix,
1405  SCIPvarGetName(consdata->z), SCIPvarGetLbGlobal(consdata->z), SCIPvarGetUbGlobal(consdata->z), zfix,
1406  SCIPconsGetName(cons));
1407  SCIPdebugPrintCons(scip, cons, NULL);
1408 
1409  /* fix x */
1410  SCIP_CALL( SCIPfixVar(scip, consdata->x, xfix, cutoff, &fixed) );
1411  if( *cutoff )
1412  return SCIP_OKAY;
1413  if( fixed )
1414  ++*nfixedvars;
1415 
1416  /* fix z */
1417  SCIP_CALL( SCIPfixVar(scip, consdata->z, zfix, cutoff, &fixed) );
1418  if( *cutoff )
1419  return SCIP_OKAY;
1420  if( fixed )
1421  ++*nfixedvars;
1422 
1423  /* delete constraint */
1424  SCIP_CALL( SCIPdelCons(scip, cons) );
1425  ++*ndelconss;
1426  }
1427  }
1428 
1429  return SCIP_OKAY;
1430 }
1431 
1432 /** given a variable and an interval, tightens the local bounds of this variable to the given interval */
1433 static
1435  SCIP* scip, /**< SCIP data structure */
1436  SCIP_VAR* var, /**< variable which bounds to tighten */
1437  SCIP_INTERVAL bounds, /**< new bounds */
1438  SCIP_Bool force, /**< force tightening even if below bound strengthening tolerance */
1439  SCIP_CONS* cons, /**< constraint that is propagated */
1440  SCIP_RESULT* result, /**< pointer to store the result of the propagation call */
1441  int* nchgbds, /**< buffer where to add the number of changed bounds */
1442  int* nfixedvars, /**< buffer where to add the number of fixed variables, can be equal to nchgbds */
1443  int* naddconss /**< buffer where to add the number of added constraints, can be NULL if force is FALSE */
1444  )
1445 {
1446  SCIP_Bool infeas;
1447  SCIP_Bool tightened;
1448 
1449  assert(scip != NULL);
1450  assert(var != NULL);
1451  assert(cons != NULL);
1452  assert(result != NULL);
1453  assert(nchgbds != NULL);
1454  assert(nfixedvars != NULL);
1455 
1456  *result = SCIP_DIDNOTFIND;
1457 
1458  if( SCIPisInfinity(scip, SCIPintervalGetInf(bounds)) || SCIPisInfinity(scip, -SCIPintervalGetSup(bounds)) )
1459  {
1460  /* domain outside [-infty, +infty] -> declare as infeasible */
1461  *result = SCIP_CUTOFF;
1462  return SCIP_OKAY;
1463  }
1464 
1465  /* if variable is not multiaggregated (or aggregated to a multiaggregated), then try SCIPfixVar or SCIPtightenVarLb/Ub
1466  * otherwise, if bound tightening is forced, add a linear constraint
1467  * otherwise, forget about the bound tightening
1468  */
1470  {
1471  /* check if variable can be fixed */
1472  if( SCIPisEQ(scip, bounds.inf, bounds.sup) )
1473  {
1474  if( !SCIPisEQ(scip, SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var)) )
1475  {
1476  /* if variable not fixed yet, then do so now */
1477  SCIP_Real fixval;
1478 
1479  if( bounds.inf != bounds.sup ) /*lint !e777*/
1480  fixval = (bounds.inf + bounds.sup) / 2.0;
1481  else
1482  fixval = bounds.inf;
1483  SCIP_CALL( SCIPfixVar(scip, var, fixval, &infeas, &tightened) );
1484 
1485  if( infeas )
1486  {
1487  SCIPdebugMsg(scip, "found <%s> infeasible due to fixing variable <%s>\n", SCIPconsGetName(cons), SCIPvarGetName(var));
1488  *result = SCIP_CUTOFF;
1489  return SCIP_OKAY;
1490  }
1491  if( tightened )
1492  {
1493  SCIPdebugMsg(scip, "fixed variable <%s> in constraint <%s> to %g\n", SCIPvarGetName(var), SCIPconsGetName(cons), SCIPvarGetLbLocal(var));
1494  ++*nfixedvars;
1495  *result = SCIP_REDUCEDDOM;
1496  }
1497  }
1498  else
1499  {
1500  /* only check if new fixing value is consistent with variable bounds, otherwise cutoff */
1501  if( SCIPisLT(scip, bounds.sup, SCIPvarGetUbLocal(var)) || SCIPisGT(scip, bounds.inf, SCIPvarGetLbLocal(var)) )
1502  {
1503  SCIPdebugMsg(scip, "found <%s> infeasible due to fixing fixed variable <%s>[%.15g,%.15g] to [%.15g,%.15g]\n",
1504  SCIPconsGetName(cons), SCIPvarGetName(var), SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var), bounds.inf, bounds.sup);
1505  *result = SCIP_CUTOFF;
1506  return SCIP_OKAY;
1507  }
1508  }
1509 
1510  return SCIP_OKAY;
1511  }
1512 
1513  /* check if lower bound can be tightened */
1514  if( SCIPintervalGetInf(bounds) > SCIPvarGetLbLocal(var) )
1515  {
1516  assert(!SCIPisInfinity(scip, -SCIPintervalGetInf(bounds)));
1517  SCIP_CALL( SCIPtightenVarLb(scip, var, SCIPintervalGetInf(bounds), force, &infeas, &tightened) );
1518  if( infeas )
1519  {
1520  SCIPdebugMsg(scip, "found %s infeasible due to domain propagation for variable %s in constraint %s\n", SCIPconsGetName(cons), SCIPvarGetName(var), SCIPconsGetName(cons));
1521  *result = SCIP_CUTOFF;
1522  return SCIP_OKAY;
1523  }
1524  if( tightened )
1525  {
1526  SCIPdebugMsg(scip, "tightened lower bound of variable %s in constraint %s to %g\n", SCIPvarGetName(var), SCIPconsGetName(cons), SCIPvarGetLbLocal(var));
1527  ++*nchgbds;
1528  *result = SCIP_REDUCEDDOM;
1529  }
1530  }
1531 
1532  /* check if upper bound can be tightened */
1533  if( SCIPintervalGetSup(bounds) < SCIPvarGetUbLocal(var) )
1534  {
1535  assert(!SCIPisInfinity(scip, SCIPintervalGetSup(bounds)));
1536  SCIP_CALL( SCIPtightenVarUb(scip, var, SCIPintervalGetSup(bounds), force, &infeas, &tightened) );
1537  if( infeas )
1538  {
1539  SCIPdebugMsg(scip, "found %s infeasible due to domain propagation for linear variable %s in constraint %s\n", SCIPconsGetName(cons), SCIPvarGetName(var), SCIPconsGetName(cons));
1540  *result = SCIP_CUTOFF;
1541  return SCIP_OKAY;
1542  }
1543  if( tightened )
1544  {
1545  SCIPdebugMsg(scip, "tightened upper bound of variable %s in constraint %s to %g\n", SCIPvarGetName(var), SCIPconsGetName(cons), SCIPvarGetUbLocal(var));
1546  ++*nchgbds;
1547  *result = SCIP_REDUCEDDOM;
1548  }
1549  }
1550  }
1551  else if( force && (SCIPisLT(scip, SCIPvarGetLbLocal(var), bounds.inf) || SCIPisGT(scip, SCIPvarGetUbLocal(var), bounds.sup)) )
1552  {
1553  /* add a linear constraint bounds.inf <= x <= bounds.sup */
1554  SCIP_CONS* auxcons;
1555  SCIP_Bool local;
1556  SCIP_Real one;
1557 
1558  assert(naddconss != NULL);
1559 
1560  /* we add constraint as local constraint if we are during probing or if we are during solve and not at the root node */
1561  local = SCIPinProbing(scip) || (SCIPgetStage(scip) == SCIP_STAGE_SOLVING && (SCIPnodeGetDepth(SCIPgetCurrentNode(scip)) > 0));
1562 
1563  one = 1.0;
1564  SCIP_CALL( SCIPcreateConsLinear(scip, &auxcons, SCIPconsGetName(cons), 1, &var, &one, bounds.inf, bounds.sup,
1566  SCIPconsIsChecked(cons), SCIPconsIsPropagated(cons), local,
1567  FALSE, FALSE, TRUE, FALSE) );
1568 
1569  if( local )
1570  {
1571  SCIP_CALL( SCIPaddConsLocal(scip, auxcons, NULL) );
1572  }
1573  else
1574  {
1575  SCIP_CALL( SCIPaddCons(scip, auxcons) );
1576  }
1577  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
1578 
1579  ++*naddconss;
1580  *result = SCIP_CONSADDED;
1581  }
1582 
1583  return SCIP_OKAY;
1584 }
1585 
1586 /** computes bounds on z in a absolute power constraints for given bounds on x */
1587 static
1588 void computeBoundsZ(
1589  SCIP* scip, /**< SCIP data structure */
1590  SCIP_CONS* cons, /**< constraint */
1591  SCIP_INTERVAL xbnds, /**< bounds on x that are to be propagated */
1592  SCIP_INTERVAL* zbnds /**< buffer to store corresponding bounds on z */
1593  )
1594 {
1595  SCIP_CONSDATA* consdata;
1596  SCIP_Real bnd;
1597  SCIP_Real x;
1598 
1599  assert(scip != NULL);
1600  assert(cons != NULL);
1601  assert(zbnds != NULL);
1602  assert(!SCIPintervalIsEmpty(SCIPinfinity(scip), xbnds));
1603 
1604  consdata = SCIPconsGetData(cons);
1605  assert(consdata != NULL);
1606 
1607  SCIPintervalSetEntire(SCIPinfinity(scip), zbnds);
1608 
1609  /* apply zcoef*z <= rhs - signedpow(xbnds.inf + offset, n) */
1610  if( !SCIPisInfinity(scip, consdata->rhs) && !SCIPisInfinity(scip, -xbnds.inf) )
1611  {
1612  x = xbnds.inf - PROPVARTOL + consdata->xoffset;
1613  bnd = consdata->rhs + PROPSIDETOL - SIGN(x) * consdata->power(REALABS(x), consdata->exponent);
1614 
1615  if( consdata->zcoef > 0.0 )
1616  zbnds->sup = bnd / consdata->zcoef;
1617  else
1618  zbnds->inf = bnd / consdata->zcoef;
1619  }
1620 
1621  /* apply zcoef*z >= lhs - signedpow(xbnds.sup + offset, n) */
1622  if( !SCIPisInfinity(scip, -consdata->lhs) && !SCIPisInfinity(scip, xbnds.sup) )
1623  {
1624  x = xbnds.sup + PROPVARTOL + consdata->xoffset;
1625  bnd = consdata->lhs - PROPSIDETOL - SIGN(x) * consdata->power(REALABS(x), consdata->exponent);
1626 
1627  if( consdata->zcoef > 0.0 )
1628  zbnds->inf = bnd / consdata->zcoef;
1629  else
1630  zbnds->sup = bnd / consdata->zcoef;
1631  }
1632 
1633  SCIPdebugMsg(scip, "given x = [%.15g, %.15g], computed z = [%.15g, %.15g] via", xbnds.inf, xbnds.sup, zbnds->inf, zbnds->sup);
1634  SCIPdebugPrintCons(scip, cons, NULL);
1635 
1636  assert(!SCIPintervalIsEmpty(SCIPinfinity(scip), *zbnds));
1637 }
1638 
1639 /** computes bounds on x in a absolute power constraints for given bounds on z */
1640 static
1641 void computeBoundsX(
1642  SCIP* scip, /**< SCIP data structure */
1643  SCIP_CONS* cons, /**< constraint */
1644  SCIP_INTERVAL zbnds, /**< bounds on x that are to be propagated */
1645  SCIP_INTERVAL* xbnds /**< buffer to store corresponding bounds on z */
1646  )
1647 {
1648  SCIP_CONSDATA* consdata;
1649  SCIP_Real bnd;
1650  SCIP_Real z;
1651 
1652  assert(scip != NULL);
1653  assert(cons != NULL);
1654  assert(xbnds != NULL);
1655  assert(!SCIPintervalIsEmpty(SCIPinfinity(scip), zbnds));
1656 
1657  consdata = SCIPconsGetData(cons);
1658  assert(consdata != NULL);
1659 
1660  SCIPintervalSetEntire(SCIPinfinity(scip), xbnds);
1661 
1662  /* apply signedpow(x+offset, n) <= rhs - (zcoef * zbnds).inf */
1663  z = (consdata->zcoef > 0.0 ? zbnds.inf : zbnds.sup);
1664  if( !SCIPisInfinity(scip, consdata->rhs) && !SCIPisInfinity(scip, REALABS(z)) )
1665  {
1666  bnd = consdata->rhs + PROPSIDETOL - consdata->zcoef * z + REALABS(consdata->zcoef) * PROPVARTOL;
1667  if( consdata->exponent == 2.0 )
1668  bnd = SIGN(bnd) * sqrt(REALABS(bnd));
1669  else
1670  bnd = SIGN(bnd) * pow(REALABS(bnd), 1.0/consdata->exponent);
1671  xbnds->sup = bnd - consdata->xoffset;
1672  }
1673 
1674  /* apply signedpow(x+offset, n) >= lhs - (zcoef * zbnds).sup */
1675  z = (consdata->zcoef > 0.0 ? zbnds.sup : zbnds.inf);
1676  if( !SCIPisInfinity(scip, consdata->rhs) && !SCIPisInfinity(scip, REALABS(z)) )
1677  {
1678  bnd = consdata->lhs - PROPSIDETOL - consdata->zcoef * z - REALABS(consdata->zcoef) * PROPVARTOL;
1679  if( consdata->exponent == 2.0 )
1680  bnd = SIGN(bnd) * sqrt(REALABS(bnd));
1681  else
1682  bnd = SIGN(bnd) * pow(REALABS(bnd), 1.0/consdata->exponent);
1683  xbnds->inf = bnd - consdata->xoffset;
1684  }
1685 
1686  SCIPdebugMsg(scip, "given z = [%.15g, %.15g], computed x = [%.15g, %.15g] via", zbnds.inf, zbnds.sup, xbnds->inf, xbnds->sup);
1687  SCIPdebugPrintCons(scip, cons, NULL);
1688 
1689  assert(!SCIPintervalIsEmpty(SCIPinfinity(scip), *xbnds));
1690 }
1691 
1692 /** checks if x or z is fixed and replaces them or deletes constraint */
1693 static
1695  SCIP* scip, /**< SCIP data structure */
1696  SCIP_CONSHDLR* conshdlr, /**< constraint handler for absolute power constraints */
1697  SCIP_CONS* cons, /**< constraint */
1698  int* ndelconss, /**< counter for number of deleted constraints */
1699  int* nupgdconss, /**< counter for number of upgraded constraints */
1700  int* nchgbds, /**< counter for number of variable bound changes */
1701  int* nfixedvars, /**< counter for number of variable fixations */
1702  SCIP_RESULT* result /**< to store result if we detect infeasibility or remove constraint */
1703  )
1704 {
1705  SCIP_CONSHDLRDATA* conshdlrdata;
1706  SCIP_CONSDATA* consdata;
1707  SCIP_Real scalar;
1708  SCIP_Real constant;
1709  SCIP_Real factor;
1710  SCIP_VAR* var;
1711 
1712  assert(scip != NULL);
1713  assert(cons != NULL);
1714  assert(ndelconss != NULL);
1715  assert(nupgdconss != NULL);
1716  assert(nchgbds != NULL);
1717  assert(nfixedvars != NULL);
1718 
1719  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1720  assert(conshdlrdata != NULL);
1721 
1722  consdata = SCIPconsGetData(cons);
1723  assert(consdata != NULL);
1724 
1725  *result = SCIP_DIDNOTFIND;
1726 
1727  if( !SCIPvarIsActive(consdata->x) && SCIPvarGetStatus(consdata->x) != SCIP_VARSTATUS_MULTAGGR )
1728  {
1729  /* replace x variable */
1730 
1731  /* get relation x = scalar * var + constant */
1732  var = consdata->x;
1733  scalar = 1.0;
1734  constant = 0.0;
1735  SCIP_CALL( SCIPgetProbvarSum(scip, &var, &scalar, &constant) );
1736 
1737  if( scalar == 0.0 )
1738  {
1739  SCIP_INTERVAL xbnds;
1740  SCIP_INTERVAL zbnds;
1741  int naddconss;
1742 
1743  naddconss = 0;
1744 
1745  /* x has been fixed to constant */
1746  assert(SCIPisFeasEQ(scip, SCIPvarGetLbGlobal(consdata->x), constant));
1747  assert(SCIPisFeasEQ(scip, SCIPvarGetUbGlobal(consdata->x), constant));
1748 
1749  /* compute corresponding bounds on z */
1750  SCIPintervalSet(&xbnds, constant);
1751  computeBoundsZ(scip, cons, xbnds, &zbnds);
1752 
1753  SCIPdebugMsg(scip, "in cons <%s>: x = <%s> fixed to %g -> tighten <%s> to [%g, %g]\n", SCIPconsGetName(cons), SCIPvarGetName(consdata->x), constant, SCIPvarGetName(consdata->z), zbnds.inf, zbnds.sup);
1754 
1755  if( SCIPisEQ(scip, consdata->lhs, consdata->rhs) )
1756  {
1757  /* if sides are equal, then we should either fix z, or declare infeasibility */
1758  if( SCIPisFeasLT(scip, SCIPvarGetUbGlobal(consdata->z), zbnds.inf) || SCIPisFeasGT(scip, SCIPvarGetLbGlobal(consdata->z), zbnds.sup) )
1759  {
1760  SCIPdebugMsg(scip, "bounds inconsistent -> cutoff\n");
1761  *result = SCIP_CUTOFF;
1762  return SCIP_OKAY;
1763  }
1764  else
1765  {
1766  /* compute fixing value for z as value corresponding to fixing of x, projected onto bounds of z */
1767  SCIP_Real zfix;
1768 
1769  zfix = consdata->rhs - SIGN(constant + consdata->xoffset) * consdata->power(REALABS(constant + consdata->xoffset), consdata->exponent);
1770  zfix /= consdata->zcoef;
1771  assert(SCIPisLE(scip, zbnds.inf, zfix));
1772  assert(SCIPisGE(scip, zbnds.sup, zfix));
1773  zfix = MIN(SCIPvarGetUbGlobal(consdata->z), MAX(SCIPvarGetLbGlobal(consdata->z), zfix)); /*lint !e666*/
1774 
1775  zbnds.inf = zfix;
1776  zbnds.sup = zfix;
1777  SCIP_CALL( tightenBounds(scip, consdata->z, zbnds, TRUE, cons, result, nchgbds, nfixedvars, &naddconss) );
1778  }
1779  }
1780  else
1781  {
1782  /* tighten bounds on z accordingly */
1783  SCIP_CALL( tightenBounds(scip, consdata->z, zbnds, TRUE, cons, result, nchgbds, nfixedvars, &naddconss) );
1784  }
1785 
1786  /* delete constraint */
1787  SCIP_CALL( SCIPdelCons(scip, cons) );
1788 
1789  /* if tightenBounds added a constraint (because z was multiaggregated), then count this as constraint upgrade, otherwise as constraint deletion */
1790  if( naddconss > 0 )
1791  ++*nupgdconss;
1792  else
1793  ++*ndelconss;
1794 
1795  return SCIP_OKAY;
1796  }
1797 
1798  SCIPdebugMsg(scip, "in cons <%s>: x = <%s> replaced by %g*<%s> + %g\n", SCIPconsGetName(cons), SCIPvarGetName(consdata->x), scalar, SCIPvarGetName(var), constant);
1799 
1800  /* constraint will be divided by scalar*pow(|scalar|,exponent-1), if scalar is not 1.0 */
1801  if( scalar == 1.0 )
1802  factor = 1.0;
1803  else if( scalar > 0.0 )
1804  factor = consdata->power( scalar, consdata->exponent);
1805  else
1806  factor = -consdata->power(-scalar, consdata->exponent);
1807 
1808  /* aggregate only if this would not lead to a vanishing or infinite coefficient for z */
1809  if( !SCIPisZero(scip, consdata->zcoef / factor) && !SCIPisInfinity(scip, REALABS(consdata->zcoef / factor)) )
1810  {
1811  /* we drop here the events for both variables, because if x is replaced by a multiaggregated variable here, then we do not need to catch bound tightenings on z anymore */
1812  SCIP_CALL( dropVarEvents(scip, conshdlrdata->eventhdlr, cons) );
1813  SCIP_CALL( SCIPunlockVarCons(scip, consdata->x, cons, !SCIPisInfinity(scip, -consdata->lhs), !SCIPisInfinity(scip, consdata->rhs)) );
1814 
1815  consdata->x = var;
1816  if( SCIPvarIsActive(consdata->x) )
1817  {
1818  SCIP_CALL( SCIPmarkDoNotMultaggrVar(scip, consdata->x) );
1819  }
1820 
1821  /* add constant to offset */
1822  consdata->xoffset += constant;
1823 
1824  /* divide constraint by factor */
1825  if( scalar == 1.0 ) ;
1826  else if( scalar > 0.0 )
1827  {
1828  if( !SCIPisInfinity(scip, -consdata->lhs) )
1829  consdata->lhs /= factor;
1830  if( !SCIPisInfinity(scip, consdata->rhs) )
1831  consdata->rhs /= factor;
1832  consdata->zcoef /= factor;
1833  consdata->xoffset /= scalar;
1834  }
1835  else
1836  {
1837  SCIP_Real oldlhs;
1838 
1839  assert(scalar < 0.0);
1840  assert(factor < 0.0);
1841 
1842  oldlhs = consdata->lhs;
1843 
1844  if( !SCIPisInfinity(scip, consdata->rhs) )
1845  consdata->lhs = consdata->rhs / factor;
1846  else
1847  consdata->lhs = -SCIPinfinity(scip);
1848  if( !SCIPisInfinity(scip, -oldlhs) )
1849  consdata->rhs = oldlhs / factor;
1850  else
1851  consdata->rhs = SCIPinfinity(scip);
1852  consdata->zcoef /= factor;
1853  consdata->xoffset /= scalar;
1854  /* since we flip both constraint sides and the sign of zcoef, the events catched for z remain the same, so update necessary there */
1855  }
1856 
1857  SCIP_CALL( SCIPlockVarCons(scip, consdata->x, cons, !SCIPisInfinity(scip, -consdata->lhs), !SCIPisInfinity(scip, consdata->rhs)) );
1858  SCIP_CALL( catchVarEvents(scip, conshdlrdata->eventhdlr, cons) );
1859 
1860  SCIPdebugPrintCons(scip, cons, NULL);
1861 
1862  /* rerun constraint comparison */
1863  conshdlrdata->comparedpairwise = FALSE;
1864  }
1865  else
1866  {
1867  SCIPwarningMessage(scip, "Skip resolving aggregation of variable <%s> in abspower constraint <%s> to avoid zcoef = %g\n",
1868  SCIPvarGetName(consdata->x), SCIPconsGetName(cons), consdata->zcoef / factor);
1869  }
1870  }
1871 
1872  if( !SCIPvarIsActive(consdata->z) && SCIPvarGetStatus(consdata->z) != SCIP_VARSTATUS_MULTAGGR )
1873  {
1874  /* replace z variable */
1875 
1876  /* get relation z = scalar * var + constant */
1877  var = consdata->z;
1878  scalar = 1.0;
1879  constant = 0.0;
1880  SCIP_CALL( SCIPgetProbvarSum(scip, &var, &scalar, &constant) );
1881 
1882  if( scalar == 0.0 )
1883  {
1884  SCIP_INTERVAL xbnds;
1885  SCIP_INTERVAL zbnds;
1886  int naddconss;
1887 
1888  naddconss = 0;
1889 
1890  /* z has been fixed to constant */
1891  assert(SCIPisFeasEQ(scip, SCIPvarGetLbGlobal(consdata->z), constant));
1892  assert(SCIPisFeasEQ(scip, SCIPvarGetUbGlobal(consdata->z), constant));
1893 
1894  /* compute corresponding bounds on x */
1895  SCIPintervalSet(&zbnds, constant);
1896  computeBoundsX(scip, cons, zbnds, &xbnds);
1897 
1898  SCIPdebugMsg(scip, "in cons <%s>: z = <%s> fixed to %g -> tighten <%s> to [%g, %g]\n", SCIPconsGetName(cons), SCIPvarGetName(consdata->z), constant, SCIPvarGetName(consdata->x), xbnds.inf, xbnds.sup);
1899 
1900  if( SCIPisEQ(scip, consdata->lhs, consdata->rhs) )
1901  {
1902  /* if sides are equal, then we should either fix x, or declare infeasibility */
1903  if( SCIPisFeasLT(scip, SCIPvarGetUbGlobal(consdata->x), xbnds.inf) || SCIPisFeasGT(scip, SCIPvarGetLbGlobal(consdata->x), xbnds.sup) )
1904  {
1905  SCIPdebugMsg(scip, "bounds inconsistent -> cutoff\n");
1906  *result = SCIP_CUTOFF;
1907  return SCIP_OKAY;
1908  }
1909  else
1910  {
1911  /* compute fixing value for x as value corresponding to fixing of z, projected onto bounds of x */
1912  SCIP_Real xfix;
1913 
1914  xfix = consdata->rhs - consdata->zcoef * constant;
1915  if( consdata->exponent == 2.0 )
1916  xfix = SIGN(xfix) * sqrt(REALABS(xfix)) - consdata->xoffset;
1917  else
1918  xfix = SIGN(xfix) * pow(REALABS(xfix), 1.0/consdata->exponent) - consdata->xoffset;
1919  assert(SCIPisLE(scip, xbnds.inf, xfix));
1920  assert(SCIPisGE(scip, xbnds.sup, xfix));
1921  xfix = MIN(SCIPvarGetUbGlobal(consdata->x), MAX(SCIPvarGetLbGlobal(consdata->x), xfix)); /*lint !e666*/
1922 
1923  xbnds.inf = xfix;
1924  xbnds.sup = xfix;
1925  SCIP_CALL( tightenBounds(scip, consdata->x, xbnds, TRUE, cons, result, nchgbds, nfixedvars, &naddconss) );
1926  }
1927  }
1928  else
1929  {
1930  /* tighten bounds on x accordingly */
1931  SCIP_CALL( tightenBounds(scip, consdata->x, xbnds, TRUE, cons, result, nchgbds, nfixedvars, &naddconss) );
1932  }
1933 
1934  /* delete constraint */
1935  SCIP_CALL( SCIPdelCons(scip, cons) );
1936 
1937  /* if tightenBounds added a constraint (because x was multiaggregated), then count this as constraint upgrade, otherwise as constraint deletion */
1938  if( naddconss > 0 )
1939  ++*nupgdconss;
1940  else
1941  ++*ndelconss;
1942 
1943  return SCIP_OKAY;
1944  }
1945 
1946  SCIPdebugMsg(scip, "in cons <%s>: z = <%s> replaced by %g*<%s> + %g\n", SCIPconsGetName(cons), SCIPvarGetName(consdata->z), scalar, SCIPvarGetName(var), constant);
1947 
1948  /* we drop here the events for both variables, because if z is replaced by a multiaggregated variable here, then we do not need to catch bound tightenings on x anymore */
1949  SCIP_CALL( dropVarEvents(scip, conshdlrdata->eventhdlr, cons) );
1950  if( consdata->zcoef > 0.0 )
1951  SCIP_CALL( SCIPunlockVarCons(scip, consdata->z, cons, !SCIPisInfinity(scip, -consdata->lhs), !SCIPisInfinity(scip, consdata->rhs)) );
1952  else
1953  SCIP_CALL( SCIPunlockVarCons(scip, consdata->z, cons, !SCIPisInfinity(scip, consdata->rhs), !SCIPisInfinity(scip, -consdata->lhs)) );
1954 
1955  consdata->z = var;
1956  if( SCIPvarIsActive(consdata->z) )
1957  {
1958  SCIP_CALL( SCIPmarkDoNotMultaggrVar(scip, consdata->z) );
1959  }
1960 
1961  /* substract constant from constraint sides */
1962  if( !SCIPisInfinity(scip, -consdata->lhs) )
1963  consdata->lhs -= consdata->zcoef * constant;
1964  if( !SCIPisInfinity(scip, consdata->rhs) )
1965  consdata->rhs -= consdata->zcoef * constant;
1966 
1967  /* multiply zcoef by scalar */
1968  consdata->zcoef *= scalar;
1969 
1970  if( consdata->zcoef > 0.0 )
1971  SCIP_CALL( SCIPlockVarCons(scip, consdata->z, cons, !SCIPisInfinity(scip, -consdata->lhs), !SCIPisInfinity(scip, consdata->rhs)) );
1972  else
1973  SCIP_CALL( SCIPlockVarCons(scip, consdata->z, cons, !SCIPisInfinity(scip, consdata->rhs), !SCIPisInfinity(scip, -consdata->lhs)) );
1974  SCIP_CALL( catchVarEvents(scip, conshdlrdata->eventhdlr, cons) );
1975 
1976  /* rerun constraint comparison */
1977  conshdlrdata->comparedpairwise = FALSE;
1978  }
1979 
1980  assert(SCIPvarIsActive(consdata->z) || SCIPvarGetStatus(consdata->z) == SCIP_VARSTATUS_MULTAGGR);
1981 
1982  return SCIP_OKAY;
1983 }
1984 
1985 /** computes violation of a constraint */
1986 static
1988  SCIP* scip, /**< SCIP data structure */
1989  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
1990  SCIP_CONS* cons, /**< constraint */
1991  SCIP_SOL* sol, /**< solution or NULL if LP solution should be used */
1992  SCIP_Real* viol, /**< pointer to store absolute (unscaled) constraint violation */
1993  SCIP_Bool* solviolbounds /**< buffer to store whether the solution violates bounds on x by more than feastol */
1994  )
1995 {
1996  SCIP_CONSDATA* consdata;
1997  SCIP_Real val;
1998  SCIP_Real xval;
1999  SCIP_Real zval;
2000  SCIP_Real relviol;
2001 
2002  assert(scip != NULL);
2003  assert(conshdlr != NULL);
2004  assert(cons != NULL);
2005  assert(viol != NULL);
2006  assert(solviolbounds != NULL);
2007 
2008  consdata = SCIPconsGetData(cons);
2009  assert(consdata != NULL);
2010 
2011  *solviolbounds = FALSE;
2012  xval = SCIPgetSolVal(scip, sol, consdata->x);
2013  zval = SCIPgetSolVal(scip, sol, consdata->z);
2014 
2015  if( SCIPisInfinity(scip, REALABS(xval)) )
2016  {
2017  consdata->lhsviol = (SCIPisInfinity(scip, -consdata->lhs) ? 0.0 : SCIPinfinity(scip));
2018  consdata->rhsviol = (SCIPisInfinity(scip, consdata->rhs) ? 0.0 : SCIPinfinity(scip));
2019 
2020  return SCIP_OKAY;
2021  }
2022 
2023  if( sol == NULL )
2024  {
2025  SCIP_Real lb;
2026  SCIP_Real ub;
2027 
2028  lb = SCIPvarGetLbLocal(consdata->x);
2029  ub = SCIPvarGetUbLocal(consdata->x);
2030 
2031  /* with non-initial columns, variables can shortly be a column variable before entering the LP and have value 0.0 in this case, which might be outside bounds */
2032  if( (!SCIPisInfinity(scip, -lb) && !SCIPisFeasGE(scip, xval, lb)) || (!SCIPisInfinity(scip, ub) && !SCIPisFeasLE(scip, xval, ub)) )
2033  *solviolbounds = TRUE;
2034  else
2035  xval = MAX(lb, MIN(ub, xval)); /* project x onto local box if just slightly outside (or even not outside) */
2036  }
2037 
2038  xval += consdata->xoffset;
2039 
2040  val = SIGN(xval) * consdata->power(REALABS(xval), consdata->exponent);
2041  val += consdata->zcoef * zval;
2042 
2043  *viol = 0.0;
2044  relviol = 0.0;
2045  if( val < consdata->lhs && !SCIPisInfinity(scip, -consdata->lhs) )
2046  {
2047  consdata->lhsviol = *viol = consdata->lhs - val;
2048  relviol = SCIPrelDiff(consdata->lhs, val);
2049  }
2050  else
2051  consdata->lhsviol = 0.0;
2052 
2053  if( val > consdata->rhs && !SCIPisInfinity(scip, consdata->rhs) )
2054  {
2055  consdata->rhsviol = *viol = val - consdata->rhs;
2056  relviol = SCIPrelDiff(val, consdata->rhs);
2057  }
2058  else
2059  consdata->rhsviol = 0.0;
2060 
2061  if( sol != NULL )
2062  SCIPupdateSolConsViolation(scip, sol, *viol, relviol);
2063 
2064  return SCIP_OKAY;
2065 }
2066 
2067 /** computes violation of a set of constraints */
2068 static
2070  SCIP* scip, /**< SCIP data structure */
2071  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
2072  SCIP_CONS** conss, /**< constraints */
2073  int nconss, /**< number of constraints */
2074  SCIP_SOL* sol, /**< solution or NULL if LP solution should be used */
2075  SCIP_Bool* solviolbounds, /**< buffer to store whether the solution violates bounds on x by more than feastol */
2076  SCIP_CONS** maxviolcon /**< buffer to store constraint with largest violation, or NULL if solution is feasible */
2077  )
2078 {
2079  SCIP_CONSDATA* consdata;
2080  SCIP_Real viol;
2081  SCIP_Real maxviol;
2082  SCIP_Bool solviolbounds1;
2083  int c;
2084 
2085  assert(scip != NULL);
2086  assert(conss != NULL || nconss == 0);
2087  assert(solviolbounds != NULL);
2088  assert(maxviolcon != NULL);
2089 
2090  *solviolbounds = FALSE;
2091  *maxviolcon = NULL;
2092 
2093  maxviol = 0.0;
2094 
2095  for( c = 0; c < nconss; ++c )
2096  {
2097  assert(conss != NULL);
2098  assert(conss[c] != NULL);
2099 
2100  SCIP_CALL( computeViolation(scip, conshdlr, conss[c], sol, &viol, &solviolbounds1) );
2101  *solviolbounds |= solviolbounds1;
2102 
2103  consdata = SCIPconsGetData(conss[c]);
2104  assert(consdata != NULL);
2105 
2106  viol = MAX(consdata->lhsviol, consdata->rhsviol);
2107  if( viol > maxviol && SCIPisGT(scip, viol, SCIPfeastol(scip)) )
2108  {
2109  maxviol = viol;
2110  *maxviolcon = conss[c];
2111  }
2112  }
2113 
2114  return SCIP_OKAY;
2115 }
2116 
2117 /** proposes branching point for constraint */
2118 static
2120  SCIP* scip, /**< SCIP data structure */
2121  SCIP_CONS* cons, /**< constraint which variable to get branching point for */
2122  SCIP_SOL* sol, /**< solution to branch on (NULL for LP or pseudosol) */
2123  int preferzero, /**< how much we prefer branching on -xoffset (0, 1, or 2) if sign is not fixed */
2124  SCIP_Bool branchminconverror /**< whether to minimize convexification error if sign is fixed */
2125  )
2126 {
2127  SCIP_CONSDATA* consdata;
2128  SCIP_VAR* x;
2129  SCIP_Real xref;
2130  SCIP_Real zref;
2131  SCIP_Real xlb;
2132  SCIP_Real xub;
2133 
2134  assert(scip != NULL);
2135  assert(cons != NULL);
2136 
2137  consdata = SCIPconsGetData(cons);
2138  assert(consdata != NULL);
2139 
2140  x = consdata->x;
2141  xlb = SCIPvarGetLbLocal(x);
2142  xub = SCIPvarGetUbLocal(x);
2143 
2144  /* check if sign of x is not fixed yet */
2145  if( SCIPisLT(scip, xlb, -consdata->xoffset) && SCIPisGT(scip, xub, -consdata->xoffset) )
2146  {
2147  /* if preferzero is 0, just return SCIP_INVALID
2148  * if preferzero is 1, then propose -xoffset if branching on -xoffset would cut off solution in both child nodes, otherwise return SCIP_INVALID
2149  * if preferzero is >1, then always propose -xoffset
2150  */
2151  assert(preferzero >= 0);
2152 
2153  if( preferzero == 0 )
2154  return SCIP_INVALID;
2155 
2156  if( preferzero > 1 || SCIPisInfinity(scip, -xlb) || SCIPisInfinity(scip, xub) )
2157  return -consdata->xoffset;
2158 
2159  xlb += consdata->xoffset;
2160  xub += consdata->xoffset;
2161 
2162  xref = SCIPgetSolVal(scip, sol, x) + consdata->xoffset;
2163  zref = SCIPgetSolVal(scip, sol, consdata->z);
2164  if( SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)) )
2165  {
2166  /* signpow(x,n,offset) + c*z <= 0 is violated
2167  * if we are close to or right of -offset, then branching on -offset gives a convex function on the right branch, this is good
2168  * otherwise if branching on -offset yields a violated secant cut in left branch, then current solution would be cutoff there, this is also still good
2169  */
2170  if( !SCIPisFeasNegative(scip, xref) || SCIPisFeasPositive(scip, -consdata->power(-xlb, consdata->exponent)*xref/xlb + consdata->zcoef * zref) )
2171  return -consdata->xoffset;
2172  return SCIP_INVALID;
2173  }
2174 
2175  assert(SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) );
2176  /* signpow(x,n) + c*z >= 0 is violated
2177  * if we are close to or left of zero, then branching on 0.0 gives a concave function on the left branch, this is good
2178  * otherwise if branching on 0.0 yields a violated secant cut in right branch, then current solution would be cutoff there, this is also still good
2179  */
2180  if( !SCIPisFeasPositive(scip, xref) || SCIPisFeasNegative(scip, -consdata->power(xub, consdata->exponent)*xref/xub + consdata->zcoef * zref) )
2181  return -consdata->xoffset;
2182  return SCIP_INVALID;
2183  }
2184 
2185  if( branchminconverror )
2186  {
2187  /* given x^n with xlb <= x <= xub, then the sum of the integrals between the function and its secant on the left and right branches are minimized
2188  * for branching on ( (ub^n - lb^n) / (n*(ub - lb)) ) ^ (1/(n-1))
2189  */
2190  if( SCIPisGE(scip, xlb, -consdata->xoffset) )
2191  {
2192  SCIP_Real ref;
2193  xlb = MAX(0.0, xlb + consdata->xoffset);
2194  xub = MAX(0.0, xub + consdata->xoffset);
2195 
2196  ref = (consdata->power(xub, consdata->exponent) - consdata->power(xlb, consdata->exponent)) / (consdata->exponent * (xub - xlb));
2197  ref = pow(ref, 1.0/(consdata->exponent-1.0));
2198  ref -= consdata->xoffset;
2199  assert(SCIPisGE(scip, ref, SCIPvarGetLbLocal(x)));
2200  assert(SCIPisLE(scip, ref, SCIPvarGetUbLocal(x)));
2201 
2202  return ref;
2203  }
2204  else
2205  {
2206  SCIP_Real ref;
2207 
2208  assert(SCIPisLE(scip, xub, -consdata->xoffset));
2209 
2210  xlb = MIN(0.0, xlb + consdata->xoffset);
2211  xub = MIN(0.0, xub + consdata->xoffset);
2212 
2213  ref = (consdata->power(-xlb, consdata->exponent) - consdata->power(-xub, consdata->exponent)) / (consdata->exponent * (-xlb + xub));
2214  ref = -pow(ref, 1.0/(consdata->exponent-1.0));
2215  ref -= consdata->xoffset;
2216  assert(SCIPisGE(scip, ref, SCIPvarGetLbLocal(x)));
2217  assert(SCIPisLE(scip, ref, SCIPvarGetUbLocal(x)));
2218 
2219  return ref;
2220  }
2221  }
2222 
2223  return SCIP_INVALID;
2224 }
2225 
2226 /** registers branching variable candidates
2227  * registers x for all violated absolute power constraints where x is not in convex region
2228  */
2229 static
2231  SCIP* scip, /**< SCIP data structure */
2232  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
2233  SCIP_CONS** conss, /**< constraints to check */
2234  int nconss, /**< number of constraints to check */
2235  SCIP_SOL* sol, /**< solution to enforce (NULL for the LP solution) */
2236  int* nnotify /**< counter for number of notifications performed */
2237  )
2238 {
2239  SCIP_CONSHDLRDATA* conshdlrdata;
2240  SCIP_CONSDATA* consdata;
2241  SCIP_Bool onlynonfixedsign;
2242  int c;
2243 
2244  assert(scip != NULL);
2245  assert(conshdlr != NULL);
2246  assert(conss != NULL || nconss == 0);
2247 
2248  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2249  assert(conshdlrdata != NULL);
2250 
2251  *nnotify = 0;
2252 
2253  onlynonfixedsign = conshdlrdata->preferzerobranch == 3;
2254 
2255  do
2256  {
2257  for( c = 0; c < nconss; ++c )
2258  {
2259  assert(conss[c] != NULL); /*lint !e613*/
2260 
2261  /* skip constraints that have been marked to be removed by propagateCons() */
2262  if( !SCIPconsIsEnabled(conss[c]) ) /*lint !e613*/
2263  continue;
2264 
2265  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
2266  assert(consdata != NULL);
2267 
2268  SCIPdebugMsg(scip, "cons <%s> violation: %g %g\n", SCIPconsGetName(conss[c]), consdata->lhsviol, consdata->rhsviol); /*lint !e613*/
2269 
2270  /* skip variables which sign is already fixed, if we are only interested in variables with unfixed sign here */
2271  if( onlynonfixedsign &&
2272  ( !SCIPisLT(scip, SCIPvarGetLbLocal(consdata->x), -consdata->xoffset) ||
2273  !SCIPisGT(scip, SCIPvarGetUbLocal(consdata->x), consdata->xoffset)) )
2274  continue;
2275 
2276  /* if the value of x lies in a concave range (i.e., where a secant approximation is used), then register x as branching variable */
2277  if( (SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)) && (SCIPisInfinity(scip, -SCIPvarGetLbLocal(consdata->x)) || SCIPgetSolVal(scip, sol, consdata->x) + consdata->xoffset <= -consdata->root * (SCIPvarGetLbLocal(consdata->x) + consdata->xoffset))) ||
2278  ( SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) && (SCIPisInfinity(scip, SCIPvarGetUbLocal(consdata->x)) || SCIPgetSolVal(scip, sol, consdata->x) + consdata->xoffset >= -consdata->root * (SCIPvarGetUbLocal(consdata->x) + consdata->xoffset))) )
2279  {
2280  /* domain propagation should have removed constraints with fixed x, at least for violated constraints */
2281  assert(!SCIPisRelEQ(scip, SCIPvarGetLbLocal(consdata->x), SCIPvarGetUbLocal(consdata->x)));
2282 
2283  SCIPdebugMsg(scip, "register var <%s> in cons <%s> with violation %g %g\n", SCIPvarGetName(consdata->x), SCIPconsGetName(conss[c]), consdata->lhsviol, consdata->rhsviol); /*lint !e613*/
2284  SCIP_CALL( SCIPaddExternBranchCand(scip, consdata->x, MAX(consdata->lhsviol, consdata->rhsviol), proposeBranchingPoint(scip, conss[c], sol, conshdlrdata->preferzerobranch, conshdlrdata->branchminconverror)) ); /*lint !e613*/
2285  ++*nnotify;
2286  }
2287  }
2288 
2289  if( onlynonfixedsign && *nnotify == 0 )
2290  {
2291  /* if we could not a variable in a violated constraint which sign is not already fixed, do another round where we consider all variables again */
2292  onlynonfixedsign = FALSE;
2293  continue;
2294  }
2295  break;
2296  }
2297  while( TRUE ); /*lint !e506 */
2298 
2299  return SCIP_OKAY; /*lint !e438*/
2300 }
2301 
2302 /** registers a variable from a violated constraint as branching candidate that has a large absolute value in the relaxation */
2303 static
2305  SCIP* scip, /**< SCIP data structure */
2306  SCIP_CONS** conss, /**< constraints */
2307  int nconss, /**< number of constraints */
2308  SCIP_SOL* sol, /**< solution to enforce (NULL for the LP solution) */
2309  SCIP_VAR** brvar /**< buffer to store branching variable */
2310  )
2311 {
2312  SCIP_CONSDATA* consdata;
2313  SCIP_Real val;
2314  SCIP_Real brvarval;
2315  int c;
2316 
2317  assert(scip != NULL);
2318  assert(conss != NULL || nconss == 0);
2319 
2320  *brvar = NULL;
2321  brvarval = -1.0;
2322 
2323  for( c = 0; c < nconss; ++c )
2324  {
2325  assert(conss != NULL);
2326  consdata = SCIPconsGetData(conss[c]);
2327  assert(consdata != NULL);
2328 
2329  /* skip constraints that have been marked to be removed by propagateCons() */
2330  if( !SCIPconsIsEnabled(conss[c]) )
2331  continue;
2332 
2333  if( !SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) && !SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)) )
2334  continue;
2335 
2336  val = SCIPgetSolVal(scip, sol, consdata->x) + consdata->xoffset;
2337  if( REALABS(val) > brvarval )
2338  {
2339  brvarval = ABS(val);
2340  *brvar = consdata->x;
2341  }
2342  }
2343 
2344  if( *brvar != NULL )
2345  {
2346  SCIP_CALL( SCIPaddExternBranchCand(scip, *brvar, brvarval, SCIP_INVALID) );
2347  }
2348 
2349  return SCIP_OKAY;
2350 }
2351 
2352 /* try to fix almost fixed x variable in violated constraint */
2353 static
2355  SCIP* scip, /**< SCIP data structure */
2356  SCIP_CONS** conss, /**< constraints */
2357  int nconss, /**< number of constraints */
2358  SCIP_Bool* infeasible, /**< buffer to store whether infeasibility was detected */
2359  SCIP_Bool* reduceddom /**< buffer to store whether some variable bound was tightened */
2360  )
2361 {
2362  SCIP_CONSDATA* consdata;
2363  SCIP_Real lb;
2364  SCIP_Real ub;
2365  SCIP_Bool tightened;
2366  int c;
2367 
2368  assert(scip != NULL);
2369  assert(conss != NULL);
2370  assert(infeasible != NULL);
2371  assert(reduceddom != NULL);
2372 
2373  *infeasible = FALSE;
2374  *reduceddom = FALSE;
2375 
2376  for( c = 0; c < nconss; ++c )
2377  {
2378  consdata = SCIPconsGetData(conss[c]);
2379  assert(consdata != NULL);
2380 
2381  /* if constraint not violated, then continue */
2382  if( !SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) && !SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)) )
2383  continue;
2384 
2385  lb = SCIPvarGetLbLocal(consdata->x);
2386  ub = SCIPvarGetUbLocal(consdata->x);
2387 
2388  /* if x not almost fixed, then continue */
2389  if( !SCIPisRelEQ(scip, lb, ub) )
2390  continue;
2391 
2392  /* if x fixed already, then continue */
2393  if( SCIPisEQ(scip, lb, ub) )
2394  continue;
2395 
2396  assert(!SCIPisInfinity(scip, -lb));
2397  assert(!SCIPisInfinity(scip, ub));
2398 
2399  /* try to fix variable */
2400  SCIP_CALL( SCIPtightenVarLb(scip, consdata->x, (lb+ub)/2.0, TRUE, infeasible, &tightened) );
2401  if( *infeasible )
2402  {
2403  SCIPdebugMsg(scip, "Fixing almost fixed variable <%s> lead to infeasibility.\n", SCIPvarGetName(consdata->x));
2404  return SCIP_OKAY;
2405  }
2406  if( tightened )
2407  {
2408  SCIPdebugMsg(scip, "Tightened lower bound of almost fixed variable <%s>.\n", SCIPvarGetName(consdata->x));
2409  *reduceddom = TRUE;
2410  }
2411 
2412  SCIP_CALL( SCIPtightenVarUb(scip, consdata->x, (lb+ub)/2.0, TRUE, infeasible, &tightened) );
2413  if( *infeasible )
2414  {
2415  SCIPdebugMsg(scip, "Fixing almost fixed variable <%s> lead to infeasibility.\n", SCIPvarGetName(consdata->x));
2416  return SCIP_OKAY;
2417  }
2418  if( tightened )
2419  {
2420  SCIPdebugMsg(scip, "Tightened upper bound of almost fixed variable <%s>.\n", SCIPvarGetName(consdata->x));
2421  *reduceddom = TRUE;
2422  }
2423 
2424  /* stop as soon as one variable has been fixed to start another enfo round */
2425  if( *reduceddom )
2426  break;
2427  }
2428 
2429  return SCIP_OKAY;
2430 }
2431 
2432 /** resolves a propagation on the given variable by supplying the variables needed for applying the corresponding
2433  * propagation rule (see propagateCons()):
2434  * see cons_varbound
2435  */
2436 static
2438  SCIP* scip, /**< SCIP data structure */
2439  SCIP_CONS* cons, /**< constraint that inferred the bound change */
2440  SCIP_VAR* infervar, /**< variable that was deduced */
2441  PROPRULE proprule, /**< propagation rule that deduced the bound change */
2442  SCIP_BOUNDTYPE boundtype, /**< the type of the changed bound (lower or upper bound) */
2443  SCIP_BDCHGIDX* bdchgidx /**< bound change index (time stamp of bound change), or NULL for current time */
2444  )
2445 {
2446  SCIP_CONSDATA* consdata;
2447 
2448  assert(scip != NULL);
2449  assert(cons != NULL);
2450  assert(infervar != NULL);
2451 
2452  consdata = SCIPconsGetData(cons);
2453  assert(consdata != NULL);
2454  assert(consdata->zcoef != 0.0);
2455 
2456  switch( proprule )
2457  {
2458  case PROPRULE_1:
2459  /* lhs <= sign(x+offset)|x+offset|^n + c*z: left hand side and bounds on z -> lower bound on x */
2460  assert(infervar == consdata->x);
2461  assert(boundtype == SCIP_BOUNDTYPE_LOWER);
2462  assert(!SCIPisInfinity(scip, -consdata->lhs));
2463  if( consdata->zcoef > 0.0 )
2464  {
2465  SCIP_CALL( SCIPaddConflictUb(scip, consdata->z, bdchgidx) );
2466  }
2467  else
2468  {
2469  SCIP_CALL( SCIPaddConflictLb(scip, consdata->z, bdchgidx) );
2470  }
2471  break;
2472 
2473  case PROPRULE_2:
2474  /* lhs <= sign(x+offset)|x+offset|^n + c*z: left hand side and upper bound on x -> bound on z */
2475  assert(infervar == consdata->z);
2476  assert(!SCIPisInfinity(scip, -consdata->lhs));
2477  SCIP_CALL( SCIPaddConflictUb(scip, consdata->x, bdchgidx) );
2478  break;
2479 
2480  case PROPRULE_3:
2481  /* sign(x+offset)|x+offset|^n + c*z <= rhs: right hand side and bounds on z -> upper bound on x */
2482  assert(infervar == consdata->x);
2483  assert(boundtype == SCIP_BOUNDTYPE_UPPER);
2484  assert(!SCIPisInfinity(scip, consdata->rhs));
2485  if( consdata->zcoef > 0.0 )
2486  {
2487  SCIP_CALL( SCIPaddConflictLb(scip, consdata->z, bdchgidx) );
2488  }
2489  else
2490  {
2491  SCIP_CALL( SCIPaddConflictUb(scip, consdata->z, bdchgidx) );
2492  }
2493  break;
2494 
2495  case PROPRULE_4:
2496  /* sign(x+offset)|x+offset|^n + c*z <= rhs: right hand side and lower bound on x -> bound on z */
2497  assert(infervar == consdata->z);
2498  assert(!SCIPisInfinity(scip, consdata->rhs));
2499  SCIP_CALL( SCIPaddConflictLb(scip, consdata->x, bdchgidx) );
2500  break;
2501 
2502  case PROPRULE_INVALID:
2503  default:
2504  SCIPerrorMessage("invalid inference information %d in absolute power constraint <%s>\n", proprule, SCIPconsGetName(cons));
2505  return SCIP_INVALIDDATA;
2506  }
2507 
2508  return SCIP_OKAY;
2509 }
2510 
2511 /** analyze infeasibility */
2512 static
2514  SCIP* scip, /**< SCIP data structure */
2515  SCIP_CONS* cons, /**< variable bound constraint */
2516  SCIP_VAR* infervar, /**< variable that was deduced */
2517  PROPRULE proprule, /**< propagation rule that deduced the bound change */
2518  SCIP_BOUNDTYPE boundtype /**< the type of the changed bound (lower or upper bound) */
2519  )
2520 {
2521  /* conflict analysis can only be applied in solving stage and if it turned on */
2523  return SCIP_OKAY;
2524 
2525  /* initialize conflict analysis, and add all variables of infeasible constraint to conflict candidate queue */
2527 
2528  /* add the bound which got violated */
2529  if( boundtype == SCIP_BOUNDTYPE_LOWER )
2530  {
2531  SCIP_CALL( SCIPaddConflictUb(scip, infervar, NULL) );
2532  }
2533  else
2534  {
2535  assert(boundtype == SCIP_BOUNDTYPE_UPPER);
2536  SCIP_CALL( SCIPaddConflictLb(scip, infervar, NULL) );
2537  }
2538 
2539  /* add the reason for the violated of the bound */
2540  SCIP_CALL( resolvePropagation(scip, cons, infervar, proprule, boundtype, NULL) );
2541 
2542  /* analyze the conflict */
2543  SCIP_CALL( SCIPanalyzeConflictCons(scip, cons, NULL) );
2544 
2545  return SCIP_OKAY;
2546 }
2547 
2548 /** propagation method for absolute power constraint
2549  * SCIPinferVarXbCons to allow for repropagation
2550  */
2551 static
2553  SCIP* scip, /**< SCIP data structure */
2554  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
2555  SCIP_CONS* cons, /**< variable bound constraint */
2556  SCIP_Bool canaddcons, /**< are we allowed to add a linear constraint when enforcing bounds for a multiaggregated variable? */
2557  SCIP_Bool* cutoff, /**< pointer to store whether the node can be cut off */
2558  int* nchgbds, /**< pointer to count number of bound changes */
2559  int* naddconss /**< pointer to count number of added constraints */
2560  )
2561 {
2562  SCIP_CONSDATA* consdata;
2563  SCIP_Real xlb;
2564  SCIP_Real xub;
2565  SCIP_Real zlb;
2566  SCIP_Real zub;
2567  SCIP_Real newlb;
2568  SCIP_Real newub;
2569  SCIP_Bool tightened;
2570  SCIP_Bool tightenedround;
2571  SCIP_Real minact;
2572  SCIP_Real maxact;
2573 
2574  assert(conshdlr != NULL);
2575  assert(cutoff != NULL);
2576  assert(nchgbds != NULL);
2577  assert(naddconss != NULL);
2578 
2579  consdata = SCIPconsGetData(cons);
2580  assert(consdata != NULL);
2581 
2582  SCIPdebugMsg(scip, "propagating absolute power constraint <%s>\n", SCIPconsGetName(cons));
2583 
2584  *cutoff = FALSE;
2585 
2586  /* get current bounds of variables */
2587  xlb = SCIPvarGetLbLocal(consdata->x);
2588  xub = SCIPvarGetUbLocal(consdata->x);
2589  zlb = SCIPvarGetLbLocal(consdata->z);
2590  zub = SCIPvarGetUbLocal(consdata->z);
2591 
2592  /* if some bound is not tightened, tighten bounds of variables as long as possible */
2593  tightenedround = SCIPconsIsMarkedPropagate(cons);
2594  while( tightenedround )
2595  {
2596  tightenedround = FALSE;
2597 
2598  /* propagate left hand side inequality: lhs <= (x+offset)*|x+offset|^n + c*z */
2599  if( !SCIPisInfinity(scip, -consdata->lhs) )
2600  {
2601  assert(!*cutoff);
2602 
2603  /* propagate bounds on x (if not multiaggregated):
2604  * (1) left hand side and bounds on z -> lower bound on x
2605  */
2606  if( SCIPvarIsActive(SCIPvarGetProbvar(consdata->x)) && (!SCIPisFeasEQ(scip, zlb, zub) || !SCIPisInfinity(scip, REALABS(zlb))) )
2607  {
2608  /* if z is fixed, first compute new lower bound on x without tolerances
2609  * if that is feasible, project new lower bound onto current bounds
2610  * otherwise, recompute with tolerances and continue as usual
2611  * do this only if variable is not essentially fixed to value of infinity
2612  */
2613  if( SCIPisFeasEQ(scip, zlb, zub) && !SCIPisInfinity(scip, zub) )
2614  {
2615  assert(!SCIPisInfinity(scip, -zlb));
2616 
2617  newlb = consdata->lhs - consdata->zcoef * (consdata->zcoef > 0.0 ? zub : zlb);
2618 
2619  /* invert sign(x+offset)|x+offset|^(n-1) = y -> x = sign(y)|y|^(1/n) - offset */
2620  if( consdata->exponent == 2.0 )
2621  newlb = SIGN(newlb) * sqrt(ABS(newlb));
2622  else
2623  newlb = SIGN(newlb) * pow(ABS(newlb), 1.0/consdata->exponent);
2624  newlb -= consdata->xoffset;
2625 
2626  if( SCIPisFeasGT(scip, newlb, xub) )
2627  {
2628  /* if new lower bound for x would yield cutoff, recompute with tolerances */
2629  newlb = consdata->lhs - PROPSIDETOL - consdata->zcoef * (consdata->zcoef > 0.0 ? (zub + PROPVARTOL) : (zlb - PROPVARTOL));
2630 
2631  /* invert sign(x+offset)|x+offset|^(n-1) = y -> x = sign(y)|y|^(1/n) - offset */
2632  if( consdata->exponent == 2.0 )
2633  newlb = SIGN(newlb) * sqrt(ABS(newlb));
2634  else
2635  newlb = SIGN(newlb) * pow(ABS(newlb), 1.0/consdata->exponent);
2636  newlb -= consdata->xoffset;
2637  }
2638  else
2639  {
2640  /* project new lower bound onto current bounds */
2641  newlb = MIN(newlb, xub);
2642  }
2643  }
2644  else
2645  {
2646  if( consdata->zcoef > 0.0 )
2647  {
2648  if( !SCIPisInfinity(scip, zub) )
2649  newlb = consdata->lhs - PROPSIDETOL - consdata->zcoef * (zub + PROPVARTOL);
2650  else
2651  newlb = -SCIPinfinity(scip);
2652  }
2653  else
2654  {
2655  if( !SCIPisInfinity(scip, -zlb) )
2656  newlb = consdata->lhs - PROPSIDETOL - consdata->zcoef * (zlb - PROPVARTOL);
2657  else
2658  newlb = -SCIPinfinity(scip);
2659  }
2660 
2661  if( !SCIPisInfinity(scip, -newlb) )
2662  {
2663  /* invert sign(x+offset)|x+offset|^(n-1) = y -> x = sign(y)|y|^(1/n) - offset */
2664  if( consdata->exponent == 2.0 )
2665  newlb = SIGN(newlb) * sqrt(ABS(newlb));
2666  else
2667  newlb = SIGN(newlb) * pow(ABS(newlb), 1.0/consdata->exponent);
2668  newlb -= consdata->xoffset;
2669  }
2670  }
2671 
2672  if( SCIPisInfinity(scip, newlb) )
2673  {
2674  /* we cannot fix a variable to +infinity, so let's report cutoff (there is no solution within SCIPs limitations to infinity) */
2675  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g] -> cutoff\n", SCIPvarGetName(consdata->x), xlb, xub, newlb, xub);
2676 
2677  *cutoff = TRUE;
2678 
2679  /* analyze infeasibility */
2680  SCIP_CALL( analyzeConflict(scip, cons, consdata->x, PROPRULE_1, SCIP_BOUNDTYPE_LOWER) );
2681  break;
2682  }
2683 
2684  if( !SCIPisInfinity(scip, -newlb) )
2685  {
2686  if( SCIPisLbBetter(scip, newlb, xlb, xub) )
2687  {
2688  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g]\n",
2689  SCIPvarGetName(consdata->x), xlb, xub, newlb, xub);
2690  SCIP_CALL( SCIPinferVarLbCons(scip, consdata->x, newlb, cons, (int)PROPRULE_1, FALSE, cutoff, &tightened) );
2691 
2692  if( *cutoff )
2693  {
2694  assert(SCIPisInfinity(scip, newlb) || SCIPisGT(scip, newlb, SCIPvarGetUbLocal(consdata->x)));
2695 
2696  /* analyze infeasibility */
2697  SCIP_CALL( analyzeConflict(scip, cons, consdata->x, PROPRULE_1, SCIP_BOUNDTYPE_LOWER) );
2698  break;
2699  }
2700 
2701  if( tightened )
2702  {
2703  tightenedround = TRUE;
2704  (*nchgbds)++;
2705  }
2706  xlb = SCIPvarGetLbLocal(consdata->x);
2707  }
2708  }
2709  }
2710 
2711  assert(!*cutoff);
2712 
2713  /* propagate bounds on z:
2714  * (2) left hand side and upper bound on x -> bound on z
2715  */
2716  if( SCIPvarGetStatus(consdata->z) != SCIP_VARSTATUS_MULTAGGR && !SCIPisInfinity(scip, xub) ) /* cannot change bounds of multaggr vars */
2717  {
2718  SCIP_Real newbd;
2719 
2720  /* if x is fixed, first compute new bound on z without tolerances
2721  * if that is feasible, project new bound onto current bounds
2722  * otherwise, recompute with tolerances and continue as usual
2723  */
2724  if( SCIPisFeasEQ(scip, xlb, xub) )
2725  {
2726  newbd = xub + consdata->xoffset;
2727  newbd = consdata->lhs - SIGN(newbd) * consdata->power(REALABS(newbd), consdata->exponent);
2728  newbd /= consdata->zcoef;
2729 
2730  if( SCIPisInfinity(scip, newbd) )
2731  newbd = SCIPinfinity(scip);
2732  else if( SCIPisInfinity(scip, -newbd) )
2733  newbd = -SCIPinfinity(scip);
2734 
2735  if( (consdata->zcoef > 0.0 && SCIPisFeasGT(scip, newbd, zub)) || (consdata->zcoef < 0.0 && SCIPisFeasLT(scip, newbd, zlb)) )
2736  {
2737  /* if infeasible, recompute with tolerances */
2738  newbd = xub + PROPVARTOL + consdata->xoffset;
2739  newbd = consdata->lhs - PROPSIDETOL - SIGN(newbd) * consdata->power(REALABS(newbd), consdata->exponent);
2740  newbd /= consdata->zcoef;
2741  }
2742  else
2743  {
2744  /* project onto current bounds of z */
2745  newbd = MIN(zub, MAX(zlb, newbd) );
2746  }
2747  }
2748  else
2749  {
2750  newbd = xub + PROPVARTOL + consdata->xoffset;
2751  newbd = consdata->lhs - PROPSIDETOL - SIGN(newbd) * consdata->power(REALABS(newbd), consdata->exponent);
2752  newbd /= consdata->zcoef;
2753  }
2754 
2755  if( consdata->zcoef > 0.0 )
2756  {
2757  newlb = newbd;
2758 
2759  if( SCIPisInfinity(scip, newlb) )
2760  {
2761  /* we cannot fix a variable to +infinity, so let's report cutoff (there is no solution within SCIPs limitations to infinity) */
2762  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g] -> cutoff\n", SCIPvarGetName(consdata->z), zlb, zub, newlb, zub);
2763 
2764  *cutoff = TRUE;
2765 
2766  /* analyze infeasibility */
2767  SCIP_CALL( analyzeConflict(scip, cons, consdata->z, PROPRULE_2, SCIP_BOUNDTYPE_LOWER) );
2768  break;
2769  }
2770 
2771  if( SCIPisLbBetter(scip, newlb, zlb, zub) )
2772  {
2773  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g]\n",
2774  SCIPvarGetName(consdata->z), zlb, zub, newlb, zub);
2775  SCIP_CALL( SCIPinferVarLbCons(scip, consdata->z, newlb, cons, (int)PROPRULE_2, FALSE, cutoff, &tightened) );
2776 
2777  if( *cutoff )
2778  {
2779  assert(SCIPisInfinity(scip, newlb) || SCIPisGT(scip, newlb, SCIPvarGetUbLocal(consdata->z)));
2780 
2781  /* analyze infeasibility */
2782  SCIP_CALL( analyzeConflict(scip, cons, consdata->z, PROPRULE_2, SCIP_BOUNDTYPE_LOWER) );
2783  break;
2784  }
2785 
2786  if( tightened )
2787  {
2788  tightenedround = TRUE;
2789  (*nchgbds)++;
2790  }
2791  zlb = SCIPvarGetLbLocal(consdata->z);
2792  }
2793  }
2794  else
2795  {
2796  newub = newbd;
2797 
2798  if( SCIPisInfinity(scip, -newub) )
2799  {
2800  /* we cannot fix a variable to -infinity, so let's report cutoff (there is no solution within SCIPs limitations to infinity) */
2801  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g] -> cutoff\n", SCIPvarGetName(consdata->z), zlb, zub, zlb, newub);
2802 
2803  *cutoff = TRUE;
2804 
2805  /* analyze infeasibility */
2806  SCIP_CALL( analyzeConflict(scip, cons, consdata->z, PROPRULE_2, SCIP_BOUNDTYPE_UPPER) );
2807  break;
2808  }
2809 
2810  if( SCIPisUbBetter(scip, newub, zlb, zub) )
2811  {
2812  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g]\n",
2813  SCIPvarGetName(consdata->z), zlb, zub, zlb, newub);
2814  SCIP_CALL( SCIPinferVarUbCons(scip, consdata->z, newub, cons, (int)PROPRULE_2, FALSE, cutoff, &tightened) );
2815 
2816  if( *cutoff )
2817  {
2818  assert(SCIPisInfinity(scip, -newub) || SCIPisLT(scip, newub, SCIPvarGetLbLocal(consdata->z)));
2819 
2820  /* analyze infeasibility */
2821  SCIP_CALL( analyzeConflict(scip, cons, consdata->z, PROPRULE_2, SCIP_BOUNDTYPE_UPPER) );
2822  break;
2823  }
2824 
2825  if( tightened )
2826  {
2827  tightenedround = TRUE;
2828  (*nchgbds)++;
2829  }
2830  zub = SCIPvarGetUbLocal(consdata->z);
2831  }
2832  }
2833  }
2834  }
2835 
2836  assert(!*cutoff);
2837 
2838  /* propagate right hand side inequality: sign(x+offset)|x+offset|^n + c*z <= rhs */
2839  if( !SCIPisInfinity(scip, consdata->rhs) )
2840  {
2841  /* propagate bounds on x:
2842  * (3) right hand side and bounds on z -> upper bound on x
2843  */
2844  if( SCIPvarIsActive(SCIPvarGetProbvar(consdata->x)) && (!SCIPisFeasEQ(scip, zlb, zub) || !SCIPisInfinity(scip, REALABS(zlb))) ) /* cannot change bounds of multaggr or fixed vars */
2845  {
2846  /* if z is fixed, first compute new upper bound on x without tolerances
2847  * if that is feasible, project new upper bound onto current bounds
2848  * otherwise, recompute with tolerances and continue as usual
2849  * do this only if variable is not essentially fixed to value of infinity
2850  */
2851  if( SCIPisFeasEQ(scip, zlb, zub) && !SCIPisInfinity(scip, zub) )
2852  {
2853  assert(!SCIPisInfinity(scip, -zlb));
2854 
2855  newub = consdata->rhs - consdata->zcoef * (consdata->zcoef > 0.0 ? zlb : zub);
2856 
2857  /* invert sign(x+offset)|x+offset|^(n-1) = y -> x = sign(y)|y|^(1/n) - offset */
2858  if( consdata->exponent == 2.0 )
2859  newub = SIGN(newub) * sqrt(ABS(newub));
2860  else
2861  newub = SIGN(newub) * pow(ABS(newub), 1.0/consdata->exponent);
2862  newub -= consdata->xoffset;
2863 
2864  if( SCIPisFeasLT(scip, newub, xlb) )
2865  {
2866  /* if new lower bound for x would yield cutoff, recompute with tolerances */
2867  newub = consdata->rhs + PROPSIDETOL - consdata->zcoef * (consdata->zcoef > 0.0 ? (zlb - PROPVARTOL) : (zub + PROPVARTOL));
2868 
2869  /* invert sign(x+offset)|x+offset|^(n-1) = y -> x = sign(y)|y|^(1/n) - offset */
2870  if( consdata->exponent == 2.0 )
2871  newub = SIGN(newub) * sqrt(ABS(newub));
2872  else
2873  newub = SIGN(newub) * pow(ABS(newub), 1.0/consdata->exponent);
2874  newub -= consdata->xoffset;
2875  }
2876  else
2877  {
2878  /* project new upper bound onto current bounds */
2879  newub = MAX(newub, xlb);
2880  }
2881  }
2882  else
2883  {
2884  if( consdata->zcoef > 0.0 )
2885  {
2886  if( !SCIPisInfinity(scip, -zlb) )
2887  newub = consdata->rhs + PROPSIDETOL - consdata->zcoef * (zlb - PROPVARTOL);
2888  else
2889  newub = SCIPinfinity(scip);
2890  }
2891  else
2892  {
2893  if( !SCIPisInfinity(scip, zub) )
2894  newub = consdata->rhs + PROPSIDETOL - consdata->zcoef * (zub + PROPVARTOL);
2895  else
2896  newub = SCIPinfinity(scip);
2897  }
2898  if( !SCIPisInfinity(scip, newub) )
2899  {
2900  /* invert sign(x+offset)|x+offset|^(n-1) = y -> x = sign(y)|y|^(1/n) - offset */
2901  if( consdata->exponent == 2.0 )
2902  newub = SIGN(newub) * sqrt(ABS(newub));
2903  else
2904  newub = SIGN(newub) * pow(ABS(newub), 1.0/consdata->exponent);
2905  newub -= consdata->xoffset;
2906  }
2907  }
2908 
2909  if( SCIPisInfinity(scip, -newub) )
2910  {
2911  /* we cannot fix a variable to -infinity, so let's report cutoff (there is no solution within SCIPs limitations to infinity) */
2912  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g] -> cutoff\n", SCIPvarGetName(consdata->x), xlb, xub, xlb, newub);
2913 
2914  *cutoff = TRUE;
2915 
2916  /* analyze infeasibility */
2917  SCIP_CALL( analyzeConflict(scip, cons, consdata->x, PROPRULE_3, SCIP_BOUNDTYPE_UPPER) );
2918  break;
2919  }
2920 
2921  if( !SCIPisInfinity(scip, newub) )
2922  {
2923  if( SCIPisUbBetter(scip, newub, xlb, xub) )
2924  {
2925  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g]\n",
2926  SCIPvarGetName(consdata->x), xlb, xub, xlb, newub);
2927  SCIP_CALL( SCIPinferVarUbCons(scip, consdata->x, newub, cons, (int)PROPRULE_3, FALSE, cutoff, &tightened) );
2928 
2929  if( *cutoff )
2930  {
2931  assert(SCIPisInfinity(scip, -newub) || SCIPisLT(scip, newub, SCIPvarGetLbLocal(consdata->x)));
2932 
2933  /* analyze infeasibility */
2934  SCIP_CALL( analyzeConflict(scip, cons, consdata->x, PROPRULE_3, SCIP_BOUNDTYPE_UPPER) );
2935  break;
2936  }
2937 
2938  if( tightened )
2939  {
2940  tightenedround = TRUE;
2941  (*nchgbds)++;
2942  }
2943  xub = SCIPvarGetUbLocal(consdata->x);
2944  }
2945  }
2946  }
2947 
2948  assert(!*cutoff);
2949 
2950  /* propagate bounds on z:
2951  * (4) right hand side and lower bound on x -> bound on z
2952  */
2953  if( SCIPvarGetStatus(consdata->z) != SCIP_VARSTATUS_MULTAGGR && !SCIPisInfinity(scip, -xlb) ) /* cannot change bounds of multaggr vars */
2954  {
2955  SCIP_Real newbd;
2956 
2957  /* if x is fixed, first compute new bound on z without tolerances
2958  * if that is feasible, project new bound onto current bounds
2959  * otherwise, recompute with tolerances and continue as usual
2960  */
2961  if( SCIPisFeasEQ(scip, xlb, xub) )
2962  {
2963  newbd = xlb + consdata->xoffset;
2964  newbd = consdata->rhs - SIGN(newbd) * consdata->power(REALABS(newbd), consdata->exponent);
2965  newbd /= consdata->zcoef;
2966 
2967  if( SCIPisInfinity(scip, newbd) )
2968  newbd = SCIPinfinity(scip);
2969  else if( SCIPisInfinity(scip, -newbd) )
2970  newbd = -SCIPinfinity(scip);
2971 
2972  if( (consdata->zcoef > 0.0 && SCIPisFeasLT(scip, newbd, zlb)) || (consdata->zcoef < 0.0 && SCIPisFeasGT(scip, newbd, zub)) )
2973  {
2974  /* if infeasible, recompute with tolerances */
2975  newbd = xlb - PROPVARTOL + consdata->xoffset;
2976  newbd = consdata->rhs + PROPSIDETOL - SIGN(newbd) * consdata->power(REALABS(newbd), consdata->exponent);
2977  newbd /= consdata->zcoef;
2978  }
2979  else
2980  {
2981  /* project onto current bounds of z */
2982  newbd = MIN(zub, MAX(zlb, newbd) );
2983  }
2984  }
2985  else
2986  {
2987  newbd = xlb - PROPVARTOL + consdata->xoffset;
2988  newbd = consdata->rhs + PROPSIDETOL - SIGN(newbd) * consdata->power(REALABS(newbd), consdata->exponent);
2989  newbd /= consdata->zcoef;
2990  }
2991 
2992  if( consdata->zcoef > 0.0 )
2993  {
2994  newub = newbd;
2995 
2996  if( SCIPisInfinity(scip, -newub) )
2997  {
2998  /* we cannot fix a variable to -infinity, so let's report cutoff (there is no solution within SCIPs limitations to infinity) */
2999  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g] -> cutoff\n", SCIPvarGetName(consdata->z), zlb, zub, zlb, newub);
3000 
3001  *cutoff = TRUE;
3002 
3003  /* analyze infeasibility */
3004  SCIP_CALL( analyzeConflict(scip, cons, consdata->z, PROPRULE_4, SCIP_BOUNDTYPE_UPPER) );
3005  break;
3006  }
3007 
3008  if( SCIPisUbBetter(scip, newub, zlb, zub) )
3009  {
3010  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g]\n",
3011  SCIPvarGetName(consdata->z), zlb, zub, zlb, newub);
3012  SCIP_CALL( SCIPinferVarUbCons(scip, consdata->z, newub, cons, (int)PROPRULE_4, FALSE, cutoff, &tightened) );
3013 
3014  if( *cutoff )
3015  {
3016  assert(SCIPisInfinity(scip, -newub) || SCIPisLT(scip, newub, SCIPvarGetLbLocal(consdata->z)));
3017 
3018  /* analyze infeasibility */
3019  SCIP_CALL( analyzeConflict(scip, cons, consdata->z, PROPRULE_4, SCIP_BOUNDTYPE_UPPER) );
3020  break;
3021  }
3022 
3023  if( tightened )
3024  {
3025  tightenedround = TRUE;
3026  (*nchgbds)++;
3027  }
3028  zub = SCIPvarGetUbLocal(consdata->z);
3029  }
3030  }
3031  else
3032  {
3033  newlb = newbd;
3034 
3035  if( SCIPisInfinity(scip, newlb) )
3036  {
3037  /* we cannot fix a variable to +infinity, so let's report cutoff (there is no solution within SCIPs limitations to infinity) */
3038  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g] -> cutoff\n", SCIPvarGetName(consdata->z), zlb, zub, newlb, zub);
3039 
3040  *cutoff = TRUE;
3041 
3042  /* analyze infeasibility */
3043  SCIP_CALL( analyzeConflict(scip, cons, consdata->z, PROPRULE_4, SCIP_BOUNDTYPE_LOWER) );
3044  break;
3045  }
3046 
3047  if( SCIPisLbBetter(scip, newlb, zlb, zub) )
3048  {
3049  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g]\n",
3050  SCIPvarGetName(consdata->z), zlb, zub, newlb, zub);
3051  SCIP_CALL( SCIPinferVarLbCons(scip, consdata->z, newlb, cons, (int)PROPRULE_4, FALSE, cutoff, &tightened) );
3052 
3053  if( *cutoff )
3054  {
3055  assert(SCIPisInfinity(scip, newlb) || SCIPisGT(scip, newlb, SCIPvarGetUbLocal(consdata->z)));
3056 
3057  /* analyze infeasibility */
3058  SCIP_CALL( analyzeConflict(scip, cons, consdata->z, PROPRULE_4, SCIP_BOUNDTYPE_LOWER) );
3059  break;
3060  }
3061 
3062  if( tightened )
3063  {
3064  tightenedround = TRUE;
3065  (*nchgbds)++;
3066  }
3067  zlb = SCIPvarGetLbLocal(consdata->z);
3068  }
3069  }
3070  }
3071  }
3072 
3073  assert(!*cutoff);
3074  }
3075 
3076  /* mark the constraint propagated */
3077  SCIP_CALL( SCIPunmarkConsPropagate(scip, cons) );
3078 
3079  if( *cutoff )
3080  return SCIP_OKAY;
3081 
3082  /* check for redundancy */
3083  if( !SCIPisInfinity(scip, -xlb) && !SCIPisInfinity(scip, consdata->zcoef > 0.0 ? -zlb : zub) )
3084  minact = SIGN(xlb + consdata->xoffset) * consdata->power(REALABS(xlb + consdata->xoffset), consdata->exponent) + consdata->zcoef * (consdata->zcoef > 0.0 ? zlb : zub);
3085  else
3086  minact = -SCIPinfinity(scip);
3087 
3088  if( !SCIPisInfinity(scip, xub) && !SCIPisInfinity(scip, consdata->zcoef > 0.0 ? zub : -zlb) )
3089  maxact = SIGN(xub + consdata->xoffset) * consdata->power(REALABS(xub + consdata->xoffset), consdata->exponent) + consdata->zcoef * (consdata->zcoef > 0.0 ? zub : zlb);
3090  else
3091  maxact = SCIPinfinity(scip);
3092 
3093  if( (SCIPisInfinity(scip, -consdata->lhs) || SCIPisGE(scip, minact, consdata->lhs)) &&
3094  (SCIPisInfinity(scip, consdata->rhs) || SCIPisLE(scip, maxact, consdata->rhs)) )
3095  {
3096  SCIPdebugMsg(scip, "absolute power constraint <%s> is redundant: <%s>[%.15g,%.15g], <%s>[%.15g,%.15g]\n",
3097  SCIPconsGetName(cons),
3098  SCIPvarGetName(consdata->x), SCIPvarGetLbLocal(consdata->x), SCIPvarGetUbLocal(consdata->x),
3099  SCIPvarGetName(consdata->z), SCIPvarGetLbLocal(consdata->z), SCIPvarGetUbLocal(consdata->z));
3100 
3101  SCIP_CALL( SCIPdelConsLocal(scip, cons) );
3102 
3103  return SCIP_OKAY;
3104  }
3105 
3106  /* delete constraint if x has been fixed */
3107  if( SCIPisRelEQ(scip, xlb, xub) && (SCIPvarIsActive(consdata->z) || canaddcons) )
3108  {
3109  SCIP_RESULT tightenresult;
3110  SCIP_INTERVAL xbnds;
3111  SCIP_INTERVAL zbnds;
3112 
3113  SCIPdebugMsg(scip, "x-variable in constraint <%s> is fixed: x = <%s>[%.15g,%.15g], z = <%s>[%.15g,%.15g]\n",
3114  SCIPconsGetName(cons), SCIPvarGetName(consdata->x), xlb, xub, SCIPvarGetName(consdata->z), zlb, zub);
3115 
3116  SCIPintervalSetBounds(&xbnds, MIN(xlb, xub), MAX(xlb, xub));
3117  computeBoundsZ(scip, cons, xbnds, &zbnds);
3118 
3119  /* in difference to the loop above, here we enforce a possible bound tightening on z, and may add a linear cons if z is multiaggregated */
3120  SCIP_CALL( tightenBounds(scip, consdata->z, zbnds, TRUE, cons, &tightenresult, nchgbds, nchgbds, naddconss) );
3121  if( tightenresult == SCIP_CUTOFF )
3122  *cutoff = TRUE;
3123 
3124  SCIP_CALL( SCIPdelConsLocal(scip, cons) );
3125 
3126  return SCIP_OKAY;
3127  }
3128 
3129  /* delete constraint if z has been fixed */
3130  if( SCIPisRelEQ(scip, zlb, zub) && (SCIPvarIsActive(consdata->x) || canaddcons) )
3131  {
3132  SCIP_RESULT tightenresult;
3133  SCIP_INTERVAL xbnds;
3134  SCIP_INTERVAL zbnds;
3135 
3136  SCIPdebugMsg(scip, "z-variable in constraint <%s> is fixed: x = <%s>[%.15g,%.15g], z = <%s>[%.15g,%.15g]\n",
3137  SCIPconsGetName(cons), SCIPvarGetName(consdata->x), xlb, xub, SCIPvarGetName(consdata->z), zlb, zub);
3138 
3139  SCIPintervalSetBounds(&zbnds, MIN(zlb, zub), MAX(zlb, zub));
3140  computeBoundsX(scip, cons, zbnds, &xbnds);
3141 
3142  /* in difference to the loop above, here we enforce a possible bound tightening on x, and may add a linear cons if x is multiaggregated */
3143  SCIP_CALL( tightenBounds(scip, consdata->x, xbnds, TRUE, cons, &tightenresult, nchgbds, nchgbds, naddconss) );
3144  if( tightenresult == SCIP_CUTOFF )
3145  *cutoff = TRUE;
3146 
3147  SCIP_CALL( SCIPdelConsLocal(scip, cons) );
3148 
3149  return SCIP_OKAY;
3150  }
3151 
3152  return SCIP_OKAY;
3153 }
3154 
3155 /** notifies SCIP about a variable bound lhs <= x + c*y <= rhs */
3156 static
3158  SCIP* scip, /**< SCIP data structure */
3159  SCIP_CONS* cons, /**< absolute power constraint this variable bound is derived form */
3160  SCIP_Bool addcons, /**< should the variable bound be added as constraint to SCIP? */
3161  SCIP_VAR* var, /**< variable x for which we want to add a variable bound */
3162  SCIP_VAR* vbdvar, /**< variable y which makes the bound a variable bound */
3163  SCIP_Real vbdcoef, /**< coefficient c of bounding variable vbdvar */
3164  SCIP_Real lhs, /**< left hand side of varbound constraint */
3165  SCIP_Real rhs, /**< right hand side of varbound constraint */
3166  SCIP_Bool* infeas, /**< pointer to store whether an infeasibility was detected */
3167  int* nbdchgs, /**< pointer where to add number of performed bound changes */
3168  int* naddconss /**< pointer where to add number of added constraints */
3169  )
3170 {
3171  int nbdchgs_local;
3172 
3173  assert(scip != NULL);
3174  assert(cons != NULL);
3175  assert(var != NULL);
3176  assert(vbdvar != NULL);
3177  assert(!SCIPisZero(scip, vbdcoef));
3178  assert(!SCIPisInfinity(scip, ABS(vbdcoef)));
3179  assert(infeas != NULL);
3180 
3181  *infeas = FALSE;
3182 
3183  /* make sure vbdvar is active, so we can search for it in SCIPvarGetVxbdVars() */
3184  if( !SCIPvarIsActive(vbdvar) )
3185  {
3186  SCIP_Real constant;
3187 
3188  constant = 0.0;
3189  SCIP_CALL( SCIPgetProbvarSum(scip, &vbdvar, &vbdcoef, &constant) );
3190  if( !SCIPvarIsActive(vbdvar) || (vbdcoef == 0.0) )
3191  return SCIP_OKAY;
3192 
3193  if( !SCIPisInfinity(scip, -lhs) )
3194  lhs -= constant;
3195  if( !SCIPisInfinity(scip, rhs) )
3196  rhs -= constant;
3197  }
3198 
3199  /* vbdvar should be a non-fixed binary variable */
3200  assert(SCIPvarIsIntegral(vbdvar));
3201  assert(SCIPisZero(scip, SCIPvarGetLbGlobal(vbdvar)));
3202  assert(SCIPisEQ(scip, SCIPvarGetUbGlobal(vbdvar), 1.0));
3203 
3204  SCIPdebugMsg(scip, "-> %g <= <%s> + %g*<%s> <= %g\n", lhs, SCIPvarGetName(var), vbdcoef, SCIPvarGetName(vbdvar), rhs);
3205 
3206  if( addcons && SCIPvarGetStatus(var) != SCIP_VARSTATUS_MULTAGGR )
3207  {
3208  SCIP_CONS* vbdcons;
3209  char name[SCIP_MAXSTRLEN];
3210 
3211  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_vbnd", SCIPconsGetName(cons));
3212 
3213  SCIP_CALL( SCIPcreateConsVarbound(scip, &vbdcons, name, var, vbdvar, vbdcoef, lhs, rhs,
3215  SCIP_CALL( SCIPaddCons(scip, vbdcons) );
3216  SCIP_CALL( SCIPreleaseCons(scip, &vbdcons) );
3217 
3218  ++*naddconss;
3219 
3220  return SCIP_OKAY;
3221  }
3222 
3223  if( !SCIPisInfinity(scip, -lhs) )
3224  {
3225  SCIP_CALL( SCIPaddVarVlb(scip, var, vbdvar, -vbdcoef, lhs, infeas, &nbdchgs_local) );
3226  if( *infeas )
3227  return SCIP_OKAY;
3228  *nbdchgs += nbdchgs_local;
3229  }
3230 
3231  if( !SCIPisInfinity(scip, rhs) )
3232  {
3233  SCIP_CALL( SCIPaddVarVub(scip, var, vbdvar, -vbdcoef, rhs, infeas, &nbdchgs_local) );
3234  if( *infeas )
3235  return SCIP_OKAY;
3236  *nbdchgs += nbdchgs_local;
3237  }
3238 
3239  return SCIP_OKAY;
3240 }
3241 
3242 /** propagates varbounds of variables
3243  * Let f(x) = sign(x+offset)|x+offset|^n, f^{-1}(y) = sign(y)|y|^(1/n) - offset.
3244  * Thus, constraint is lhs <= f(x) + c*z <= rhs.
3245  *
3246  * Given a variable bound constraint x <= a*y + b with y a binary variable, one obtains
3247  * y = 0 -> f(x) <= f(b) -> lhs <= f(b) + c*z
3248  * y = 1 -> f(x) <= f(a+b) -> lhs <= f(a+b) + c*z
3249  * => lhs <= f(b) + y * (f(a+b)-f(b)) + c*z
3250  *
3251  * Given a variable bound constraint x >= a*y + b with y a binary variable, one obtains analogously
3252  * f(b) + y * (f(a+b)-f(b)) + c*z <= rhs
3253  *
3254  * Given a variable bound constraint c*z <= a*y + b with y a binary variable, one obtains
3255  * y = 0 -> lhs <= f(x) + b -> x >= f^{-1}(lhs - b)
3256  * y = 1 -> lhs <= f(x) + a+b -> x >= f^{-1}(lhs - (a+b))
3257  * => x >= f^{-1}(lhs - b) + y * (f^{-1}(lhs - (a+b)) - f^{-1}(lhs - b))
3258  *
3259  * Given a variable bound constraint c*z >= a*y + b with y a binary variable, one obtains analogously
3260  * x <= f^{-1}(rhs - b) + y * (f^{-1}(rhs - (a+b)) - f^{-1}(rhs - b))
3261  */
3262 static
3264  SCIP* scip, /**< SCIP data structure */
3265  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3266  SCIP_CONS* cons, /**< absolute power constraint */
3267  SCIP_Bool* infeas, /**< pointer to store whether an infeasibility was detected */
3268  int* nbdchgs, /**< pointer where to add number of performed bound changes */
3269  int* naddconss /**< pointer where to add number of added constraints */
3270  )
3271 {
3272  SCIP_CONSHDLRDATA* conshdlrdata;
3273  SCIP_CONSDATA* consdata;
3274  SCIP_VAR* y;
3275  SCIP_Real a;
3276  SCIP_Real b;
3277  SCIP_Real fb;
3278  SCIP_Real fab;
3279  SCIP_Real vbcoef;
3280  SCIP_Real vbconst;
3281  int i;
3282 
3283  assert(scip != NULL);
3284  assert(conshdlr != NULL);
3285  assert(cons != NULL);
3286  assert(infeas != NULL);
3287  assert(nbdchgs != NULL);
3288  assert(naddconss != NULL);
3289 
3290  *infeas = FALSE;
3291 
3292  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3293  assert(conshdlrdata != NULL);
3294 
3295  consdata = SCIPconsGetData(cons);
3296  assert(consdata != NULL);
3297  assert(consdata->z != NULL);
3298 
3299  /* don't do anything if it looks like we have numerical troubles */
3300  if( SCIPisZero(scip, consdata->zcoef) )
3301  return SCIP_OKAY;
3302 
3303  if( !SCIPisInfinity(scip, -consdata->lhs) )
3304  {
3305  /* propagate varbounds x <= a*y+b onto z
3306  * lhs <= f(b) + y * (f(a+b)-f(b)) + c*z
3307  * -> c*z >= lhs-f(b) + y * (f(b)-f(a+b))
3308  */
3309  for( i = 0; i < SCIPvarGetNVubs(consdata->x); ++i )
3310  {
3311  y = SCIPvarGetVubVars(consdata->x)[i];
3312  a = SCIPvarGetVubCoefs(consdata->x)[i];
3313  b = SCIPvarGetVubConstants(consdata->x)[i];
3314 
3315  /* skip variable bound if y is not integer or its valid values are not {0,1}
3316  * @todo extend to arbitrary integer variables
3317  */
3318  if( !SCIPvarIsBinary(y) || SCIPvarGetLbGlobal(y) > 0.5 || SCIPvarGetUbGlobal(y) < 0.5 )
3319  continue;
3320 
3321  /* skip variable bound if coefficient is very small */
3322  if( SCIPisFeasZero(scip, consdata->power(a, consdata->exponent)) )
3323  continue;
3324 
3325  SCIPdebugMsg(scip, "propagate variable bound <%s> <= %g*<%s> + %g\n", SCIPvarGetName(consdata->x), a, SCIPvarGetName(y), b);
3326 
3327  fb = SIGN( b + consdata->xoffset) * consdata->power( b + consdata->xoffset, consdata->exponent); /* f( b) = sign( b) | b|^n */
3328  fab = SIGN(a+b + consdata->xoffset) * consdata->power(a+b + consdata->xoffset, consdata->exponent); /* f(a+b) = sign(a+b) |a+b|^n */
3329 
3330  vbcoef = (fb - fab) / consdata->zcoef;
3331  vbconst = (consdata->lhs - fb) / consdata->zcoef;
3332 
3333  if( consdata->zcoef > 0.0 )
3334  {
3335  /* add varbound z >= (lhs-f(b))/c + y * (f(b)-f(a+b))/c */
3336  SCIP_CALL( addVarbound(scip, cons, conshdlrdata->addvarboundcons, consdata->z, y, -vbcoef, vbconst, SCIPinfinity(scip), infeas, nbdchgs, naddconss) );
3337  }
3338  else
3339  {
3340  /* add varbound z <= (lhs-f(b))/c + y * (f(b)-f(a+b))/c */
3341  SCIP_CALL( addVarbound(scip, cons, conshdlrdata->addvarboundcons, consdata->z, y, -vbcoef, -SCIPinfinity(scip), vbconst, infeas, nbdchgs, naddconss) );
3342  }
3343  if( *infeas )
3344  return SCIP_OKAY;
3345  }
3346  }
3347 
3348  /* propagate varbounds x >= a*y+b onto z
3349  * f(b) + y * (f(a+b)-f(b)) + c*z <= rhs
3350  * -> c*z <= rhs-f(b) + y * (f(b)-f(a+b))
3351  */
3352  if( !SCIPisInfinity(scip, consdata->rhs) )
3353  {
3354  for( i = 0; i < SCIPvarGetNVlbs(consdata->x); ++i )
3355  {
3356  y = SCIPvarGetVlbVars(consdata->x)[i];
3357  a = SCIPvarGetVlbCoefs(consdata->x)[i];
3358  b = SCIPvarGetVlbConstants(consdata->x)[i];
3359 
3360  /* skip variable bound if y is not integer or its valid values are not {0,1}
3361  * @todo extend to arbitrary integer variables
3362  */
3363  if( !SCIPvarIsBinary(y) || SCIPvarGetLbGlobal(y) > 0.5 || SCIPvarGetUbGlobal(y) < 0.5 )
3364  continue;
3365 
3366  /* skip variable bound if coefficient is very small */
3367  if( SCIPisFeasZero(scip, consdata->power(a, consdata->exponent)) )
3368  continue;
3369 
3370  SCIPdebugMsg(scip, "propagate variable bound <%s> >= %g*<%s> + %g\n", SCIPvarGetName(consdata->x), a, SCIPvarGetName(y), b);
3371 
3372  fb = SIGN( b + consdata->xoffset) * consdata->power( b + consdata->xoffset, consdata->exponent); /* f( b) = sign( b) | b|^n */
3373  fab = SIGN(a+b + consdata->xoffset) * consdata->power(a+b + consdata->xoffset, consdata->exponent); /* f(a+b) = sign(a+b) |a+b|^n */
3374 
3375  vbcoef = (fb - fab) / consdata->zcoef;
3376  vbconst = (consdata->rhs - fb) / consdata->zcoef;
3377 
3378  if( consdata->zcoef > 0.0 )
3379  {
3380  /* add varbound z <= (rhs-f(b))/c + y * (f(b)-f(a+b))/c */
3381  SCIP_CALL( addVarbound(scip, cons, conshdlrdata->addvarboundcons, consdata->z, y, -vbcoef, -SCIPinfinity(scip), vbconst, infeas, nbdchgs, naddconss) );
3382  }
3383  else
3384  {
3385  /* add varbound z >= (rhs-f(b))/c + y * (f(b)-f(a+b))/c */
3386  SCIP_CALL( addVarbound(scip, cons, conshdlrdata->addvarboundcons, consdata->z, y, -vbcoef, vbconst, SCIPinfinity(scip), infeas, nbdchgs, naddconss) );
3387  }
3388  if( *infeas )
3389  return SCIP_OKAY;
3390  }
3391  }
3392 
3393  /* propagate variable upper bounds on z onto x
3394  * c*z <= a*y+b -> x >= f^{-1}(lhs - b) + y * (f^{-1}(lhs - (a+b)) - f^{-1}(lhs - b))
3395  * c*z >= a*y+b -> x <= f^{-1}(rhs - b) + y * (f^{-1}(rhs - (a+b)) - f^{-1}(rhs - b))
3396  */
3397  if( (consdata->zcoef > 0.0 && !SCIPisInfinity(scip, -consdata->lhs)) ||
3398  ( consdata->zcoef < 0.0 && !SCIPisInfinity(scip, consdata->rhs)) )
3399  for( i = 0; i < SCIPvarGetNVubs(consdata->z); ++i )
3400  {
3401  y = SCIPvarGetVubVars(consdata->z)[i];
3402  a = SCIPvarGetVubCoefs(consdata->z)[i] * consdata->zcoef;
3403  b = SCIPvarGetVubConstants(consdata->z)[i] * consdata->zcoef;
3404 
3405  SCIPdebugMsg(scip, "propagate variable bound %g*<%s> %c= %g*<%s> + %g\n", consdata->zcoef, SCIPvarGetName(consdata->z), consdata->zcoef > 0 ? '<' : '>', a, SCIPvarGetName(y), b);
3406 
3407  /* skip variable bound if y is not integer or its valid values are not {0,1}
3408  * @todo extend to arbitrary integer variables
3409  */
3410  if( !SCIPvarIsBinary(y) || SCIPvarGetLbGlobal(y) > 0.5 || SCIPvarGetUbGlobal(y) < 0.5 )
3411  continue;
3412 
3413  if( consdata->zcoef > 0.0 )
3414  {
3415  fb = consdata->lhs - b;
3416  fb = SIGN(fb) * pow(ABS(fb), 1.0/consdata->exponent);
3417  fab = consdata->lhs - (a+b);
3418  fab = SIGN(fab) * pow(ABS(fab), 1.0/consdata->exponent);
3419  SCIP_CALL( addVarbound(scip, cons, conshdlrdata->addvarboundcons, consdata->x, y, fb - fab, fb - consdata->xoffset, SCIPinfinity(scip), infeas, nbdchgs, naddconss) );
3420  }
3421  else
3422  {
3423  fb = consdata->rhs - b;
3424  fb = SIGN(fb) * pow(ABS(fb), 1.0/consdata->exponent);
3425  fab = consdata->rhs - (a+b);
3426  fab = SIGN(fab) * pow(ABS(fab), 1.0/consdata->exponent);
3427  SCIP_CALL( addVarbound(scip, cons, conshdlrdata->addvarboundcons, consdata->x, y, fb - fab, -SCIPinfinity(scip), fb - consdata->xoffset, infeas, nbdchgs, naddconss) );
3428  }
3429  if( *infeas )
3430  return SCIP_OKAY;
3431  }
3432 
3433  /* propagate variable lower bounds on z onto x
3434  * c*z <= a*y+b -> x >= f^{-1}(lhs - b) + y * (f^{-1}(lhs - (a+b)) - f^{-1}(lhs - b))
3435  * c*z >= a*y+b -> x <= f^{-1}(rhs - b) + y * (f^{-1}(rhs - (a+b)) - f^{-1}(rhs - b))
3436  */
3437  if( (consdata->zcoef < 0.0 && !SCIPisInfinity(scip, -consdata->lhs)) ||
3438  ( consdata->zcoef > 0.0 && !SCIPisInfinity(scip, consdata->rhs)) )
3439  for( i = 0; i < SCIPvarGetNVlbs(consdata->z); ++i )
3440  {
3441  y = SCIPvarGetVlbVars(consdata->z)[i];
3442  a = SCIPvarGetVlbCoefs(consdata->z)[i] * consdata->zcoef;
3443  b = SCIPvarGetVlbConstants(consdata->z)[i] * consdata->zcoef;
3444 
3445  SCIPdebugMsg(scip, "propagate variable bound %g*<%s> %c= %g*<%s> + %g\n", consdata->zcoef, SCIPvarGetName(consdata->z), consdata->zcoef > 0 ? '>' : '<', a, SCIPvarGetName(y), b);
3446 
3447  /* skip variable bound if y is not integer or its valid values are not {0,1}
3448  * @todo extend to arbitrary integer variables
3449  */
3450  if( !SCIPvarIsBinary(y) || SCIPvarGetLbGlobal(y) > 0.5 || SCIPvarGetUbGlobal(y) < 0.5 )
3451  continue;
3452 
3453  if( consdata->zcoef > 0.0 )
3454  {
3455  fb = consdata->rhs - b;
3456  fb = SIGN(fb) * pow(ABS(fb), 1.0/consdata->exponent);
3457  fab = consdata->rhs - (a+b);
3458  fab = SIGN(fab) * pow(ABS(fab), 1.0/consdata->exponent);
3459  SCIP_CALL( addVarbound(scip, cons, conshdlrdata->addvarboundcons, consdata->x, y, fb - fab, -SCIPinfinity(scip), fb - consdata->xoffset, infeas, nbdchgs, naddconss) );
3460  }
3461  else
3462  {
3463  fb = consdata->lhs - b;
3464  fb = SIGN(fb) * pow(ABS(fb), 1.0/consdata->exponent);
3465  fab = consdata->lhs - (a+b);
3466  fab = SIGN(fab) * pow(ABS(fab), 1.0/consdata->exponent);
3467  SCIP_CALL( addVarbound(scip, cons, conshdlrdata->addvarboundcons, consdata->x, y, fb - fab, fb - consdata->xoffset, SCIPinfinity(scip), infeas, nbdchgs, naddconss) );
3468  }
3469  if( *infeas )
3470  return SCIP_OKAY;
3471  }
3472 
3473  return SCIP_OKAY;
3474 }
3475 
3476 /** computes linear underestimator for (x+offset)^n + c*z <= rhs by linearization in x
3477  *
3478  * the generated cut is xmul * n * (refpoint+offset)^(n-1) * x + c*z <= rhs + ((n-1)*refpoint-offset) * (refpoint+offset)^(n-1)
3479  */
3480 static
3482  SCIP* scip, /**< SCIP data structure */
3483  SCIP_ROWPREP** rowprep, /**< buffer to store rowprep */
3484  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3485  SCIP_Real refpoint, /**< base point for linearization */
3486  SCIP_Real exponent, /**< exponent n in sign(x)abs(x)^n */
3487  SCIP_Real xoffset, /**< offset of x */
3488  SCIP_Real xmult, /**< multiplier for coefficient of x */
3489  SCIP_Real zcoef, /**< coefficient of z */
3490  SCIP_Real rhs, /**< right hand side */
3491  SCIP_VAR* x, /**< variable x */
3492  SCIP_VAR* z, /**< variable z */
3493  SCIP_Bool islocal /**< whether the cut is valid only locally */
3494  )
3495 {
3496  SCIP_CONSHDLRDATA* conshdlrdata;
3497  SCIP_Real tmp;
3498 
3499  assert(scip != NULL);
3500  assert(rowprep != NULL);
3501  assert(!SCIPisFeasNegative(scip, refpoint+xoffset));
3502  assert(!SCIPisInfinity(scip, refpoint));
3503 
3504  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3505  assert(conshdlrdata != NULL);
3506 
3507  if( refpoint < -xoffset )
3508  refpoint = -xoffset;
3509 
3510  tmp = exponent == 2.0 ? refpoint+xoffset : pow(refpoint+xoffset, exponent-1);
3511  if( SCIPisInfinity(scip, tmp) )
3512  {
3513  SCIPdebugMsg(scip, "skip linearization cut because (refpoint+offset)^(exponent-1) > infinity\n");
3514  *rowprep = NULL;
3515  return SCIP_OKAY;
3516  }
3517 
3518  rhs += ((exponent-1)*refpoint-xoffset)*tmp; /* now rhs is the rhs of the cut */
3519  /* do not change the right hand side to a value > infinity (this would trigger an assertion in lp.c) */
3520  if( SCIPisInfinity(scip, rhs) )
3521  {
3522  SCIPdebugMsg(scip, "skip linearization cut because its rhs would be > infinity\n");
3523  *rowprep = NULL;
3524  return SCIP_OKAY;
3525  }
3526 
3527  SCIP_CALL( SCIPcreateRowprep(scip, rowprep, SCIP_SIDETYPE_RIGHT, islocal) );
3528  (void) SCIPsnprintf((*rowprep)->name, (int)sizeof((*rowprep)->name), "signpowlinearizecut_%u", ++(conshdlrdata->ncuts));
3529  SCIPaddRowprepSide(*rowprep, rhs);
3530  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, x, xmult*exponent*tmp) );
3531  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, z, zcoef) );
3532 
3533  return SCIP_OKAY;
3534 }
3535 
3536 /** computes linear underestimator for (x+xoffset)^n + c*z <= rhs by linearization in x
3537  *
3538  * the generated cut is xmul * n * (refpoint+offset)^(n-1) * x + c*z <= rhs + ((n-1)*refpoint-offset) * (refpoint+offset)^(n-1)
3539  * where refpoint is computed by projecting (xref, zref) onto the graph of (x+offset)^n w.r.t. euclidean norm
3540  *
3541  * Thus, the projection is computed by minimizing 1/2(x-xref)^2 + 1/2(((x+offset)^n-rhs)/(-c) - zref)^2.
3542  * I.e., we aim to find a root of
3543  * g(x) = x - xref + n/c (x+offset)^(n-1) (zref - rhs/c) + n/c^2 (x+offset)^(2n-1)
3544  * We do this numerically by executing up to five newton iterations. It is
3545  * g'(x) = 1 + n(n-1)/c (x+offset)^(n-2) (zref - rhs/c) + n(2n-1)/c^2 (x+offset)^(2n-2)
3546  */
3547 static
3549  SCIP* scip, /**< SCIP data structure */
3550  SCIP_ROWPREP** rowprep, /**< buffer to store rowprep */
3551  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3552  SCIP_Real xref, /**< reference point for x */
3553  SCIP_Real zref, /**< reference point for z */
3554  SCIP_Real xmin, /**< minimal value x is allowed to take */
3555  SCIP_Real exponent, /**< exponent n in sign(x+offset)abs(x+offset)^n */
3556  SCIP_Real xoffset, /**< offset of x */
3557  SCIP_Real xmult, /**< multiplier for coefficient of x */
3558  SCIP_Real zcoef, /**< coefficient of z */
3559  SCIP_Real rhs, /**< right hand side */
3560  SCIP_VAR* x, /**< variable x */
3561  SCIP_VAR* z, /**< variable z */
3562  SCIP_Bool islocal /**< whether the cut is valid only locally */
3563  )
3564 {
3565  SCIP_Real tmp;
3566  SCIP_Real xproj;
3567  SCIP_Real gval;
3568  SCIP_Real gderiv;
3569  int iter;
3570 
3571  assert(scip != NULL);
3572  assert(!SCIPisFeasNegative(scip, xref+xoffset));
3573  assert(!SCIPisInfinity(scip, xref));
3574 
3575  if( xref < xmin )
3576  xref = xmin;
3577 
3578  xproj = xref;
3579  iter = 0;
3580  if( exponent == 2.0 )
3581  {
3582  do
3583  {
3584  tmp = (xproj+xoffset) * (xproj+xoffset);
3585  gval = xproj - xref + 2*(xproj+xoffset) / zcoef * ((tmp-rhs)/zcoef + zref);
3586  if( !SCIPisFeasPositive(scip, ABS(gval)) )
3587  break;
3588 
3589  gderiv = 1 + 6 * tmp / (zcoef*zcoef) + 2 / zcoef * (zref - rhs/zcoef);
3590  xproj -= gval / gderiv;
3591  }
3592  while( ++iter <= 5 );
3593  }
3594  else
3595  {
3596  do
3597  {
3598  tmp = pow(xproj + xoffset, exponent-1);
3599  gval = xproj - xref + exponent / zcoef * (pow(xproj+xoffset, 2*exponent-1)/zcoef + tmp * (zref-rhs/zcoef));
3600  if( !SCIPisFeasPositive(scip, ABS(gval)) )
3601  break;
3602 
3603  gderiv = 1 + exponent / zcoef * ( (2*exponent-1)*tmp*tmp/zcoef + (exponent-1)*pow(xproj+xoffset, exponent-2) * (zref-rhs/zcoef) );
3604  xproj -= gval / gderiv;
3605  }
3606  while( ++iter <= 5 );
3607  }
3608 
3609  if( xproj < xmin )
3610  xproj = xmin;
3611 
3612  SCIP_CALL( generateLinearizationCut(scip, rowprep, conshdlr, xproj, exponent, xoffset, xmult, zcoef, rhs, x, z, islocal) );
3613 
3614  return SCIP_OKAY;
3615 }
3616 
3617 /** computes secant underestimator for sign(x+offset)abs(x+offset)^n + c*z <= rhs
3618  *
3619  * the generated cut is slope*xmult*x + c*z <= rhs + (-xlb-offset)^n + slope*xlb,
3620  * where slope = (sign(xub+offset)*abs(xub+offset)^n + (-xlb-offset)^n) / (xub - xlb).
3621  *
3622  * the cut is not generated if the given solution (or the LP solution) would not be cutoff
3623  */
3624 static
3626  SCIP* scip, /**< SCIP data structure */
3627  SCIP_ROWPREP** rowprep, /**< buffer to store rowprep */
3628  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3629  SCIP_SOL* sol, /**< point we want to cut off, or NULL for LP solution */
3630  SCIP_Real xlb, /**< lower bound of x */
3631  SCIP_Real xub, /**< upper bound of x */
3632  SCIP_Real exponent, /**< exponent n in sign(x+offset)abs(x+offset)^n */
3633  SCIP_Real xoffset, /**< offset of x */
3634  DECL_MYPOW ((*mypow)), /**< function to use for computing power */
3635  SCIP_Real xmult, /**< multiplier for coefficient of x */
3636  SCIP_Real zcoef, /**< coefficient of z */
3637  SCIP_Real rhs, /**< right hand side */
3638  SCIP_VAR* x, /**< variable x */
3639  SCIP_VAR* z /**< variable z */
3640  )
3641 {
3642  SCIP_CONSHDLRDATA* conshdlrdata;
3643  SCIP_Real slope, tmp, val;
3644 
3645  assert(scip != NULL);
3646  assert(SCIPisLE(scip, xlb, xub));
3647  assert(!SCIPisPositive(scip, xlb+xoffset));
3648 
3649  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3650  assert(conshdlrdata != NULL);
3651 
3652  /* ignore constraints with fixed x (should be removed soon) */
3653  if( SCIPisRelEQ(scip, xlb, xub) )
3654  {
3655  SCIPdebugMsg(scip, "skip secant cut because <%s> is fixed [%.15g,%.15g]\n", SCIPvarGetName(x), SCIPvarGetLbLocal(x), SCIPvarGetUbLocal(x));
3656  return SCIP_OKAY;
3657  }
3658 
3659  if( xlb > -xoffset )
3660  xlb = -xoffset;
3661 
3662  tmp = mypow(-xlb-xoffset, exponent);
3663  slope = SIGN(xub+xoffset) * mypow(ABS(xub+xoffset), exponent) + tmp;
3664  slope /= xub - xlb;
3665 
3666  /* check if cut would violated solution, check that slope is not above value of infinity */
3667  val = -tmp + slope * (xmult * SCIPgetSolVal(scip, sol, x) - xlb) + zcoef * SCIPgetSolVal(scip, sol, z) - rhs;
3668  if( !SCIPisFeasPositive(scip, val) || SCIPisInfinity(scip, REALABS(slope)) )
3669  {
3670  *rowprep = NULL;
3671  return SCIP_OKAY;
3672  }
3673 
3674  SCIP_CALL( SCIPcreateRowprep(scip, rowprep, SCIP_SIDETYPE_RIGHT, SCIPnodeGetDepth(SCIPgetCurrentNode(scip)) > 0 /* local */) );
3675  (void) SCIPsnprintf((*rowprep)->name, SCIP_MAXSTRLEN, "signpowsecantcut_%u", ++(conshdlrdata->nsecantcuts));
3676 
3677  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, x, xmult*slope) );
3678  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, z, zcoef) );
3679  SCIPaddRowprepSide(*rowprep, rhs + tmp + slope*xlb);
3680 
3681  return SCIP_OKAY;
3682 }
3683 
3684 /** computes secant underestimator for sign(x+xoffset)abs(x+xoffset)^n + c*z <= rhs
3685  *
3686  * The generated cut is slope*xmult*x + c*z <= rhs + (-xlb-xoffset)^n + slope*xlb,
3687  * where slope = (sign(xub+xoffset)*abs(xub+xoffset)^n + (-xlb-xoffset)^n) / (xub - xlb).
3688  */
3689 static
3691  SCIP* scip, /**< SCIP data structure */
3692  SCIP_ROWPREP** rowprep, /**< buffer to store rowprep */
3693  SCIP_Real xlb, /**< lower bound of x */
3694  SCIP_Real xub, /**< upper bound of x */
3695  SCIP_Real exponent, /**< exponent n in sign(x)abs(x)^n */
3696  SCIP_Real xoffset, /**< offset of x */
3697  DECL_MYPOW ((*mypow)), /**< function to use for computing power */
3698  SCIP_Real xmult, /**< multiplier for coefficient of x */
3699  SCIP_Real zcoef, /**< coefficient of z */
3700  SCIP_Real rhs, /**< right hand side */
3701  SCIP_VAR* x, /**< variable x */
3702  SCIP_VAR* z /**< variable z */
3703  )
3704 {
3705  SCIP_Real slope, tmp;
3706 
3707  assert(scip != NULL);
3708  assert(rowprep != NULL);
3709  assert(SCIPisLE(scip, xlb, xub));
3710  assert(!SCIPisPositive(scip, xlb + xoffset));
3711 
3712  /* ignore constraints with fixed x (should be removed soon) */
3713  if( SCIPisRelEQ(scip, xlb, xub) )
3714  return SCIP_OKAY;
3715 
3716  if( xlb > -xoffset )
3717  xlb = -xoffset;
3718 
3719  tmp = mypow(-xlb-xoffset, exponent);
3720  slope = SIGN(xub+xoffset) * mypow(ABS(xub+xoffset), exponent) + tmp;
3721  slope /= xub - xlb;
3722 
3723  if( SCIPisInfinity(scip, REALABS(slope)) )
3724  return SCIP_OKAY;
3725 
3726  SCIP_CALL( SCIPcreateRowprep(scip, rowprep, SCIP_SIDETYPE_RIGHT, SCIPnodeGetDepth(SCIPgetCurrentNode(scip)) > 0 /* local */) );
3727  (void)SCIPmemccpy((*rowprep)->name, "signpowcut", '\0', 11);
3728  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, x, xmult*slope) );
3729  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, z, zcoef) );
3730  SCIPaddRowprepSide(*rowprep, rhs + tmp + slope*xlb);
3731 
3732  return SCIP_OKAY;
3733 }
3734 
3735 /** generates a cut
3736  * based on Liberti and Pantelides, Convex Envelopes of Monomials of Odd Degree, J. Global Optimization 25, 157-168, 2003, and previous publications
3737  */
3738 static
3740  SCIP* scip, /**< SCIP data structure */
3741  SCIP_CONS* cons, /**< constraint */
3742  SCIP_SIDETYPE violside, /**< side to separate */
3743  SCIP_SOL* sol, /**< solution to separate, or NULL if LP solution should be used */
3744  SCIP_ROW** row, /**< storage for cut */
3745  SCIP_Bool onlyinbounds, /**< whether linearization is allowed only in variable bounds */
3746  SCIP_Real minviol /**< a minimal violation in sol we hope to achieve */
3747  )
3748 {
3749  SCIP_CONSHDLRDATA* conshdlrdata;
3750  SCIP_CONSDATA* consdata;
3751  SCIP_ROWPREP* rowprep = NULL;
3752  SCIP_Real c;
3753  SCIP_Real xlb;
3754  SCIP_Real xglb;
3755  SCIP_Real xub;
3756  SCIP_Real xval;
3757  SCIP_Real xoffset;
3758  SCIP_Real xmult;
3759  SCIP_Real zcoef;
3760  SCIP_Real rhs;
3761 
3762  assert(scip != NULL);
3763  assert(cons != NULL);
3764  assert(row != NULL);
3765 
3766  conshdlrdata = SCIPconshdlrGetData(SCIPconsGetHdlr(cons));
3767  assert(conshdlrdata != NULL);
3768 
3769  consdata = SCIPconsGetData(cons);
3770  assert(consdata != NULL);
3771 
3772  assert(SCIPisGT(scip, violside == SCIP_SIDETYPE_LEFT ? consdata->lhsviol : consdata->rhsviol, SCIPfeastol(scip)));
3773 
3774  *row = NULL;
3775 
3776  SCIPdebugMsg(scip, "generate cut for constraint <%s> with violated side %d\n", SCIPconsGetName(cons), violside);
3777  SCIPdebugPrintCons(scip, cons, NULL);
3778  SCIPdebugMsg(scip, "xlb = %g xub = %g xval = %g zval = %.15g\n", SCIPvarGetLbLocal(consdata->x), SCIPvarGetUbLocal(consdata->x), SCIPgetSolVal(scip, sol, consdata->x), SCIPgetSolVal(scip, sol, consdata->z));
3779 
3780  if( violside == SCIP_SIDETYPE_RIGHT )
3781  {
3782  xglb = SCIPvarGetLbGlobal(consdata->x);
3783  xlb = SCIPvarGetLbLocal(consdata->x);
3784  xub = SCIPvarGetUbLocal(consdata->x);
3785  xval = SCIPgetSolVal(scip, sol, consdata->x);
3786  xoffset = consdata->xoffset;
3787  xmult = 1.0;
3788  zcoef = consdata->zcoef;
3789  rhs = consdata->rhs;
3790  }
3791  else
3792  {
3793  xglb = -SCIPvarGetUbGlobal(consdata->x);
3794  xlb = -SCIPvarGetUbLocal(consdata->x);
3795  xub = -SCIPvarGetLbLocal(consdata->x);
3796  xval = -SCIPgetSolVal(scip, sol, consdata->x);
3797  xoffset = -consdata->xoffset;
3798  xmult = -1.0;
3799  zcoef = -consdata->zcoef;
3800  rhs = -consdata->lhs;
3801  }
3802  /* move reference point onto local domain, if clearly (>eps) outside */
3803  if( SCIPisLT(scip, xval, xlb) )
3804  xval = xlb;
3805  else if( SCIPisGT(scip, xval, xub) )
3806  xval = xub;
3807 
3808  if( SCIPisInfinity(scip, REALABS(xval)) )
3809  {
3810  SCIPdebugMsg(scip, "skip separation since x is at infinity\n");
3811  return SCIP_OKAY;
3812  }
3813 
3814  if( !SCIPisNegative(scip, xlb+xoffset) )
3815  {
3816  /* [xlb, xub] completely in positive orthant -> function is convex on whole domain */
3817  SCIP_Bool islocal;
3818 
3819  islocal = (!SCIPconsIsGlobal(cons) || SCIPisNegative(scip, xglb+xoffset)) && SCIPnodeGetDepth(SCIPgetCurrentNode(scip)) > 0;
3820  if( conshdlrdata->projectrefpoint && !onlyinbounds )
3821  {
3822  SCIP_CALL( generateLinearizationCutProject(scip, &rowprep, SCIPconsGetHdlr(cons), xval, SCIPgetSolVal(scip, sol, consdata->z), -xoffset, consdata->exponent,
3823  xoffset, xmult, zcoef, rhs, consdata->x, consdata->z, islocal) );
3824  }
3825  else if( !onlyinbounds )
3826  {
3827  SCIP_CALL( generateLinearizationCut(scip, &rowprep, SCIPconsGetHdlr(cons), xval, consdata->exponent, xoffset, xmult, zcoef, rhs,
3828  consdata->x, consdata->z, islocal) );
3829  }
3830  else
3831  {
3832  SCIP_CALL( generateLinearizationCut(scip, &rowprep, SCIPconsGetHdlr(cons), 2.0*xval > xlb + xub ? xub : xlb, consdata->exponent, xoffset, xmult, zcoef, rhs,
3833  consdata->x, consdata->z, islocal) );
3834  }
3835  }
3836  else if( !SCIPisPositive(scip, xub+xoffset) )
3837  {
3838  /* [xlb, xub] completely in negative orthant -> function is concave on whole domain */
3839  if( SCIPisInfinity(scip, -xlb) )
3840  return SCIP_OKAY;
3841  SCIP_CALL( generateSecantCut(scip, &rowprep, SCIPconsGetHdlr(cons), sol, xlb, xub, consdata->exponent, xoffset, consdata->power, xmult, zcoef, rhs, consdata->x, consdata->z) );
3842  }
3843  else if( (c = - consdata->root * (xlb+xoffset) - xoffset) > xub )
3844  {
3845  /* c is right of xub -> use secant */
3846  if( SCIPisInfinity(scip, -xlb) || SCIPisInfinity(scip, xub) )
3847  return SCIP_OKAY;
3848  SCIP_CALL( generateSecantCut(scip, &rowprep, SCIPconsGetHdlr(cons), sol, xlb, xub, consdata->exponent, xoffset, consdata->power, xmult, zcoef, rhs, consdata->x, consdata->z) );
3849  }
3850  else if( xval >= c )
3851  {
3852  /* xval is right of c -> use linearization */
3853  if( conshdlrdata->projectrefpoint && !onlyinbounds )
3854  SCIP_CALL( generateLinearizationCutProject(scip, &rowprep, SCIPconsGetHdlr(cons), xval, SCIPgetSolVal(scip, sol, consdata->z), c, consdata->exponent,
3855  xoffset, xmult, zcoef, rhs, consdata->x, consdata->z, SCIPnodeGetDepth(SCIPgetCurrentNode(scip)) > 0) );
3856  else if( !onlyinbounds )
3857  SCIP_CALL( generateLinearizationCut(scip, &rowprep, SCIPconsGetHdlr(cons), xval, consdata->exponent, xoffset, xmult, zcoef, rhs,
3858  consdata->x, consdata->z, xval+xoffset < - consdata->root * (xglb+xoffset) && SCIPnodeGetDepth(SCIPgetCurrentNode(scip)) > 0) );
3859  else
3860  SCIP_CALL( generateLinearizationCut(scip, &rowprep, SCIPconsGetHdlr(cons), xub, consdata->exponent, xoffset, xmult, zcoef, rhs,
3861  consdata->x, consdata->z, xub+xoffset < - consdata->root * (xglb+xoffset) && SCIPnodeGetDepth(SCIPgetCurrentNode(scip)) > 0) );
3862  }
3863  else
3864  {
3865  /* xval between xlb and c -> use secant */
3866  if( SCIPisInfinity(scip, -xlb) || SCIPisInfinity(scip, c) )
3867  return SCIP_OKAY;
3868  SCIP_CALL( generateSecantCut(scip, &rowprep, SCIPconsGetHdlr(cons), sol, xlb, c, consdata->exponent, xoffset, consdata->power, xmult, zcoef, rhs, consdata->x, consdata->z) );
3869  }
3870 
3871  /* check and improve numerics */
3872  if( rowprep != NULL )
3873  {
3874  SCIP_Real coefrange;
3875 
3876  SCIPdebug( SCIPprintRowprep(scip, rowprep, NULL) );
3877 
3878  /* we should not need SCIPmergeRowprep() with only 2 vars in the row */
3879  assert(rowprep->nvars <= 2);
3880 
3881  SCIP_CALL( SCIPcleanupRowprep(scip, rowprep, sol, conshdlrdata->cutmaxrange, minviol, &coefrange, NULL) );
3882 
3883  if( coefrange >= conshdlrdata->cutmaxrange )
3884  {
3885  SCIPdebugMsg(scip, "skip cut for constraint <%s> because of very large range: %g\n", SCIPconsGetName(cons), coefrange);
3886  }
3887  else if( SCIPisInfinity(scip, REALABS(rowprep->side)) )
3888  {
3889  SCIPdebugMsg(scip, "skip cut for constraint <%s> because of very large side: %g\n", SCIPconsGetName(cons), rowprep->side);
3890  }
3891  else if( rowprep->nvars > 0 && SCIPisInfinity(scip, REALABS(rowprep->coefs[0])) )
3892  {
3893  SCIPdebugMsg(scip, "skip cut for constraint <%s> because of very large coef: %g\n", SCIPconsGetName(cons), rowprep->coefs[0]);
3894  }
3895  else
3896  {
3897  SCIP_CALL( SCIPgetRowprepRowConshdlr(scip, row, rowprep, SCIPconsGetHdlr(cons)) );
3898  }
3899 
3900  SCIPfreeRowprep(scip, &rowprep);
3901  }
3902 
3903  return SCIP_OKAY;
3904 }
3905 
3906 /** tries to separate solution or LP solution by a linear cut
3907  * assumes that constraint violations have been computed
3908  */
3909 static
3911  SCIP* scip, /**< SCIP data structure */
3912  SCIP_CONSHDLR* conshdlr, /**< quadratic constraints handler */
3913  SCIP_CONS** conss, /**< constraints */
3914  int nconss, /**< number of constraints */
3915  int nusefulconss, /**< number of constraints that seem to be useful */
3916  SCIP_SOL* sol, /**< solution to separate, or NULL if LP solution should be used */
3917  SCIP_Real minefficacy, /**< minimal efficacy of a cut if it should be added to the LP */
3918  SCIP_Bool inenforcement, /**< whether we are in constraint enforcement */
3919  SCIP_Bool onlyinbounds, /**< whether linearization is allowed only in variable bounds */
3920  SCIP_Bool* success, /**< result of separation: separated point (TRUE) or not (FALSE) */
3921  SCIP_Bool* cutoff, /**< whether a cutoff has been detected */
3922  SCIP_Real* bestefficacy /**< buffer to store best efficacy of a cut that was added to the LP, if found; or NULL if not of interest */
3923  )
3924 {
3925  SCIP_CONSHDLRDATA* conshdlrdata;
3926  SCIP_CONSDATA* consdata;
3927  SCIP_SIDETYPE side;
3928  SCIP_Real efficacy;
3929  int c;
3930  SCIP_ROW* row;
3931 
3932  assert(scip != NULL);
3933  assert(conshdlr != NULL);
3934  assert(conss != NULL || nconss == 0);
3935  assert(success != NULL);
3936  assert(cutoff != NULL);
3937 
3938  *success = FALSE;
3939  *cutoff = FALSE;
3940 
3941  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3942  assert(conshdlrdata != NULL);
3943 
3944  if( bestefficacy != NULL )
3945  *bestefficacy = 0.0;
3946 
3947  for( c = 0, side = SCIP_SIDETYPE_LEFT; c < nconss && ! (*cutoff); c = (side == SCIP_SIDETYPE_RIGHT ? c+1 : c), side = (side == SCIP_SIDETYPE_LEFT ? SCIP_SIDETYPE_RIGHT : SCIP_SIDETYPE_LEFT) )
3948  {
3949  assert(conss[c] != NULL); /*lint !e613*/
3950 
3951  /* skip constraints that are not enabled, e.g., because they were already marked for deletion at this node */
3952  if( !SCIPconsIsEnabled(conss[c]) ) /*lint !e613*/
3953  continue;
3954 
3955  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
3956  assert(consdata != NULL);
3957 
3958  if( SCIPisGT(scip, side == SCIP_SIDETYPE_LEFT ? consdata->lhsviol : consdata->rhsviol, SCIPfeastol(scip)) )
3959  {
3960  /* try to generate a cut */
3961  SCIP_CALL( generateCut(scip, conss[c], side, sol, &row, onlyinbounds, minefficacy) ); /*lint !e613*/
3962  if( row == NULL ) /* failed to generate cut */
3963  continue;
3964 
3965  /* if cut is violated sufficiently, then add,
3966  * unless it corresponds to a bound change that is too weak (<eps) to be added
3967  */
3968  efficacy = -SCIPgetRowSolFeasibility(scip, row, sol);
3969  if( SCIPisGT(scip, efficacy, minefficacy) && SCIPisCutApplicable(scip, row) )
3970  {
3971  SCIP_Bool infeasible;
3972 
3973  SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
3974  if ( infeasible )
3975  *cutoff = TRUE;
3976  else
3977  *success = TRUE;
3978  if( bestefficacy != NULL && efficacy > *bestefficacy )
3979  *bestefficacy = efficacy;
3980 
3981  /* notify indicator constraint handler about this cut */
3982  if( conshdlrdata->conshdlrindicator != NULL && !SCIProwIsLocal(row) )
3983  {
3984  SCIP_CALL( SCIPaddRowIndicator(scip, conshdlrdata->conshdlrindicator, row) );
3985  }
3986 
3987  /* mark row as not removable from LP for current node, if in enforcement */
3988  if( inenforcement && !conshdlrdata->enfocutsremovable )
3989  SCIPmarkRowNotRemovableLocal(scip, row);
3990  }
3991 
3992  SCIP_CALL( SCIPreleaseRow (scip, &row) );
3993  }
3994 
3995  /* enforce only useful constraints
3996  * others are only checked and enforced if we are still feasible or have not found a separating cut yet
3997  */
3998  if( c >= nusefulconss && *success )
3999  break;
4000  }
4001 
4002  return SCIP_OKAY;
4003 }
4004 
4005 /** adds linearizations cuts for convex constraints w.r.t. a given reference point to cutpool and sepastore
4006  * if separatedlpsol is not NULL, then a cut that separates the LP solution is added to the sepastore and is forced to enter the LP
4007  * if separatedlpsol is not NULL, but cut does not separate the LP solution, then it is added to the cutpool only
4008  * if separatedlpsol is NULL, then cut is added to cutpool only
4009  */
4010 static
4012  SCIP* scip, /**< SCIP data structure */
4013  SCIP_CONSHDLR* conshdlr, /**< quadratic constraints handler */
4014  SCIP_CONS** conss, /**< constraints */
4015  int nconss, /**< number of constraints */
4016  SCIP_SOL* ref, /**< reference point where to linearize, or NULL for LP solution */
4017  SCIP_Bool* separatedlpsol, /**< buffer to store whether a cut that separates the current LP solution was found and added to LP, or NULL if adding to cutpool only */
4018  SCIP_Real minefficacy /**< minimal efficacy of a cut when checking for separation of LP solution */
4019  )
4020 {
4021  SCIP_CONSDATA* consdata;
4022  SCIP_Bool addedtolp;
4023  SCIP_ROW* row;
4024  int c;
4025 
4026  assert(scip != NULL);
4027  assert(conshdlr != NULL);
4028  assert(conss != NULL || nconss == 0);
4029 
4030  if( separatedlpsol != NULL )
4031  *separatedlpsol = FALSE;
4032 
4033  for( c = 0; c < nconss; ++c )
4034  {
4035  assert(conss[c] != NULL); /*lint !e613*/
4036 
4037  if( SCIPconsIsLocal(conss[c]) ) /*lint !e613*/
4038  continue;
4039 
4040  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
4041  assert(consdata != NULL);
4042 
4043  if( !SCIPisGT(scip, SCIPvarGetUbGlobal(consdata->x), -consdata->xoffset) && !SCIPisInfinity(scip, -consdata->lhs) )
4044  {
4045  /* constraint function is concave for x+offset <= 0.0, so can linearize w.r.t. lhs */
4046  consdata->lhsviol = 1.0;
4047  consdata->rhsviol = 0.0;
4048  SCIP_CALL( generateCut(scip, conss[c], SCIP_SIDETYPE_LEFT, ref, &row, FALSE, minefficacy) ); /*lint !e613*/
4049  }
4050  else if( !SCIPisLT(scip, SCIPvarGetLbGlobal(consdata->x), -consdata->xoffset) && !SCIPisInfinity(scip, -consdata->rhs) )
4051  {
4052  /* constraint function is convex for x+offset >= 0.0, so can linearize w.r.t. rhs */
4053  consdata->lhsviol = 0.0;
4054  consdata->rhsviol = 1.0;
4055  SCIP_CALL( generateCut(scip, conss[c], SCIP_SIDETYPE_RIGHT, ref, &row, FALSE, minefficacy) ); /*lint !e613*/
4056  }
4057  else
4058  {
4059  /* sign not fixed or nothing to linearize */
4060  continue;
4061  }
4062 
4063  if( row == NULL )
4064  continue;
4065 
4066  addedtolp = FALSE;
4067 
4068  /* if caller wants, then check if cut separates LP solution and add to sepastore if so */
4069  if( separatedlpsol != NULL )
4070  {
4071  if( -SCIPgetRowLPFeasibility(scip, row) >= minefficacy )
4072  {
4073  SCIP_Bool infeasible;
4074 
4075  *separatedlpsol = TRUE;
4076  addedtolp = TRUE;
4077  SCIP_CALL( SCIPaddRow(scip, row, TRUE, &infeasible) );
4078  assert( ! infeasible );
4079  }
4080  }
4081 
4082  if( !addedtolp && !SCIProwIsLocal(row) )
4083  {
4084  SCIP_CALL( SCIPaddPoolCut(scip, row) );
4085  }
4086 
4087  SCIP_CALL( SCIPreleaseRow(scip, &row) );
4088  }
4089 
4090  return SCIP_OKAY;
4091 }
4092 
4093 /** processes the event that a new primal solution has been found */
4094 static
4095 SCIP_DECL_EVENTEXEC(processNewSolutionEvent)
4097  SCIP_CONSHDLRDATA* conshdlrdata;
4098  SCIP_CONSHDLR* conshdlr;
4099  SCIP_CONS** conss;
4100  int nconss;
4101  SCIP_SOL* sol;
4102 
4103  assert(scip != NULL);
4104  assert(event != NULL);
4105  assert(eventdata != NULL);
4106  assert(eventhdlr != NULL);
4107 
4108  assert((SCIPeventGetType(event) & SCIP_EVENTTYPE_SOLFOUND) != 0);
4109 
4110  conshdlr = (SCIP_CONSHDLR*)eventdata;
4111 
4112  nconss = SCIPconshdlrGetNConss(conshdlr);
4113 
4114  if( nconss == 0 )
4115  return SCIP_OKAY;
4116 
4117  sol = SCIPeventGetSol(event);
4118  assert(sol != NULL);
4119 
4120  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4121  assert(conshdlrdata != NULL);
4122 
4123  /* we are only interested in solution coming from some heuristic other than trysol, but not from the tree
4124  * the reason for ignoring trysol solutions is that they may come from an NLP solve in sepalp, where we already added linearizations,
4125  * or are from the tree, but postprocessed via proposeFeasibleSolution
4126  */
4127  if( SCIPsolGetHeur(sol) == NULL || SCIPsolGetHeur(sol) == conshdlrdata->trysolheur )
4128  return SCIP_OKAY;
4129 
4130  conss = SCIPconshdlrGetConss(conshdlr);
4131  assert(conss != NULL);
4132 
4133  SCIPdebugMsg(scip, "catched new sol event %" SCIP_EVENTTYPE_FORMAT " from heur <%s>; have %d conss\n", SCIPeventGetType(event), SCIPheurGetName(SCIPsolGetHeur(sol)), nconss);
4134 
4135  SCIP_CALL( addLinearizationCuts(scip, conshdlr, conss, nconss, sol, NULL, 0.0) );
4136 
4137  return SCIP_OKAY;
4138 }
4139 
4140 /** given a solution, try to make absolute power constraints feasible by shifting the linear variable z and pass this solution to the trysol heuristic */
4141 static
4143  SCIP* scip, /**< SCIP data structure */
4144  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
4145  SCIP_CONS** conss, /**< constraints to process */
4146  int nconss, /**< number of constraints */
4147  SCIP_SOL* sol /**< solution to process */
4148  )
4149 {
4150  SCIP_CONSDATA* consdata;
4151  SCIP_SOL* newsol;
4152  SCIP_Real xtermval;
4153  SCIP_Real zval;
4154  SCIP_Real viol;
4155  SCIP_Bool solviolbounds;
4156  int c;
4157 
4158  assert(scip != NULL);
4159  assert(conshdlr != NULL);
4160  assert(conss != NULL || nconss == 0);
4161 
4162  /* don't propose new solutions if not in presolve or solving */
4164  return SCIP_OKAY;
4165 
4166  if( sol != NULL )
4167  {
4168  SCIP_CALL( SCIPcreateSolCopy(scip, &newsol, sol) );
4169  }
4170  else
4171  {
4172  SCIP_CALL( SCIPcreateLPSol(scip, &newsol, NULL) );
4173  }
4174  SCIP_CALL( SCIPunlinkSol(scip, newsol) );
4175 
4176  for( c = 0; c < nconss; ++c )
4177  {
4178  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
4179  assert(consdata != NULL);
4180  assert(consdata->z != NULL);
4181  assert(consdata->zcoef != 0.0);
4182 
4183  /* recompute violation w.r.t. current solution */
4184  SCIP_CALL( computeViolation(scip, conshdlr, conss[c], newsol, &viol, &solviolbounds) ); /*lint !e613*/
4185  assert(!solviolbounds);
4186 
4187  /* do nothing if constraint is satisfied */
4188  if( !SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) && !SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)) )
4189  continue;
4190 
4191  /* if violation is at infinity, give up */
4192  if( SCIPisInfinity(scip, MAX(consdata->lhsviol, consdata->rhsviol)) )
4193  break;
4194 
4195  /* @todo could also adjust x while keeping z fixed */
4196 
4197  /* if variable is multiaggregated, then cannot set its solution value, so give up */
4198  if( SCIPvarGetStatus(consdata->z) == SCIP_VARSTATUS_MULTAGGR )
4199  break;
4200 
4201  /* compute value of x-term */
4202  xtermval = SCIPgetSolVal(scip, newsol, consdata->x);
4203  xtermval += consdata->xoffset;
4204  xtermval = SIGN(xtermval) * consdata->power(ABS(xtermval), consdata->exponent);
4205 
4206  /* if left hand side is violated, try to set z such that lhs is active */
4207  if( SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) )
4208  {
4209  assert(!SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip))); /* should only have one side violated (otherwise some variable is at infinity) */
4210 
4211  zval = (consdata->lhs - xtermval)/consdata->zcoef;
4212  /* bad luck: z would get value outside of its domain */
4213  if( SCIPisInfinity(scip, REALABS(zval)) || SCIPisFeasLT(scip, zval, SCIPvarGetLbGlobal(consdata->z)) || SCIPisFeasGT(scip, zval, SCIPvarGetUbGlobal(consdata->z)) )
4214  break;
4215  SCIP_CALL( SCIPsetSolVal(scip, newsol, consdata->z, zval) );
4216  }
4217 
4218  /* if right hand side is violated, try to set z such that rhs is active */
4219  if( SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)) )
4220  {
4221  zval = (consdata->rhs - xtermval)/consdata->zcoef;
4222  /* bad luck: z would get value outside of its domain */
4223  if( SCIPisInfinity(scip, REALABS(zval)) || SCIPisFeasLT(scip, zval, SCIPvarGetLbGlobal(consdata->z)) || SCIPisFeasGT(scip, zval, SCIPvarGetUbGlobal(consdata->z)) )
4224  break;
4225  SCIP_CALL( SCIPsetSolVal(scip, newsol, consdata->z, zval) );
4226  }
4227  }
4228 
4229  /* if we have a solution that should satisfy all absolute power constraints and has a better objective than the current upper bound, then pass it to the trysol heuristic */
4230  if( c == nconss )
4231  {
4232  SCIP_CONSHDLRDATA* conshdlrdata;
4233 
4234  SCIPdebugMsg(scip, "pass solution with objective val %g to trysol heuristic\n", SCIPgetSolTransObj(scip, newsol));
4235 
4236  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4237  assert(conshdlrdata != NULL);
4238  assert(conshdlrdata->trysolheur != NULL);
4239 
4240  SCIP_CALL( SCIPheurPassSolTrySol(scip, conshdlrdata->trysolheur, newsol) );
4241  }
4242 
4243  SCIP_CALL( SCIPfreeSol(scip, &newsol) );
4244 
4245  return SCIP_OKAY;
4246 }
4247 
4248 /** create a nonlinear row representation of the constraint and stores them in consdata */
4249 static
4251  SCIP* scip, /**< SCIP data structure */
4252  SCIP_CONS* cons /**< absolute power constraint */
4253  )
4254 {
4255  SCIP_CONSDATA* consdata;
4256  SCIP_EXPRTREE* exprtree;
4257  SCIP_EXPRCURV curv;
4258  SCIP_QUADELEM quadelem;
4259  SCIP_VAR* linvars[2];
4260  SCIP_Real lincoefs[2];
4261  SCIP_VAR* quadvar;
4262  SCIP_Real constant;
4263  SCIP_Bool expisint;
4264  int sign;
4265  int nlinvars;
4266  int nquadvars;
4267  int nquadelems;
4268  int n;
4269 
4270  assert(scip != NULL);
4271  assert(cons != NULL);
4272 
4273  consdata = SCIPconsGetData(cons);
4274  assert(consdata != NULL);
4275 
4276  if( consdata->nlrow != NULL )
4277  {
4278  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
4279  }
4280 
4281  nlinvars = 0;
4282  nquadvars = 0;
4283  nquadelems = 0;
4284  exprtree = NULL;
4285  constant = 0.0;
4286 
4287  /* check if sign of x is fixed, determine curvature of abspower function */
4288  if( !SCIPisNegative(scip, SCIPvarGetLbGlobal(consdata->x)+consdata->xoffset) )
4289  {
4290  sign = 1;
4291  curv = SCIP_EXPRCURV_CONVEX;
4292  }
4293  else if( !SCIPisPositive(scip, SCIPvarGetUbGlobal(consdata->x)+consdata->xoffset) )
4294  {
4295  sign = -1;
4296  curv = SCIP_EXPRCURV_CONCAVE;
4297  }
4298  else
4299  {
4300  sign = 0;
4301  curv = SCIP_EXPRCURV_UNKNOWN;
4302  }
4303 
4304  /* check if exponent is integral */
4305  expisint = SCIPisIntegral(scip, consdata->exponent);
4306  n = (int)SCIPround(scip, consdata->exponent);
4307 
4308  /* create quadelem or expression tree for nonlinear part sign(x+offset)abs(x+offset)^n */
4309  if( sign != 0 || (expisint && (n % 2 == 1)) )
4310  {
4311  /* sign is fixes or exponent is odd integer */
4312  if( expisint && n == 2 )
4313  {
4314  /* sign of x is clear and exponent is 2.0 -> generate quadratic, linear, and constant term for +/- (x+offset)^n */
4315  assert(sign == -1 || sign == 1);
4316  nquadelems = 1;
4317  quadelem.idx1 = 0;
4318  quadelem.idx2 = 0;
4319  quadelem.coef = (SCIP_Real)sign;
4320  nquadvars = 1;
4321  quadvar = consdata->x;
4322 
4323  if( consdata->xoffset != 0.0 )
4324  {
4325  linvars[0] = consdata->x;
4326  lincoefs[0] = sign * 2.0 * consdata->xoffset;
4327  nlinvars = 1;
4328  constant = sign * consdata->xoffset * consdata->xoffset;
4329  }
4330  }
4331  else
4332  {
4333  /* exponent is odd or sign of x is clear, generate expression tree for +/- (+/-(x+offset))^exponent */
4334  SCIP_EXPR* expr;
4335  SCIP_EXPR* expr2;
4336 
4337  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, 0) ); /* x */
4338  if( consdata->xoffset != 0.0 )
4339  {
4340  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_CONST, consdata->xoffset) );
4341  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_PLUS, expr, expr2) ); /* x + offset */
4342  }
4343  if( sign == -1 && !expisint )
4344  {
4345  /* if exponent is not integer and x is negative, then negate */
4346  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_CONST, -1.0) ); /* -1 */
4347  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_MUL, expr, expr2) ); /* -(x+offset) */
4348  }
4349  /* use intpower for integer exponent and realpower for fractional exponent */
4350  if( expisint )
4351  {
4352  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_INTPOWER, expr, n) ); /* (x+offset)^n */
4353  }
4354  else
4355  {
4356  assert(sign == 1 || sign == -1);
4357  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_REALPOWER, expr, consdata->exponent) ); /* abs(x+offset)^exponent */
4358  }
4359  /* if exponent is even integer, then negate result; if it's an odd integer, then intpower already takes care of correct sign */
4360  if( sign == -1 && !(expisint && n % 2 == 1) )
4361  {
4362  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_CONST, -1.0) ); /* -1 */
4363  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_MUL, expr, expr2) ); /* -abs(x+offset)^exponent */
4364  }
4365  SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &exprtree, expr, 1, 0, NULL) );
4366  }
4367  }
4368  else
4369  {
4370  /* exponent is not odd integer and sign of x is not fixed -> generate expression tree for signpower(x+offset, n) */
4371  SCIP_EXPR* expr;
4372  SCIP_EXPR* expr2;
4373 
4374  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, 0) ); /* x */
4375  if( consdata->xoffset != 0.0 )
4376  {
4377  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_CONST, consdata->xoffset) );
4378  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_PLUS, expr, expr2) ); /* x + offset */
4379  }
4380  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_SIGNPOWER, expr, (SCIP_Real)consdata->exponent) ); /* signpower(x+offset, n) */
4381 
4382  SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &exprtree, expr, 1, 0, NULL) );
4383  }
4384  assert(exprtree != NULL || nquadelems > 0);
4385 
4386  /* tell expression tree, which is its variable */
4387  if( exprtree != NULL )
4388  {
4389  SCIP_CALL( SCIPexprtreeSetVars(exprtree, 1, &consdata->x) );
4390  }
4391 
4392  assert(nlinvars < 2);
4393  linvars[nlinvars] = consdata->z;
4394  lincoefs[nlinvars] = consdata->zcoef;
4395  ++nlinvars;
4396 
4397  /* create nlrow */
4398  SCIP_CALL( SCIPcreateNlRow(scip, &consdata->nlrow, SCIPconsGetName(cons), constant,
4399  nlinvars, linvars, lincoefs,
4400  nquadvars, &quadvar, nquadelems, &quadelem,
4401  exprtree, consdata->lhs, consdata->rhs,
4402  curv
4403  ) );
4404 
4405  if( exprtree != NULL )
4406  {
4407  SCIP_CALL( SCIPexprtreeFree(&exprtree) );
4408  }
4409 
4410  return SCIP_OKAY;
4411 }
4412 
4413 /** upgrades a quadratic constraint where the quadratic part is only a single square term and the quadratic variable sign is fixed to a signpower constraint */
4414 static
4415 SCIP_DECL_QUADCONSUPGD(quadconsUpgdAbspower)
4416 { /*lint --e{715}*/
4417  SCIP_QUADVARTERM quadvarterm;
4418  SCIP_VAR* x;
4419  SCIP_VAR* z;
4420  SCIP_Real xoffset;
4421  SCIP_Real zcoef;
4422  SCIP_Real signpowcoef;
4423  SCIP_Real lhs;
4424  SCIP_Real rhs;
4425 
4426  *nupgdconss = 0;
4427 
4428  /* need at least one linear variable */
4429  if( SCIPgetNLinearVarsQuadratic(scip, cons) == 0 )
4430  return SCIP_OKAY;
4431 
4432  /* consider only quadratic constraints with a single square term */
4433  if( SCIPgetNQuadVarTermsQuadratic(scip, cons) != 1 )
4434  return SCIP_OKAY;
4435  assert(SCIPgetNBilinTermsQuadratic(scip, cons) == 0);
4436 
4437  quadvarterm = SCIPgetQuadVarTermsQuadratic(scip, cons)[0];
4438  if( SCIPisZero(scip, quadvarterm.sqrcoef) )
4439  return SCIP_OKAY;
4440 
4441  /* don't upgrade if upgrade would scale the constraint down (divide by |sqrcoef|)
4442  * @todo we could still allow this if we were keeping the scaling factor around for the feasibility check
4443  */
4444  if( REALABS(quadvarterm.sqrcoef) > 1.0 )
4445  return SCIP_OKAY;
4446 
4447  x = quadvarterm.var;
4448  xoffset = quadvarterm.lincoef / (2.0 * quadvarterm.sqrcoef);
4449 
4450  /* check that x has fixed sign */
4451  if( SCIPisNegative(scip, SCIPvarGetLbGlobal(x) + xoffset) && SCIPisPositive(scip, SCIPvarGetUbGlobal(x) + xoffset) )
4452  return SCIP_OKAY;
4453 
4454  /* check whether upgdconss array has enough space to store 1 or 2 constraints */
4455  if( SCIPgetNLinearVarsQuadratic(scip, cons) > 1 )
4456  *nupgdconss = -2;
4457  else
4458  *nupgdconss = -1;
4459  if( -*nupgdconss > upgdconsssize )
4460  return SCIP_OKAY;
4461 
4462  *nupgdconss = 0;
4463 
4464  SCIPdebugMsg(scip, "upgrade quadratic constraint <%s> to absolute power, x = [%g,%g], offset = %g\n", SCIPconsGetName(cons), SCIPvarGetLbGlobal(x), SCIPvarGetUbGlobal(x), xoffset);
4465  SCIPdebugPrintCons(scip, cons, NULL);
4466 
4467  lhs = SCIPgetLhsQuadratic(scip, cons);
4468  rhs = SCIPgetRhsQuadratic(scip, cons);
4469 
4470  /* get z and its coefficient */
4471  if( SCIPgetNLinearVarsQuadratic(scip, cons) > 1 )
4472  {
4473  /* create auxiliary variable and constraint for linear part, since we can handle only at most one variable in cons_signpower */
4474  char name[SCIP_MAXSTRLEN];
4475  SCIP_VAR* auxvar;
4476 
4477  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_linpart", SCIPconsGetName(cons));
4478  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS,
4480  SCIP_CALL( SCIPaddVar(scip, auxvar) );
4481 
4482  SCIP_CALL( SCIPcreateConsLinear(scip, &upgdconss[0], name, SCIPgetNLinearVarsQuadratic(scip, cons),
4484  SCIPisInfinity(scip, -lhs) ? -SCIPinfinity(scip) : 0.0,
4485  SCIPisInfinity(scip, rhs) ? SCIPinfinity(scip) : 0.0,
4489  SCIPconsIsStickingAtNode(cons)) );
4490  SCIP_CALL( SCIPaddCoefLinear(scip, upgdconss[*nupgdconss], auxvar, -1.0) );
4491 
4492  z = auxvar;
4493  zcoef = 1.0;
4494 
4495  ++*nupgdconss;
4496 
4497  /* compute and set value of auxvar in debug solution */
4498 #ifdef WITH_DEBUG_SOLUTION
4499  if( SCIPdebugIsMainscip(scip) )
4500  {
4501  SCIP_Real debugval;
4502  SCIP_Real debugvarval;
4503  int i;
4504 
4505  debugval = 0.0;
4506  for( i = 0; i < SCIPgetNLinearVarsQuadratic(scip, cons); ++i )
4507  {
4508  SCIP_CALL( SCIPdebugGetSolVal(scip, SCIPgetLinearVarsQuadratic(scip, cons)[i], &debugvarval) );
4509  debugval += SCIPgetCoefsLinearVarsQuadratic(scip, cons)[i] * debugvarval;
4510  }
4511 
4512  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, debugval) );
4513  }
4514 #endif
4515 
4516  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
4517  }
4518  else
4519  {
4520  assert(SCIPgetNLinearVarsQuadratic(scip, cons) == 1);
4521  z = SCIPgetLinearVarsQuadratic(scip, cons)[0];
4522  zcoef = SCIPgetCoefsLinearVarsQuadratic(scip, cons)[0];
4523  }
4524 
4525  /* we now have lhs <= sqrcoef * (x + offset)^2 - sqrcoef * offset^2 + zcoef * z <= rhs */
4526 
4527  /* move sqrcoef * offset^2 into lhs and rhs */
4528  if( !SCIPisInfinity(scip, -lhs) )
4529  lhs += quadvarterm.sqrcoef * xoffset * xoffset;
4530  if( !SCIPisInfinity(scip, rhs) )
4531  rhs += quadvarterm.sqrcoef * xoffset * xoffset;
4532 
4533  /* divide by sqrcoef if x+offset > 0 and by -sqrcoef if < 0 */
4534  signpowcoef = quadvarterm.sqrcoef;
4535  if( SCIPisNegative(scip, SCIPvarGetLbGlobal(x) + xoffset) )
4536  signpowcoef = -signpowcoef;
4537  if( signpowcoef > 0.0 )
4538  {
4539  if( !SCIPisInfinity(scip, -lhs) )
4540  lhs /= signpowcoef;
4541  if( !SCIPisInfinity(scip, rhs) )
4542  rhs /= signpowcoef;
4543  }
4544  else
4545  {
4546  SCIP_Real newrhs;
4547 
4548  if( !SCIPisInfinity(scip, -lhs) )
4549  newrhs = lhs / signpowcoef;
4550  else
4551  newrhs = SCIPinfinity(scip);
4552  if( !SCIPisInfinity(scip, rhs) )
4553  lhs = rhs / signpowcoef;
4554  else
4555  lhs = -SCIPinfinity(scip);
4556  rhs = newrhs;
4557  }
4558  zcoef /= signpowcoef;
4559 
4560  /* create the absolute power constraint */
4561  SCIP_CALL( SCIPcreateConsAbspower(scip, &upgdconss[*nupgdconss], SCIPconsGetName(cons), x, z, 2.0,
4562  xoffset, zcoef, lhs, rhs,
4566  SCIPconsIsStickingAtNode(cons)) );
4567  SCIPdebugPrintCons(scip, upgdconss[*nupgdconss], NULL);
4568  ++*nupgdconss;
4569 
4570  return SCIP_OKAY;
4571 }
4572 
4573 /** tries to upgrade a nonlinear constraint into a absolute power constraint */
4574 static
4575 SCIP_DECL_NONLINCONSUPGD(nonlinconsUpgdAbspower)
4577  SCIP_EXPRGRAPH* exprgraph;
4578  SCIP_EXPRGRAPHNODE* node;
4579  SCIP_EXPRGRAPHNODE* child;
4580  SCIP_Real exponent;
4581  SCIP_VAR* x;
4582  SCIP_VAR* z;
4583  SCIP_Real signpowcoef;
4584  SCIP_Real zcoef;
4585  SCIP_Real xoffset;
4586  SCIP_Real constant;
4587  SCIP_Real lhs;
4588  SCIP_Real rhs;
4589 
4590  assert(nupgdconss != NULL);
4591  assert(upgdconss != NULL);
4592 
4593  *nupgdconss = 0;
4594 
4595  /* absolute power needs at least one linear variable (constraint is trivial, otherwise) */
4596  if( SCIPgetNLinearVarsNonlinear(scip, cons) == 0 )
4597  return SCIP_OKAY;
4598 
4599  node = SCIPgetExprgraphNodeNonlinear(scip, cons);
4600 
4601  /* no interest in linear constraints */
4602  if( node == NULL )
4603  return SCIP_OKAY;
4604 
4605  /* need exactly one argument */
4606  if( SCIPexprgraphGetNodeNChildren(node) != 1 )
4607  return SCIP_OKAY;
4608 
4609  constant = 0.0;
4610  signpowcoef = 1.0; /* coefficient of sign(x)abs(x)^n term, to be reformulated away... */
4611 
4612  child = SCIPexprgraphGetNodeChildren(node)[0];
4613 
4614  /* check if node expression fits to absolute power constraint */
4615  switch( SCIPexprgraphGetNodeOperator(node) )
4616  {
4617  case SCIP_EXPR_REALPOWER:
4618  /* realpower with exponent > 1.0 can always be signpower, since it assumes that argument is >= 0.0 */
4619  exponent = SCIPexprgraphGetNodeRealPowerExponent(node);
4620  if( exponent <= 1.0 )
4621  return SCIP_OKAY;
4622 
4623  assert(SCIPexprgraphGetNodeBounds(child).inf >= 0.0);
4624  break;
4625 
4626  case SCIP_EXPR_INTPOWER:
4627  {
4628  /* check if exponent > 1.0 and either odd or even with child having fixed sign */
4629  SCIP_INTERVAL childbounds;
4630 
4632  if( exponent <= 1.0 )
4633  return SCIP_OKAY;
4634 
4635  childbounds = SCIPexprgraphGetNodeBounds(child);
4636  if( (int)exponent % 2 == 0 && childbounds.inf < 0.0 && childbounds.sup > 0.0 )
4637  return SCIP_OKAY;
4638 
4639  /* use x^exponent = -sign(x) |x|^exponent if exponent is even and x always negative */
4640  if( (int)exponent % 2 == 0 && childbounds.inf < 0.0 )
4641  signpowcoef = -1.0;
4642 
4643  break;
4644  }
4645 
4646  case SCIP_EXPR_SQUARE:
4647  {
4648  /* check if child has fixed sign */
4649  SCIP_INTERVAL childbounds;
4650 
4651  childbounds = SCIPexprgraphGetNodeBounds(child);
4652  if( childbounds.inf < 0.0 && childbounds.sup > 0.0 )
4653  return SCIP_OKAY;
4654 
4655  /* use x^2 = -sign(x) |x|^2 if x is always negative */
4656  if( childbounds.inf < 0.0 )
4657  signpowcoef = -1.0;
4658 
4659  exponent = 2.0;
4660  break;
4661  }
4662 
4663  case SCIP_EXPR_SIGNPOWER:
4664  /* check if exponent > 1.0 */
4665  exponent = SCIPexprgraphGetNodeSignPowerExponent(node);
4666  if( exponent <= 1.0 )
4667  return SCIP_OKAY;
4668  break;
4669 
4670  case SCIP_EXPR_POLYNOMIAL:
4671  {
4672  SCIP_EXPRDATA_MONOMIAL* monomial;
4673  SCIP_INTERVAL childbounds;
4674 
4675  /* check if only one univariate monomial with exponent > 1.0 */
4676 
4677  /* if sum of univariate monomials, then this should have been taken care of by exprgraphnodeReformSignpower */
4679  return SCIP_OKAY;
4680  assert(SCIPexprgraphGetNodePolynomialNMonomials(node) == 1); /* assume simplified, i.e., no constant polynomial */
4681 
4682  monomial = SCIPexprgraphGetNodePolynomialMonomials(node)[0];
4683  assert(SCIPexprGetMonomialNFactors(monomial) == 1); /* since we have only one children and assume simplified */
4684 
4685  exponent = SCIPexprGetMonomialExponents(monomial)[0];
4686  if( exponent <= 1.0 )
4687  return SCIP_OKAY;
4688 
4689  /* if exponent is even integer and child has mixed sign, then cannot do
4690  * if exponent is even integer and child is always negative, then can do via multiplication by -1.0 */
4691  childbounds = SCIPexprgraphGetNodeBounds(child);
4692  if( SCIPisIntegral(scip, exponent) && ((int)SCIPround(scip, exponent) % 2 == 0) && childbounds.inf < 0.0 )
4693  {
4694  if( childbounds.sup > 0.0 )
4695  return SCIP_OKAY;
4696  signpowcoef = -1.0;
4697  }
4698 
4699  constant = SCIPexprgraphGetNodePolynomialConstant(node);
4700  signpowcoef *= SCIPexprGetMonomialCoef(monomial);
4701 
4702  break;
4703  }
4704 
4705  default:
4706  return SCIP_OKAY;
4707  } /*lint !e788*/
4708  assert(SCIPexprgraphGetNodeNChildren(node) == 1);
4709 
4710  /* check magnitue of coefficient of z in signpower constraint */
4711  zcoef = 1.0;
4712  if( SCIPgetNLinearVarsNonlinear(scip, cons) == 1 )
4713  zcoef = SCIPgetLinearCoefsNonlinear(scip, cons)[0];
4714  zcoef /= signpowcoef;
4716  {
4717  zcoef /= pow(REALABS(SCIPexprgraphGetNodeLinearCoefs(child)[0]), exponent);
4718  }
4719  if( SCIPisZero(scip, zcoef) )
4720  {
4721  SCIPdebugMsg(scip, "skip upgrade to signpower since |zcoef| = %g would be zero\n", zcoef);
4722  return SCIP_OKAY;
4723  }
4724 
4725  /* count how many constraints we need to add (use negative numbers, for convenience):
4726  * one constraint for absolute power,
4727  * plus one if we need to replace the linear part by single variable,
4728  * plus one if we need to replace the argument of absolute power by a single variable
4729  */
4730  *nupgdconss = -1;
4731 
4733  {
4734  /* if node has known curvature and we would add auxiliary var for child, then don't upgrade
4735  * it's not really necessary, but may introduce more numerical troubles
4736  * @todo maybe still do if child is linear?
4737  */
4739  {
4740  *nupgdconss = 0;
4741  return SCIP_OKAY;
4742  }
4743 
4744  --*nupgdconss;
4745  }
4746 
4747  if( SCIPgetNLinearVarsNonlinear(scip, cons) > 1 )
4748  --*nupgdconss;
4749 
4750  /* request larger upgdconss array */
4751  if( upgdconsssize < -*nupgdconss )
4752  return SCIP_OKAY;
4753 
4754  SCIPdebugMsg(scip, "upgrading constraint <%s>\n", SCIPconsGetName(cons));
4755 
4756  /* start counting at zero again */
4757  *nupgdconss = 0;
4758 
4759  exprgraph = SCIPgetExprgraphNonlinear(scip, SCIPconsGetHdlr(cons));
4760 
4761  lhs = SCIPgetLhsNonlinear(scip, cons);
4762  rhs = SCIPgetRhsNonlinear(scip, cons);
4763 
4764  /* get x and it's offset */
4766  {
4767  x = (SCIP_VAR*)SCIPexprgraphGetNodeVar(exprgraph, child);
4768  xoffset = 0.0;
4769  }
4771  {
4772  SCIP_Real xcoef;
4773 
4775  xcoef = SCIPexprgraphGetNodeLinearCoefs(child)[0];
4776  assert(!SCIPisZero(scip, xcoef));
4777 
4778  signpowcoef *= (xcoef < 0.0 ? -1.0 : 1.0) * pow(REALABS(xcoef), exponent);
4779  xoffset = SCIPexprgraphGetNodeLinearConstant(child) / xcoef;
4780  }
4781  else
4782  {
4783  /* reformulate by adding auxiliary variable and constraint for child */
4784  char name[SCIP_MAXSTRLEN];
4785  SCIP_INTERVAL bounds;
4786  SCIP_VAR* auxvar;
4787  SCIP_Real minusone;
4788 
4789  bounds = SCIPexprgraphGetNodeBounds(child);
4790  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_powerarg", SCIPconsGetName(cons));
4791 
4792  SCIPdebugMsg(scip, "add auxiliary variable and constraint %s for node %p(%d,%d)\n", name, (void*)child, SCIPexprgraphGetNodeDepth(child), SCIPexprgraphGetNodePosition(child));
4793 
4794  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, SCIPintervalGetInf(bounds), SCIPintervalGetSup(bounds), 0.0,
4796  SCIP_CALL( SCIPaddVar(scip, auxvar) );
4797 
4798  /* create new constraint child == auxvar
4799  * since signpower is monotonic, we need only child <= auxvar or child >= auxvar, if not both sides are finite, and depending on signpowcoef
4800  * i.e., we need child - auxvar <= 0.0 if rhs is finite and signpowcoef > 0.0 or lhs is finite and signpowcoef < 0.0
4801  * and we need 0.0 <= child - auxvar if lhs is finite and signpowcoef > 0.0 or rhs is finite and signpowcoef < 0.0
4802  */
4803  minusone = -1.0;
4804  assert(upgdconsssize > *nupgdconss);
4805  SCIP_CALL( SCIPcreateConsNonlinear2(scip, &upgdconss[*nupgdconss], name, 1, &auxvar, &minusone, child,
4806  ((signpowcoef > 0.0 && !SCIPisInfinity(scip, -lhs)) || (signpowcoef < 0.0 && !SCIPisInfinity(scip, rhs))) ? 0.0 : -SCIPinfinity(scip),
4807  ((signpowcoef > 0.0 && !SCIPisInfinity(scip, rhs)) || (signpowcoef < 0.0 && !SCIPisInfinity(scip, -lhs))) ? 0.0 : SCIPinfinity(scip),
4808  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
4809  ++*nupgdconss;
4810 
4811  /* use auxvar to setup absolute power constraint */
4812  x = auxvar;
4813  xoffset = 0.0;
4814 
4815  /* compute and set value of auxvar in debug solution, if debugging is enabled */
4816  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, SCIPexprgraphGetNodeVal(child)) ); /*lint !e506 !e774*/
4817 
4818  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
4819  }
4820 
4821  /* get z and its coefficient */
4822  if( SCIPgetNLinearVarsNonlinear(scip, cons) > 1 )
4823  {
4824  /* create auxiliary variable and constraint for linear part, since we can handle only at most one variable in cons_signpower */
4825  char name[SCIP_MAXSTRLEN];
4826  SCIP_VAR* auxvar;
4827 
4828  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_linpart", SCIPconsGetName(cons));
4829  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS,
4831  SCIP_CALL( SCIPaddVar(scip, auxvar) );
4832 
4833  assert(upgdconsssize > *nupgdconss);
4834  SCIP_CALL( SCIPcreateConsLinear(scip, &upgdconss[*nupgdconss], name, SCIPgetNLinearVarsNonlinear(scip, cons),
4836  SCIPisInfinity(scip, -lhs) ? -SCIPinfinity(scip) : 0.0,
4837  SCIPisInfinity(scip, rhs) ? SCIPinfinity(scip) : 0.0,
4841  SCIPconsIsStickingAtNode(cons)) );
4842  SCIP_CALL( SCIPaddCoefLinear(scip, upgdconss[*nupgdconss], auxvar, -1.0) );
4843 
4844  z = auxvar;
4845  zcoef = 1.0;
4846 
4847  ++*nupgdconss;
4848 
4849  /* compute and set value of auxvar in debug solution */
4850 #ifdef WITH_DEBUG_SOLUTION
4851  if( SCIPdebugIsMainscip(scip) )
4852  {
4853  SCIP_Real debugval;
4854  SCIP_Real debugvarval;
4855  int i;
4856 
4857  debugval = 0.0;
4858  for( i = 0; i < SCIPgetNLinearVarsNonlinear(scip, cons); ++i )
4859  {
4860  SCIP_CALL( SCIPdebugGetSolVal(scip, SCIPgetLinearVarsNonlinear(scip, cons)[i], &debugvarval) );
4861  debugval += SCIPgetLinearCoefsNonlinear(scip, cons)[i] * debugvarval;
4862  }
4863 
4864  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, debugval) );
4865  }
4866 #endif
4867 
4868  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
4869  }
4870  else
4871  {
4872  assert(SCIPgetNLinearVarsNonlinear(scip, cons) == 1);
4873  z = SCIPgetLinearVarsNonlinear(scip, cons)[0];
4874  zcoef = SCIPgetLinearCoefsNonlinear(scip, cons)[0];
4875  }
4876 
4877  if( constant != 0.0 )
4878  {
4879  if( !SCIPisInfinity(scip, -lhs) )
4880  lhs -= constant;
4881  if( !SCIPisInfinity(scip, rhs) )
4882  rhs -= constant;
4883  }
4884 
4885  /* divide absolute power constraint by signpowcoef */
4886  if( signpowcoef != 1.0 )
4887  {
4888  zcoef /= signpowcoef;
4889  if( signpowcoef < 0.0 )
4890  {
4891  SCIP_Real newrhs;
4892 
4893  newrhs = SCIPisInfinity(scip, -lhs) ? SCIPinfinity(scip) : lhs/signpowcoef;
4894  lhs = SCIPisInfinity(scip, rhs) ? -SCIPinfinity(scip) : rhs/signpowcoef;
4895  rhs = newrhs;
4896  }
4897  else
4898  {
4899  if( !SCIPisInfinity(scip, -lhs) )
4900  lhs /= signpowcoef;
4901  if( !SCIPisInfinity(scip, rhs) )
4902  rhs /= signpowcoef;
4903  }
4904  }
4905 
4906  /* finally setup a absolute power constraint */
4907 
4908  assert(*nupgdconss < upgdconsssize);
4909  SCIP_CALL( SCIPcreateConsAbspower(scip, &upgdconss[*nupgdconss], SCIPconsGetName(cons),
4910  x, z, exponent, xoffset, zcoef, lhs, rhs,
4914  SCIPconsIsStickingAtNode(cons)) );
4915  ++*nupgdconss;
4916 
4917  return SCIP_OKAY;
4918 }
4919 
4920 /** tries to reformulate a expression graph node via introducing a absolute power constraint
4921  * if node fits to absolute power and has indefinte curvature and has no nonlinear parents and has siblings, then replace by auxvar and absolute power constraint
4922  * if it still has nonlinear parents, then we wait to see if reformulation code move node into auxiliary constraint,
4923  * so we do not add unnessary auxiliary variables for something like an x^2 in an exp(x^2)
4924  * if it has no siblings, then we let the upgrading for nonlinear constraints take care of it,
4925  * since it may be able to upgrade the constraint as a whole and can take the constraint sides into account too (may need only <=/>= auxcons)
4926  */
4927 static
4928 SCIP_DECL_EXPRGRAPHNODEREFORM(exprgraphnodeReformAbspower)
4930  SCIP_EXPRGRAPHNODE* child;
4931  char name[SCIP_MAXSTRLEN];
4932  SCIP_CONS* cons;
4933  SCIP_Real exponent;
4934  SCIP_VAR* x;
4935  SCIP_VAR* z;
4936  SCIP_Real signpowcoef;
4937  SCIP_Real xoffset;
4938  SCIP_Real constant;
4939 
4940  assert(scip != NULL);
4941  assert(exprgraph != NULL);
4942  assert(node != NULL);
4943  assert(naddcons != NULL);
4944  assert(reformnode != NULL);
4945 
4946  *reformnode = NULL;
4947 
4949  return SCIP_OKAY;
4950 
4951  constant = 0.0;
4952  signpowcoef = 1.0; /* coefficient of sign(x)abs(x)^n term, to be move in from of z... */
4953 
4954  /* check if node expression fits to absolute power constraint */
4955  switch( SCIPexprgraphGetNodeOperator(node) )
4956  {
4957  case SCIP_EXPR_REALPOWER:
4958  /* realpower with exponent > 1.0 can always be absolute power, since it assumes that argument is >= 0.0
4959  * @todo we should also ensure that argument is >= 0.0
4960  */
4961  exponent = SCIPexprgraphGetNodeRealPowerExponent(node);
4962  if( exponent <= 1.0 )
4963  return SCIP_OKAY;
4964 
4965  child = SCIPexprgraphGetNodeChildren(node)[0];
4966  assert(SCIPexprgraphGetNodeBounds(child).inf >= 0.0);
4967  break;
4968 
4969  case SCIP_EXPR_INTPOWER:
4970  {
4971  /* check if exponent > 1.0 and either odd or even with child having fixed sign */
4972  SCIP_INTERVAL childbounds;
4973 
4975  if( exponent <= 1.0 )
4976  return SCIP_OKAY;
4977 
4978  child = SCIPexprgraphGetNodeChildren(node)[0];
4979  childbounds = SCIPexprgraphGetNodeBounds(child);
4980  if( (int)exponent % 2 == 0 && childbounds.inf < 0.0 && childbounds.sup > 0.0 )
4981  return SCIP_OKAY;
4982 
4983  /* use x^exponent = -sign(x) |x|^exponent if exponent is even and x always negative */
4984  if( (int)exponent % 2 == 0 && childbounds.inf < 0.0 )
4985  signpowcoef = -1.0;
4986 
4987  break;
4988  }
4989 
4990  case SCIP_EXPR_SQUARE:
4991  {
4992  /* check if child has fixed sign */
4993  SCIP_INTERVAL childbounds;
4994 
4995  child = SCIPexprgraphGetNodeChildren(node)[0];
4996  childbounds = SCIPexprgraphGetNodeBounds(child);
4997  if( childbounds.inf < 0.0 && childbounds.sup > 0.0 )
4998  return SCIP_OKAY;
4999 
5000  /* use x^2 = -sign(x) |x|^2 if x is always negative */
5001  if( childbounds.inf < 0.0 )
5002  signpowcoef = -1.0;
5003 
5004  exponent = 2.0;
5005  break;
5006  }
5007 
5008  case SCIP_EXPR_SIGNPOWER:
5009  /* check if exponent > 1.0 */
5010  exponent = SCIPexprgraphGetNodeSignPowerExponent(node);
5011  if( exponent <= 1.0 )
5012  return SCIP_OKAY;
5013  child = SCIPexprgraphGetNodeChildren(node)[0];
5014  break;
5015 
5016  case SCIP_EXPR_POLYNOMIAL:
5017  {
5018  SCIP_EXPRDATA_MONOMIAL* monomial;
5019  SCIP_INTERVAL childbounds;
5020 
5021  /* check if only one monomial */
5023  return SCIP_OKAY;
5024  assert(SCIPexprgraphGetNodePolynomialNMonomials(node) == 1); /* assume simplified, i.e., no constant polynomial */
5025 
5026  monomial = SCIPexprgraphGetNodePolynomialMonomials(node)[0];
5027 
5028  assert(SCIPexprgraphGetNodeNChildren(node) >= 1);
5029  if( SCIPexprgraphGetNodeNChildren(node) == 1 )
5030  {
5031  /* handle univariate case with exponent > 1.0 */
5032  assert(SCIPexprGetMonomialNFactors(monomial) == 1); /* since we have only one child and assume simplified */
5033 
5034  exponent = SCIPexprGetMonomialExponents(monomial)[0];
5035  if( exponent <= 1.0 )
5036  return SCIP_OKAY;
5037 
5038  /* if exponent is even integer and child has mixed sign, then cannot do
5039  * if exponent is even integer and child is always negative, then can do via multiplication by -1.0 */
5040  child = SCIPexprgraphGetNodeChildren(node)[0];
5041  childbounds = SCIPexprgraphGetNodeBounds(child);
5042  if( SCIPisIntegral(scip, exponent) && ((int)SCIPround(scip, exponent) % 2 == 0) && childbounds.inf < 0.0 )
5043  {
5044  if( childbounds.sup > 0.0 )
5045  return SCIP_OKAY;
5046  signpowcoef = -1.0;
5047  }
5048 
5049  constant = SCIPexprgraphGetNodePolynomialConstant(node);
5050  signpowcoef *= SCIPexprGetMonomialCoef(monomial);
5051  }
5052  else if( SCIPexprgraphGetNodeNChildren(node) == 2 )
5053  {
5054  /* look for abs(x)^p * x or x * abs(x)^p */
5055  SCIP_EXPRGRAPHNODE* child1;
5056  SCIP_EXPRGRAPHNODE* child2;
5057 
5058  assert(SCIPexprGetMonomialNFactors(monomial) == 2); /* since we have two children and assume simplified */
5059 
5062 
5063  if( SCIPexprgraphGetNodeOperator(child1) == SCIP_EXPR_ABS && SCIPexprgraphGetNodeChildren(child1)[0] == child2 )
5064  {
5065  /* have abs(x)^p1 * x^p2; need p1 > 0, p2 = 1 */
5066  if( SCIPexprGetMonomialExponents(monomial)[0] <= 0.0 || SCIPexprGetMonomialExponents(monomial)[1] != 1.0 )
5067  return SCIP_OKAY;
5068  exponent = 1.0 + SCIPexprGetMonomialExponents(monomial)[0];
5069  child = child2;
5070  }
5071  else if( SCIPexprgraphGetNodeOperator(child2) == SCIP_EXPR_ABS && SCIPexprgraphGetNodeChildren(child2)[0] == child1 )
5072  {
5073  /* have x^p1 * abs(x)^p2; need p1 = 1, p2 > 0 */
5074  if( SCIPexprGetMonomialExponents(monomial)[0] != 1.0 || SCIPexprGetMonomialExponents(monomial)[1] <= 0.0 )
5075  return SCIP_OKAY;
5076  exponent = 1.0 + SCIPexprGetMonomialExponents(monomial)[1];
5077  child = child1;
5078  }
5079  else
5080  return SCIP_OKAY;
5081 
5082  constant = SCIPexprgraphGetNodePolynomialConstant(node);
5083  signpowcoef = SCIPexprGetMonomialCoef(monomial);
5084  }
5085  else
5086  return SCIP_OKAY;
5087 
5088  break;
5089  }
5090 
5091  case SCIP_EXPR_MUL:
5092  {
5093  /* look for abs(x) * x and x * abs(x) */
5094  SCIP_EXPRGRAPHNODE* child1;
5095  SCIP_EXPRGRAPHNODE* child2;
5096 
5097  child1 = SCIPexprgraphGetNodeChildren(node)[0];
5098  child2 = SCIPexprgraphGetNodeChildren(node)[1];
5099 
5100  if( SCIPexprgraphGetNodeNChildren(child1) == 1 && SCIPexprgraphGetNodeOperator(child1) == SCIP_EXPR_ABS && SCIPexprgraphGetNodeChildren(child1)[0] == child2 )
5101  {
5102  /* abs(x) * x */
5103  exponent = 2.0;
5104  child = child2;
5105  break;
5106  }
5107 
5108  if( SCIPexprgraphGetNodeNChildren(child2) == 1 && SCIPexprgraphGetNodeOperator(child2) == SCIP_EXPR_ABS && SCIPexprgraphGetNodeChildren(child2)[0] == child1 )
5109  {
5110  /* x * abs(x) */
5111  exponent = 2.0;
5112  child = child1;
5113  break;
5114  }
5115 
5116  /* the cases below should not happening, because node will be of polynomial type after simplifty */
5117 #ifdef SCIP_DISABLED_CODE
5120  {
5121  /* abs(x)^p * x */
5123  child = child2;
5124  break;
5125  }
5126 
5129  {
5130  /* x * abs(x)^p */
5132  child = child1;
5133  break;
5134  }
5135 #endif
5136  return SCIP_OKAY;
5137  }
5138 
5139  default:
5140  return SCIP_OKAY;
5141  } /*lint !e788*/
5142 
5144  return SCIP_OKAY;
5145  if( SCIPexprgraphGetNodeNChildren(node) == 1 && !SCIPexprgraphHasNodeSibling(node) ) /* why did I require siblings??? */
5146  return SCIP_OKAY;
5147 
5148  SCIPdebugMsg(scip, "reformulate node %p via signpower\n", (void*)node);
5149 
5150  /* get x and its offset */
5152  {
5153  x = (SCIP_VAR*)SCIPexprgraphGetNodeVar(exprgraph, child);
5154  xoffset = 0.0;
5155  }
5157  {
5158  SCIP_Real xcoef;
5159 
5161  xcoef = SCIPexprgraphGetNodeLinearCoefs(child)[0];
5162  assert(!SCIPisZero(scip, xcoef));
5163 
5164  signpowcoef *= (xcoef < 0.0 ? -1.0 : 1.0) * pow(REALABS(xcoef), exponent);
5165  xoffset = SCIPexprgraphGetNodeLinearConstant(child) / xcoef;
5166  }
5167  else
5168  {
5169  /* reformulate by adding auxiliary variable and constraint for child */
5170  SCIP_INTERVAL bounds;
5171  SCIP_VAR* auxvar;
5172  SCIP_Real minusone;
5173 
5174  bounds = SCIPexprgraphGetNodeBounds(child);
5175  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "nlreform%dsp", *naddcons);
5176 
5177  SCIPdebugMsg(scip, "add auxiliary variable and constraint %s for node %p(%d,%d)\n", name, (void*)child, SCIPexprgraphGetNodeDepth(child), SCIPexprgraphGetNodePosition(child));
5178 
5179  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, SCIPintervalGetInf(bounds), SCIPintervalGetSup(bounds), 0.0, SCIP_VARTYPE_CONTINUOUS,
5180  TRUE, TRUE, NULL, NULL, NULL, NULL, NULL) );
5181  SCIP_CALL( SCIPaddVar(scip, auxvar) );
5182 
5183  /* create new constraint child == auxvar */
5184  minusone = -1.0;
5185  SCIP_CALL( SCIPcreateConsNonlinear2(scip, &cons, name, 1, &auxvar, &minusone, child, 0.0, 0.0,
5186  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
5187  SCIP_CALL( SCIPaddCons(scip, cons) );
5188  ++*naddcons;
5189 
5190  /* use auxvar to setup signpower constraint */
5191  x = auxvar;
5192  xoffset = 0.0;
5193 
5194  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, SCIPexprgraphGetNodeVal(child)) ); /*lint !e506 !e774*/
5195 
5196  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
5197  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
5198  }
5199 
5200  /* create auxiliary variable z and add to expression graph */
5201  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "nlreform%dsp", *naddcons);
5202  SCIP_CALL( SCIPcreateVar(scip, &z, name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS,
5203  TRUE, TRUE, NULL, NULL, NULL, NULL, NULL) );
5204  SCIP_CALL( SCIPaddVar(scip, z) );
5205  SCIP_CALL( SCIPexprgraphAddVars(exprgraph, 1, (void**)&z, reformnode) );
5206 
5207  /* setup a absolute power constraint */
5208  if( REALABS(signpowcoef) * SCIPfeastol(scip) < 1.0 )
5209  {
5210  /* if signpowcoef is not huge (<10^6), then put it into absolute power constraint */
5211  SCIP_CALL( SCIPcreateConsAbspower(scip, &cons, name,
5212  x, z, exponent, xoffset, -1.0/signpowcoef, -constant/signpowcoef, -constant/signpowcoef,
5213  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
5214  SCIP_CALL( SCIPaddCons(scip, cons) );
5215  SCIPdebugPrintCons(scip, cons, NULL);
5216  ++*naddcons;
5217 
5218  /* compute value of z and reformnode and set in debug solution and expression graph, resp. */
5219 #ifdef WITH_DEBUG_SOLUTION
5220  if( SCIPdebugIsMainscip(scip) )
5221  {
5222  SCIP_Real xval;
5223  SCIP_Real zval;
5224 
5225  SCIP_CALL( SCIPdebugGetSolVal(scip, x, &xval) );
5226  zval = signpowcoef * SIGN(xval + xoffset) * pow(REALABS(xval + xoffset), exponent) + constant;
5227 
5228  SCIP_CALL( SCIPdebugAddSolVal(scip, z, zval) );
5229  SCIPexprgraphSetVarNodeValue(*reformnode, zval);
5230  }
5231 #endif
5232  }
5233  else
5234  {
5235  /* if signpowcoef is huge, then avoid very small coefficient of z
5236  * instead create additional node on top of current reformnode */
5237  SCIP_EXPRGRAPHNODE* linnode;
5238 
5239  SCIP_CALL( SCIPcreateConsAbspower(scip, &cons, name,
5240  x, z, exponent, xoffset, -1.0, 0.0, 0.0,
5241  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
5242  SCIP_CALL( SCIPaddCons(scip, cons) );
5243  SCIPdebugPrintCons(scip, cons, NULL);
5244  ++*naddcons;
5245 
5246  /* compute value of z and reformnode and set in debug solution and expression graph, resp. */
5247 #ifdef WITH_DEBUG_SOLUTION
5248  if( SCIPdebugIsMainscip(scip) )
5249  {
5250  SCIP_Real xval;
5251  SCIP_Real zval;
5252 
5253  SCIP_CALL( SCIPdebugGetSolVal(scip, x, &xval) );
5254  zval = SIGN(xval + xoffset) * pow(REALABS(xval + xoffset), exponent);
5255 
5256  SCIP_CALL( SCIPdebugAddSolVal(scip, z, zval) );
5257  SCIPexprgraphSetVarNodeValue(*reformnode, zval);
5258  }
5259 #endif
5260 
5261  SCIP_CALL( SCIPexprgraphCreateNodeLinear(SCIPblkmem(scip), &linnode, 1, &signpowcoef, constant) );
5262  SCIP_CALL( SCIPexprgraphAddNode(exprgraph, linnode, -1, 1, reformnode) );
5263 
5264  *reformnode = linnode;
5265  }
5266 
5267  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
5268  SCIP_CALL( SCIPreleaseVar(scip, &z) );
5269 
5270  return SCIP_OKAY;
5271 }
5272 
5273 /** helper function to enforce constraints */
5274 static
5276  SCIP* scip, /**< SCIP data structure */
5277  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
5278  SCIP_CONS** conss, /**< constraints to process */
5279  int nconss, /**< number of constraints */
5280  int nusefulconss, /**< number of useful (non-obsolete) constraints to process */
5281  SCIP_SOL* sol, /**< solution to enforce (NULL for the LP solution) */
5282  SCIP_Bool solinfeasible, /**< was the solution already declared infeasible by a constraint handler? */
5283  SCIP_RESULT* result /**< pointer to store the result of the enforcing call */
5284  )
5285 {
5286  SCIP_CONSHDLRDATA* conshdlrdata;
5287  SCIP_CONS* maxviolcons;
5288  SCIP_CONSDATA* consdata;
5289  SCIP_Bool success;
5290  SCIP_Bool cutoff;
5291  SCIP_Bool solviolbounds;
5292  SCIP_Real sepaefficacy;
5293  SCIP_Real maxviol;
5294  int nnotify;
5295  int c;
5296 
5297  assert(scip != NULL);
5298  assert(conshdlr != NULL);
5299  assert(conss != NULL || nconss == 0);
5300  assert(result != NULL);
5301 
5302  conshdlrdata = SCIPconshdlrGetData(conshdlr);
5303  assert(conshdlrdata != NULL);
5304 
5305  SCIP_CALL( computeViolations(scip, conshdlr, conss, nconss, sol, &solviolbounds, &maxviolcons) );
5306 
5307  if( maxviolcons == NULL )
5308  {
5309  *result = SCIP_FEASIBLE;
5310  return SCIP_OKAY;
5311  }
5312 
5313  *result = SCIP_INFEASIBLE;
5314 
5315  if( solviolbounds )
5316  {
5317  /* if LP solution violates variable bounds, then this should be because a row was added that
5318  * introduced this variable newly to the LP, in which case it gets value 0.0; the row should
5319  * have been added to resolve an infeasibility, so solinfeasible should be TRUE
5320  * see also issue #627
5321  */
5322  assert(solinfeasible);
5323  /* however, if solinfeasible is actually not TRUE, then better cut off the node to avoid that SCIP
5324  * stops because infeasible cannot be resolved */
5325  /*lint --e{774} */
5326  if( !solinfeasible )
5327  *result = SCIP_CUTOFF;
5328  return SCIP_OKAY;
5329  }
5330 
5331  /* if we are above the 100'th enforcement round for this node, something is strange
5332  * (maybe the LP does not think that the cuts we add are violated, or we do ECP on a high-dimensional convex function)
5333  * in this case, check if some limit is hit or SCIP should stop for some other reason and terminate enforcement by creating a dummy node
5334  * (in optimized more, returning SCIP_INFEASIBLE in *result would be sufficient, but in debug mode this would give an assert in scip.c)
5335  * the reason to wait for 100 rounds is to avoid calls to SCIPisStopped in normal runs, which may be expensive
5336  * we only increment nenforounds until 101 to avoid an overflow
5337  */