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 scip.zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15
16 /**@file nlhdlr_bilinear.c
17 * @ingroup DEFPLUGINS_NLHDLR
18 * @brief bilinear nonlinear handler
19 * @author Benjamin Mueller
20 */
21
22 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
23
24 #include <string.h>
25
26 #include "scip/nlhdlr_bilinear.h"
27 #include "scip/cons_nonlinear.h"
28 #include "scip/expr_product.h"
29 #include "scip/expr_var.h"
30
31 /* fundamental nonlinear handler properties */
32 #define NLHDLR_NAME "bilinear"
33 #define NLHDLR_DESC "bilinear handler for expressions"
34 #define NLHDLR_DETECTPRIORITY -10 /**< it is important that the nlhdlr runs after the default nldhlr */
35 #define NLHDLR_ENFOPRIORITY -10
36
37 #define MIN_INTERIORITY 0.01 /**< minimum interiority for a reference point for applying separation */
38 #define MIN_ABSBOUNDSIZE 0.1 /**< minimum size of variable bounds for applying separation */
39
40 /* properties of the bilinear nlhdlr statistics table */
41 #define TABLE_NAME_BILINEAR "nlhdlr_bilinear"
42 #define TABLE_DESC_BILINEAR "bilinear nlhdlr statistics table"
43 #define TABLE_POSITION_BILINEAR 14800 /**< the position of the statistics table */
44 #define TABLE_EARLIEST_STAGE_BILINEAR SCIP_STAGE_INITSOLVE /**< output of the statistics table is only printed from this stage onwards */
45
46
47 /*
48 * Data structures
49 */
50
51 /** nonlinear handler expression data */
52 struct SCIP_NlhdlrExprData
53 {
54 SCIP_Real underineqs[6]; /**< inequalities for underestimation */
55 int nunderineqs; /**< total number of inequalities for underestimation */
56 SCIP_Real overineqs[6]; /**< inequalities for overestimation */
57 int noverineqs; /**< total number of inequalities for overestimation */
58 SCIP_Longint lastnodeid; /**< id of the last node that has been used for separation */
59 int nseparoundslastnode; /**< number of separation calls of the last node */
60 };
61
62 /** nonlinear handler data */
63 struct SCIP_NlhdlrData
64 {
65 SCIP_EXPR** exprs; /**< expressions that have been detected by the nlhdlr */
66 int nexprs; /**< total number of expression that have been detected */
67 int exprsize; /**< size of exprs array */
68 SCIP_HASHMAP* exprmap; /**< hashmap to store the position of each expression in the exprs array */
69
70 /* parameter */
71 SCIP_Bool useinteval; /**< whether to use the interval evaluation callback of the nlhdlr */
72 SCIP_Bool usereverseprop; /**< whether to use the reverse propagation callback of the nlhdlr */
73 int maxseparoundsroot; /**< maximum number of separation rounds in the root node */
74 int maxseparounds; /**< maximum number of separation rounds in a local node */
75 int maxsepadepth; /**< maximum depth to apply separation */
76 };
77
78 /*
79 * Local methods
80 */
81
82 /** helper function to compute the violation of an inequality of the form xcoef * x <= ycoef * y + constant for two
83 * corner points of the domain [lbx,ubx] x [lby,uby]
84 */
85 static
86 void getIneqViol(
87 SCIP_VAR* x, /**< first variable */
88 SCIP_VAR* y, /**< second variable */
89 SCIP_Real xcoef, /**< x-coefficient */
90 SCIP_Real ycoef, /**< y-coefficient */
91 SCIP_Real constant, /**< constant */
92 SCIP_Real* viol1, /**< buffer to store the violation of the first corner point */
93 SCIP_Real* viol2 /**< buffer to store the violation of the second corner point */
94 )
95 {
96 SCIP_Real norm;
97 assert(viol1 != NULL);
98 assert(viol2 != NULL);
99
100 norm = SQRT(SQR(xcoef) + SQR(ycoef));
101
102 /* inequality can be used for underestimating xy if and only if xcoef * ycoef > 0 */
103 if( xcoef * ycoef >= 0 )
104 {
105 /* violation for top-left and bottom-right corner */
106 *viol1 = MAX(0, (xcoef * SCIPvarGetLbLocal(x) - ycoef * SCIPvarGetUbLocal(y) - constant) / norm); /*lint !e666*/
107 *viol2 = MAX(0, (xcoef * SCIPvarGetUbLocal(x) - ycoef * SCIPvarGetLbLocal(y) - constant) / norm); /*lint !e666*/
108 }
109 else
110 {
111 /* violation for top-right and bottom-left corner */
112 *viol1 = MAX(0, (xcoef * SCIPvarGetUbLocal(x) - ycoef * SCIPvarGetUbLocal(y) - constant) / norm); /*lint !e666*/
113 *viol2 = MAX(0, (xcoef * SCIPvarGetLbLocal(x) - ycoef * SCIPvarGetLbLocal(y) - constant) / norm); /*lint !e666*/
114 }
115 }
116
117 /** auxiliary function to decide whether to use inequalities for a strong relaxation of bilinear terms or not */
118 static
119 SCIP_Bool useBilinIneqs(
120 SCIP* scip, /**< SCIP data structure */
121 SCIP_VAR* x, /**< x variable */
122 SCIP_VAR* y, /**< y variable */
123 SCIP_Real refx, /**< reference point for x */
124 SCIP_Real refy /**< reference point for y */
125 )
126 {
127 SCIP_Real lbx;
128 SCIP_Real ubx;
129 SCIP_Real lby;
130 SCIP_Real uby;
131 SCIP_Real interiorityx;
132 SCIP_Real interiorityy;
133 SCIP_Real interiority;
134
135 assert(x != NULL);
136 assert(y != NULL);
137 assert(x != y);
138
139 /* get variable bounds */
140 lbx = SCIPvarGetLbLocal(x);
141 ubx = SCIPvarGetUbLocal(x);
142 lby = SCIPvarGetLbLocal(y);
143 uby = SCIPvarGetUbLocal(y);
144
145 /* compute interiority */
146 interiorityx = MIN(refx-lbx, ubx-refx) / MAX(ubx-lbx, SCIPepsilon(scip)); /*lint !e666*/
147 interiorityy = MIN(refy-lby, uby-refy) / MAX(uby-lby, SCIPepsilon(scip)); /*lint !e666*/
148 interiority = 2.0*MIN(interiorityx, interiorityy);
149
150 return ubx - lbx >= MIN_ABSBOUNDSIZE && uby - lby >= MIN_ABSBOUNDSIZE && interiority >= MIN_INTERIORITY;
151 }
152
153 /** helper function to update the best relaxation for a bilinear term when using valid linear inequalities */
154 static
155 void updateBilinearRelaxation(
156 SCIP* scip, /**< SCIP data structure */
157 SCIP_VAR* RESTRICT x, /**< first variable */
158 SCIP_VAR* RESTRICT y, /**< second variable */
159 SCIP_Real bilincoef, /**< coefficient of the bilinear term */
160 SCIP_SIDETYPE violside, /**< side of quadratic constraint that is violated */
161 SCIP_Real refx, /**< reference point for the x variable */
162 SCIP_Real refy, /**< reference point for the y variable */
163 SCIP_Real* RESTRICT ineqs, /**< coefficients of each linear inequality; stored as triple (xcoef,ycoef,constant) */
164 int nineqs, /**< total number of inequalities */
165 SCIP_Real mccormickval, /**< value of the McCormick relaxation at the reference point */
166 SCIP_Real* RESTRICT bestcoefx, /**< pointer to update the x coefficient */
167 SCIP_Real* RESTRICT bestcoefy, /**< pointer to update the y coefficient */
168 SCIP_Real* RESTRICT bestconst, /**< pointer to update the constant */
169 SCIP_Real* RESTRICT bestval, /**< value of the best relaxation that have been found so far */
170 SCIP_Bool* success /**< buffer to store whether we found a better relaxation */
171 )
172 {
173 SCIP_Real constshift[2] = {0.0, 0.0};
174 SCIP_Real constant;
175 SCIP_Real xcoef;
176 SCIP_Real ycoef;
177 SCIP_Real lbx;
178 SCIP_Real ubx;
179 SCIP_Real lby;
180 SCIP_Real uby;
181 SCIP_Bool update;
182 SCIP_Bool overestimate;
183 int i;
184
185 assert(x != y);
186 assert(!SCIPisZero(scip, bilincoef));
187 assert(nineqs >= 0 && nineqs <= 2);
188 assert(bestcoefx != NULL);
189 assert(bestcoefy != NULL);
190 assert(bestconst != NULL);
191 assert(bestval != NULL);
192
193 /* no inequalities available */
194 if( nineqs == 0 )
195 return;
196 assert(ineqs != NULL);
197
198 lbx = SCIPvarGetLbLocal(x);
199 ubx = SCIPvarGetUbLocal(x);
200 lby = SCIPvarGetLbLocal(y);
201 uby = SCIPvarGetUbLocal(y);
202 overestimate = (violside == SCIP_SIDETYPE_LEFT);
203
204 /* check cases for which we can't compute a tighter relaxation */
205 if( SCIPisFeasLE(scip, refx, lbx) || SCIPisFeasGE(scip, refx, ubx)
206 || SCIPisFeasLE(scip, refy, lby) || SCIPisFeasGE(scip, refy, uby) )
207 return;
208
209 /* due to the feasibility tolerances of the LP and NLP solver, it might possible that the reference point is
210 * violating the linear inequalities; to ensure that we compute a valid underestimate, we relax the linear
211 * inequality by changing its constant part
212 */
213 for( i = 0; i < nineqs; ++i )
214 {
215 constshift[i] = MAX(0.0, ineqs[3*i] * refx - ineqs[3*i+1] * refy - ineqs[3*i+2]);
216 SCIPdebugMsg(scip, "constant shift of inequality %d = %.16f\n", i, constshift[i]);
217 }
218
219 /* try to use both inequalities */
220 if( nineqs == 2 )
221 {
222 SCIPcomputeBilinEnvelope2(scip, bilincoef, lbx, ubx, refx, lby, uby, refy, overestimate, ineqs[0], ineqs[1],
223 ineqs[2] + constshift[0], ineqs[3], ineqs[4], ineqs[5] + constshift[1], &xcoef, &ycoef, &constant, &update);
224
225 if( update )
226 {
227 SCIP_Real val = xcoef * refx + ycoef * refy + constant;
228 SCIP_Real relimpr = 1.0 - (REALABS(val - bilincoef * refx * refy) + 1e-4) / (REALABS(*bestval - bilincoef * refx * refy) + 1e-4);
229 SCIP_Real absimpr = REALABS(val - (*bestval));
230
231 /* update relaxation if possible */
232 if( relimpr > 0.05 && absimpr > 1e-3 && ((overestimate && SCIPisRelLT(scip, val, *bestval))
233 || (!overestimate && SCIPisRelGT(scip, val, *bestval))) )
234 {
235 *bestcoefx = xcoef;
236 *bestcoefy = ycoef;
237 *bestconst = constant;
238 *bestval = val;
239 *success = TRUE;
240 }
241 }
242 }
243
244 /* use inequalities individually */
245 for( i = 0; i < nineqs; ++i )
246 {
247 SCIPcomputeBilinEnvelope1(scip, bilincoef, lbx, ubx, refx, lby, uby, refy, overestimate, ineqs[3*i], ineqs[3*i+1],
248 ineqs[3*i+2] + constshift[i], &xcoef, &ycoef, &constant, &update);
249
250 if( update )
251 {
252 SCIP_Real val = xcoef * refx + ycoef * refy + constant;
253 SCIP_Real relimpr = 1.0 - (REALABS(val - bilincoef * refx * refy) + 1e-4)
254 / (REALABS(mccormickval - bilincoef * refx * refy) + 1e-4);
255 SCIP_Real absimpr = REALABS(val - (*bestval));
256
257 /* update relaxation if possible */
258 if( relimpr > 0.05 && absimpr > 1e-3 && ((overestimate && SCIPisRelLT(scip, val, *bestval))
259 || (!overestimate && SCIPisRelGT(scip, val, *bestval))) )
260 {
261 *bestcoefx = xcoef;
262 *bestcoefy = ycoef;
263 *bestconst = constant;
264 *bestval = val;
265 *success = TRUE;
266 }
267 }
268 }
269 }
270
271 /** helper function to determine whether a given point satisfy given inequalities */
272 static
273 SCIP_Bool isPointFeasible(
274 SCIP* scip, /**< SCIP data structure */
275 SCIP_Real x, /**< x-coordinate */
276 SCIP_Real y, /**< y-coordinate */
277 SCIP_Real lbx, /**< lower bound of x */
278 SCIP_Real ubx, /**< upper bound of x */
279 SCIP_Real lby, /**< lower bound of y */
280 SCIP_Real uby, /**< upper bound of y */
281 SCIP_Real* ineqs, /**< inequalities of the form coefx x <= coefy y + constant */
282 int nineqs /**< total number of inequalities */
283 )
284 {
285 int i;
286
287 assert(ineqs != NULL);
288 assert(nineqs > 0);
289
290 /* check whether point satisfies the bounds */
291 if( SCIPisLT(scip, x, lbx) || SCIPisGT(scip, x, ubx)
292 || SCIPisLT(scip, y, lby) || SCIPisGT(scip, y, uby) )
293 return FALSE;
294
295 /* check whether point satisfy the linear inequalities */
296 for( i = 0; i < nineqs; ++i )
297 {
298 SCIP_Real coefx = ineqs[3*i];
299 SCIP_Real coefy = ineqs[3*i+1];
300 SCIP_Real constant = ineqs[3*i+2];
301
302 /* TODO check with an absolute comparison? */
303 if( SCIPisGT(scip, coefx*x - coefy*y - constant, 0.0) )
304 return FALSE;
305 }
306
307 return TRUE;
308 }
309
310 /** helper function for computing all vertices of the polytope described by the linear inequalities and the local
311 * extrema of the bilinear term along each inequality
312 *
313 * @note there are at most 22 points where the min/max can be achieved (given that there are at most 4 inequalities)
314 * - corners of [lbx,ubx]x[lby,uby] (4)
315 * - two intersection points for each inequality with the box (8)
316 * - global maximum / minimum on each inequality (4)
317 * - intersection between two inequalities (6)
318 */
319 static
320 void getFeasiblePointsBilinear(
321 SCIP* scip, /**< SCIP data structure */
322 SCIP_CONSHDLR* conshdlr, /**< constraint handler, if levelset == TRUE, otherwise can be NULL */
323 SCIP_EXPR* expr, /**< product expression */
324 SCIP_INTERVAL exprbounds, /**< bounds on product expression, only used if levelset == TRUE */
325 SCIP_Real* underineqs, /**< inequalities for underestimation */
326 int nunderineqs, /**< total number of inequalities for underestimation */
327 SCIP_Real* overineqs, /**< inequalities for overestimation */
328 int noverineqs, /**< total number of inequalities for overestimation */
329 SCIP_Bool levelset, /**< should the level set be considered? */
330 SCIP_Real* xs, /**< array to store x-coordinates of computed points */
331 SCIP_Real* ys, /**< array to store y-coordinates of computed points */
332 int* npoints /**< buffer to store the total number of computed points */
333 )
334 {
335 SCIP_EXPR* child1;
336 SCIP_EXPR* child2;
337 SCIP_Real ineqs[12];
338 SCIP_INTERVAL boundsx;
339 SCIP_INTERVAL boundsy;
340 SCIP_Real lbx;
341 SCIP_Real ubx;
342 SCIP_Real lby;
343 SCIP_Real uby;
344 int nineqs = 0;
345 int i;
346
347 assert(scip != NULL);
348 assert(conshdlr != NULL || !levelset);
349 assert(expr != NULL);
350 assert(xs != NULL);
351 assert(ys != NULL);
352 assert(SCIPexprGetNChildren(expr) == 2);
353 assert(noverineqs + nunderineqs > 0);
354 assert(noverineqs + nunderineqs <= 4);
355
356 *npoints = 0;
357
358 /* collect inequalities */
359 for( i = 0; i < noverineqs; ++i )
360 {
361 ineqs[3*nineqs] = overineqs[3*i];
362 ineqs[3*nineqs+1] = overineqs[3*i+1];
363 ineqs[3*nineqs+2] = overineqs[3*i+2];
364 ++nineqs;
365 }
366 for( i = 0; i < nunderineqs; ++i )
367 {
368 ineqs[3*nineqs] = underineqs[3*i];
369 ineqs[3*nineqs+1] = underineqs[3*i+1];
370 ineqs[3*nineqs+2] = underineqs[3*i+2];
371 ++nineqs;
372 }
373 assert(nineqs == noverineqs + nunderineqs);
374
375 /* collect children */
376 child1 = SCIPexprGetChildren(expr)[0];
377 child2 = SCIPexprGetChildren(expr)[1];
378 assert(child1 != NULL && child2 != NULL);
379 assert(child1 != child2);
380
381 /* collect bounds of children */
382 if( !levelset )
383 {
384 /* if called from inteval, then use activity */
385 boundsx = SCIPexprGetActivity(child1);
386 boundsy = SCIPexprGetActivity(child2);
387 }
388 else
389 {
390 /* if called from reverseprop, then use bounds */
391 boundsx = SCIPgetExprBoundsNonlinear(scip, child1);
392 boundsy = SCIPgetExprBoundsNonlinear(scip, child2);
393
394 /* if children bounds are empty, then returning with *npoints==0 is the way to go */
395 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, boundsx)
396 || SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, boundsy) )
397 return;
398 }
399 lbx = boundsx.inf;
400 ubx = boundsx.sup;
401 lby = boundsy.inf;
402 uby = boundsy.sup;
403
404 /* corner points that satisfy all inequalities */
405 for( i = 0; i < 4; ++i )
406 {
407 SCIP_Real cx = i < 2 ? lbx : ubx;
408 SCIP_Real cy = (i % 2) == 0 ? lby : uby;
409
410 if( isPointFeasible(scip, cx, cy, lbx, ubx, lby, uby, ineqs, nineqs) )
411 {
412 xs[*npoints] = cx;
413 ys[*npoints] = cy;
414 ++(*npoints);
415 }
416 }
417
418 /* intersection point of inequalities with [lbx,ubx] x [lby,uby] and extremum of xy on each inequality */
419 for( i = 0; i < nineqs; ++i )
420 {
421 SCIP_Real coefx = ineqs[3*i];
422 SCIP_Real coefy = ineqs[3*i+1];
423 SCIP_Real constant = ineqs[3*i+2];
424 SCIP_Real px[5] = {lbx, ubx, (coefy*lby + constant)/coefx, (coefy*uby + constant)/coefx, 0.0};
425 SCIP_Real py[5] = {(coefx*lbx - constant)/coefy, (coefx*ubx - constant)/coefy, lby, uby, 0.0};
426 int j;
427
428 /* the last entry corresponds to the extremum of xy on the line */
429 py[4] = (-constant) / (2.0 * coefy);
430 px[4] = constant / (2.0 * coefx);
431
432 for( j = 0; j < 5; ++j )
433 {
434 if( isPointFeasible(scip, px[j], py[j], lbx, ubx, lby, uby, ineqs, nineqs) )
435 {
436 xs[*npoints] = px[j];
437 ys[*npoints] = py[j];
438 ++(*npoints);
439 }
440 }
441 }
442
443 /* intersection point between two inequalities */
444 for( i = 0; i < nineqs - 1; ++i )
445 {
446 SCIP_Real coefx1 = ineqs[3*i];
447 SCIP_Real coefy1 = ineqs[3*i+1];
448 SCIP_Real constant1 = ineqs[3*i+2];
449 int j;
450
451 for( j = i + 1; j < nineqs; ++j )
452 {
453 SCIP_Real coefx2 = ineqs[3*j];
454 SCIP_Real coefy2 = ineqs[3*j+1];
455 SCIP_Real constant2 = ineqs[3*j+2];
456 SCIP_Real px;
457 SCIP_Real py;
458
459 /* no intersection point -> skip */
460 if( SCIPisZero(scip, coefx2*coefy1 - coefx1 * coefy2) )
461 continue;
462
463 py = (constant2 * coefx1 - constant1 * coefx2)/ (coefx2 * coefy1 - coefx1 * coefy2);
464 px = (coefy1 * py + constant1) / coefx1;
465 assert(SCIPisRelEQ(scip, px, (coefy2 * py + constant2) / coefx2));
466
467 if( isPointFeasible(scip, px, py, lbx, ubx, lby, uby, ineqs, nineqs) )
468 {
469 xs[*npoints] = px;
470 ys[*npoints] = py;
471 ++(*npoints);
472 }
473 }
474 }
475
476 assert(*npoints <= 22);
477
478 /* consider the intersection of the level set with
479 *
480 * 1. the boundary of the box
481 * 2. the linear inequalities
482 *
483 * this adds at most for 4 (level set curves) * 4 (inequalities) * 2 (intersection points) for all linear
484 * inequalities and 4 (level set curves) * 2 (intersection points) with the boundary of the box
485 */
486 if( !levelset )
487 return;
488
489 /* compute intersection of level sets with the boundary */
490 for( i = 0; i < 2; ++i )
491 {
492 SCIP_Real vals[4] = {lbx, ubx, lby, uby};
493 SCIP_Real val;
494 int k;
495
496 /* fix auxiliary variable to its lower or upper bound and consider the coefficient of the product */
497 val = (i == 0) ? exprbounds.inf : exprbounds.sup;
498 val /= SCIPgetCoefExprProduct(expr);
499
500 for( k = 0; k < 4; ++k )
501 {
502 if( !SCIPisZero(scip, vals[k]) )
503 {
504 SCIP_Real res = val / vals[k];
505
506 assert(SCIPisRelGE(scip, SCIPgetCoefExprProduct(expr)*res*vals[k], exprbounds.inf));
507 assert(SCIPisRelLE(scip, SCIPgetCoefExprProduct(expr)*res*vals[k], exprbounds.sup));
508
509 /* fix x to lbx or ubx */
510 if( k < 2 && isPointFeasible(scip, vals[k], res, lbx, ubx, lby, uby, ineqs, nineqs) )
511 {
512 xs[*npoints] = vals[k];
513 ys[*npoints] = res;
514 ++(*npoints);
515 }
516 /* fix y to lby or uby */
517 else if( k >= 2 && isPointFeasible(scip, res, vals[k], lbx, ubx, lby, uby, ineqs, nineqs) )
518 {
519 xs[*npoints] = res;
520 ys[*npoints] = vals[k];
521 ++(*npoints);
522 }
523 }
524 }
525 }
526
527 /* compute intersection points of level sets with the linear inequalities */
528 for( i = 0; i < nineqs; ++i )
529 {
530 SCIP_INTERVAL result;
531 SCIP_Real coefx = ineqs[3*i];
532 SCIP_Real coefy = ineqs[3*i+1];
533 SCIP_Real constant = ineqs[3*i+2];
534 SCIP_INTERVAL sqrcoef;
535 SCIP_INTERVAL lincoef;
536 SCIP_Real px;
537 SCIP_Real py;
538 int k;
539
540 /* solve system of coefx x = coefy y + constant and X = xy which is the same as computing the solutions of
541 *
542 * (coefy / coefx) y^2 + (constant / coefx) y = inf(X) or sup(X)
543 */
544 SCIPintervalSet(&sqrcoef, coefy / coefx);
545 SCIPintervalSet(&lincoef, constant / coefx);
546
547 for( k = 0; k < 2; ++k )
548 {
549 SCIP_INTERVAL rhs;
550 SCIP_INTERVAL ybnds;
551
552 /* set right-hand side */
553 if( k == 0 )
554 SCIPintervalSet(&rhs, exprbounds.inf);
555 else
556 SCIPintervalSet(&rhs, exprbounds.sup);
557
558 SCIPintervalSetBounds(&ybnds, lby, uby);
559 SCIPintervalSolveUnivariateQuadExpression(SCIP_INTERVAL_INFINITY, &result, sqrcoef, lincoef, rhs, ybnds);
560
561 /* interval is empty -> no solution available */
562 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, result) )
563 continue;
564
565 /* compute and check point */
566 py = SCIPintervalGetInf(result);
567 px = (coefy * py + constant) / coefx;
568
569 if( isPointFeasible(scip, px, py, lbx, ubx, lby, uby, ineqs, nineqs) )
570 {
571 xs[*npoints] = px;
572 ys[*npoints] = py;
573 ++(*npoints);
574 }
575
576 /* check for a second solution */
577 if( SCIPintervalGetInf(result) != SCIPintervalGetSup(result) ) /*lint !e777*/
578 {
579 py = SCIPintervalGetSup(result);
580 px = (coefy * py + constant) / coefx;
581
582 if( isPointFeasible(scip, px, py, lbx, ubx, lby, uby, ineqs, nineqs) )
583 {
584 xs[*npoints] = px;
585 ys[*npoints] = py;
586 ++(*npoints);
587 }
588 }
589 }
590 }
591
592 assert(*npoints <= 62);
593 }
594
595 /** computes interval for a bilinear term when using at least one inequality */
596 static
597 SCIP_INTERVAL intevalBilinear(
598 SCIP* scip, /**< SCIP data structure */
599 SCIP_EXPR* expr, /**< product expression */
600 SCIP_Real* underineqs, /**< inequalities for underestimation */
601 int nunderineqs, /**< total number of inequalities for underestimation */
602 SCIP_Real* overineqs, /**< inequalities for overestimation */
603 int noverineqs /**< total number of inequalities for overestimation */
604 )
605 {
606 SCIP_INTERVAL interval = {0., 0.};
607 SCIP_Real xs[22];
608 SCIP_Real ys[22];
609 SCIP_Real inf;
610 SCIP_Real sup;
611 int npoints;
612 int i;
613
614 assert(scip != NULL);
615 assert(expr != NULL);
616 assert(SCIPexprGetNChildren(expr) == 2);
617 assert(noverineqs + nunderineqs <= 4);
618
619 /* no inequalities available -> skip computation */
620 if( noverineqs == 0 && nunderineqs == 0 )
621 {
622 SCIPintervalSetEntire(SCIP_INTERVAL_INFINITY, &interval);
623 return interval;
624 }
625
626 /* x or y has empty interval -> empty */
627 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(SCIPexprGetChildren(expr)[0])) ||
628 SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(SCIPexprGetChildren(expr)[1])) )
629 {
630 SCIPintervalSetEmpty(&interval);
631 return interval;
632 }
633
634 /* compute all feasible points (since we use levelset == FALSE, the value of interval doesn't matter) */
635 getFeasiblePointsBilinear(scip, NULL, expr, interval, underineqs, nunderineqs, overineqs,
636 noverineqs, FALSE, xs, ys, &npoints);
637
638 /* no feasible point left -> return an empty interval */
639 if( npoints == 0 )
640 {
641 SCIPintervalSetEmpty(&interval);
642 return interval;
643 }
644
645 /* compute the minimum and maximum over all computed points */
646 inf = xs[0] * ys[0];
647 sup = inf;
648 for( i = 1; i < npoints; ++i )
649 {
650 inf = MIN(inf, xs[i] * ys[i]);
651 sup = MAX(sup, xs[i] * ys[i]);
652 }
653 assert(inf <= sup);
654
655 /* adjust infinite values */
656 inf = MAX(inf, -SCIP_INTERVAL_INFINITY);
657 sup = MIN(sup, SCIP_INTERVAL_INFINITY);
658
659 /* multiply resulting interval with coefficient of the product expression */
660 SCIPintervalSetBounds(&interval, inf, sup);
661 if( SCIPgetCoefExprProduct(expr) != 1.0 )
662 SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &interval, interval, SCIPgetCoefExprProduct(expr));
663
664 return interval;
665 }
666
667 /** uses inequalities for bilinear terms to get stronger bounds during reverse propagation */
668 static
669 void reversePropBilinear(
670 SCIP* scip, /**< SCIP data structure */
671 SCIP_CONSHDLR* conshdlr, /**< constraint handler */
672 SCIP_EXPR* expr, /**< product expression */
673 SCIP_INTERVAL exprbounds, /**< bounds on product expression */
674 SCIP_Real* underineqs, /**< inequalities for underestimation */
675 int nunderineqs, /**< total number of inequalities for underestimation */
676 SCIP_Real* overineqs, /**< inequalities for overestimation */
677 int noverineqs, /**< total number of inequalities for overestimation */
678 SCIP_INTERVAL* intervalx, /**< buffer to store the new interval for x */
679 SCIP_INTERVAL* intervaly /**< buffer to store the new interval for y */
680 )
681 {
682 SCIP_Real xs[62];
683 SCIP_Real ys[62];
684 SCIP_Real exprinf;
685 SCIP_Real exprsup;
686 SCIP_Bool first = TRUE;
687 int npoints;
688 int i;
689
690 assert(scip != NULL);
691 assert(conshdlr != NULL);
692 assert(expr != NULL);
693 assert(intervalx != NULL);
694 assert(intervaly != NULL);
695 assert(SCIPexprGetNChildren(expr) == 2);
696
697 assert(noverineqs + nunderineqs > 0);
698
699 /* set intervals to be empty */
700 SCIPintervalSetEmpty(intervalx);
701 SCIPintervalSetEmpty(intervaly);
702
703 /* compute feasible points */
704 getFeasiblePointsBilinear(scip, conshdlr, expr, exprbounds, underineqs, nunderineqs, overineqs,
705 noverineqs, TRUE, xs, ys, &npoints);
706
707 /* no feasible points left -> problem is infeasible */
708 if( npoints == 0 )
709 return;
710
711 /* get bounds of the product expression */
712 exprinf = exprbounds.inf;
713 exprsup = exprbounds.sup;
714
715 /* update intervals with the computed points */
716 for( i = 0; i < npoints; ++i )
717 {
718 SCIP_Real val = SCIPgetCoefExprProduct(expr) * xs[i] * ys[i];
719
720 #ifndef NDEBUG
721 {
722 SCIP_Real lbx = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]).inf;
723 SCIP_Real ubx = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]).sup;
724 SCIP_Real lby = SCIPexprGetActivity(SCIPexprGetChildren(expr)[1]).inf;
725 SCIP_Real uby = SCIPexprGetActivity(SCIPexprGetChildren(expr)[1]).sup;
726
727 assert(nunderineqs == 0 || isPointFeasible(scip, xs[i], ys[i], lbx, ubx, lby, uby, underineqs, nunderineqs));
728 assert(noverineqs == 0 || isPointFeasible(scip, xs[i], ys[i], lbx, ubx, lby, uby, overineqs, noverineqs));
729 }
730 #endif
731
732 /* only accept points for which the value of x*y is in the interval of the product expression
733 *
734 * NOTE: in order to consider all relevant points, we are a bit conservative here and relax the interval of
735 * the expression by SCIPfeastol()
736 */
737 if( SCIPisRelGE(scip, val, exprinf - SCIPfeastol(scip)) && SCIPisRelLE(scip, val, exprsup + SCIPfeastol(scip)) )
738 {
739 if( first )
740 {
741 SCIPintervalSet(intervalx, xs[i]);
742 SCIPintervalSet(intervaly, ys[i]);
743 first = FALSE;
744 }
745 else
746 {
747 (*intervalx).inf = MIN((*intervalx).inf, xs[i]);
748 (*intervalx).sup = MAX((*intervalx).sup, xs[i]);
749 (*intervaly).inf = MIN((*intervaly).inf, ys[i]);
750 (*intervaly).sup = MAX((*intervaly).sup, ys[i]);
751 }
752
753 SCIPdebugMsg(scip, "consider points (%g,%g)=%g for reverse propagation\n", xs[i], ys[i], val);
754 }
755 }
756 }
757
758 /** helper function to compute the convex envelope of a bilinear term when two linear inequalities are given; we
759 * use the same notation and formulas as in Locatelli 2016
760 */
761 static
762 void computeBilinEnvelope2(
763 SCIP* scip, /**< SCIP data structure */
764 SCIP_Real x, /**< reference point for x */
765 SCIP_Real y, /**< reference point for y */
766 SCIP_Real mi, /**< coefficient of x in the first linear inequality */
767 SCIP_Real qi, /**< constant in the first linear inequality */
768 SCIP_Real mj, /**< coefficient of x in the second linear inequality */
769 SCIP_Real qj, /**< constant in the second linear inequality */
770 SCIP_Real* RESTRICT xi, /**< buffer to store x coordinate of the first point */
771 SCIP_Real* RESTRICT yi, /**< buffer to store y coordinate of the first point */
772 SCIP_Real* RESTRICT xj, /**< buffer to store x coordinate of the second point */
773 SCIP_Real* RESTRICT yj, /**< buffer to store y coordinate of the second point */
774 SCIP_Real* RESTRICT xcoef, /**< buffer to store the x coefficient of the envelope */
775 SCIP_Real* RESTRICT ycoef, /**< buffer to store the y coefficient of the envelope */
776 SCIP_Real* RESTRICT constant /**< buffer to store the constant of the envelope */
777 )
778 {
779 SCIP_Real QUAD(xiq);
780 SCIP_Real QUAD(yiq);
781 SCIP_Real QUAD(xjq);
782 SCIP_Real QUAD(yjq);
783 SCIP_Real QUAD(xcoefq);
784 SCIP_Real QUAD(ycoefq);
785 SCIP_Real QUAD(constantq);
786 SCIP_Real QUAD(tmpq);
787
788 assert(xi != NULL);
789 assert(yi != NULL);
790 assert(xj != NULL);
791 assert(yj != NULL);
792 assert(xcoef != NULL);
793 assert(ycoef != NULL);
794 assert(constant != NULL);
795
796 if( SCIPisEQ(scip, mi, mj) )
797 {
798 /* xi = (x + mi * y - qi) / (2.0*mi) */
799 SCIPquadprecProdDD(xiq, mi, y);
800 SCIPquadprecSumQD(xiq, xiq, x);
801 SCIPquadprecSumQD(xiq, xiq, -qi);
802 SCIPquadprecDivQD(xiq, xiq, 2.0 * mi);
803 assert(EPSEQ((x + mi * y - qi) / (2.0*mi), QUAD_TO_DBL(xiq), 1e-3));
804
805 /* yi = mi*(*xi) + qi */
806 SCIPquadprecProdQD(yiq, xiq, mi);
807 SCIPquadprecSumQD(yiq, yiq, qi);
808 assert(EPSEQ(mi*QUAD_TO_DBL(xiq) + qi, QUAD_TO_DBL(yiq), 1e-3));
809
810 /* xj = (*xi) + (qi - qj)/ (2.0*mi) */
811 SCIPquadprecSumDD(xjq, qi, -qj);
812 SCIPquadprecDivQD(xjq, xjq, 2.0 * mi);
813 SCIPquadprecSumQQ(xjq, xjq, xiq);
814 assert(EPSEQ(QUAD_TO_DBL(xiq) + (qi - qj)/ (2.0*mi), QUAD_TO_DBL(xjq), 1e-3));
815
816 /* yj = mj * (*xj) + qj */
817 SCIPquadprecProdQD(yjq, xjq, mj);
818 SCIPquadprecSumQD(yjq, yjq, qj);
819 assert(EPSEQ(mj * QUAD_TO_DBL(xjq) + qj, QUAD_TO_DBL(yjq), 1e-3));
820
821 /* ycoef = (*xi) + (qi - qj) / (4.0*mi) note that this is wrong in Locatelli 2016 */
822 SCIPquadprecSumDD(ycoefq, qi, -qj);
823 SCIPquadprecDivQD(ycoefq, ycoefq, 4.0 * mi);
824 SCIPquadprecSumQQ(ycoefq, ycoefq, xiq);
825 assert(EPSEQ(QUAD_TO_DBL(xiq) + (qi - qj) / (4.0*mi), QUAD_TO_DBL(ycoefq), 1e-3));
826
827 /* xcoef = 2.0*mi*(*xi) - mi * (*ycoef) + qi */
828 SCIPquadprecProdQD(xcoefq, xiq, 2.0 * mi);
829 SCIPquadprecProdQD(tmpq, ycoefq, -mi);
830 SCIPquadprecSumQQ(xcoefq, xcoefq, tmpq);
831 SCIPquadprecSumQD(xcoefq, xcoefq, qi);
832 assert(EPSEQ(2.0*mi*QUAD_TO_DBL(xiq) - mi * QUAD_TO_DBL(ycoefq) + qi, QUAD_TO_DBL(xcoefq), 1e-3));
833
834 /* constant = -mj*SQR(*xj) - (*ycoef) * qj */
835 SCIPquadprecSquareQ(constantq, xjq);
836 SCIPquadprecProdQD(constantq, constantq, -mj);
837 SCIPquadprecProdQD(tmpq, ycoefq, -qj);
838 SCIPquadprecSumQQ(constantq, constantq, tmpq);
839 /* assert(EPSEQ(-mj*SQR(QUAD_TO_DBL(xjq)) - QUAD_TO_DBL(ycoefq) * qj, QUAD_TO_DBL(constantq), 1e-3)); */
840
841 *xi = QUAD_TO_DBL(xiq);
842 *yi = QUAD_TO_DBL(yiq);
843 *xj = QUAD_TO_DBL(xjq);
844 *yj = QUAD_TO_DBL(yjq);
845 *ycoef = QUAD_TO_DBL(ycoefq);
846 *xcoef = QUAD_TO_DBL(xcoefq);
847 *constant = QUAD_TO_DBL(constantq);
848 }
849 else if( mi > 0.0 )
850 {
851 assert(mj > 0.0);
852
853 /* xi = (y + SQRT(mi*mj)*x - qi) / (REALABS(mi) + SQRT(mi*mj)) */
854 SCIPquadprecProdDD(xiq, mi, mj);
855 SCIPquadprecSqrtQ(xiq, xiq);
856 SCIPquadprecProdQD(xiq, xiq, x);
857 SCIPquadprecSumQD(xiq, xiq, y);
858 SCIPquadprecSumQD(xiq, xiq, -qi); /* (y + SQRT(mi*mj)*x - qi) */
859 SCIPquadprecProdDD(tmpq, mi, mj);
860 SCIPquadprecSqrtQ(tmpq, tmpq);
861 SCIPquadprecSumQD(tmpq, tmpq, REALABS(mi)); /* REALABS(mi) + SQRT(mi*mj) */
862 SCIPquadprecDivQQ(xiq, xiq, tmpq);
863 assert(EPSEQ((y + SQRT(mi*mj)*x - qi) / (REALABS(mi) + SQRT(mi*mj)), QUAD_TO_DBL(xiq), 1e-3));
864
865 /* yi = mi*(*xi) + qi */
866 SCIPquadprecProdQD(yiq, xiq, mi);
867 SCIPquadprecSumQD(yiq, yiq, qi);
868 assert(EPSEQ(mi*(QUAD_TO_DBL(xiq)) + qi, QUAD_TO_DBL(yiq), 1e-3));
869
870 /* xj = (y + SQRT(mi*mj)*x - qj) / (REALABS(mj) + SQRT(mi*mj)) */
871 SCIPquadprecProdDD(xjq, mi, mj);
872 SCIPquadprecSqrtQ(xjq, xjq);
873 SCIPquadprecProdQD(xjq, xjq, x);
874 SCIPquadprecSumQD(xjq, xjq, y);
875 SCIPquadprecSumQD(xjq, xjq, -qj); /* (y + SQRT(mi*mj)*x - qj) */
876 SCIPquadprecProdDD(tmpq, mi, mj);
877 SCIPquadprecSqrtQ(tmpq, tmpq);
878 SCIPquadprecSumQD(tmpq, tmpq, REALABS(mj)); /* REALABS(mj) + SQRT(mi*mj) */
879 SCIPquadprecDivQQ(xjq, xjq, tmpq);
880 assert(EPSEQ((y + SQRT(mi*mj)*x - qj) / (REALABS(mj) + SQRT(mi*mj)), QUAD_TO_DBL(xjq), 1e-3));
881
882 /* yj = mj*(*xj) + qj */
883 SCIPquadprecProdQD(yjq, xjq, mj);
884 SCIPquadprecSumQD(yjq, yjq, qj);
885 assert(EPSEQ(mj*QUAD_TO_DBL(xjq) + qj, QUAD_TO_DBL(yjq), 1e-3));
886
887 /* ycoef = (2.0*mj*(*xj) + qj - 2.0*mi*(*xi) - qi) / (mj - mi) */
888 SCIPquadprecProdQD(ycoefq, xjq, 2.0 * mj);
889 SCIPquadprecSumQD(ycoefq, ycoefq, qj);
890 SCIPquadprecProdQD(tmpq, xiq, -2.0 * mi);
891 SCIPquadprecSumQQ(ycoefq, ycoefq, tmpq);
892 SCIPquadprecSumQD(ycoefq, ycoefq, -qi);
893 SCIPquadprecSumDD(tmpq, mj, -mi);
894 SCIPquadprecDivQQ(ycoefq, ycoefq, tmpq);
895 assert(EPSEQ((2.0*mj*QUAD_TO_DBL(xjq) + qj - 2.0*mi*QUAD_TO_DBL(xiq) - qi) / (mj - mi), QUAD_TO_DBL(ycoefq), 1e-3));
896
897 /* xcoef = 2.0*mj*(*xj) + qj - mj*(*ycoef) */
898 SCIPquadprecProdQD(xcoefq, xjq, 2.0 * mj);
899 SCIPquadprecSumQD(xcoefq, xcoefq, qj);
900 SCIPquadprecProdQD(tmpq, ycoefq, -mj);
901 SCIPquadprecSumQQ(xcoefq, xcoefq, tmpq);
902 assert(EPSEQ(2.0*mj*QUAD_TO_DBL(xjq) + qj - mj*QUAD_TO_DBL(ycoefq), QUAD_TO_DBL(xcoefq), 1e-3));
903
904 /* constant = -mj*SQR(*xj) - (*ycoef) * qj */
905 SCIPquadprecSquareQ(constantq, xjq);
906 SCIPquadprecProdQD(constantq, constantq, -mj);
907 SCIPquadprecProdQD(tmpq, ycoefq, -qj);
908 SCIPquadprecSumQQ(constantq, constantq, tmpq);
909 /* assert(EPSEQ(-mj*SQR(QUAD_TO_DBL(xjq)) - QUAD_TO_DBL(ycoefq) * qj, QUAD_TO_DBL(constantq), 1e-3)); */
910
911 *xi = QUAD_TO_DBL(xiq);
912 *yi = QUAD_TO_DBL(yiq);
913 *xj = QUAD_TO_DBL(xjq);
914 *yj = QUAD_TO_DBL(yjq);
915 *ycoef = QUAD_TO_DBL(ycoefq);
916 *xcoef = QUAD_TO_DBL(xcoefq);
917 *constant = QUAD_TO_DBL(constantq);
918 }
919 else
920 {
921 assert(mi < 0.0 && mj < 0.0);
922
923 /* apply variable transformation x = -x in case for overestimation */
924 computeBilinEnvelope2(scip, -x, y, -mi, qi, -mj, qj, xi, yi, xj, yj, xcoef, ycoef, constant);
925
926 /* revert transformation; multiply cut by -1 and change -x by x */
927 *xi = -(*xi);
928 *xj = -(*xj);
929 *ycoef = -(*ycoef);
930 *constant = -(*constant);
931 }
932 }
933
934 /** output method of statistics table to output file stream 'file' */
935 static
936 SCIP_DECL_TABLEOUTPUT(tableOutputBilinear)
937 { /*lint --e{715}*/
938 SCIP_NLHDLR* nlhdlr;
939 SCIP_NLHDLRDATA* nlhdlrdata;
940 SCIP_CONSHDLR* conshdlr;
941 SCIP_HASHMAP* hashmap;
942 SCIP_EXPRITER* it;
943 int resfound = 0;
944 int restotal = 0;
945 int c;
946
947 conshdlr = SCIPfindConshdlr(scip, "nonlinear");
948 assert(conshdlr != NULL);
949 nlhdlr = SCIPfindNlhdlrNonlinear(conshdlr, NLHDLR_NAME);
950 assert(nlhdlr != NULL);
951 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
952 assert(nlhdlrdata != NULL);
953
954 /* allocate memory */
955 SCIP_CALL( SCIPhashmapCreate(&hashmap, SCIPblkmem(scip), nlhdlrdata->nexprs) );
956 SCIP_CALL( SCIPcreateExpriter(scip, &it) );
957
958 for( c = 0; c < nlhdlrdata->nexprs; ++c )
959 {
960 assert(!SCIPhashmapExists(hashmap, nlhdlrdata->exprs[c]));
961 SCIP_CALL( SCIPhashmapInsertInt(hashmap, nlhdlrdata->exprs[c], 0) );
962 }
963
964 /* count in how many constraints each expression is contained */
965 for( c = 0; c < SCIPconshdlrGetNConss(conshdlr); ++c )
966 {
967 SCIP_CONS* cons = SCIPconshdlrGetConss(conshdlr)[c];
968 SCIP_EXPR* expr;
969
970 SCIP_CALL( SCIPexpriterInit(it, SCIPgetExprNonlinear(cons), SCIP_EXPRITER_DFS, FALSE) );
971
972 for( expr = SCIPexpriterGetCurrent(it); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) ) /*lint !e441*/
973 {
974 if( SCIPhashmapExists(hashmap, expr) )
975 {
976 int nuses = SCIPhashmapGetImageInt(hashmap, expr);
977 SCIP_CALL( SCIPhashmapSetImageInt(hashmap, expr, nuses + 1) );
978 }
979 }
980 }
981
982 /* compute success ratio */
983 for( c = 0; c < nlhdlrdata->nexprs; ++c )
984 {
985 SCIP_NLHDLREXPRDATA* nlhdlrexprdata;
986 int nuses;
987
988 nuses = SCIPhashmapGetImageInt(hashmap, nlhdlrdata->exprs[c]);
989 assert(nuses > 0);
990
991 nlhdlrexprdata = SCIPgetNlhdlrExprDataNonlinear(nlhdlr, nlhdlrdata->exprs[c]);
992 assert(nlhdlrexprdata != NULL);
993
994 if( nlhdlrexprdata->nunderineqs > 0 || nlhdlrexprdata->noverineqs > 0 )
995 resfound += nuses;
996 restotal += nuses;
997 }
998
999 /* print statistics */
1000 SCIPinfoMessage(scip, file, "Bilinear Nlhdlr : %10s %10s\n", "#found", "#total");
1001 SCIPinfoMessage(scip, file, " %-17s:", "");
1002 SCIPinfoMessage(scip, file, " %10d", resfound);
1003 SCIPinfoMessage(scip, file, " %10d", restotal);
1004 SCIPinfoMessage(scip, file, "\n");
1005
1006 /* free memory */
1007 SCIPfreeExpriter(&it);
1008 SCIPhashmapFree(&hashmap);
1009
1010 return SCIP_OKAY;
1011 }
1012
1013
1014 /*
1015 * Callback methods of nonlinear handler
1016 */
1017
1018 /** nonlinear handler copy callback */
1019 static
1020 SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrBilinear)
1021 { /*lint --e{715}*/
1022 assert(targetscip != NULL);
1023 assert(sourcenlhdlr != NULL);
1024 assert(strcmp(SCIPnlhdlrGetName(sourcenlhdlr), NLHDLR_NAME) == 0);
1025
1026 SCIP_CALL( SCIPincludeNlhdlrBilinear(targetscip) );
1027
1028 return SCIP_OKAY;
1029 }
1030
1031 /** callback to free data of handler */
1032 static
1033 SCIP_DECL_NLHDLRFREEHDLRDATA(nlhdlrFreehdlrdataBilinear)
1034 { /*lint --e{715}*/
1035 assert(nlhdlrdata != NULL);
1036 assert((*nlhdlrdata)->nexprs == 0);
1037
1038 if( (*nlhdlrdata)->exprmap != NULL )
1039 {
1040 assert(SCIPhashmapGetNElements((*nlhdlrdata)->exprmap) == 0);
1041 SCIPhashmapFree(&(*nlhdlrdata)->exprmap);
1042 }
1043
1044 SCIPfreeBlockMemoryArrayNull(scip, &(*nlhdlrdata)->exprs, (*nlhdlrdata)->exprsize);
1045 SCIPfreeBlockMemory(scip, nlhdlrdata);
1046
1047 return SCIP_OKAY;
1048 }
1049
1050 /** callback to free expression specific data */
1051 static
1052 SCIP_DECL_NLHDLRFREEEXPRDATA(nlhdlrFreeExprDataBilinear)
1053 { /*lint --e{715}*/
1054 SCIP_NLHDLRDATA* nlhdlrdata;
1055 int pos;
1056
1057 assert(expr != NULL);
1058
1059 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1060 assert(nlhdlrdata != NULL);
1061 assert(nlhdlrdata->nexprs > 0);
1062 assert(nlhdlrdata->exprs != NULL);
1063 assert(nlhdlrdata->exprmap != NULL);
1064 assert(SCIPhashmapExists(nlhdlrdata->exprmap, (void*)expr));
1065
1066 pos = SCIPhashmapGetImageInt(nlhdlrdata->exprmap, (void*)expr);
1067 assert(pos >= 0 && pos < nlhdlrdata->nexprs);
1068 assert(nlhdlrdata->exprs[pos] == expr);
1069
1070 /* move the last expression to the free position */
1071 if( nlhdlrdata->nexprs > 0 && pos != nlhdlrdata->nexprs - 1 )
1072 {
1073 SCIP_EXPR* lastexpr = nlhdlrdata->exprs[nlhdlrdata->nexprs - 1];
1074 assert(expr != lastexpr);
1075 assert(SCIPhashmapExists(nlhdlrdata->exprmap, (void*)lastexpr));
1076
1077 nlhdlrdata->exprs[pos] = lastexpr;
1078 nlhdlrdata->exprs[nlhdlrdata->nexprs - 1] = NULL;
1079 SCIP_CALL( SCIPhashmapSetImageInt(nlhdlrdata->exprmap, (void*)lastexpr, pos) );
1080 }
1081
1082 /* remove expression from the nonlinear handler data */
1083 SCIP_CALL( SCIPhashmapRemove(nlhdlrdata->exprmap, (void*)expr) );
1084 SCIP_CALL( SCIPreleaseExpr(scip, &expr) );
1085 --nlhdlrdata->nexprs;
1086
1087 /* free nonlinear handler expression data */
1088 SCIPfreeBlockMemoryNull(scip, nlhdlrexprdata);
1089
1090 return SCIP_OKAY;
1091 }
1092
1093 /** callback to be called in initialization */
1094 #define nlhdlrInitBilinear NULL
1095
1096 /** callback to be called in deinitialization */
1097 static
1098 SCIP_DECL_NLHDLREXIT(nlhdlrExitBilinear)
1099 { /*lint --e{715}*/
1100 assert(SCIPnlhdlrGetData(nlhdlr) != NULL);
1101 assert(SCIPnlhdlrGetData(nlhdlr)->nexprs == 0);
1102
1103 return SCIP_OKAY;
1104 }
1105
1106 /** callback to detect structure in expression tree */
1107 static
1108 SCIP_DECL_NLHDLRDETECT(nlhdlrDetectBilinear)
1109 { /*lint --e{715}*/
1110 SCIP_NLHDLRDATA* nlhdlrdata;
1111
1112 assert(expr != NULL);
1113 assert(participating != NULL);
1114
1115 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1116 assert(nlhdlrdata);
1117
1118 /* only during solving will we have the extra inequalities that we rely on so much here */
1119 if( SCIPgetStage(scip) < SCIP_STAGE_INITSOLVE )
1120 return SCIP_OKAY;
1121
1122 /* check for product expressions with two children */
1123 if( SCIPisExprProduct(scip, expr) && SCIPexprGetNChildren(expr) == 2
1124 && (nlhdlrdata->exprmap == NULL || !SCIPhashmapExists(nlhdlrdata->exprmap, (void*)expr)) )
1125 {
1126 SCIP_EXPR** children;
1127 SCIP_Bool valid;
1128 int c;
1129
1130 children = SCIPexprGetChildren(expr);
1131 assert(children != NULL);
1132
1133 /* detection is only successful if both children will have auxiliary variable or are variables
1134 * that are not binary variables */
1135 valid = TRUE;
1136 for( c = 0; c < 2; ++c )
1137 {
1138 assert(children[c] != NULL);
1139 if( SCIPgetExprNAuxvarUsesNonlinear(children[c]) == 0 &&
1140 (!SCIPisExprVar(scip, children[c]) || SCIPvarIsBinary(SCIPgetVarExprVar(children[c]))) )
1141 {
1142 valid = FALSE;
1143 break;
1144 }
1145 }
1146
1147 if( valid )
1148 {
1149 /* create expression data for the nonlinear handler */
1150 SCIP_CALL( SCIPallocClearBlockMemory(scip, nlhdlrexprdata) );
1151 (*nlhdlrexprdata)->lastnodeid = -1;
1152
1153 /* ensure that there is enough memory to store the detected expression */
1154 if( nlhdlrdata->exprsize < nlhdlrdata->nexprs + 1 )
1155 {
1156 int newsize = SCIPcalcMemGrowSize(scip, nlhdlrdata->nexprs + 1);
1157 assert(newsize > nlhdlrdata->exprsize);
1158
1159 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &nlhdlrdata->exprs, nlhdlrdata->exprsize, newsize) );
1160 nlhdlrdata->exprsize = newsize;
1161 }
1162
1163 /* create expression map, if not done so far */
1164 if( nlhdlrdata->exprmap == NULL )
1165 {
1166 SCIP_CALL( SCIPhashmapCreate(&nlhdlrdata->exprmap, SCIPblkmem(scip), SCIPgetNVars(scip)) );
1167 }
1168
1169 #ifndef NDEBUG
1170 {
1171 int i;
1172
1173 for( i = 0; i < nlhdlrdata->nexprs; ++i )
1174 assert(nlhdlrdata->exprs[i] != expr);
1175 }
1176 #endif
1177
1178 /* add expression to nlhdlrdata and capture it */
1179 nlhdlrdata->exprs[nlhdlrdata->nexprs] = expr;
1180 SCIPcaptureExpr(expr);
1181 SCIP_CALL( SCIPhashmapInsertInt(nlhdlrdata->exprmap, (void*)expr, nlhdlrdata->nexprs) );
1182 ++nlhdlrdata->nexprs;
1183
1184 /* tell children that we will use their auxvar and use its activity for both estimate and domain propagation */
1185 SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, children[0], TRUE, nlhdlrdata->useinteval
1186 || nlhdlrdata->usereverseprop, TRUE, TRUE) );
1187 SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, children[1], TRUE, nlhdlrdata->useinteval
1188 || nlhdlrdata->usereverseprop, TRUE, TRUE) );
1189 }
1190 }
1191
1192 if( *nlhdlrexprdata != NULL )
1193 {
1194 /* we want to join separation and domain propagation, if not disabled by parameter */
1195 *participating = SCIP_NLHDLR_METHOD_SEPABOTH;
1196 if( nlhdlrdata->useinteval || nlhdlrdata->usereverseprop )
1197 *participating |= SCIP_NLHDLR_METHOD_ACTIVITY;
1198 }
1199
1200 #ifdef SCIP_DEBUG
1201 if( *participating )
1202 {
1203 SCIPdebugMsg(scip, "detected expr ");
1204 SCIPprintExpr(scip, expr, NULL);
1205 SCIPinfoMessage(scip, NULL, " participating: %d\n", *participating);
1206 }
1207 #endif
1208
1209 return SCIP_OKAY;
1210 }
1211
1212 /** auxiliary evaluation callback of nonlinear handler */
1213 static
1214 SCIP_DECL_NLHDLREVALAUX(nlhdlrEvalauxBilinear)
1215 { /*lint --e{715}*/
1216 SCIP_VAR* var1;
1217 SCIP_VAR* var2;
1218 SCIP_Real coef;
1219
1220 assert(SCIPisExprProduct(scip, expr));
1221 assert(SCIPexprGetNChildren(expr) == 2);
1222
1223 var1 = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[0]);
1224 assert(var1 != NULL);
1225 var2 = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[1]);
1226 assert(var2 != NULL);
1227 coef = SCIPgetCoefExprProduct(expr);
1228
1229 *auxvalue = coef * SCIPgetSolVal(scip, sol, var1) * SCIPgetSolVal(scip, sol, var2);
1230
1231 return SCIP_OKAY;
1232 }
1233
1234 /** separation initialization method of a nonlinear handler (called during CONSINITLP) */
1235 #define nlhdlrInitSepaBilinear NULL
1236
1237 /** separation deinitialization method of a nonlinear handler (called during CONSEXITSOL) */
1238 #define nlhdlrExitSepaBilinear NULL
1239
1240 /** nonlinear handler separation callback */
1241 #define nlhdlrEnfoBilinear NULL
1242
1243 /** nonlinear handler under/overestimation callback */
1244 static
1245 SCIP_DECL_NLHDLRESTIMATE(nlhdlrEstimateBilinear)
1246 { /*lint --e{715}*/
1247 SCIP_NLHDLRDATA* nlhdlrdata;
1248 SCIP_VAR* x;
1249 SCIP_VAR* y;
1250 SCIP_VAR* auxvar;
1251 SCIP_Real lincoefx = 0.0;
1252 SCIP_Real lincoefy = 0.0;
1253 SCIP_Real linconstant = 0.0;
1254 SCIP_Real refpointx;
1255 SCIP_Real refpointy;
1256 SCIP_Real violation;
1257 SCIP_Longint nodeid;
1258 SCIP_Bool mccsuccess = TRUE;
1259 SCIP_ROWPREP* rowprep;
1260
1261 assert(rowpreps != NULL);
1262
1263 *success = FALSE;
1264 *addedbranchscores = FALSE;
1265
1266 /* check whether an inequality is available */
1267 if( nlhdlrexprdata->noverineqs == 0 && nlhdlrexprdata->nunderineqs == 0 )
1268 return SCIP_OKAY;
1269
1270 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1271 assert(nlhdlrdata != NULL);
1272
1273 nodeid = SCIPnodeGetNumber(SCIPgetCurrentNode(scip));
1274
1275 /* update last node */
1276 if( nlhdlrexprdata->lastnodeid != nodeid )
1277 {
1278 nlhdlrexprdata->lastnodeid = nodeid;
1279 nlhdlrexprdata->nseparoundslastnode = 0;
1280 }
1281
1282 /* update separation round */
1283 ++nlhdlrexprdata->nseparoundslastnode;
1284
1285 /* check working limits */
1286 if( (SCIPgetDepth(scip) == 0 && nlhdlrexprdata->nseparoundslastnode > nlhdlrdata->maxseparoundsroot)
1287 || (SCIPgetDepth(scip) > 0 && nlhdlrexprdata->nseparoundslastnode > nlhdlrdata->maxseparounds)
1288 || SCIPgetDepth(scip) > nlhdlrdata->maxsepadepth )
1289 return SCIP_OKAY;
1290
1291 /* collect variables */
1292 x = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[0]);
1293 assert(x != NULL);
1294 y = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[1]);
1295 assert(y != NULL);
1296 auxvar = SCIPgetExprAuxVarNonlinear(expr);
1297 assert(auxvar != NULL);
1298
1299 /* get and adjust the reference points */
1300 refpointx = MIN(MAX(SCIPgetSolVal(scip, sol, x), SCIPvarGetLbLocal(x)),SCIPvarGetUbLocal(x)); /*lint !e666*/
1301 refpointy = MIN(MAX(SCIPgetSolVal(scip, sol, y), SCIPvarGetLbLocal(y)),SCIPvarGetUbLocal(y)); /*lint !e666*/
1302 assert(SCIPisLE(scip, refpointx, SCIPvarGetUbLocal(x)) && SCIPisGE(scip, refpointx, SCIPvarGetLbLocal(x)));
1303 assert(SCIPisLE(scip, refpointy, SCIPvarGetUbLocal(y)) && SCIPisGE(scip, refpointy, SCIPvarGetLbLocal(y)));
1304
1305 /* use McCormick inequalities to decide whether we want to separate or not */
1306 SCIPaddBilinMcCormick(scip, SCIPgetCoefExprProduct(expr), SCIPvarGetLbLocal(x), SCIPvarGetUbLocal(x), refpointx,
1307 SCIPvarGetLbLocal(y), SCIPvarGetUbLocal(y), refpointy, overestimate, &lincoefx, &lincoefy, &linconstant,
1308 &mccsuccess);
1309
1310 /* too large values in McCormick inequalities -> skip */
1311 if( !mccsuccess )
1312 return SCIP_OKAY;
1313
1314 /* compute violation for the McCormick relaxation */
1315 violation = lincoefx * refpointx + lincoefy * refpointy + linconstant - SCIPgetSolVal(scip, sol, auxvar);
1316 if( overestimate )
1317 violation = -violation;
1318
1319 /* only use a tighter relaxations if McCormick does not separate the reference point */
1320 if( SCIPisFeasLE(scip, violation, 0.0) && useBilinIneqs(scip, x, y, refpointx, refpointy) )
1321 {
1322 SCIP_Bool useoverestineq = SCIPgetCoefExprProduct(expr) > 0.0 ? overestimate : !overestimate;
1323 SCIP_Real mccormickval = lincoefx * refpointx + lincoefy * refpointy + linconstant;
1324 SCIP_Real* ineqs;
1325 SCIP_Real bestval;
1326 int nineqs;
1327
1328 /* McCormick relaxation is too weak */
1329 bestval = mccormickval;
1330
1331 /* get the inequalities that might lead to a tighter relaxation */
1332 if( useoverestineq )
1333 {
1334 ineqs = nlhdlrexprdata->overineqs;
1335 nineqs = nlhdlrexprdata->noverineqs;
1336 }
1337 else
1338 {
1339 ineqs = nlhdlrexprdata->underineqs;
1340 nineqs = nlhdlrexprdata->nunderineqs;
1341 }
1342
1343 /* use linear inequalities to update relaxation */
1344 updateBilinearRelaxation(scip, x, y, SCIPgetCoefExprProduct(expr),
1345 overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT,
1346 refpointx, refpointy, ineqs, nineqs, mccormickval,
1347 &lincoefx, &lincoefy, &linconstant, &bestval,
1348 success);
1349
1350 #ifndef NDEBUG
1351 /* check whether cut is really valid */
1352 if( *success )
1353 {
1354 assert(!overestimate || SCIPisLE(scip, bestval, mccormickval));
1355 assert(overestimate || SCIPisGE(scip, bestval, mccormickval));
1356 }
1357 #endif
1358 }
1359
1360 if( *success )
1361 {
1362 SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) );
1363 SCIProwprepAddConstant(rowprep, linconstant);
1364 SCIP_CALL( SCIPensureRowprepSize(scip, rowprep, 2) );
1365 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, x, lincoefx) );
1366 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, y, lincoefy) );
1367 SCIP_CALL( SCIPsetPtrarrayVal(scip, rowpreps, 0, rowprep) );
1368 }
1369
1370 return SCIP_OKAY;
1371 }
1372
1373 /** nonlinear handler interval evaluation callback */
1374 static
1375 SCIP_DECL_NLHDLRINTEVAL(nlhdlrIntevalBilinear)
1376 { /*lint --e{715}*/
1377 SCIP_NLHDLRDATA* nlhdlrdata;
1378 assert(nlhdlrexprdata != NULL);
1379
1380 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1381 assert(nlhdlrdata != NULL);
1382
1383 if( nlhdlrdata->useinteval && nlhdlrexprdata->nunderineqs + nlhdlrexprdata->noverineqs > 0 )
1384 {
1385 SCIP_INTERVAL tmp = intevalBilinear(scip, expr, nlhdlrexprdata->underineqs, nlhdlrexprdata->nunderineqs,
1386 nlhdlrexprdata->overineqs, nlhdlrexprdata->noverineqs);
1387
1388 /* intersect intervals if we have learned a tighter interval */
1389 if( SCIPisGT(scip, tmp.inf, (*interval).inf) || SCIPisLT(scip, tmp.sup, (*interval).sup) )
1390 SCIPintervalIntersect(interval, *interval, tmp);
1391 }
1392
1393 return SCIP_OKAY;
1394 }
1395
1396 /** nonlinear handler callback for reverse propagation */
1397 static
1398 SCIP_DECL_NLHDLRREVERSEPROP(nlhdlrReversepropBilinear)
1399 { /*lint --e{715}*/
1400 SCIP_NLHDLRDATA* nlhdlrdata;
1401
1402 assert(nlhdlrexprdata != NULL);
1403
1404 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1405 assert(nlhdlrdata != NULL);
1406
1407 if( nlhdlrdata->usereverseprop && nlhdlrexprdata->nunderineqs + nlhdlrexprdata->noverineqs > 0 )
1408 {
1409 SCIP_EXPR* childx;
1410 SCIP_EXPR* childy;
1411 SCIP_INTERVAL intervalx;
1412 SCIP_INTERVAL intervaly;
1413
1414 assert(SCIPexprGetNChildren(expr) == 2);
1415 childx = SCIPexprGetChildren(expr)[0];
1416 childy = SCIPexprGetChildren(expr)[1];
1417 assert(childx != NULL && childy != NULL);
1418
1419 SCIPintervalSetEntire(SCIP_INTERVAL_INFINITY, &intervalx);
1420 SCIPintervalSetEntire(SCIP_INTERVAL_INFINITY, &intervaly);
1421
1422 /* compute bounds on x and y */
1423 reversePropBilinear(scip, conshdlr, expr, bounds, nlhdlrexprdata->underineqs, nlhdlrexprdata->nunderineqs,
1424 nlhdlrexprdata->overineqs, nlhdlrexprdata->noverineqs, &intervalx, &intervaly);
1425
1426 /* tighten bounds of x */
1427 SCIPdebugMsg(scip, "try to tighten bounds of x: [%g,%g] -> [%g,%g]\n",
1428 SCIPgetExprBoundsNonlinear(scip, childx).inf, SCIPgetExprBoundsNonlinear(scip, childx).sup,
1429 intervalx.inf, intervalx.sup);
1430
1431 SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, SCIPexprGetChildren(expr)[0], intervalx, infeasible,
1432 nreductions) );
1433
1434 if( !(*infeasible) )
1435 {
1436 /* tighten bounds of y */
1437 SCIPdebugMsg(scip, "try to tighten bounds of y: [%g,%g] -> [%g,%g]\n",
1438 SCIPgetExprBoundsNonlinear(scip, childy).inf, SCIPgetExprBoundsNonlinear(scip, childy).sup,
1439 intervaly.inf, intervaly.sup);
1440 SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, SCIPexprGetChildren(expr)[1], intervaly,
1441 infeasible, nreductions) );
1442 }
1443 }
1444
1445 return SCIP_OKAY;
1446 }
1447
1448 /*
1449 * nonlinear handler specific interface methods
1450 */
1451
1452 /** includes bilinear nonlinear handler in nonlinear constraint handler */
1453 SCIP_RETCODE SCIPincludeNlhdlrBilinear(
1454 SCIP* scip /**< SCIP data structure */
1455 )
1456 {
1457 SCIP_NLHDLRDATA* nlhdlrdata;
1458 SCIP_NLHDLR* nlhdlr;
1459
1460 assert(scip != NULL);
1461
1462 /* create nonlinear handler specific data */
|
(1) Event cond_false: |
Condition "(nlhdlrdata = BMSallocBlockMemory_call(SCIPblkmem(scip), 48UL /* sizeof (*nlhdlrdata) */, "/sapmnt/home1/d029903/my_work/SCIPSoPlex_coverity/scipoptsuite-v800-rc07/scipoptsuite-8.0.0/scip/src/scip/nlhdlr_bilinear.c", 1463)) == NULL", taking false branch. |
|
(2) Event cond_false: |
Condition "(_restat_ = (((nlhdlrdata = BMSallocBlockMemory_call(SCIPblkmem(scip), 48UL /* sizeof (*nlhdlrdata) */, "/sapmnt/home1/d029903/my_work/SCIPSoPlex_coverity/scipoptsuite-v800-rc07/scipoptsuite-8.0.0/scip/src/scip/nlhdlr_bilinear.c", 1463)) == NULL) ? SCIP_NOMEMORY : SCIP_OKAY)) != SCIP_OKAY", taking false branch. |
|
(3) Event if_end: |
End of if statement. |
1463 SCIP_CALL( SCIPallocBlockMemory(scip, &nlhdlrdata) );
1464 BMSclearMemory(nlhdlrdata);
1465
|
(4) Event alloc_arg: |
"SCIPincludeNlhdlrNonlinear" allocates memory that is stored into "nlhdlr". [details] |
|
(5) Event cond_true: |
Condition "(_restat_ = SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, "bilinear", "bilinear handler for expressions", -10, -10, nlhdlrDetectBilinear, nlhdlrEvalauxBilinear, nlhdlrdata)) != SCIP_OKAY", taking true branch. |
|
(6) Event leaked_storage: |
Variable "nlhdlr" going out of scope leaks the storage it points to. |
1466 SCIP_CALL( SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, NLHDLR_NAME, NLHDLR_DESC, NLHDLR_DETECTPRIORITY,
1467 NLHDLR_ENFOPRIORITY, nlhdlrDetectBilinear, nlhdlrEvalauxBilinear, nlhdlrdata) );
1468 assert(nlhdlr != NULL);
1469
1470 SCIPnlhdlrSetCopyHdlr(nlhdlr, nlhdlrCopyhdlrBilinear);
1471 SCIPnlhdlrSetFreeHdlrData(nlhdlr, nlhdlrFreehdlrdataBilinear);
1472 SCIPnlhdlrSetFreeExprData(nlhdlr, nlhdlrFreeExprDataBilinear);
1473 SCIPnlhdlrSetInitExit(nlhdlr, nlhdlrInitBilinear, nlhdlrExitBilinear);
1474 SCIPnlhdlrSetSepa(nlhdlr, nlhdlrInitSepaBilinear, nlhdlrEnfoBilinear, nlhdlrEstimateBilinear, nlhdlrExitSepaBilinear);
1475 SCIPnlhdlrSetProp(nlhdlr, nlhdlrIntevalBilinear, nlhdlrReversepropBilinear);
1476
1477 /* parameters */
1478 SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/useinteval",
1479 "whether to use the interval evaluation callback of the nlhdlr",
1480 &nlhdlrdata->useinteval, FALSE, TRUE, NULL, NULL) );
1481
1482 SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/usereverseprop",
1483 "whether to use the reverse propagation callback of the nlhdlr",
1484 &nlhdlrdata->usereverseprop, FALSE, TRUE, NULL, NULL) );
1485
1486 SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/maxseparoundsroot",
1487 "maximum number of separation rounds in the root node",
1488 &nlhdlrdata->maxseparoundsroot, FALSE, 10, 0, INT_MAX, NULL, NULL) );
1489
1490 SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/maxseparounds",
1491 "maximum number of separation rounds in a local node",
1492 &nlhdlrdata->maxseparounds, FALSE, 1, 0, INT_MAX, NULL, NULL) );
1493
1494 SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/maxsepadepth",
1495 "maximum depth to apply separation",
1496 &nlhdlrdata->maxsepadepth, FALSE, INT_MAX, 0, INT_MAX, NULL, NULL) );
1497
1498 /* statistic table */
1499 assert(SCIPfindTable(scip, TABLE_NAME_BILINEAR) == NULL);
1500 SCIP_CALL( SCIPincludeTable(scip, TABLE_NAME_BILINEAR, TABLE_DESC_BILINEAR, FALSE,
1501 NULL, NULL, NULL, NULL, NULL, NULL, tableOutputBilinear,
1502 NULL, TABLE_POSITION_BILINEAR, TABLE_EARLIEST_STAGE_BILINEAR) );
1503
1504 return SCIP_OKAY;
1505 }
1506
1507 /** returns an array of expressions that have been detected by the bilinear nonlinear handler */
1508 SCIP_EXPR** SCIPgetNlhdlrBilinearExprs(
1509 SCIP_NLHDLR* nlhdlr /**< nonlinear handler */
1510 )
1511 {
1512 SCIP_NLHDLRDATA* nlhdlrdata;
1513
1514 assert(nlhdlr != NULL);
1515 assert(strcmp(SCIPnlhdlrGetName(nlhdlr), NLHDLR_NAME) == 0);
1516
1517 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1518 assert(nlhdlrdata);
1519
1520 return nlhdlrdata->exprs;
1521 }
1522
1523 /** returns the total number of expressions that have been detected by the bilinear nonlinear handler */
1524 int SCIPgetNlhdlrBilinearNExprs(
1525 SCIP_NLHDLR* nlhdlr /**< nonlinear handler */
1526 )
1527 {
1528 SCIP_NLHDLRDATA* nlhdlrdata;
1529
1530 assert(nlhdlr != NULL);
1531 assert(strcmp(SCIPnlhdlrGetName(nlhdlr), NLHDLR_NAME) == 0);
1532
1533 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1534 assert(nlhdlrdata);
1535
1536 return nlhdlrdata->nexprs;
1537 }
1538
1539 /** adds a globally valid inequality of the form \f$\text{xcoef}\cdot x \leq \text{ycoef} \cdot y + \text{constant}\f$ to a product expression of the form \f$x\cdot y\f$ */
1540 SCIP_RETCODE SCIPaddNlhdlrBilinearIneq(
1541 SCIP* scip, /**< SCIP data structure */
1542 SCIP_NLHDLR* nlhdlr, /**< nonlinear handler */
1543 SCIP_EXPR* expr, /**< product expression */
1544 SCIP_Real xcoef, /**< x coefficient */
1545 SCIP_Real ycoef, /**< y coefficient */
1546 SCIP_Real constant, /**< constant part */
1547 SCIP_Bool* success /**< buffer to store whether inequality has been accepted */
1548 )
1549 {
1550 SCIP_NLHDLREXPRDATA* nlhdlrexprdata;
1551 SCIP_VAR* x;
1552 SCIP_VAR* y;
1553 SCIP_Real* ineqs;
1554 SCIP_Real viol1;
1555 SCIP_Real viol2;
1556 SCIP_Bool underestimate;
1557 int nineqs;
1558 int i;
1559
1560 assert(scip != NULL);
1561 assert(nlhdlr != NULL);
1562 assert(strcmp(SCIPnlhdlrGetName(nlhdlr), NLHDLR_NAME) == 0);
1563 assert(expr != NULL);
1564 assert(SCIPexprGetNChildren(expr) == 2);
1565 assert(xcoef != SCIP_INVALID); /*lint !e777 */
1566 assert(ycoef != SCIP_INVALID); /*lint !e777 */
1567 assert(constant != SCIP_INVALID); /*lint !e777 */
1568 assert(success != NULL);
1569
1570 *success = FALSE;
1571
1572 /* find nonlinear handler expression handler data */
1573 nlhdlrexprdata = SCIPgetNlhdlrExprDataNonlinear(nlhdlr, expr);
1574
1575 if( nlhdlrexprdata == NULL )
1576 {
1577 SCIPwarningMessage(scip, "nonlinear expression data has not been found. "
1578 "Skip SCIPaddConsExprExprProductBilinearIneq()\n");
1579 return SCIP_OKAY;
1580 }
1581
1582 /* ignore inequalities that only yield to a (possible) bound tightening */
1583 if( SCIPisFeasZero(scip, xcoef) || SCIPisFeasZero(scip, ycoef) )
1584 return SCIP_OKAY;
1585
1586 /* collect variables */
1587 x = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[0]);
1588 y = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[1]);
1589 assert(x != NULL);
1590 assert(y != NULL);
1591 assert(x != y);
1592
1593 /* normalize inequality s.t. xcoef in {-1,1} */
1594 if( !SCIPisEQ(scip, REALABS(xcoef), 1.0) )
1595 {
1596 constant /= REALABS(xcoef);
1597 ycoef /= REALABS(xcoef);
1598 xcoef /= REALABS(xcoef);
1599 }
1600
1601 /* coefficients of the inequality determine whether the inequality can be used for under- or overestimation */
1602 underestimate = xcoef * ycoef > 0;
1603
1604 SCIPdebugMsg(scip, "add inequality for a bilinear term: %g %s <= %g %s + %g (underestimate=%u)\n", xcoef,
1605 SCIPvarGetName(x), ycoef, SCIPvarGetName(y), constant, underestimate);
1606
1607 /* compute violation of the inequality of the important corner points */
1608 getIneqViol(x, y, xcoef, ycoef, constant, &viol1, &viol2);
1609 SCIPdebugMsg(scip, "violations of inequality = (%g,%g)\n", viol1, viol2);
1610
1611 /* inequality does not cutoff one of the important corner points -> skip */
1612 if( SCIPisFeasLE(scip, MAX(viol1, viol2), 0.0) )
1613 return SCIP_OKAY;
1614
1615 if( underestimate )
1616 {
1617 ineqs = nlhdlrexprdata->underineqs;
1618 nineqs = nlhdlrexprdata->nunderineqs;
1619 }
1620 else
1621 {
1622 ineqs = nlhdlrexprdata->overineqs;
1623 nineqs = nlhdlrexprdata->noverineqs;
1624 }
1625
1626 /* check for a duplicate */
1627 for( i = 0; i < nineqs; ++i )
1628 {
1629 if( SCIPisFeasEQ(scip, xcoef, ineqs[3*i]) && SCIPisFeasEQ(scip, ycoef, ineqs[3*i+1])
1630 && SCIPisFeasEQ(scip, constant, ineqs[3*i+2]) )
1631 {
1632 SCIPdebugMsg(scip, "inequality already found -> skip\n");
1633 return SCIP_OKAY;
1634 }
1635 }
1636
1637 {
1638 SCIP_Real viols1[2] = {0.0, 0.0};
1639 SCIP_Real viols2[2] = {0.0, 0.0};
1640
1641 /* compute violations of existing inequalities */
1642 for( i = 0; i < nineqs; ++i )
1643 {
1644 getIneqViol(x, y, ineqs[3*i], ineqs[3*i+1], ineqs[3*i+2], &viols1[i], &viols2[i]);
1645
1646 /* check whether an existing inequality is dominating the candidate */
1647 if( SCIPisGE(scip, viols1[i], viol1) && SCIPisGE(scip, viols2[i], viol2) )
1648 {
1649 SCIPdebugMsg(scip, "inequality is dominated by %d -> skip\n", i);
1650 return SCIP_OKAY;
1651 }
1652
1653 /* replace inequality if candidate is dominating it */
1654 if( SCIPisLT(scip, viols1[i], viol1) && SCIPisLT(scip, viols2[i], viol2) )
1655 {
1656 SCIPdebugMsg(scip, "inequality dominates %d -> replace\n", i);
1657 ineqs[3*i] = xcoef;
1658 ineqs[3*i+1] = ycoef;
1659 ineqs[3*i+2] = constant;
1660 *success = TRUE;
1661 }
1662 }
1663
1664 /* inequality is not dominated by other inequalities -> add if we have less than 2 inequalities */
1665 if( nineqs < 2 )
1666 {
1667 ineqs[3*nineqs] = xcoef;
1668 ineqs[3*nineqs + 1] = ycoef;
1669 ineqs[3*nineqs + 2] = constant;
1670 *success = TRUE;
1671 SCIPdebugMsg(scip, "add inequality\n");
1672
1673 /* increase number of inequalities */
1674 if( underestimate )
1675 ++(nlhdlrexprdata->nunderineqs);
1676 else
1677 ++(nlhdlrexprdata->noverineqs);
1678 }
1679 }
1680
1681 if( *success )
1682 {
1683 /* With the added inequalities, we can potentially compute tighter activities for the expression,
1684 * so constraints that contain this expression should be propagated again.
1685 * We don't have a direct expression to constraint mapping, though. This call marks all expr-constraints
1686 * which include any of the variables that this expression depends on for propagation.
1687 */
1688 SCIP_CALL( SCIPmarkExprPropagateNonlinear(scip, expr) );
1689 }
1690
1691 return SCIP_OKAY;
1692 }
1693
1694 /** computes coefficients of linearization of a bilinear term in a reference point */
1695 void SCIPaddBilinLinearization(
1696 SCIP* scip, /**< SCIP data structure */
1697 SCIP_Real bilincoef, /**< coefficient of bilinear term */
1698 SCIP_Real refpointx, /**< point where to linearize first variable */
1699 SCIP_Real refpointy, /**< point where to linearize second variable */
1700 SCIP_Real* lincoefx, /**< buffer to add coefficient of first variable in linearization */
1701 SCIP_Real* lincoefy, /**< buffer to add coefficient of second variable in linearization */
1702 SCIP_Real* linconstant, /**< buffer to add constant of linearization */
1703 SCIP_Bool* success /**< buffer to set to FALSE if linearization has failed due to large numbers */
1704 )
1705 {
1706 SCIP_Real constant;
1707
1708 assert(scip != NULL);
1709 assert(lincoefx != NULL);
1710 assert(lincoefy != NULL);
1711 assert(linconstant != NULL);
1712 assert(success != NULL);
1713
1714 if( bilincoef == 0.0 )
1715 return;
1716
1717 if( SCIPisInfinity(scip, REALABS(refpointx)) || SCIPisInfinity(scip, REALABS(refpointy)) )
1718 {
1719 *success = FALSE;
1720 return;
1721 }
1722
1723 /* bilincoef * x * y -> bilincoef * (refpointx * refpointy + refpointy * (x - refpointx) + refpointx * (y - refpointy))
1724 * = -bilincoef * refpointx * refpointy + bilincoef * refpointy * x + bilincoef * refpointx * y
1725 */
1726
1727 constant = -bilincoef * refpointx * refpointy;
1728
1729 if( SCIPisInfinity(scip, REALABS(bilincoef * refpointx)) || SCIPisInfinity(scip, REALABS(bilincoef * refpointy))
1730 || SCIPisInfinity(scip, REALABS(constant)) )
1731 {
1732 *success = FALSE;
1733 return;
1734 }
1735
1736 *lincoefx += bilincoef * refpointy;
1737 *lincoefy += bilincoef * refpointx;
1738 *linconstant += constant;
1739 }
1740
1741 /** computes coefficients of McCormick under- or overestimation of a bilinear term */
1742 void SCIPaddBilinMcCormick(
1743 SCIP* scip, /**< SCIP data structure */
1744 SCIP_Real bilincoef, /**< coefficient of bilinear term */
1745 SCIP_Real lbx, /**< lower bound on first variable */
1746 SCIP_Real ubx, /**< upper bound on first variable */
1747 SCIP_Real refpointx, /**< reference point for first variable */
1748 SCIP_Real lby, /**< lower bound on second variable */
1749 SCIP_Real uby, /**< upper bound on second variable */
1750 SCIP_Real refpointy, /**< reference point for second variable */
1751 SCIP_Bool overestimate, /**< whether to compute an overestimator instead of an underestimator */
1752 SCIP_Real* lincoefx, /**< buffer to add coefficient of first variable in linearization */
1753 SCIP_Real* lincoefy, /**< buffer to add coefficient of second variable in linearization */
1754 SCIP_Real* linconstant, /**< buffer to add constant of linearization */
1755 SCIP_Bool* success /**< buffer to set to FALSE if linearization has failed due to large numbers */
1756 )
1757 {
1758 SCIP_Real constant;
1759 SCIP_Real coefx;
1760 SCIP_Real coefy;
1761
1762 assert(scip != NULL);
1763 assert(!SCIPisInfinity(scip, lbx));
1764 assert(!SCIPisInfinity(scip, -ubx));
1765 assert(!SCIPisInfinity(scip, lby));
1766 assert(!SCIPisInfinity(scip, -uby));
1767 assert(SCIPisInfinity(scip, -lbx) || SCIPisLE(scip, lbx, ubx));
1768 assert(SCIPisInfinity(scip, -lby) || SCIPisLE(scip, lby, uby));
1769 assert(SCIPisInfinity(scip, -lbx) || SCIPisLE(scip, lbx, refpointx));
1770 assert(SCIPisInfinity(scip, -lby) || SCIPisLE(scip, lby, refpointy));
1771 assert(SCIPisInfinity(scip, ubx) || SCIPisGE(scip, ubx, refpointx));
1772 assert(SCIPisInfinity(scip, uby) || SCIPisGE(scip, uby, refpointy));
1773 assert(lincoefx != NULL);
1774 assert(lincoefy != NULL);
1775 assert(linconstant != NULL);
1776 assert(success != NULL);
1777
1778 if( bilincoef == 0.0 )
1779 return;
1780
1781 if( overestimate )
1782 bilincoef = -bilincoef;
1783
1784 if( SCIPisRelEQ(scip, lbx, ubx) && SCIPisRelEQ(scip, lby, uby) )
1785 {
1786 /* both x and y are mostly fixed */
1787 SCIP_Real cand1;
1788 SCIP_Real cand2;
1789 SCIP_Real cand3;
1790 SCIP_Real cand4;
1791
1792 coefx = 0.0;
1793 coefy = 0.0;
1794
1795 /* estimate x * y by constant */
1796 cand1 = lbx * lby;
1797 cand2 = lbx * uby;
1798 cand3 = ubx * lby;
1799 cand4 = ubx * uby;
1800
1801 /* take most conservative value for underestimator */
1802 if( bilincoef < 0.0 )
1803 constant = bilincoef * MAX( MAX(cand1, cand2), MAX(cand3, cand4) );
1804 else
1805 constant = bilincoef * MIN( MIN(cand1, cand2), MIN(cand3, cand4) );
1806 }
1807 else if( bilincoef > 0.0 )
1808 {
1809 /* either x or y is not fixed and coef > 0.0 */
1810 if( !SCIPisInfinity(scip, -lbx) && !SCIPisInfinity(scip, -lby) &&
1811 (SCIPisInfinity(scip, ubx) || SCIPisInfinity(scip, uby)
1812 || (uby - refpointy) * (ubx - refpointx) >= (refpointy - lby) * (refpointx - lbx)) )
1813 {
1814 if( SCIPisRelEQ(scip, lbx, ubx) )
1815 {
1816 /* x*y = lbx * y + (x-lbx) * y >= lbx * y + (x-lbx) * lby >= lbx * y + min{(ubx-lbx) * lby, 0 * lby} */
1817 coefx = 0.0;
1818 coefy = bilincoef * lbx;
1819 constant = bilincoef * (lby < 0.0 ? (ubx-lbx) * lby : 0.0);
1820 }
1821 else if( SCIPisRelEQ(scip, lby, uby) )
1822 {
1823 /* x*y = lby * x + (y-lby) * x >= lby * x + (y-lby) * lbx >= lby * x + min{(uby-lby) * lbx, 0 * lbx} */
1824 coefx = bilincoef * lby;
1825 coefy = 0.0;
1826 constant = bilincoef * (lbx < 0.0 ? (uby-lby) * lbx : 0.0);
1827 }
1828 else
1829 {
1830 coefx = bilincoef * lby;
1831 coefy = bilincoef * lbx;
1832 constant = -bilincoef * lbx * lby;
1833 }
1834 }
1835 else if( !SCIPisInfinity(scip, ubx) && !SCIPisInfinity(scip, uby) )
1836 {
1837 if( SCIPisRelEQ(scip, lbx, ubx) )
1838 {
1839 /* x*y = ubx * y + (x-ubx) * y >= ubx * y + (x-ubx) * uby >= ubx * y + min{(lbx-ubx) * uby, 0 * uby} */
1840 coefx = 0.0;
1841 coefy = bilincoef * ubx;
1842 constant = bilincoef * (uby > 0.0 ? (lbx - ubx) * uby : 0.0);
1843 }
1844 else if( SCIPisRelEQ(scip, lby, uby) )
1845 {
1846 /* x*y = uby * x + (y-uby) * x >= uby * x + (y-uby) * ubx >= uby * x + min{(lby-uby) * ubx, 0 * ubx} */
1847 coefx = bilincoef * uby;
1848 coefy = 0.0;
1849 constant = bilincoef * (ubx > 0.0 ? (lby - uby) * ubx : 0.0);
1850 }
1851 else
1852 {
1853 coefx = bilincoef * uby;
1854 coefy = bilincoef * ubx;
1855 constant = -bilincoef * ubx * uby;
1856 }
1857 }
1858 else
1859 {
1860 *success = FALSE;
1861 return;
1862 }
1863 }
1864 else
1865 {
1866 /* either x or y is not fixed and coef < 0.0 */
1867 if( !SCIPisInfinity(scip, ubx) && !SCIPisInfinity(scip, -lby) &&
1868 (SCIPisInfinity(scip, -lbx) || SCIPisInfinity(scip, uby)
1869 || (ubx - lbx) * (refpointy - lby) <= (uby - lby) * (refpointx - lbx)) )
1870 {
1871 if( SCIPisRelEQ(scip, lbx, ubx) )
1872 {
1873 /* x*y = ubx * y + (x-ubx) * y <= ubx * y + (x-ubx) * lby <= ubx * y + max{(lbx-ubx) * lby, 0 * lby} */
1874 coefx = 0.0;
1875 coefy = bilincoef * ubx;
1876 constant = bilincoef * (lby < 0.0 ? (lbx - ubx) * lby : 0.0);
1877 }
1878 else if( SCIPisRelEQ(scip, lby, uby) )
1879 {
1880 /* x*y = lby * x + (y-lby) * x <= lby * x + (y-lby) * ubx <= lby * x + max{(uby-lby) * ubx, 0 * ubx} */
1881 coefx = bilincoef * lby;
1882 coefy = 0.0;
1883 constant = bilincoef * (ubx > 0.0 ? (uby - lby) * ubx : 0.0);
1884 }
1885 else
1886 {
1887 coefx = bilincoef * lby;
1888 coefy = bilincoef * ubx;
1889 constant = -bilincoef * ubx * lby;
1890 }
1891 }
1892 else if( !SCIPisInfinity(scip, -lbx) && !SCIPisInfinity(scip, uby) )
1893 {
1894 if( SCIPisRelEQ(scip, lbx, ubx) )
1895 {
1896 /* x*y = lbx * y + (x-lbx) * y <= lbx * y + (x-lbx) * uby <= lbx * y + max{(ubx-lbx) * uby, 0 * uby} */
1897 coefx = 0.0;
1898 coefy = bilincoef * lbx;
1899 constant = bilincoef * (uby > 0.0 ? (ubx - lbx) * uby : 0.0);
1900 }
1901 else if( SCIPisRelEQ(scip, lby, uby) )
1902 {
1903 /* x*y = uby * x + (y-uby) * x <= uby * x + (y-uby) * lbx <= uby * x + max{(lby-uby) * lbx, 0 * lbx} */
1904 coefx = bilincoef * uby;
1905 coefy = 0.0;
1906 constant = bilincoef * (lbx < 0.0 ? (lby - uby) * lbx : 0.0);
1907 }
1908 else
1909 {
1910 coefx = bilincoef * uby;
1911 coefy = bilincoef * lbx;
1912 constant = -bilincoef * lbx * uby;
1913 }
1914 }
1915 else
1916 {
1917 *success = FALSE;
1918 return;
1919 }
1920 }
1921
1922 if( SCIPisInfinity(scip, REALABS(coefx)) || SCIPisInfinity(scip, REALABS(coefy))
1923 || SCIPisInfinity(scip, REALABS(constant)) )
1924 {
1925 *success = FALSE;
1926 return;
1927 }
1928
1929 if( overestimate )
1930 {
1931 coefx = -coefx;
1932 coefy = -coefy;
1933 constant = -constant;
1934 }
1935
1936 SCIPdebugMsg(scip, "%.15g * x[%.15g,%.15g] * y[%.15g,%.15g] %c= %.15g * x %+.15g * y %+.15g\n", bilincoef, lbx, ubx,
1937 lby, uby, overestimate ? '<' : '>', coefx, coefy, constant);
1938
1939 *lincoefx += coefx;
1940 *lincoefy += coefy;
1941 *linconstant += constant;
1942 }
1943
1944 /** computes coefficients of linearization of a bilinear term in a reference point when given a linear inequality
1945 * involving only the variables of the bilinear term
1946 *
1947 * @note the formulas are extracted from "Convex envelopes of bivariate functions through the solution of KKT systems"
1948 * by Marco Locatelli
1949 */
1950 void SCIPcomputeBilinEnvelope1(
1951 SCIP* scip, /**< SCIP data structure */
1952 SCIP_Real bilincoef, /**< coefficient of bilinear term */
1953 SCIP_Real lbx, /**< lower bound on first variable */
1954 SCIP_Real ubx, /**< upper bound on first variable */
1955 SCIP_Real refpointx, /**< reference point for first variable */
1956 SCIP_Real lby, /**< lower bound on second variable */
1957 SCIP_Real uby, /**< upper bound on second variable */
1958 SCIP_Real refpointy, /**< reference point for second variable */
1959 SCIP_Bool overestimate, /**< whether to compute an overestimator instead of an underestimator */
1960 SCIP_Real xcoef, /**< x coefficient of linear inequality; must be in {-1,0,1} */
1961 SCIP_Real ycoef, /**< y coefficient of linear inequality */
1962 SCIP_Real constant, /**< constant of linear inequality */
1963 SCIP_Real* RESTRICT lincoefx, /**< buffer to store coefficient of first variable in linearization */
1964 SCIP_Real* RESTRICT lincoefy, /**< buffer to store coefficient of second variable in linearization */
1965 SCIP_Real* RESTRICT linconstant, /**< buffer to store constant of linearization */
1966 SCIP_Bool* RESTRICT success /**< buffer to store whether linearization was successful */
1967 )
1968 {
1969 SCIP_Real xs[2] = {lbx, ubx};
1970 SCIP_Real ys[2] = {lby, uby};
1971 SCIP_Real minx;
1972 SCIP_Real maxx;
1973 SCIP_Real miny;
1974 SCIP_Real maxy;
1975 SCIP_Real QUAD(lincoefyq);
1976 SCIP_Real QUAD(lincoefxq);
1977 SCIP_Real QUAD(linconstantq);
1978 SCIP_Real QUAD(denomq);
1979 SCIP_Real QUAD(mjq);
1980 SCIP_Real QUAD(qjq);
1981 SCIP_Real QUAD(xjq);
1982 SCIP_Real QUAD(yjq);
1983 SCIP_Real QUAD(tmpq);
1984 SCIP_Real vx;
1985 SCIP_Real vy;
1986 int n;
1987 int i;
1988
1989 assert(scip != NULL);
1990 assert(!SCIPisInfinity(scip, lbx));
1991 assert(!SCIPisInfinity(scip, -ubx));
1992 assert(!SCIPisInfinity(scip, lby));
1993 assert(!SCIPisInfinity(scip, -uby));
1994 assert(SCIPisLE(scip, lbx, ubx));
1995 assert(SCIPisLE(scip, lby, uby));
1996 assert(SCIPisLE(scip, lbx, refpointx));
1997 assert(SCIPisGE(scip, ubx, refpointx));
1998 assert(SCIPisLE(scip, lby, refpointy));
1999 assert(SCIPisGE(scip, uby, refpointy));
2000 assert(lincoefx != NULL);
2001 assert(lincoefy != NULL);
2002 assert(linconstant != NULL);
2003 assert(success != NULL);
2004 assert(xcoef == 0.0 || xcoef == -1.0 || xcoef == 1.0); /*lint !e777*/
2005 assert(ycoef != SCIP_INVALID && ycoef != 0.0); /*lint !e777*/
2006 assert(constant != SCIP_INVALID); /*lint !e777*/
2007
2008 *success = FALSE;
2009 *lincoefx = SCIP_INVALID;
2010 *lincoefy = SCIP_INVALID;
2011 *linconstant = SCIP_INVALID;
2012
2013 /* reference point does not satisfy linear inequality */
2014 if( SCIPisFeasGT(scip, xcoef * refpointx - ycoef * refpointy - constant, 0.0) )
2015 return;
2016
2017 /* compute minimal and maximal bounds on x and y for accepting the reference point */
2018 minx = lbx + 0.01 * (ubx-lbx);
2019 maxx = ubx - 0.01 * (ubx-lbx);
2020 miny = lby + 0.01 * (uby-lby);
2021 maxy = uby - 0.01 * (uby-lby);
2022
2023 /* check whether the reference point is in [minx,maxx]x[miny,maxy] */
2024 if( SCIPisLE(scip, refpointx, minx) || SCIPisGE(scip, refpointx, maxx)
2025 || SCIPisLE(scip, refpointy, miny) || SCIPisGE(scip, refpointy, maxy) )
2026 return;
2027
2028 /* always consider xy without the bilinear coefficient */
2029 if( bilincoef < 0.0 )
2030 overestimate = !overestimate;
2031
2032 /* we use same notation as in "Convex envelopes of bivariate functions through the solution of KKT systems", 2016 */
2033 /* mj = xcoef / ycoef */
2034 SCIPquadprecDivDD(mjq, xcoef, ycoef);
2035
2036 /* qj = -constant / ycoef */
2037 SCIPquadprecDivDD(qjq, -constant, ycoef);
2038
2039 /* mj > 0 => underestimate; mj < 0 => overestimate */
2040 if( SCIPisNegative(scip, QUAD_TO_DBL(mjq)) != overestimate )
2041 return;
2042
2043 /* get the corner point that satisfies the linear inequality xcoef*x <= ycoef*y + constant */
2044 if( !overestimate )
2045 {
2046 ys[0] = uby;
2047 ys[1] = lby;
2048 }
2049
2050 vx = SCIP_INVALID;
2051 vy = SCIP_INVALID;
2052 n = 0;
2053 for( i = 0; i < 2; ++i )
2054 {
2055 SCIP_Real activity = xcoef * xs[i] - ycoef * ys[i] - constant;
2056 if( SCIPisLE(scip, activity, 0.0) )
2057 {
2058 /* corner point is satisfies inequality */
2059 vx = xs[i];
2060 vy = ys[i];
2061 }
2062 else if( SCIPisFeasGT(scip, activity, 0.0) )
2063 /* corner point is clearly cut off */
2064 ++n;
2065 }
2066
2067 /* skip if no corner point satisfies the inequality or if no corner point is cut off
2068 * (that is, all corner points satisfy the inequality almost [1e-9..1e-6]) */
2069 if( n != 1 || vx == SCIP_INVALID || vy == SCIP_INVALID ) /*lint !e777*/
2070 return;
2071
2072 /* denom = mj*(refpointx - vx) + vy - refpointy */
2073 SCIPquadprecSumDD(denomq, refpointx, -vx); /* refpoint - vx */
2074 SCIPquadprecProdQQ(denomq, denomq, mjq); /* mj * (refpoint - vx) */
2075 SCIPquadprecSumQD(denomq, denomq, vy); /* mj * (refpoint - vx) + vy */
2076 SCIPquadprecSumQD(denomq, denomq, -refpointy); /* mj * (refpoint - vx) + vy - refpointy */
2077
2078 if( SCIPisZero(scip, QUAD_TO_DBL(denomq)) )
2079 return;
2080
2081 /* (xj,yj) is the projection onto the line xcoef*x = ycoef*y + constant */
2082 /* xj = (refpointx*(vy - qj) - vx*(refpointy - qj)) / denom */
2083 SCIPquadprecProdQD(xjq, qjq, -1.0); /* - qj */
2084 SCIPquadprecSumQD(xjq, xjq, vy); /* vy - qj */
2085 SCIPquadprecProdQD(xjq, xjq, refpointx); /* refpointx * (vy - qj) */
2086 SCIPquadprecProdQD(tmpq, qjq, -1.0); /* - qj */
2087 SCIPquadprecSumQD(tmpq, tmpq, refpointy); /* refpointy - qj */
2088 SCIPquadprecProdQD(tmpq, tmpq, -vx); /* - vx * (refpointy - qj) */
2089 SCIPquadprecSumQQ(xjq, xjq, tmpq); /* refpointx * (vy - qj) - vx * (refpointy - qj) */
2090 SCIPquadprecDivQQ(xjq, xjq, denomq); /* (refpointx * (vy - qj) - vx * (refpointy - qj)) / denom */
2091
2092 /* yj = mj * xj + qj */
2093 SCIPquadprecProdQQ(yjq, mjq, xjq);
2094 SCIPquadprecSumQQ(yjq, yjq, qjq);
2095
2096 assert(SCIPisFeasEQ(scip, xcoef*QUAD_TO_DBL(xjq) - ycoef*QUAD_TO_DBL(yjq) - constant, 0.0));
2097
2098 /* check whether the projection is in [minx,maxx] x [miny,maxy]; this avoids numerical difficulties when the
2099 * projection is close to the variable bounds
2100 */
2101 if( SCIPisLE(scip, QUAD_TO_DBL(xjq), minx) || SCIPisGE(scip, QUAD_TO_DBL(xjq), maxx)
2102 || SCIPisLE(scip, QUAD_TO_DBL(yjq), miny) || SCIPisGE(scip, QUAD_TO_DBL(yjq), maxy) )
2103 return;
2104
2105 assert(vy - QUAD_TO_DBL(mjq)*vx - QUAD_TO_DBL(qjq) != 0.0);
2106
2107 /* lincoefy = (mj*SQR(xj) - 2.0*mj*vx*xj - qj*vx + vx*vy) / (vy - mj*vx - qj) */
2108 SCIPquadprecSquareQ(lincoefyq, xjq); /* xj^2 */
2109 SCIPquadprecProdQQ(lincoefyq, lincoefyq, mjq); /* mj * xj^2 */
2110 SCIPquadprecProdQQ(tmpq, mjq, xjq); /* mj * xj */
2111 SCIPquadprecProdQD(tmpq, tmpq, -2.0 * vx); /* -2 * vx * mj * xj */
2112 SCIPquadprecSumQQ(lincoefyq, lincoefyq, tmpq); /* mj * xj^2 -2 * vx * mj * xj */
2113 SCIPquadprecProdQD(tmpq, qjq, -vx); /* -qj * vx */
2114 SCIPquadprecSumQQ(lincoefyq, lincoefyq, tmpq); /* mj * xj^2 -2 * vx * mj * xj -qj * vx */
2115 SCIPquadprecProdDD(tmpq, vx, vy); /* vx * vy */
2116 SCIPquadprecSumQQ(lincoefyq, lincoefyq, tmpq); /* mj * xj^2 -2 * vx * mj * xj -qj * vx + vx * vy */
2117 SCIPquadprecProdQD(tmpq, mjq, vx); /* mj * vx */
2118 SCIPquadprecSumQD(tmpq, tmpq, -vy); /* -vy + mj * vx */
2119 SCIPquadprecSumQQ(tmpq, tmpq, qjq); /* -vy + mj * vx + qj */
2120 QUAD_SCALE(tmpq, -1.0); /* vy - mj * vx - qj */
2121 SCIPquadprecDivQQ(lincoefyq, lincoefyq, tmpq); /* (mj * xj^2 -2 * vx * mj * xj -qj * vx + vx * vy) / (vy - mj * vx - qj) */
2122
2123 /* lincoefx = 2.0*mj*xj + qj - mj*(*lincoefy) */
2124 SCIPquadprecProdQQ(lincoefxq, mjq, xjq); /* mj * xj */
2125 QUAD_SCALE(lincoefxq, 2.0); /* 2 * mj * xj */
2126 SCIPquadprecSumQQ(lincoefxq, lincoefxq, qjq); /* 2 * mj * xj + qj */
2127 SCIPquadprecProdQQ(tmpq, mjq, lincoefyq); /* mj * lincoefy */
2128 QUAD_SCALE(tmpq, -1.0); /* - mj * lincoefy */
2129 SCIPquadprecSumQQ(lincoefxq, lincoefxq, tmpq); /* 2 * mj * xj + qj - mj * lincoefy */
2130
2131 /* linconstant = -mj*SQR(xj) - (*lincoefy)*qj */
2132 SCIPquadprecSquareQ(linconstantq, xjq); /* xj^2 */
2133 SCIPquadprecProdQQ(linconstantq, linconstantq, mjq); /* mj * xj^2 */
2134 QUAD_SCALE(linconstantq, -1.0); /* - mj * xj^2 */
2135 SCIPquadprecProdQQ(tmpq, lincoefyq, qjq); /* lincoefy * qj */
2136 QUAD_SCALE(tmpq, -1.0); /* - lincoefy * qj */
2137 SCIPquadprecSumQQ(linconstantq, linconstantq, tmpq); /* - mj * xj^2 - lincoefy * qj */
2138
2139 /* consider the bilinear coefficient */
2140 SCIPquadprecProdQD(lincoefxq, lincoefxq, bilincoef);
2141 SCIPquadprecProdQD(lincoefyq, lincoefyq, bilincoef);
2142 SCIPquadprecProdQD(linconstantq, linconstantq, bilincoef);
2143 *lincoefx = QUAD_TO_DBL(lincoefxq);
2144 *lincoefy = QUAD_TO_DBL(lincoefyq);
2145 *linconstant = QUAD_TO_DBL(linconstantq);
2146
2147 /* cut needs to be tight at (vx,vy) and (xj,yj); otherwise we consider the cut to be numerically bad */
2148 *success = SCIPisFeasEQ(scip, (*lincoefx)*vx + (*lincoefy)*vy + (*linconstant), bilincoef*vx*vy)
2149 && SCIPisFeasEQ(scip, (*lincoefx)*QUAD_TO_DBL(xjq) + (*lincoefy)*QUAD_TO_DBL(yjq) + (*linconstant),
2150 bilincoef*QUAD_TO_DBL(xjq)*QUAD_TO_DBL(yjq));
2151
2152 #ifndef NDEBUG
2153 {
2154 SCIP_Real activity = (*lincoefx)*refpointx + (*lincoefy)*refpointy + (*linconstant);
2155
2156 /* cut needs to under- or overestimate the bilinear term at the reference point */
2157 if( bilincoef < 0.0 )
2158 overestimate = !overestimate;
2159
2160 if( overestimate )
2161 assert(SCIPisFeasGE(scip, activity, bilincoef*refpointx*refpointy));
2162 else
2163 assert(SCIPisFeasLE(scip, activity, bilincoef*refpointx*refpointy));
2164 }
2165 #endif
2166 }
2167
2168 /** computes coefficients of linearization of a bilinear term in a reference point when given two linear inequalities
2169 * involving only the variables of the bilinear term
2170 *
2171 * @note the formulas are extracted from "Convex envelopes of bivariate functions through the solution of KKT systems"
2172 * by Marco Locatelli
2173 */
2174 void SCIPcomputeBilinEnvelope2(
2175 SCIP* scip, /**< SCIP data structure */
2176 SCIP_Real bilincoef, /**< coefficient of bilinear term */
2177 SCIP_Real lbx, /**< lower bound on first variable */
2178 SCIP_Real ubx, /**< upper bound on first variable */
2179 SCIP_Real refpointx, /**< reference point for first variable */
2180 SCIP_Real lby, /**< lower bound on second variable */
2181 SCIP_Real uby, /**< upper bound on second variable */
2182 SCIP_Real refpointy, /**< reference point for second variable */
2183 SCIP_Bool overestimate, /**< whether to compute an overestimator instead of an underestimator */
2184 SCIP_Real xcoef1, /**< x coefficient of linear inequality; must be in {-1,0,1} */
2185 SCIP_Real ycoef1, /**< y coefficient of linear inequality */
2186 SCIP_Real constant1, /**< constant of linear inequality */
2187 SCIP_Real xcoef2, /**< x coefficient of linear inequality; must be in {-1,0,1} */
2188 SCIP_Real ycoef2, /**< y coefficient of linear inequality */
2189 SCIP_Real constant2, /**< constant of linear inequality */
2190 SCIP_Real* RESTRICT lincoefx, /**< buffer to store coefficient of first variable in linearization */
2191 SCIP_Real* RESTRICT lincoefy, /**< buffer to store coefficient of second variable in linearization */
2192 SCIP_Real* RESTRICT linconstant, /**< buffer to store constant of linearization */
2193 SCIP_Bool* RESTRICT success /**< buffer to store whether linearization was successful */
2194 )
2195 {
2196 SCIP_Real mi, mj, qi, qj, xi, xj, yi, yj;
2197 SCIP_Real xcoef, ycoef, constant;
2198 SCIP_Real minx, maxx, miny, maxy;
2199
2200 assert(scip != NULL);
2201 assert(!SCIPisInfinity(scip, lbx));
2202 assert(!SCIPisInfinity(scip, -ubx));
2203 assert(!SCIPisInfinity(scip, lby));
2204 assert(!SCIPisInfinity(scip, -uby));
2205 assert(SCIPisLE(scip, lbx, ubx));
2206 assert(SCIPisLE(scip, lby, uby));
2207 assert(SCIPisLE(scip, lbx, refpointx));
2208 assert(SCIPisGE(scip, ubx, refpointx));
2209 assert(SCIPisLE(scip, lby, refpointy));
2210 assert(SCIPisGE(scip, uby, refpointy));
2211 assert(lincoefx != NULL);
2212 assert(lincoefy != NULL);
2213 assert(linconstant != NULL);
2214 assert(success != NULL);
2215 assert(xcoef1 != 0.0 && xcoef1 != SCIP_INVALID); /*lint !e777*/
2216 assert(ycoef1 != SCIP_INVALID && ycoef1 != 0.0); /*lint !e777*/
2217 assert(constant1 != SCIP_INVALID); /*lint !e777*/
2218 assert(xcoef2 != 0.0 && xcoef2 != SCIP_INVALID); /*lint !e777*/
2219 assert(ycoef2 != SCIP_INVALID && ycoef2 != 0.0); /*lint !e777*/
2220 assert(constant2 != SCIP_INVALID); /*lint !e777*/
2221
2222 *success = FALSE;
2223 *lincoefx = SCIP_INVALID;
2224 *lincoefy = SCIP_INVALID;
2225 *linconstant = SCIP_INVALID;
2226
2227 /* reference point does not satisfy linear inequalities */
2228 if( SCIPisFeasGT(scip, xcoef1 * refpointx - ycoef1 * refpointy - constant1, 0.0)
2229 || SCIPisFeasGT(scip, xcoef2 * refpointx - ycoef2 * refpointy - constant2, 0.0) )
2230 return;
2231
2232 /* compute minimal and maximal bounds on x and y for accepting the reference point */
2233 minx = lbx + 0.01 * (ubx-lbx);
2234 maxx = ubx - 0.01 * (ubx-lbx);
2235 miny = lby + 0.01 * (uby-lby);
2236 maxy = uby - 0.01 * (uby-lby);
2237
2238 /* check the reference point is in the interior of the domain */
2239 if( SCIPisLE(scip, refpointx, minx) || SCIPisGE(scip, refpointx, maxx)
2240 || SCIPisLE(scip, refpointy, miny) || SCIPisFeasGE(scip, refpointy, maxy) )
2241 return;
2242
2243 /* the sign of the x-coefficients of the two inequalities must be different; otherwise the convex or concave
2244 * envelope can be computed via SCIPcomputeBilinEnvelope1 for each inequality separately
2245 */
2246 if( (xcoef1 > 0) == (xcoef2 > 0) )
2247 return;
2248
2249 /* always consider xy without the bilinear coefficient */
2250 if( bilincoef < 0.0 )
2251 overestimate = !overestimate;
2252
2253 /* we use same notation as in "Convex envelopes of bivariate functions through the solution of KKT systems", 2016 */
2254 mi = xcoef1 / ycoef1;
2255 qi = -constant1 / ycoef1;
2256 mj = xcoef2 / ycoef2;
2257 qj = -constant2 / ycoef2;
2258
2259 /* mi, mj > 0 => underestimate; mi, mj < 0 => overestimate */
2260 if( SCIPisNegative(scip, mi) != overestimate || SCIPisNegative(scip, mj) != overestimate )
2261 return;
2262
2263 /* compute cut according to Locatelli 2016 */
2264 computeBilinEnvelope2(scip, refpointx, refpointy, mi, qi, mj, qj, &xi, &yi, &xj, &yj, &xcoef, &ycoef, &constant);
2265 assert(SCIPisRelEQ(scip, mi*xi + qi, yi));
2266 assert(SCIPisRelEQ(scip, mj*xj + qj, yj));
2267
2268 /* it might happen that (xi,yi) = (xj,yj) if the two lines intersect */
2269 if( SCIPisEQ(scip, xi, xj) && SCIPisEQ(scip, yi, yj) )
2270 return;
2271
2272 /* check whether projected points are in the interior */
2273 if( SCIPisLE(scip, xi, minx) || SCIPisGE(scip, xi, maxx) || SCIPisLE(scip, yi, miny) || SCIPisGE(scip, yi, maxy) )
2274 return;
2275 if( SCIPisLE(scip, xj, minx) || SCIPisGE(scip, xj, maxx) || SCIPisLE(scip, yj, miny) || SCIPisGE(scip, yj, maxy) )
2276 return;
2277
2278 *lincoefx = bilincoef * xcoef;
2279 *lincoefy = bilincoef * ycoef;
2280 *linconstant = bilincoef * constant;
2281
2282 /* cut needs to be tight at (vx,vy) and (xj,yj) */
2283 *success = SCIPisFeasEQ(scip, (*lincoefx)*xi + (*lincoefy)*yi + (*linconstant), bilincoef*xi*yi)
2284 && SCIPisFeasEQ(scip, (*lincoefx)*xj + (*lincoefy)*yj + (*linconstant), bilincoef*xj*yj);
2285
2286 #ifndef NDEBUG
2287 {
2288 SCIP_Real activity = (*lincoefx)*refpointx + (*lincoefy)*refpointy + (*linconstant);
2289
2290 /* cut needs to under- or overestimate the bilinear term at the reference point */
2291 if( bilincoef < 0.0 )
2292 overestimate = !overestimate;
2293
2294 if( overestimate )
2295 assert(SCIPisFeasGE(scip, activity, bilincoef*refpointx*refpointy));
2296 else
2297 assert(SCIPisFeasLE(scip, activity, bilincoef*refpointx*refpointy));
2298 }
2299 #endif
2300 }
2301