1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 2 /* */ 3 /* This file is part of the program and library */ 4 /* SCIP --- Solving Constraint Integer Programs */ 5 /* */ 6 /* Copyright (c) 2002-2023 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 35 #include "scip/nlhdlr_bilinear.h" 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 */ 61 struct 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 */ 72 struct 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 */ 94 static 95 void getIneqViol( 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 */ 127 static 128 SCIP_Bool useBilinIneqs( 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 */ 163 static 164 void updateBilinearRelaxation( 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 */ 281 static 282 SCIP_Bool isPointFeasible( 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 */ 328 static 329 void getFeasiblePointsBilinear( 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 */ 406 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, boundsx) 407 || SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, boundsy) ) 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 */ 577 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, result) ) 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 */ 611 static 612 SCIP_INTERVAL intevalBilinear( 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 { 637 SCIPintervalSetEntire(SCIP_INTERVAL_INFINITY, &interval); 638 return interval; 639 } 640 641 /* x or y has empty interval -> empty */ 642 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(SCIPexprGetChildren(expr)[0])) || 643 SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(SCIPexprGetChildren(expr)[1])) ) 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 ) 679 SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &interval, interval, SCIPgetCoefExprProduct(expr)); 680 681 return interval; 682 } 683 684 /** uses inequalities for bilinear terms to get stronger bounds during reverse propagation */ 685 static 686 void reversePropBilinear( 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 { 739 SCIP_Real lbx = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]).inf; 740 SCIP_Real ubx = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]).sup; 741 SCIP_Real lby = SCIPexprGetActivity(SCIPexprGetChildren(expr)[1]).inf; 742 SCIP_Real uby = SCIPexprGetActivity(SCIPexprGetChildren(expr)[1]).sup; 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 */ 778 static 779 void computeBilinEnvelope2( 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' */ 952 static 953 SCIP_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) ); 973 SCIP_CALL( SCIPcreateExpriter(scip, &it) ); 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 987 SCIP_CALL( SCIPexpriterInit(it, SCIPgetExprNonlinear(cons), SCIP_EXPRITER_DFS, FALSE) ); 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 */ 1036 static 1037 SCIP_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 */ 1049 static 1050 SCIP_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 */ 1068 static 1069 SCIP_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 */ 1114 static 1115 SCIP_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 */ 1124 static 1125 SCIP_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 */ 1136 if( SCIPgetStage(scip) < SCIP_STAGE_INITSOLVE ) 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 */ 1230 static 1231 SCIP_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 1240 var1 = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[0]); 1241 assert(var1 != NULL); 1242 var2 = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[1]); 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 */ 1261 static 1262 SCIP_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 1290 nodeid = SCIPnodeGetNumber(SCIPgetCurrentNode(scip)); 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 */ 1309 x = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[0]); 1310 assert(x != NULL); 1311 y = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[1]); 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 */ 1323 SCIPaddBilinMcCormick(scip, SCIPgetCoefExprProduct(expr), SCIPvarGetLbLocal(x), SCIPvarGetUbLocal(x), refpointx, 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 */ 1361 updateBilinearRelaxation(scip, x, y, SCIPgetCoefExprProduct(expr), 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 { 1379 SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) ); 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 */ 1391 static 1392 SCIP_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 */ 1414 static 1415 SCIP_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 1436 SCIPintervalSetEntire(SCIP_INTERVAL_INFINITY, &intervalx); 1437 SCIPintervalSetEntire(SCIP_INTERVAL_INFINITY, &intervaly); 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", 1445 SCIPgetExprBoundsNonlinear(scip, childx).inf, SCIPgetExprBoundsNonlinear(scip, childx).sup, 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", 1455 SCIPgetExprBoundsNonlinear(scip, childy).inf, SCIPgetExprBoundsNonlinear(scip, childy).sup, 1456 intervaly.inf, intervaly.sup); 1457 SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, SCIPexprGetChildren(expr)[1], intervaly, 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 */ 1470 SCIP_RETCODE SCIPincludeNlhdlrBilinear( 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 1484 SCIP_CALL( SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, NLHDLR_NAME, NLHDLR_DESC, NLHDLR_DETECTPRIORITY, 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); 1492 SCIPnlhdlrSetSepa(nlhdlr, nlhdlrInitSepaBilinear, nlhdlrEnfoBilinear, nlhdlrEstimateBilinear, nlhdlrExitSepaBilinear); 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 */ 1517 assert(SCIPfindTable(scip, TABLE_NAME_BILINEAR) == NULL); 1518 SCIP_CALL( SCIPincludeTable(scip, TABLE_NAME_BILINEAR, TABLE_DESC_BILINEAR, FALSE, 1519 NULL, NULL, NULL, NULL, NULL, NULL, tableOutputBilinear, 1520 NULL, TABLE_POSITION_BILINEAR, TABLE_EARLIEST_STAGE_BILINEAR) ); 1521 /**! [SnippetIncludeNlhdlrBilinear] */ 1522 1523 return SCIP_OKAY; 1524 } 1525 1526 /** returns an array of expressions that have been detected by the bilinear nonlinear handler */ 1527 SCIP_EXPR** SCIPgetExprsBilinear( 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 */ 1543 int SCIPgetNExprsBilinear( 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$ */ 1559 SCIP_RETCODE SCIPaddIneqBilinear( 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 */ 1606 x = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[0]); 1607 y = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[1]); 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 */ 1708 SCIP_CALL( SCIPmarkExprPropagateNonlinear(scip, expr) ); 1709 } 1710 1711 return SCIP_OKAY; 1712 } 1713 1714 /** computes coefficients of linearization of a bilinear term in a reference point */ 1715 void SCIPaddBilinLinearization( 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 */ 1762 void SCIPaddBilinMcCormick( 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 */ 1970 void SCIPcomputeBilinEnvelope1( 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 */ 2194 void SCIPcomputeBilinEnvelope2( 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 } 2321