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_quotient.c 26 * @ingroup DEFPLUGINS_NLHDLR 27 * @brief quotient nonlinear handler 28 * @author Benjamin Mueller 29 * @author Fabian Wegscheider 30 * 31 * @todo implement INITSEPA 32 * @todo use the convex envelope for x/y described in Tawarmalani and Sahinidis (2002) if y has a finite upper bound 33 */ 34 35 #include <string.h> 36 37 #include "scip/nlhdlr_quotient.h" 38 #include "scip/cons_nonlinear.h" 39 #include "scip/pub_misc_rowprep.h" 40 #include "scip/nlhdlr.h" 41 #include "scip/nlhdlr_bilinear.h" 42 43 /* fundamental nonlinear handler properties */ 44 #define NLHDLR_NAME "quotient" 45 #define NLHDLR_DESC "nonlinear handler for quotient expressions" 46 #define NLHDLR_DETECTPRIORITY 20 47 #define NLHDLR_ENFOPRIORITY 20 48 49 /** translate from one value of infinity to another 50 * 51 * if val is ≥ infty1, then give infty2, else give val 52 */ 53 #define infty2infty(infty1, infty2, val) ((val) >= (infty1) ? (infty2) : (val)) 54 55 /*lint -e666*/ 56 57 /* 58 * Data structures 59 */ 60 61 /** nonlinear handler expression data */ 62 struct SCIP_NlhdlrExprData 63 { 64 SCIP_EXPR* numexpr; /**< expression of the numerator */ 65 SCIP_Real numcoef; /**< coefficient of the numerator */ 66 SCIP_Real numconst; /**< constant of the numerator */ 67 SCIP_EXPR* denomexpr; /**< expression of the denominator */ 68 SCIP_Real denomcoef; /**< coefficient of the denominator */ 69 SCIP_Real denomconst; /**< constant of the denominator */ 70 SCIP_Real constant; /**< constant */ 71 }; 72 73 /* 74 * Local methods 75 */ 76 77 /** helper method to create nonlinear handler expression data */ 78 static 79 SCIP_RETCODE exprdataCreate( 80 SCIP* scip, /**< SCIP data structure */ 81 SCIP_NLHDLREXPRDATA** nlhdlrexprdata, /**< nonlinear handler expression data */ 82 SCIP_EXPR* numexpr, /**< expression of the numerator */ 83 SCIP_Real numcoef, /**< coefficient of the numerator */ 84 SCIP_Real numconst, /**< constant of the numerator */ 85 SCIP_EXPR* denomexpr, /**< expression of the denominator */ 86 SCIP_Real denomcoef, /**< coefficient of the denominator */ 87 SCIP_Real denomconst, /**< constant of the denominator */ 88 SCIP_Real constant /**< constant */ 89 ) 90 { 91 assert(nlhdlrexprdata != NULL); 92 assert(numexpr != NULL); 93 assert(denomexpr != NULL); 94 assert(!SCIPisZero(scip, numcoef)); 95 assert(!SCIPisZero(scip, denomcoef)); 96 97 /* allocate memory */ 98 SCIP_CALL( SCIPallocBlockMemory(scip, nlhdlrexprdata) ); 99 100 /* store values */ 101 (*nlhdlrexprdata)->numexpr = numexpr; 102 (*nlhdlrexprdata)->numcoef = numcoef; 103 (*nlhdlrexprdata)->numconst = numconst; 104 (*nlhdlrexprdata)->denomexpr = denomexpr; 105 (*nlhdlrexprdata)->denomcoef = denomcoef; 106 (*nlhdlrexprdata)->denomconst = denomconst; 107 (*nlhdlrexprdata)->constant = constant; 108 109 /* capture expressions */ 110 SCIPcaptureExpr(numexpr); 111 SCIPcaptureExpr(denomexpr); 112 113 return SCIP_OKAY; 114 } 115 116 /** helper method to free nonlinear handler expression data */ 117 static 118 SCIP_RETCODE exprdataFree( 119 SCIP* scip, /**< SCIP data structure */ 120 SCIP_NLHDLREXPRDATA** nlhdlrexprdata /**< nonlinear handler expression data */ 121 ) 122 { 123 assert(nlhdlrexprdata != NULL); 124 assert(*nlhdlrexprdata != NULL); 125 assert((*nlhdlrexprdata)->numexpr != NULL); 126 assert((*nlhdlrexprdata)->denomexpr != NULL); 127 128 /* release expressions */ 129 SCIP_CALL( SCIPreleaseExpr(scip, &(*nlhdlrexprdata)->denomexpr) ); 130 SCIP_CALL( SCIPreleaseExpr(scip, &(*nlhdlrexprdata)->numexpr) ); 131 132 /* free expression data of nonlinear handler */ 133 SCIPfreeBlockMemory(scip, nlhdlrexprdata); 134 135 return SCIP_OKAY; 136 } 137 138 /** helper method to transform an expression g(x) into a*f(x) + b */ 139 static 140 void transformExpr( 141 SCIP* scip, /**< SCIP data structure */ 142 SCIP_EXPR* expr, /**< expression */ 143 SCIP_EXPR** target, /**< pointer to store the expression f(x) */ 144 SCIP_Real* coef, /**< pointer to store the coefficient */ 145 SCIP_Real* constant /**< pointer to store the constant */ 146 ) 147 { 148 assert(expr != NULL); 149 assert(target != NULL); 150 assert(coef != NULL); 151 assert(constant != NULL); 152 153 /* expression is a sum with one child */ 154 if( SCIPisExprSum(scip, expr) && SCIPexprGetNChildren(expr) == 1 ) 155 { 156 *target = SCIPexprGetChildren(expr)[0]; 157 *coef = SCIPgetCoefsExprSum(expr)[0]; 158 *constant = SCIPgetConstantExprSum(expr); 159 } 160 else /* otherwise return 1 * f(x) + 0 */ 161 { 162 *target = expr; 163 *coef = 1.0; 164 *constant = 0.0; 165 } 166 } 167 168 /** helper method to detect an expression of the form (a*x + b) / (c*y + d) + e 169 * 170 * Due to the expansion of products, there are two types of expressions that can be detected: 171 * 172 * 1. prod(f(x), pow(g(y),-1)) 173 * 2. sum(prod(f(x),pow(g(y),-1)), pow(g(y),-1)) 174 * 175 * @todo At the moment quotients like xy / z are not detected, because they are turned into a product expression 176 * with three children, i.e., x * y * (1 / z). 177 */ 178 static 179 SCIP_RETCODE detectExpr( 180 SCIP* scip, /**< SCIP data structure */ 181 SCIP_EXPR* expr, /**< expression */ 182 SCIP_NLHDLREXPRDATA** nlhdlrexprdata, /**< pointer to store nonlinear handler expression data */ 183 SCIP_Bool* success /**< pointer to store whether nonlinear handler should be called for this expression */ 184 ) 185 { 186 SCIP_EXPR** children; 187 SCIP_EXPR* denomexpr = NULL; 188 SCIP_EXPR* numexpr = NULL; 189 SCIP_EXPR* xexpr = NULL; 190 SCIP_EXPR* yexpr = NULL; 191 SCIP_Real a, b, c, d, e; 192 SCIP_Real nomfac = 1.0; 193 SCIP_Real numconst = 0.0; 194 195 assert(scip != NULL); 196 assert(expr != NULL); 197 198 *success = FALSE; 199 a = 0.0; 200 b = 0.0; 201 c = 0.0; 202 d = 0.0; 203 e = 0.0; 204 205 /* possible structures only have two children */ 206 if( SCIPexprGetNChildren(expr) != 2 ) 207 return SCIP_OKAY; 208 209 /* expression must be either a product or a sum */ 210 if( !SCIPisExprProduct(scip, expr) && !SCIPisExprSum(scip, expr) ) 211 return SCIP_OKAY; 212 213 children = SCIPexprGetChildren(expr); 214 assert(children != NULL); 215 216 /* case: prod(f(x), pow(g(y),-1)) */ 217 if( SCIPisExprProduct(scip, expr) ) 218 { 219 if( SCIPisExprPower(scip, children[0]) && SCIPgetExponentExprPow(children[0]) == -1.0 ) /*lint !e777*/ 220 { 221 denomexpr = SCIPexprGetChildren(children[0])[0]; 222 numexpr = children[1]; 223 } 224 else if( SCIPisExprPower(scip, children[1]) && SCIPgetExponentExprPow(children[1]) == -1.0 ) /*lint !e777*/ 225 { 226 denomexpr = SCIPexprGetChildren(children[1])[0]; 227 numexpr = children[0]; 228 } 229 230 /* remember to scale the numerator by the coefficient stored in the product expression */ 231 nomfac = SCIPgetCoefExprProduct(expr); 232 } 233 /* case: sum(prod(f(x),pow(g(y),-1)), pow(g(y),-1)) */ 234 else 235 { 236 SCIP_Real* sumcoefs; 237 238 assert(SCIPisExprSum(scip, expr)); 239 sumcoefs = SCIPgetCoefsExprSum(expr); 240 241 /* children[0] is 1/g(y) and children[1] is a product of f(x) and 1/g(y) */ 242 if( SCIPisExprPower(scip, children[0]) && SCIPgetExponentExprPow(children[0]) == -1.0 243 && SCIPisExprProduct(scip, children[1]) && SCIPexprGetNChildren(children[1]) == 2 ) /* lint !e777 */ 244 { 245 SCIP_Real prodcoef = SCIPgetCoefExprProduct(children[1]); 246 247 if( children[0] == SCIPexprGetChildren(children[1])[0] ) 248 { 249 denomexpr = SCIPexprGetChildren(children[0])[0]; 250 numexpr = SCIPexprGetChildren(children[1])[1]; 251 } 252 else if( children[0] == SCIPexprGetChildren(children[1])[1] ) 253 { 254 denomexpr = SCIPexprGetChildren(children[0])[0]; 255 numexpr = SCIPexprGetChildren(children[1])[0]; 256 } 257 258 /* remember scalar and constant for numerator */ 259 nomfac = sumcoefs[1] * prodcoef; 260 numconst = sumcoefs[0]; 261 } 262 /* children[1] is 1/g(y) and children[0] is a product of f(x) and 1/g(y) */ 263 else if( SCIPisExprPower(scip, children[1]) && SCIPgetExponentExprPow(children[1]) == -1.0 264 && SCIPisExprProduct(scip, children[0]) && SCIPexprGetNChildren(children[0]) == 2 ) /* lint !e777 */ 265 { 266 SCIP_Real prodcoef = SCIPgetCoefExprProduct(children[0]); 267 268 if( children[1] == SCIPexprGetChildren(children[0])[0] ) 269 { 270 denomexpr = SCIPexprGetChildren(children[1])[0]; 271 numexpr = SCIPexprGetChildren(children[0])[1]; 272 } 273 else if( children[1] == SCIPexprGetChildren(children[0])[1] ) 274 { 275 denomexpr = SCIPexprGetChildren(children[1])[0]; 276 numexpr = SCIPexprGetChildren(children[0])[0]; 277 } 278 279 /* remember scalar and constant for numerator */ 280 nomfac = sumcoefs[0] * prodcoef; 281 numconst = sumcoefs[1]; 282 } 283 284 /* remember the constant of the sum expression */ 285 e = SCIPgetConstantExprSum(expr); 286 } 287 288 if( denomexpr != NULL && numexpr != NULL ) 289 { 290 /* transform numerator and denominator to detect structures like (a * f(x) + b) / (c * f(x) + d) */ 291 transformExpr(scip, numexpr, &xexpr, &a, &b); 292 transformExpr(scip, denomexpr, &yexpr, &c, &d); 293 294 SCIPdebugMsg(scip, "detected numerator (%g * %p + %g) and denominator (%g * %p + %g)\n", a, (void*)xexpr, b, 295 c, (void*)yexpr, d); 296 297 /* detection is only be successful if the expression of the numerator an denominator are the same 298 * (so boundtightening can be stronger than default) or we are going to provide estimators (there will be an auxvar) 299 */ 300 *success = (xexpr == yexpr) || (SCIPgetExprNAuxvarUsesNonlinear(expr) > 0); 301 302 #ifdef SCIP_DEBUG 303 SCIPinfoMessage(scip, NULL, "Expression for numerator: "); 304 SCIP_CALL( SCIPprintExpr(scip, xexpr, NULL) ); 305 SCIPinfoMessage(scip, NULL, "\nExpression for denominator: "); 306 SCIP_CALL( SCIPprintExpr(scip, yexpr, NULL) ); 307 SCIPinfoMessage(scip, NULL, "\n"); 308 #endif 309 } 310 311 /* register usage of xexpr and yexpr 312 * create nonlinear handler expression data 313 */ 314 if( *success ) 315 { 316 assert(xexpr != NULL); 317 assert(xexpr != NULL); 318 assert(a != 0.0); 319 assert(c != 0.0); 320 321 assert(SCIPgetExprNAuxvarUsesNonlinear(expr) > 0 || xexpr == yexpr); 322 323 /* request auxiliary variables for xexpr and yexpr if we will estimate 324 * mark that the bounds of the expression are important to construct the estimators 325 * (TODO check the curvature of the univariate quotient, as bounds may actually not be used) 326 * if univariate, then we also do inteval and reverseprop, so mark that the activities will be used for inteval 327 */ 328 SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, xexpr, 329 SCIPgetExprNAuxvarUsesNonlinear(expr) > 0, 330 xexpr == yexpr, 331 SCIPgetExprNAuxvarUsesNonlinear(expr) > 0, 332 SCIPgetExprNAuxvarUsesNonlinear(expr) > 0) ); 333 334 if( xexpr != yexpr && SCIPgetExprNAuxvarUsesNonlinear(expr) > 0 ) 335 { 336 SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, yexpr, TRUE, FALSE, TRUE, TRUE) ); 337 } 338 339 a = nomfac * a; 340 b = nomfac * b + numconst; 341 342 SCIPdebug( SCIP_CALL( SCIPprintExpr(scip, expr, NULL) ); ) 343 SCIPdebug( SCIPinfoMessage(scip, NULL, "\n") ); 344 SCIPdebugMsg(scip, "detected quotient expression (%g * %p + %g) / (%g * %p + %g) + %g\n", a, (void*)xexpr, 345 b, c, (void*)yexpr, d, e); 346 SCIP_CALL( exprdataCreate(scip, nlhdlrexprdata, xexpr, a, b, yexpr, c, d, e) ); 347 } 348 349 return SCIP_OKAY; 350 } 351 352 /** helper method to compute interval for (a x + b) / (c x + d) + e */ 353 static 354 SCIP_INTERVAL intEvalQuotient( 355 SCIP* scip, /**< SCIP data structure */ 356 SCIP_INTERVAL bnds, /**< bounds on x */ 357 SCIP_Real a, /**< coefficient in numerator */ 358 SCIP_Real b, /**< constant in numerator */ 359 SCIP_Real c, /**< coefficient in denominator */ 360 SCIP_Real d, /**< constant in denominator */ 361 SCIP_Real e /**< constant */ 362 ) 363 { 364 SCIP_INTERVAL result; 365 SCIP_INTERVAL denominterval; 366 SCIP_INTERVAL numinterval; 367 int i; 368 369 assert(scip != NULL); 370 371 /* return empty interval if the domain of x is empty */ 372 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, bnds) ) 373 { 374 SCIPintervalSetEmpty(&result); 375 return result; 376 } 377 378 /* compute bounds for denominator */ 379 SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &denominterval, bnds, c); 380 SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &denominterval, denominterval, d); 381 382 /* there is no useful interval if 0 is in the interior of the interval of the denominator */ 383 if( SCIPintervalGetInf(denominterval) < 0.0 && SCIPintervalGetSup(denominterval) > 0.0 ) 384 { 385 SCIPintervalSetEntire(SCIP_INTERVAL_INFINITY, &result); 386 return result; 387 } 388 389 /* a d = b c implies that f(x) = b / d + e, i.e., f is constant */ 390 if( a*d - b*c == 0.0 ) 391 { 392 SCIPintervalSet(&result, b / d + e); 393 return result; 394 } 395 396 /* 397 * evaluate for [x.inf,x.inf] and [x.sup,x.sup] independently 398 */ 399 SCIPintervalSetEmpty(&result); 400 401 for( i = 0; i < 2; ++i ) 402 { 403 SCIP_INTERVAL quotinterval; 404 SCIP_Real val = (i == 0) ? bnds.inf : bnds.sup; 405 406 /* set the resulting interval to a / c if the bounds is infinite */ 407 if( SCIPisInfinity(scip, REALABS(val)) ) 408 { 409 SCIPintervalSet("interval, a); 410 SCIPintervalDivScalar(SCIP_INTERVAL_INFINITY, "interval, quotinterval, c); 411 } 412 else 413 { 414 /* a x' + b */ 415 SCIPintervalSet(&numinterval, val); 416 SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &numinterval, numinterval, a); 417 SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &numinterval, numinterval, b); 418 419 /* c x' + d */ 420 SCIPintervalSet(&denominterval, val); 421 SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &denominterval, denominterval, c); 422 SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &denominterval, denominterval, d); 423 424 /* (a x' + b) / (c x' + d) + e */ 425 SCIPintervalDiv(SCIP_INTERVAL_INFINITY, "interval, numinterval, denominterval); 426 SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, "interval, quotinterval, e); 427 } 428 429 /* unify with the resulting interval */ 430 SCIPintervalUnify(&result, result, quotinterval); 431 } 432 433 return result; 434 } 435 436 /** helper method to compute reverse propagation for (a x + b) / (c x + d) + e */ 437 static 438 SCIP_INTERVAL reversepropQuotient( 439 SCIP_INTERVAL bnds, /**< bounds on (a x + b) / (c x + d) + e */ 440 SCIP_Real a, /**< coefficient in numerator */ 441 SCIP_Real b, /**< constant in numerator */ 442 SCIP_Real c, /**< coefficient in denominator */ 443 SCIP_Real d, /**< constant in denominator */ 444 SCIP_Real e /**< constant */ 445 ) 446 { 447 SCIP_INTERVAL result; 448 int i; 449 450 SCIPintervalSetEmpty(&result); 451 452 /* return empty interval if the domain of the expression is empty */ 453 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, bnds) ) 454 return result; 455 456 /* substract constant from bounds of the expression */ 457 SCIPintervalSubScalar(SCIP_INTERVAL_INFINITY, &bnds, bnds, e); 458 459 /* if the expression is constant or the limit lies inside the domain, nothing can be propagated */ 460 if( a*d - b*c == 0.0 || (bnds.inf < a / c && bnds.sup > a / c) ) 461 { 462 SCIPintervalSetEntire(SCIP_INTERVAL_INFINITY, &result); 463 return result; 464 } 465 466 /* compute bounds for [x.inf,x.inf] and [x.sup,x.sup] independently */ 467 for( i = 0; i < 2; ++i ) 468 { 469 SCIP_INTERVAL denominator; 470 SCIP_INTERVAL numerator; 471 SCIP_INTERVAL quotient; 472 SCIP_Real val = (i == 0) ? bnds.inf : bnds.sup; 473 474 /* (d * x' - b) */ 475 SCIPintervalSet(&numerator, d); 476 SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &numerator, numerator, val); 477 SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &numerator, numerator, -b); 478 479 /* (a - c * x') */ 480 SCIPintervalSet(&denominator, -c); 481 SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &denominator, denominator, val); 482 SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &denominator, denominator, a); 483 484 /* (d * x' - b) / (a - c * x') */ 485 SCIPintervalDiv(SCIP_INTERVAL_INFINITY, "ient, numerator, denominator); 486 487 /* unify with the resulting interval */ 488 SCIPintervalUnify(&result, result, quotient); 489 } 490 491 return result; 492 } 493 494 /** adds data to given rowprep; the generated estimator is always locally valid 495 * 496 * @note the constant is moved to the left- or right-hand side 497 * @note other than the name of this function may indicate, it does not create a rowprep 498 */ 499 static 500 SCIP_RETCODE createRowprep( 501 SCIP* scip, /**< SCIP data structure */ 502 SCIP_ROWPREP* rowprep, /**< a rowprep where to store the estimator */ 503 SCIP_VAR** vars, /**< variables */ 504 SCIP_Real* coefs, /**< coefficients */ 505 SCIP_Real constant, /**< constant */ 506 int nlinvars /**< total number of variables */ 507 ) 508 { 509 assert(scip != NULL); 510 assert(rowprep != NULL); 511 assert(coefs != NULL); 512 assert(vars != NULL); 513 514 /* create rowprep */ 515 SCIProwprepAddSide(rowprep, -constant); 516 SCIP_CALL( SCIPensureRowprepSize(scip, rowprep, nlinvars + 1) ); 517 518 /* add coefficients */ 519 SCIP_CALL( SCIPaddRowprepTerms(scip, rowprep, nlinvars, vars, coefs) ); 520 521 return SCIP_OKAY; 522 } 523 524 /** computes an estimator at a given point for the univariate case (ax + b) / (cx + d) + e 525 * 526 * Depending on the reference point, the estimator is a tangent or a secant on the graph. 527 * It depends on whether we are under- or overestimating, whether we are on the left or 528 * on the right side of the singularity at -d/c, and whether it is the monotone increasing 529 * (ad - bc > 0) or decreasing part (ad - bc < 0). Together, there are 8 cases: 530 * 531 * - mon. incr. + overestimate + left hand side --> secant 532 * - mon. incr. + overestimate + right hand side --> tangent 533 * - mon. incr. + understimate + left hand side --> tangent 534 * - mon. incr. + understimate + right hand side --> secant 535 * - mon. decr. + overestimate + left hand side --> tangent 536 * - mon. decr. + overestimate + right hand side --> secant 537 * - mon. decr. + understimate + left hand side --> secant 538 * - mon. decr. + understimate + right hand side --> tangent 539 */ 540 static 541 SCIP_RETCODE estimateUnivariate( 542 SCIP* scip, /**< SCIP data structure */ 543 SCIP_Real lbx, /**< local lower bound of x */ 544 SCIP_Real ubx, /**< local upper bound of x */ 545 SCIP_Real gllbx, /**< global lower bound of x */ 546 SCIP_Real glubx, /**< global upper bound of x */ 547 SCIP_Real solx, /**< solution value of x */ 548 SCIP_Real a, /**< coefficient in numerator */ 549 SCIP_Real b, /**< constant in numerator */ 550 SCIP_Real c, /**< coefficient in denominator */ 551 SCIP_Real d, /**< constant in denominator */ 552 SCIP_Real e, /**< constant */ 553 SCIP_Real* coef, /**< pointer to store the coefficient */ 554 SCIP_Real* constant, /**< pointer to store the constant */ 555 SCIP_Bool overestimate, /**< whether the expression should be overestimated */ 556 SCIP_Bool* local, /**< pointer to store whether the estimate is locally valid */ 557 SCIP_Bool* branchinguseful, /**< pointer to store whether branching on the expression would improve the estimator */ 558 SCIP_Bool* success /**< buffer to store whether separation was successful */ 559 ) 560 { 561 SCIP_Real singularity; 562 SCIP_Bool isinleftpart; 563 SCIP_Bool monincreasing; 564 565 assert(lbx <= solx && solx <= ubx); 566 assert(coef != NULL); 567 assert(constant != NULL); 568 assert(local != NULL); 569 assert(branchinguseful != NULL); 570 assert(success != NULL); 571 572 *branchinguseful = TRUE; 573 *success = FALSE; 574 *coef = 0.0; 575 *constant = 0.0; 576 singularity = -d / c; 577 578 /* estimate is globally valid if local and global bounds are equal */ 579 *local = gllbx != lbx || glubx != ubx; /*lint !e777*/ 580 581 /* if 0 is in the denom interval, estimation is not possible */ 582 if( SCIPisLE(scip, lbx, singularity) && SCIPisGE(scip, ubx, singularity) ) 583 return SCIP_OKAY; 584 585 isinleftpart = (ubx < singularity); 586 monincreasing = (a * d - b * c > 0.0); 587 588 /* this encodes the 8 cases explained above */ 589 if( monincreasing == (overestimate == isinleftpart) ) 590 { 591 SCIP_Real lbeval; 592 SCIP_Real ubeval; 593 594 /* if one of the bounds is infinite, secant cannot be computed */ 595 if( SCIPisInfinity(scip, -lbx) || SCIPisInfinity(scip, ubx) ) 596 return SCIP_OKAY; 597 598 lbeval = (a * lbx + b) / (c * lbx + d) + e; 599 ubeval = (a * ubx + b) / (c * ubx + d) + e; 600 601 /* compute coefficient and constant of linear estimator */ 602 *coef = (ubeval - lbeval) / (ubx - lbx); 603 *constant = ubeval - (*coef) * ubx; 604 } 605 else 606 { 607 SCIP_Real soleval; 608 609 soleval = (a * solx + b) / (c * solx + d) + e; 610 611 /* compute coefficient and constant of linear estimator */ 612 *coef = (a * d - b * c) / SQR(d + c * solx); 613 *constant = soleval - (*coef) * solx; 614 615 /* gradient cuts are globally valid if the singularity is not in [gllbx,glubx] */ 616 *local = SCIPisLE(scip, gllbx, singularity) && SCIPisGE(scip, glubx, singularity); 617 618 /* branching will not improve the convexification via tangent cuts */ 619 *branchinguseful = FALSE; 620 } 621 622 /* avoid huge values in the cut */ 623 if( SCIPisHugeValue(scip, REALABS(*coef)) || SCIPisHugeValue(scip, REALABS(*constant)) ) 624 return SCIP_OKAY; 625 626 *success = TRUE; 627 628 return SCIP_OKAY; 629 } 630 631 /** helper method to compute estimator for the univariate case; the estimator is stored in a given rowprep */ 632 static 633 SCIP_RETCODE estimateUnivariateQuotient( 634 SCIP* scip, /**< SCIP data structure */ 635 SCIP_SOL* sol, /**< solution point (or NULL for the LP solution) */ 636 SCIP_EXPR* xexpr, /**< argument expression */ 637 SCIP_Real a, /**< coefficient in numerator */ 638 SCIP_Real b, /**< constant in numerator */ 639 SCIP_Real c, /**< coefficient in denominator */ 640 SCIP_Real d, /**< constant in denominator */ 641 SCIP_Real e, /**< constant */ 642 SCIP_Bool overestimate, /**< whether the expression should be overestimated */ 643 SCIP_ROWPREP* rowprep, /**< a rowprep where to store the estimator */ 644 SCIP_Bool* branchinguseful, /**< pointer to store whether branching on the expression would improve the estimator */ 645 SCIP_Bool* success /**< buffer to store whether separation was successful */ 646 ) 647 { 648 SCIP_VAR* x; 649 SCIP_Real constant; 650 SCIP_Real coef; 651 SCIP_Real gllbx; 652 SCIP_Real glubx; 653 SCIP_Real lbx; 654 SCIP_Real ubx; 655 SCIP_Real solx; 656 SCIP_Bool local; 657 SCIP_INTERVAL bnd; 658 659 assert(rowprep != NULL); 660 assert(branchinguseful != NULL); 661 assert(success != NULL); 662 663 x = SCIPgetExprAuxVarNonlinear(xexpr); 664 665 /* get local bounds on xexpr */ 666 SCIPintervalSetBounds(&bnd, 667 -infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, -SCIPvarGetLbGlobal(x)), 668 infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, SCIPvarGetUbGlobal(x))); 669 SCIP_CALL( SCIPevalExprActivity(scip, xexpr) ); 670 SCIPintervalIntersectEps(&bnd, SCIPepsilon(scip), bnd, SCIPexprGetActivity(xexpr)); 671 lbx = bnd.inf; 672 ubx = bnd.sup; 673 674 /* check whether variable has been fixed or has empty interval */ 675 if( SCIPisEQ(scip, lbx, ubx) || ubx < lbx ) 676 { 677 *success = FALSE; 678 return SCIP_OKAY; 679 } 680 681 /* get global variable bounds */ 682 gllbx = SCIPvarGetLbGlobal(x); 683 glubx = SCIPvarGetUbGlobal(x); 684 685 /* get and adjust solution value */ 686 solx = SCIPgetSolVal(scip, sol, x); 687 solx = MIN(MAX(solx, lbx), ubx); 688 689 /* compute an estimator */ 690 SCIP_CALL( estimateUnivariate(scip, lbx, ubx, gllbx, glubx, solx, a, b, c, d, e, &coef, &constant, overestimate, &local, branchinguseful, success) ); 691 692 /* add estimator to rowprep, if successful */ 693 if( *success ) 694 { 695 (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "quot_%s_%lld", SCIPvarGetName(x), SCIPgetNLPs(scip)); 696 SCIP_CALL( createRowprep(scip, rowprep, &x, &coef, constant, 1) ); 697 SCIProwprepSetLocal(rowprep, local); 698 } 699 700 return SCIP_OKAY; 701 } 702 703 /** helper method to compute a gradient cut for 704 * \f[ 705 * h^c(x,y) := \frac{1}{y} \left(\frac{x + \sqrt{\text{lbx}\cdot\text{ubx}}}{\sqrt{\text{lbx}} + \sqrt{\text{ubx}}}\right)^2 706 * \f] 707 * at a given reference point 708 * 709 * See Zamora and Grossmann (1988) for more details. 710 */ 711 static 712 void hcGradCut( 713 SCIP_Real lbx, /**< lower bound of x */ 714 SCIP_Real ubx, /**< upper bound of x */ 715 SCIP_Real solx, /**< solution value of x */ 716 SCIP_Real soly, /**< solution value of y */ 717 SCIP_Real* coefx, /**< pointer to store the coefficient of x */ 718 SCIP_Real* coefy, /**< pointer to store the coefficient of y */ 719 SCIP_Real* constant /**< pointer to store the constant */ 720 ) 721 { 722 SCIP_Real tmp1; 723 SCIP_Real tmp2; 724 725 assert(lbx >= 0.0); 726 assert(lbx <= ubx); 727 assert(soly > 0.0); 728 assert(coefx != NULL); 729 assert(coefy != NULL); 730 assert(constant != NULL); 731 732 tmp1 = SQRT(lbx * ubx) + solx; 733 tmp2 = SQR(SQRT(lbx) + SQRT(ubx)) * soly; /*lint !e666*/ 734 assert(tmp2 > 0.0); 735 736 *coefx = 2.0 * tmp1 / tmp2; 737 *coefy = -SQR(tmp1) / (tmp2 * soly); 738 *constant = 2.0 * SQRT(lbx * ubx) * tmp1 / tmp2; 739 } 740 741 /** computes an over- or underestimator at a given point for the bivariate case x/y ≤/≥ z 742 * 743 * There are the following cases for y > 0: 744 * 745 * 1. lbx < 0 < ubx: 746 * Rewrite x / y = z as x = y * z and use McCormick to compute a valid inequality of the form 747 * x = y * z ≤ a * y + b * z + c. Note that b > 0 because of y > 0. The inequality is then transformed 748 * to x / b - a/b * y - c/b ≤ z, which results in a valid underestimator for x / y over the set 749 * {(x,y) | lbz ≤ x / y ≤ ubz}. Note that overestimating/underestimating the bilinear term with McCormick 750 * results in an underestimator/overestimator for x / y. 751 * 752 * 2. lbx ≥ 0 or ubx ≤ 0: 753 * - overestimation: use \f$z \leq \frac{1}{\text{lby}\cdot\text{uby}} \min(\text{uby}\cdot x - \text{lbx}\cdot y + \text{lbx}\cdot\text{lby}, \text{lby}\cdot x - \text{ubx}\cdot y + \text{ubx}\cdot\text{uby})\f$ 754 * - underestimation: use \f$z \geq x/y \geq \frac{1}{y} \frac{x + \sqrt{\text{lbx}\cdot\text{ubx}}}{\sqrt{\text{lbx} + \sqrt{\text{ubx}}}}\f$ and build gradient cut 755 * 756 * If y < 0, swap and negate its bounds and compute the respective opposite estimator (and negate it). 757 * 758 * If 0 is in the interval of y, nothing is possible. 759 */ 760 static 761 SCIP_RETCODE estimateBivariate( 762 SCIP* scip, /**< SCIP data structure */ 763 SCIP_Real lbx, /**< lower bound of x */ 764 SCIP_Real ubx, /**< upper bound of x */ 765 SCIP_Real lby, /**< lower bound of y */ 766 SCIP_Real uby, /**< upper bound of y */ 767 SCIP_Real lbz, /**< lower bound of z */ 768 SCIP_Real ubz, /**< lower bound of z */ 769 SCIP_Real solx, /**< reference point for x */ 770 SCIP_Real soly, /**< reference point for y */ 771 SCIP_Real solz, /**< reference point for z */ 772 SCIP_Bool overestimate, /**< whether the expression should be overestimated */ 773 SCIP_Real* coefx, /**< pointer to store the x coefficient */ 774 SCIP_Real* coefy, /**< pointer to store the y coefficient */ 775 SCIP_Real* constant, /**< pointer to store the constant */ 776 SCIP_Bool* branchingusefulx, /**< pointer to store whether branching on x would improve the estimator */ 777 SCIP_Bool* branchingusefuly, /**< pointer to store whether branching on y would improve the estimator */ 778 SCIP_Bool* success /**< buffer to store whether computing the estimator was successful */ 779 ) 780 { 781 SCIP_Bool negatedx = FALSE; 782 SCIP_Bool negatedy = FALSE; 783 784 assert(lbx <= solx && solx <= ubx); 785 assert(lby <= soly && soly <= uby); 786 assert(lbz <= solz && solz <= ubz); 787 assert(coefx != NULL); 788 assert(coefy != NULL); 789 assert(constant != NULL); 790 assert(branchingusefulx != NULL); 791 assert(branchingusefuly != NULL); 792 assert(success != NULL); 793 794 *branchingusefulx = TRUE; 795 *branchingusefuly = TRUE; 796 *success = TRUE; 797 *coefx = 0.0; 798 *coefy = 0.0; 799 *constant = 0.0; 800 801 /* if 0 is in [lby,uby], then it is not possible to compute an estimator */ 802 if( SCIPisLE(scip, lby, 0.0) && SCIPisGE(scip, uby, 0.0) ) 803 { 804 *success = FALSE; 805 return SCIP_OKAY; 806 } 807 808 /* negate bounds of y if it is not positive */ 809 if( uby < 0.0 ) 810 { 811 SCIP_Real tmp = uby; 812 813 uby = -lby; 814 lby = -tmp; 815 soly = -soly; 816 negatedy = TRUE; 817 overestimate = !overestimate; 818 } 819 820 /* case 1: 0 is in the interior of [lbx,ubx] */ 821 if( lbx < 0.0 && 0.0 < ubx ) 822 { 823 SCIP_Real mccoefy = 0.0; 824 SCIP_Real mccoefaux = 0.0; 825 SCIP_Real mcconst = 0.0; 826 827 /* as explained in the description of this method, overestimating/underestimating the bilinear term results in an 828 * underestimator/overestimator for x / y 829 */ 830 SCIPaddBilinMcCormick(scip, 1.0, lbz, ubz, solz, lby, uby, soly, !overestimate, &mccoefaux, &mccoefy, &mcconst, 831 success); 832 assert(mccoefaux >= 0.0); 833 834 if( !(*success) ) 835 return SCIP_OKAY; 836 837 /* resulting estimator is x/b - a/b * y - c/b, where a*y + b*z + c is the estimator for y*z */ 838 *coefx = 1.0 / mccoefaux; 839 *coefy = -mccoefy / mccoefaux; 840 *constant = -mcconst / mccoefaux; 841 } 842 /* case 2: 0 is not in the interior of [lbx,ubx] */ 843 else 844 { 845 /* negate bounds of x if it is negative */ 846 if( ubx <= 0.0 ) 847 { 848 SCIP_Real tmp = ubx; 849 850 ubx = -lbx; 851 lbx = -tmp; 852 solx = -solx; 853 negatedx = TRUE; 854 overestimate = !overestimate; 855 } 856 857 /* case 2a */ 858 if( overestimate ) 859 { 860 /* check where the minimum is attained */ 861 if( uby * solx - lbx * soly + lbx * lby <= lby * solx - ubx * soly + ubx * uby ) 862 { 863 *coefx = 1.0 / lby; 864 *coefy = -lbx / (lby * uby); 865 *constant = lbx / uby; 866 } 867 else 868 { 869 *coefx = 1.0 / uby; 870 *coefy = -ubx / (lby * uby); 871 *constant = ubx / lby; 872 } 873 } 874 /* case 2b */ 875 else 876 { 877 /* compute gradient cut for h^c(x,y) at (solx,soly) */ 878 hcGradCut(lbx, ubx, solx, soly, coefx, coefy, constant); 879 880 /* estimator is independent of the bounds of y */ 881 *branchingusefuly = FALSE; 882 } 883 } 884 885 /* reverse negations of x and y in the resulting estimator */ 886 if( negatedx ) 887 *coefx = -(*coefx); 888 if( negatedy ) 889 *coefy = -(*coefy); 890 891 /* if exactly one variable has been negated, then we have computed an underestimate/overestimate for the negated 892 * expression, which results in an overestimate/underestimate for the original expression 893 */ 894 if( negatedx != negatedy ) 895 { 896 *coefx = -(*coefx); 897 *coefy = -(*coefy); 898 *constant = -(*constant); 899 } 900 901 /* avoid huge values in the estimator */ 902 if( SCIPisHugeValue(scip, REALABS(*coefx)) || SCIPisHugeValue(scip, REALABS(*coefy)) 903 || SCIPisHugeValue(scip, REALABS(*constant)) ) 904 { 905 *success = FALSE; 906 return SCIP_OKAY; 907 } 908 909 return SCIP_OKAY; 910 } 911 912 /** construct an estimator for a quotient expression of the form (ax + b) / (cy + d) + e 913 * 914 * The resulting estimator is stored in a rowprep. 915 * 916 * The method first computes an estimator for x' / y' with x := ax + b and y := cy + d 917 * and then transforms this estimator to one for the quotient (ax + b) / (cy + d) + e. 918 */ 919 static 920 SCIP_RETCODE estimateBivariateQuotient( 921 SCIP* scip, /**< SCIP data structure */ 922 SCIP_EXPR* xexpr, /**< numerator expression */ 923 SCIP_EXPR* yexpr, /**< denominator expression */ 924 SCIP_VAR* auxvar, /**< auxiliary variable */ 925 SCIP_SOL* sol, /**< solution point (or NULL for the LP solution) */ 926 SCIP_Real a, /**< coefficient of numerator */ 927 SCIP_Real b, /**< constant of numerator */ 928 SCIP_Real c, /**< coefficient of denominator */ 929 SCIP_Real d, /**< constant of denominator */ 930 SCIP_Real e, /**< constant term */ 931 SCIP_Bool overestimate, /**< whether the expression should be overestimated */ 932 SCIP_ROWPREP* rowprep, /**< a rowprep where to store the estimator */ 933 SCIP_Bool* branchingusefulx, /**< pointer to store whether branching on x would improve the estimator */ 934 SCIP_Bool* branchingusefuly, /**< pointer to store whether branching on y would improve the estimator */ 935 SCIP_Bool* success /**< buffer to store whether separation was successful */ 936 ) 937 { 938 SCIP_VAR* vars[2]; 939 SCIP_Real coefs[2] = {0.0, 0.0}; 940 SCIP_Real constant = 0.0; 941 SCIP_Real solx; 942 SCIP_Real soly; 943 SCIP_Real solz; 944 SCIP_Real lbx; 945 SCIP_Real ubx; 946 SCIP_Real lby; 947 SCIP_Real uby; 948 SCIP_Real lbz; 949 SCIP_Real ubz; 950 SCIP_INTERVAL bnd; 951 952 assert(xexpr != NULL); 953 assert(yexpr != NULL); 954 assert(xexpr != yexpr); 955 assert(auxvar != NULL); 956 assert(rowprep != NULL); 957 assert(branchingusefulx != NULL); 958 assert(branchingusefuly != NULL); 959 assert(success != NULL); 960 961 vars[0] = SCIPgetExprAuxVarNonlinear(xexpr); 962 vars[1] = SCIPgetExprAuxVarNonlinear(yexpr); 963 964 /* get bounds for x, y, and z */ 965 SCIPintervalSetBounds(&bnd, 966 -infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, -SCIPvarGetLbGlobal(vars[0])), 967 infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, SCIPvarGetUbGlobal(vars[0]))); 968 SCIP_CALL( SCIPevalExprActivity(scip, xexpr) ); 969 SCIPintervalIntersectEps(&bnd, SCIPepsilon(scip), bnd, SCIPexprGetActivity(xexpr)); 970 lbx = bnd.inf; 971 ubx = bnd.sup; 972 973 SCIPintervalSetBounds(&bnd, 974 -infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, -SCIPvarGetLbGlobal(vars[1])), 975 infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, SCIPvarGetUbGlobal(vars[1]))); 976 SCIP_CALL( SCIPevalExprActivity(scip, yexpr) ); 977 SCIPintervalIntersectEps(&bnd, SCIPepsilon(scip), bnd, SCIPexprGetActivity(yexpr)); 978 lby = bnd.inf; 979 uby = bnd.sup; 980 981 lbz = SCIPvarGetLbLocal(auxvar); 982 ubz = SCIPvarGetUbLocal(auxvar); 983 984 /* check whether one of the variables has been fixed or has empty domain */ 985 if( SCIPisEQ(scip, lbx, ubx) || SCIPisEQ(scip, lby, uby) || ubx < lbx || uby < lby ) 986 { 987 *success = FALSE; 988 return SCIP_OKAY; 989 } 990 991 /* get and adjust solution values */ 992 solx = SCIPgetSolVal(scip, sol, vars[0]); 993 soly = SCIPgetSolVal(scip, sol, vars[1]); 994 solz = SCIPgetSolVal(scip, sol, auxvar); 995 solx = MIN(MAX(solx, lbx), ubx); 996 soly = MIN(MAX(soly, lby), uby); 997 solz = MIN(MAX(solz, lbz), ubz); 998 999 /* compute an estimator */ 1000 SCIP_CALL( estimateBivariate(scip, 1001 MIN(a * lbx, a * ubx) + b, MAX(a * lbx, a * ubx) + b, /* bounds of x' */ 1002 MIN(c * lby, c * uby) + d, MAX(c * lby, c * uby) + d, /* bounds of y' */ 1003 lbz, ubz, a * solx + b, c * soly + d, solz, overestimate, &coefs[0], &coefs[1], &constant, 1004 branchingusefulx, branchingusefuly, success) ); 1005 1006 /* add estimator to rowprep, if successful */ 1007 if( *success ) 1008 { 1009 /* transform estimator Ax' + By'+ C = A(ax + b) + B (cy + d) + C = (Aa) x + (Bc) y + (C + Ab + Bd); 1010 * add the constant e separately 1011 */ 1012 constant += coefs[0] * b + coefs[1] * d + e; 1013 coefs[0] *= a; 1014 coefs[1] *= c; 1015 1016 /* prepare rowprep */ 1017 (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "quot_%s_%s_%lld", SCIPvarGetName(vars[0]), SCIPvarGetName(vars[1]), 1018 SCIPgetNLPs(scip)); 1019 SCIP_CALL( createRowprep(scip, rowprep, vars, coefs, constant, 2) ); 1020 } 1021 1022 return SCIP_OKAY; 1023 } 1024 1025 /* 1026 * Callback methods of nonlinear handler 1027 */ 1028 1029 /** nonlinear handler copy callback */ 1030 static 1031 SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrQuotient) 1032 { /*lint --e{715}*/ 1033 assert(targetscip != NULL); 1034 assert(sourcenlhdlr != NULL); 1035 assert(strcmp(SCIPnlhdlrGetName(sourcenlhdlr), NLHDLR_NAME) == 0); 1036 1037 SCIP_CALL( SCIPincludeNlhdlrQuotient(targetscip) ); 1038 1039 return SCIP_OKAY; 1040 } 1041 1042 1043 /** callback to free expression specific data */ 1044 static 1045 SCIP_DECL_NLHDLRFREEEXPRDATA(nlhdlrFreeExprDataQuotient) 1046 { /*lint --e{715}*/ 1047 assert(nlhdlrexprdata != NULL); 1048 assert(*nlhdlrexprdata != NULL); 1049 1050 /* free expression data of nonlinear handler */ 1051 SCIP_CALL( exprdataFree(scip, nlhdlrexprdata) ); 1052 1053 return SCIP_OKAY; 1054 } 1055 1056 1057 /** callback to detect structure in expression tree */ 1058 static 1059 SCIP_DECL_NLHDLRDETECT(nlhdlrDetectQuotient) 1060 { /*lint --e{715}*/ 1061 SCIP_Bool success; 1062 1063 assert(nlhdlrexprdata != NULL); 1064 1065 /* call detection routine */ 1066 SCIP_CALL( detectExpr(scip, expr, nlhdlrexprdata, &success) ); 1067 1068 if( success ) 1069 { 1070 if( SCIPgetExprNAuxvarUsesNonlinear(expr) > 0 ) 1071 *participating = SCIP_NLHDLR_METHOD_SEPABOTH; 1072 1073 if( (*nlhdlrexprdata)->numexpr == (*nlhdlrexprdata)->denomexpr ) 1074 { 1075 /* if univariate, then we also do inteval and reverseprop */ 1076 *participating |= SCIP_NLHDLR_METHOD_ACTIVITY; 1077 1078 /* if univariate, then all our methods are enforcing */ 1079 *enforcing |= *participating; 1080 } 1081 } 1082 1083 return SCIP_OKAY; 1084 } 1085 1086 1087 /** auxiliary evaluation callback of nonlinear handler */ 1088 static 1089 SCIP_DECL_NLHDLREVALAUX(nlhdlrEvalauxQuotient) 1090 { /*lint --e{715}*/ 1091 SCIP_VAR* auxvarx; 1092 SCIP_VAR* auxvary; 1093 SCIP_Real solvalx; 1094 SCIP_Real solvaly; 1095 SCIP_Real nomval; 1096 SCIP_Real denomval; 1097 1098 assert(expr != NULL); 1099 assert(auxvalue != NULL); 1100 1101 /**! [SnippetNlhdlrEvalauxQuotient] */ 1102 /* get auxiliary variables */ 1103 auxvarx = SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->numexpr); 1104 auxvary = SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->denomexpr); 1105 assert(auxvarx != NULL); 1106 assert(auxvary != NULL); 1107 1108 /* get solution values of the auxiliary variables */ 1109 solvalx = SCIPgetSolVal(scip, sol, auxvarx); 1110 solvaly = SCIPgetSolVal(scip, sol, auxvary); 1111 1112 /* evaluate expression w.r.t. the values of the auxiliary variables */ 1113 nomval = nlhdlrexprdata->numcoef * solvalx + nlhdlrexprdata->numconst; 1114 denomval = nlhdlrexprdata->denomcoef * solvaly + nlhdlrexprdata->denomconst; 1115 1116 /* return SCIP_INVALID if the denominator evaluates to zero */ 1117 *auxvalue = (denomval != 0.0) ? nlhdlrexprdata->constant + nomval / denomval : SCIP_INVALID; 1118 /**! [SnippetNlhdlrEvalauxQuotient] */ 1119 1120 return SCIP_OKAY; 1121 } 1122 1123 1124 /** nonlinear handler under/overestimation callback 1125 * 1126 * @todo which of the paramters did I not use, but have to be taken into consideration? 1127 */ 1128 static 1129 SCIP_DECL_NLHDLRESTIMATE(nlhdlrEstimateQuotient) 1130 { /*lint --e{715}*/ 1131 SCIP_Bool branchingusefulx = FALSE; 1132 SCIP_Bool branchingusefuly = FALSE; 1133 SCIP_ROWPREP* rowprep; 1134 1135 assert(nlhdlr != NULL); 1136 assert(expr != NULL); 1137 assert(nlhdlrexprdata != NULL); 1138 assert(rowpreps != NULL); 1139 1140 /** ![SnippetNlhdlrEstimateQuotient] */ 1141 *addedbranchscores = FALSE; 1142 *success = FALSE; 1143 1144 SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) ); 1145 1146 if( nlhdlrexprdata->numexpr == nlhdlrexprdata->denomexpr ) 1147 { 1148 /* univariate case */ 1149 SCIP_CALL( estimateUnivariateQuotient(scip, sol, nlhdlrexprdata->numexpr, nlhdlrexprdata->numcoef, nlhdlrexprdata->numconst, 1150 nlhdlrexprdata->denomcoef, nlhdlrexprdata->denomconst, nlhdlrexprdata->constant, overestimate, rowprep, 1151 &branchingusefulx, success) ); 1152 } 1153 else 1154 { 1155 /* bivariate case */ 1156 SCIP_CALL( estimateBivariateQuotient(scip, nlhdlrexprdata->numexpr, nlhdlrexprdata->denomexpr, SCIPgetExprAuxVarNonlinear(expr), sol, 1157 nlhdlrexprdata->numcoef, nlhdlrexprdata->numconst, nlhdlrexprdata->denomcoef, nlhdlrexprdata->denomconst, 1158 nlhdlrexprdata->constant, overestimate, rowprep, 1159 &branchingusefulx, &branchingusefuly, success) ); 1160 } 1161 1162 if( *success ) 1163 { 1164 SCIP_CALL( SCIPsetPtrarrayVal(scip, rowpreps, 0, rowprep) ); 1165 } 1166 else 1167 { 1168 SCIPfreeRowprep(scip, &rowprep); 1169 } 1170 1171 /* add branching scores if requested */ 1172 if( addbranchscores ) 1173 { 1174 SCIP_EXPR* exprs[2]; 1175 SCIP_Real violation; 1176 int nexprs = 0; 1177 1178 if( branchingusefulx ) 1179 exprs[nexprs++] = nlhdlrexprdata->numexpr; 1180 if( branchingusefuly ) 1181 exprs[nexprs++] = nlhdlrexprdata->denomexpr; 1182 1183 /* compute violation w.r.t. the auxiliary variable(s) */ 1184 #ifndef BRSCORE_ABSVIOL 1185 SCIP_CALL( SCIPgetExprRelAuxViolationNonlinear(scip, expr, auxvalue, sol, &violation, NULL, NULL) ); 1186 #else 1187 SCIP_CALL( SCIPgetExprAbsAuxViolationNonlinear(scip, expr, auxvalue, sol, &violation, NULL, NULL) ); 1188 #endif 1189 assert(violation > 0.0); /* there should be a violation if we were called to enforce */ 1190 1191 SCIP_CALL( SCIPaddExprsViolScoreNonlinear(scip, exprs, nexprs, violation, sol, addedbranchscores) ); 1192 } 1193 /** ![SnippetNlhdlrEstimateQuotient] */ 1194 1195 return SCIP_OKAY; 1196 } 1197 1198 /** nonlinear handler solution linearization callback */ 1199 static 1200 SCIP_DECL_NLHDLRSOLLINEARIZE(nlhdlrSollinearizeQuotient) 1201 { /*lint --e{715}*/ 1202 SCIP_VAR* x; 1203 SCIP_Real lbx; 1204 SCIP_Real ubx; 1205 SCIP_Real solx; 1206 int c; 1207 1208 assert(nlhdlr != NULL); 1209 assert(expr != NULL); 1210 assert(nlhdlrexprdata != NULL); 1211 1212 /* in the bivariate case, we do not get globally valid estimators */ 1213 if( nlhdlrexprdata->numexpr != nlhdlrexprdata->denomexpr ) 1214 return SCIP_OKAY; 1215 1216 x = SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->numexpr); 1217 lbx = SCIPvarGetLbGlobal(x); 1218 ubx = SCIPvarGetUbGlobal(x); 1219 solx = SCIPgetSolVal(scip, sol, x); 1220 solx = MAX(MIN(solx, ubx), lbx); 1221 1222 for( c = (overestimate ? 0 : 1); c < (underestimate ? 2 : 1); ++c ) /* c == 0: overestimate, c == 1: underestimate */ 1223 { 1224 SCIP_ROWPREP* rowprep; 1225 SCIP_Bool success = FALSE; 1226 SCIP_Bool local = TRUE; 1227 SCIP_Bool branchinguseful; 1228 SCIP_Real coef; 1229 SCIP_Real constant; 1230 1231 /* compute estimator, will be secant or gradient */ 1232 SCIP_CALL( estimateUnivariate(scip, lbx, ubx, lbx, ubx, solx, 1233 nlhdlrexprdata->numcoef, nlhdlrexprdata->numconst, nlhdlrexprdata->denomcoef, nlhdlrexprdata->denomconst, nlhdlrexprdata->constant, 1234 &coef, &constant, c == 0, &local, &branchinguseful, &success) ); 1235 1236 /* skip if not successful, only locally valid, or just a secant (branchinguseful is TRUE) */ 1237 if( !success || local || branchinguseful ) 1238 continue; 1239 1240 SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, c == 0 ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, FALSE) ); 1241 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, x, coef) ); 1242 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPgetExprAuxVarNonlinear(expr), -1.0) ); 1243 SCIProwprepAddConstant(rowprep, constant); 1244 (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "quot_%s_sol%d", SCIPvarGetName(x), SCIPsolGetIndex(sol)); 1245 1246 SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, sol, SCIPgetHugeValue(scip), &success) ); 1247 1248 /* if cleanup succeeded and rowprep is still global, add to cutpool */ 1249 if( success && !SCIProwprepIsLocal(rowprep) ) 1250 { 1251 SCIP_ROW* row; 1252 1253 SCIP_CALL( SCIPgetRowprepRowCons(scip, &row, rowprep, cons) ); 1254 SCIP_CALL( SCIPaddPoolCut(scip, row) ); 1255 SCIP_CALL( SCIPreleaseRow(scip, &row) ); 1256 } 1257 1258 SCIPfreeRowprep(scip, &rowprep); 1259 } 1260 1261 return SCIP_OKAY; 1262 } 1263 1264 /** nonlinear handler interval evaluation callback */ 1265 static 1266 SCIP_DECL_NLHDLRINTEVAL(nlhdlrIntevalQuotient) 1267 { /*lint --e{715}*/ 1268 SCIP_INTERVAL bnds; 1269 1270 assert(nlhdlrexprdata != NULL); 1271 assert(nlhdlrexprdata->numexpr != NULL); 1272 assert(nlhdlrexprdata->denomexpr != NULL); 1273 1274 /* it is not possible to compute tighter intervals if both expressions are different 1275 * we should not be called in this case, as we haven't said we would participate in this activity in detect 1276 */ 1277 assert(nlhdlrexprdata->numexpr == nlhdlrexprdata->denomexpr); 1278 1279 /**! [SnippetNlhdlrIntevalQuotient] */ 1280 /* get activity of the numerator (= denominator) expression */ 1281 bnds = SCIPexprGetActivity(nlhdlrexprdata->numexpr); 1282 1283 /* call interval evaluation for the univariate quotient expression */ 1284 *interval = intEvalQuotient(scip, bnds, nlhdlrexprdata->numcoef, nlhdlrexprdata->numconst, 1285 nlhdlrexprdata->denomcoef, nlhdlrexprdata->denomconst, nlhdlrexprdata->constant); 1286 /**! [SnippetNlhdlrIntevalQuotient] */ 1287 1288 return SCIP_OKAY; 1289 } 1290 1291 1292 /** nonlinear handler callback for reverse propagation */ 1293 static 1294 SCIP_DECL_NLHDLRREVERSEPROP(nlhdlrReversepropQuotient) 1295 { /*lint --e{715}*/ 1296 SCIP_INTERVAL result; 1297 1298 assert(nlhdlrexprdata != NULL); 1299 assert(nlhdlrexprdata->numexpr != NULL); 1300 assert(nlhdlrexprdata->denomexpr != NULL); 1301 1302 /* it is not possible to compute tighter intervals if both expressions are different 1303 * we should not be called in this case, as we haven't said we would participate in this activity in detect 1304 */ 1305 assert(nlhdlrexprdata->numexpr == nlhdlrexprdata->denomexpr); 1306 1307 SCIPdebugMsg(scip, "call reverse propagation for expression (%g %p + %g) / (%g %p + %g) + %g bounds [%g,%g]\n", 1308 nlhdlrexprdata->numcoef, (void*)nlhdlrexprdata->numexpr, nlhdlrexprdata->numconst, 1309 nlhdlrexprdata->denomcoef, (void*)nlhdlrexprdata->denomexpr, nlhdlrexprdata->denomconst, 1310 nlhdlrexprdata->constant, bounds.inf, bounds.sup); 1311 1312 /* call reverse propagation */ 1313 /**! [SnippetNlhdlrReversepropQuotient] */ 1314 result = reversepropQuotient(bounds, nlhdlrexprdata->numcoef, nlhdlrexprdata->numconst, 1315 nlhdlrexprdata->denomcoef, nlhdlrexprdata->denomconst, nlhdlrexprdata->constant); 1316 1317 SCIPdebugMsg(scip, "try to tighten bounds of %p: [%g,%g] -> [%g,%g]\n", 1318 (void*)nlhdlrexprdata->numexpr, SCIPgetExprBoundsNonlinear(scip, nlhdlrexprdata->numexpr).inf, 1319 SCIPgetExprBoundsNonlinear(scip, nlhdlrexprdata->numexpr).sup, result.inf, result.sup); 1320 1321 /* tighten bounds of the expression */ 1322 SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, nlhdlrexprdata->numexpr, result, infeasible, nreductions) ); 1323 /**! [SnippetNlhdlrReversepropQuotient] */ 1324 1325 return SCIP_OKAY; 1326 } 1327 1328 1329 /* 1330 * nonlinear handler specific interface methods 1331 */ 1332 1333 /** includes quotient nonlinear handler in nonlinear constraint handler */ 1334 SCIP_RETCODE SCIPincludeNlhdlrQuotient( 1335 SCIP* scip /**< SCIP data structure */ 1336 ) 1337 { 1338 SCIP_NLHDLRDATA* nlhdlrdata; 1339 SCIP_NLHDLR* nlhdlr; 1340 1341 assert(scip != NULL); 1342 1343 /* create nonlinear handler data */ 1344 nlhdlrdata = NULL; 1345 1346 SCIP_CALL( SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, NLHDLR_NAME, NLHDLR_DESC, NLHDLR_DETECTPRIORITY, 1347 NLHDLR_ENFOPRIORITY, nlhdlrDetectQuotient, nlhdlrEvalauxQuotient, nlhdlrdata) ); 1348 assert(nlhdlr != NULL); 1349 1350 SCIPnlhdlrSetCopyHdlr(nlhdlr, nlhdlrCopyhdlrQuotient); 1351 SCIPnlhdlrSetFreeExprData(nlhdlr, nlhdlrFreeExprDataQuotient); 1352 SCIPnlhdlrSetSepa(nlhdlr, NULL, NULL, nlhdlrEstimateQuotient, NULL); 1353 SCIPnlhdlrSetSollinearize(nlhdlr, nlhdlrSollinearizeQuotient); 1354 SCIPnlhdlrSetProp(nlhdlr, nlhdlrIntevalQuotient, nlhdlrReversepropQuotient); 1355 1356 return SCIP_OKAY; 1357 } 1358