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