Scippy

SCIP

Solving Constraint Integer Programs

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