1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 2 /* */ 3 /* This file is part of the program and library */ 4 /* SCIP --- Solving Constraint Integer Programs */ 5 /* */ 6 /* Copyright 2002-2022 Zuse Institute Berlin */ 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_signomial.c 26 * @ingroup DEFPLUGINS_NLHDLR 27 * @brief signomial nonlinear handler 28 * @author Liding Xu 29 */ 30 31 #ifdef SCIP_SIGCUT_DEBUG 32 #ifndef SCIP_SIG_DEBUG 33 #define SCIP_SIG_DEBUG 34 #endif 35 #endif 36 37 #ifdef SCIP_SIG_DEBUG 38 #define SCIP_DEBUG 39 #endif 40 41 #include "scip/nlhdlr_signomial.h" 42 #include "scip/cons_nonlinear.h" 43 #include "scip/pub_misc_rowprep.h" 44 #include "scip/pub_nlhdlr.h" 45 #include "scip/scip_expr.h" 46 #include "scip/expr_var.h" 47 #include "scip/expr_pow.h" 48 49 /* fundamental nonlinear handler properties */ 50 #define NLHDLR_NAME "signomial" 51 #define NLHDLR_DESC "handler for signomial expressions" 52 #define NLHDLR_DETECTPRIORITY 30 53 #define NLHDLR_ENFOPRIORITY 30 54 55 /* handler specific parameters */ 56 #define NLHDLR_MAXNUNDERVARS 14 /**< maximum number of variables when underestimating a concave power function (maximum: 14) */ 57 #define NLHDLR_MINCUTSCALE 1e-5 /**< minimum scale factor when scaling a cut (minimum: 1e-6) */ 58 59 60 /** nonlinear handler expression data 61 * 62 * A signomial expression admits the form \f$ cx^a = y \f$, where \f$ y \f$ is an auxiliary variable representing this 63 * expression and all \f$ x \f$ are positive. 64 * 65 * The natural formulation of the expression is defined as \f$ x^a = t = y/c \f$, where \f$ t \f$ is a 66 * non-existant slack variable denoting \f$ y/c \f$. The variables in \f$x\f$ with positive exponents form positive 67 * variables \f$ u \f$, and the associated exponents form positive exponents \f$ f \f$. The variables in \f$ x \f$ with 68 * negative exponents and \f$ t \f$ form negative variables \f$ v \f$, and the associated exponents form negative 69 * exponents \f$ g \f$. Let \f$ s = \max(|f|,|g|) \f$ be a normalization constant, where \f$ |\cdot| \f$ denotes the L1 norm. Apply a scaling step: 70 * Dividing the entries of \f$ f \f$ by \f$ s \f$, and dividing the entries of \f$ g \f$ by \f$ s \f$ as well. Then \f$ x^a = t \f$ has 71 * a reformulation \f$ u^f = v^g \f$, where \f$ u^f, v^g \$ are two concave power functions. 72 */ 73 struct SCIP_NlhdlrExprData 74 { 75 SCIP_Real coef; /**< coefficient \f$c\f$ */ 76 SCIP_EXPR** factors; /**< expression factors representing \f$x\f$ */ 77 int nfactors; /**< number of factors */ 78 int nvars; /**< number of variables \f$(x,y)\f$ */ 79 SCIP_Real* exponents; /**< exponents \f$e\f$ */ 80 int nposvars; /**< number of positive variables \f$u\f$ */ 81 int nnegvars; /**< number of negative variables \f$v\f$ */ 82 SCIP_Bool* signs; /**< indicators for sign of variables after reformulation, TRUE for positive, FALSE for negative */ 83 SCIP_Real* refexponents; /**< exponents of \f$(x,t)\f$ after reformulation (always positive) */ 84 SCIP_Bool isstorecapture; /**< are all variables already got? */ 85 86 /* working parameters will be modified after getting all variables */ 87 SCIP_VAR** vars; /**< variables \f$(x,y)\f$ */ 88 SCIP_INTERVAL* intervals; /**< intervals storing lower and upper bounds of variables \f$(x,y)\f$ */ 89 SCIP_Real* box; /**< the upper/lower bounds of variables, used in SCIPcomputeFacetVertexPolyhedralNonlinear() */ 90 SCIP_Real* xstar; /**< the values of variables, used in SCIPcomputeFacetVertexPolyhedralNonlinear() */ 91 SCIP_Real* facetcoefs; /**< the coefficients of variables, returned by SCIPcomputeFacetVertexPolyhedralNonlinear() */ 92 }; 93 94 /** nonlinear handler data */ 95 struct SCIP_NlhdlrData 96 { 97 /* parameters */ 98 int maxnundervars; /**< maximum number of variables in underestimating a concave power function */ 99 SCIP_Real mincutscale; /**< minimum scale factor when scaling a cut */ 100 }; 101 102 /** data struct to be passed on to vertexpoly-evalfunction (see SCIPcomputeFacetVertexPolyhedralNonlinear) */ 103 typedef struct 104 { 105 SCIP_NLHDLREXPRDATA* nlhdlrexprdata; /**< expression data */ 106 SCIP_Bool sign; /**< the sign of variables in the reformulated constraint for vertexpoly-evalfunction */ 107 int nsignvars; /**< the number of variables in the reformulated constraint for vertexpoly-evalfunction */ 108 SCIP* scip; /**< SCIP data structure */ 109 } VERTEXPOLYFUN_EVALDATA; 110 111 112 /* 113 * Local methods 114 */ 115 116 #ifdef SCIP_SIGCUT_DEBUG 117 118 /** print the information on a signomial term */ 119 static 120 SCIP_RETCODE printSignomial( 121 SCIP* scip, /**< SCIP data structure */ 122 SCIP_EXPR* expr, /**< expression */ 123 SCIP_NLHDLREXPRDATA* nlhdlrexprdata /**< expression data */ 124 ) 125 { 126 assert(expr != NULL); 127 128 /* print overall statistics and the expression */ 129 SCIPdebugMsg(scip, " #all variables: %d, #positive exponent variables: %d, #negative exponent variables: %d, auxvar: %s \n expr: ", 130 nlhdlrexprdata->nvars, nlhdlrexprdata->nposvars, nlhdlrexprdata->nnegvars, SCIPvarGetName(SCIPgetExprAuxVarNonlinear(expr)) ); 131 SCIPprintExpr(scip, expr, NULL); 132 133 /* if variables are detected, we can print more information */ 134 if( !nlhdlrexprdata->isstorecapture ) 135 { 136 SCIPinfoMessage(scip, NULL, "\n"); 137 return SCIP_OKAY; 138 } 139 140 /* print the natural formulation of the expression */ 141 SCIPinfoMessage(scip, NULL, ". natural formulation c x^a = y: %.2f", nlhdlrexprdata->coef); 142 for( int i = 0; i < nlhdlrexprdata->nvars - 1; i++ ) 143 { 144 SCIPinfoMessage(scip, NULL, "%s^%.2f", SCIPvarGetName(nlhdlrexprdata->vars[i]), nlhdlrexprdata->exponents[i]); 145 } 146 SCIPinfoMessage(scip, NULL, " = %s", SCIPvarGetName(nlhdlrexprdata->vars[nlhdlrexprdata->nvars - 1])); 147 148 /* print the reformulation of the expression */ 149 SCIPinfoMessage(scip, NULL, ". Reformulation u^f = v^g: "); 150 if( nlhdlrexprdata->nposvars == 0 ) 151 { 152 SCIPinfoMessage(scip, NULL, "1"); 153 } 154 else 155 { 156 for( int i = 0; i < nlhdlrexprdata->nvars; i++ ) 157 if( nlhdlrexprdata->signs[i] ) 158 { 159 SCIPinfoMessage(scip, NULL, "%s^%.2f", SCIPvarGetName(nlhdlrexprdata->vars[i]), nlhdlrexprdata->refexponents[i]); 160 } 161 } 162 SCIPinfoMessage(scip, NULL, " = "); 163 if( nlhdlrexprdata->nnegvars == 0 ) 164 { 165 SCIPinfoMessage(scip, NULL, "1"); 166 } 167 else 168 { 169 for( int i = 0; i < nlhdlrexprdata->nvars; i++ ) 170 { 171 if( !nlhdlrexprdata->signs[i] ) 172 { 173 if( i == nlhdlrexprdata->nvars - 1 ) 174 { 175 SCIPinfoMessage(scip, NULL, "(%s/%.2f)^%.2f", SCIPvarGetName(nlhdlrexprdata->vars[i]), nlhdlrexprdata->coef, nlhdlrexprdata->refexponents[i]); 176 } 177 else 178 { 179 SCIPinfoMessage(scip, NULL, "%s^%.2f", SCIPvarGetName(nlhdlrexprdata->vars[i]), nlhdlrexprdata->refexponents[i]); 180 } 181 } 182 } 183 } 184 SCIPinfoMessage(scip, NULL, "\n"); 185 186 return SCIP_OKAY; 187 } 188 189 #endif 190 191 /** free the memory of expression data */ 192 static 193 void freeExprDataMem( 194 SCIP* scip, /**< SCIP data structure */ 195 SCIP_NLHDLREXPRDATA** nlhdlrexprdata, /**< expression data */ 196 SCIP_Bool ispartial /**< free the partially allocated memory or the fully allocated memory? */ 197 ) 198 { 199 200 SCIPfreeBlockMemoryArrayNull(scip, &(*nlhdlrexprdata)->factors, (*nlhdlrexprdata)->nfactors); 201 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->exponents, (*nlhdlrexprdata)->nfactors); 202 if( !ispartial ) 203 { 204 int nvars = (*nlhdlrexprdata)->nvars; 205 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->signs, nvars); 206 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->refexponents, nvars); 207 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->vars, nvars); 208 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->intervals, nvars); 209 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->box, 2*nvars); /*lint !e647*/ 210 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->xstar, nvars); 211 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->facetcoefs, nvars); 212 } 213 SCIPfreeBlockMemory(scip, nlhdlrexprdata); 214 *nlhdlrexprdata = NULL; 215 } 216 217 /** reforms a rowprep to a standard form for nonlinear handlers 218 * 219 * The input rowprep is of the form \f$ a_u u + a_v v + b \f$. 220 * If in the overestimate mode, then we relax \f$ t \le x^a \f$, i.e., \f$ 0 \le u^f - v^g \f$. This implies that \f$ t \f$ is in \f$ v = (v',t) \f$. 221 * Therefore, the valid inequality is \f$ 0 \le a_u u + a_v v + b \f$. As overestimate mode requires that \f$ t \f$ is in the left hand side, 222 * the coefficients of \f$ t \f$ must be made negative while keeping the sign of the inequality, we can show that \f$ a_t \le 0 \f$, so it suffices 223 * to divide the both sides by \f$ -a_t \ge 0\f$, which yields \f$ t \le (a_u u + a_v' v' + b) / -a_t \f$. 224 * A rowprep in standard form only contains an estimator of the expression and no auxvar. 225 * If in the underestimate mode, then we relax \f$ x^a \le t \f$, i.e., \f$ u^f - v^g \le 0 \f$. This implies that \f$ t \f$ is in \f$ v = (v',t) \f$. 226 * Therefore, the valid inequality is \f$ a_u u + a_v v + b \le 0 \f$. As overestimate mode requires that \f$ t \f$ is in the left hand side, 227 * the coefficients of \f$ t \f$ must be made negative while keeping the sign of the inequality, we can show that \f$ a_t \le 0 \f$, so it suffices 228 * to divide the both sides by \f$ -a_t \ge 0 \f$, which yields \f$ (a_u u + a_v' v' + b) / -a_t \le t \f$. 229 * A rowprep in standard form only contains an estimator of the expression and no auxvar. 230 */ 231 static 232 void reformRowprep( 233 SCIP* scip, /**< SCIP data structure */ 234 SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< expression data */ 235 SCIP_ROWPREP* rowprep, /**< cut to be reformulated */ 236 SCIP_Real mincutscale, /**< min scaling factor for the cut in rowprep */ 237 SCIP_Bool* success /**< pointer to store whether the reformulating was successful */ 238 ) 239 { 240 int i; 241 int nvars; 242 SCIP_Real coefauxvar; 243 SCIP_Real scale; 244 SCIP_Real* coefs; 245 SCIP_VAR* auxvar; 246 SCIP_VAR** vars; 247 248 assert(rowprep != NULL); 249 assert(nlhdlrexprdata != NULL); 250 251 coefs = SCIProwprepGetCoefs(rowprep); 252 vars = SCIProwprepGetVars(rowprep); 253 nvars = SCIProwprepGetNVars(rowprep); 254 255 /* find auxvar's cut coefficient and set it to zero */ 256 auxvar = nlhdlrexprdata->vars[nlhdlrexprdata->nfactors]; 257 coefauxvar = 1.0; 258 for( i = 0; i < nvars; i++ ) 259 { 260 if( vars[i] == auxvar ) 261 { 262 coefauxvar = coefs[i]; 263 coefs[i] = 0.0; 264 break; 265 } 266 } 267 268 if( SCIPisZero(scip, coefauxvar) || coefauxvar >= 0.0 ) 269 { 270 *success = FALSE; 271 } 272 else 273 { 274 /* the reformation scales the cut so that coefficients and constant are divided by the absolute value of coefauxvar */ 275 assert(coefauxvar < 0.0); 276 scale = -1.0 / coefauxvar; 277 for( i = 0; i < nvars; i++ ) 278 { 279 if( vars[i] == auxvar ) 280 continue; 281 coefs[i] *= scale; 282 } 283 /* set the side to scale*side by adding -side and adding scale*side */ 284 SCIProwprepAddSide(rowprep, SCIProwprepGetSide(rowprep) * (-1.0 + scale)); 285 286 *success = SCIPisGT(scip, scale, mincutscale); 287 } 288 } 289 290 /** store and capture variables associated with the expression and its subexpressions */ 291 static 292 SCIP_RETCODE storeCaptureVars( 293 SCIP* scip, /**< SCIP data structure */ 294 SCIP_EXPR* expr, /**< expression */ 295 SCIP_NLHDLREXPRDATA* nlhdlrexprdata /**< expression data */ 296 ) 297 { 298 int c; 299 300 assert(!nlhdlrexprdata->isstorecapture); 301 302 /* get and capture variables \f$x\f$ */ 303 for( c = 0; c < nlhdlrexprdata->nfactors; ++c ) 304 { 305 nlhdlrexprdata->vars[c] = SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->factors[c]); 306 assert(nlhdlrexprdata->vars[c] != NULL); 307 308 SCIP_CALL( SCIPcaptureVar(scip, nlhdlrexprdata->vars[c]) ); 309 } 310 311 /* get and capture variable \f$y\f$ */ 312 nlhdlrexprdata->vars[c] = SCIPgetExprAuxVarNonlinear(expr); 313 assert(nlhdlrexprdata->vars[c] != NULL); 314 SCIP_CALL( SCIPcaptureVar(scip, nlhdlrexprdata->vars[c]) ); 315 316 nlhdlrexprdata->isstorecapture = TRUE; 317 318 return SCIP_OKAY; 319 } 320 321 /** get bounds of variables x,t and check whether they are box constrained signomial variables */ 322 static 323 void checkSignomialBounds( 324 SCIP* scip, /**< SCIP data structure */ 325 SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< expression data */ 326 SCIP_Bool* isboxsignomial /**< buffer to store whether variables are box constrained signomial variables */ 327 ) 328 { 329 int c; 330 SCIP_Real powinf; 331 SCIP_Real powsup; 332 SCIP_Real productinf = 1; 333 SCIP_Real productsup = 1; 334 335 assert(nlhdlrexprdata->isstorecapture); 336 337 *isboxsignomial = FALSE; 338 339 /* get bounds of x */ 340 for( c = 0; c < nlhdlrexprdata->nfactors; c++ ) 341 { 342 SCIP_Real inf = SCIPvarGetLbLocal(nlhdlrexprdata->vars[c]); 343 SCIP_Real sup = SCIPvarGetUbLocal(nlhdlrexprdata->vars[c]); 344 345 /* if the bounds of the variable are not positive and finite, or (bounds equal) then the expression is not a signomial */ 346 if( !SCIPisPositive(scip, inf) || !SCIPisPositive(scip, sup) || SCIPisInfinity(scip, sup) || SCIPisEQ(scip, inf, sup) ) 347 return; 348 349 nlhdlrexprdata->intervals[c].inf = inf; 350 nlhdlrexprdata->intervals[c].sup = MAX(sup, inf + 0.1); 351 powinf = pow(inf, nlhdlrexprdata->exponents[c]); 352 powsup = pow(sup, nlhdlrexprdata->exponents[c]); 353 productinf *= MIN(powinf, powsup); 354 productsup *= MAX(powinf, powsup); 355 } 356 357 /* compute bounds of t by bounds of x */ 358 nlhdlrexprdata->intervals[c].inf = productinf; 359 nlhdlrexprdata->intervals[c].sup = productsup; 360 361 *isboxsignomial = TRUE; 362 } 363 364 /** evaluate expression at solution w.r.t. auxiliary variables */ 365 static 366 SCIP_DECL_VERTEXPOLYFUN(nlhdlrExprEvalPower) 367 { 368 int i; 369 int j; 370 SCIP_Real val; 371 VERTEXPOLYFUN_EVALDATA* evaldata; 372 SCIP_NLHDLREXPRDATA* nlhdlrexprdata; 373 374 assert(args != NULL); 375 376 evaldata = (VERTEXPOLYFUN_EVALDATA *)funcdata; 377 378 assert(evaldata != NULL); 379 assert(nargs == evaldata->nsignvars); 380 381 nlhdlrexprdata = evaldata->nlhdlrexprdata; 382 383 #ifdef SCIP_MORE_DEBUG 384 SCIPdebugMsg(evaldata->scip, "eval vertexpolyfun\n"); 385 #endif 386 val = 1.0; 387 for( i = 0, j = 0; i < nlhdlrexprdata->nvars; ++i ) 388 { 389 if( nlhdlrexprdata->signs[i] != evaldata->sign ) 390 continue; 391 /* the reformulated exponent of args[j] is found */ 392 val *= pow(args[j], nlhdlrexprdata->refexponents[i]); 393 j++; 394 } 395 396 return val; 397 } 398 399 /** determine whether a power function \f$ w^h \f$ is special and add an overunderestimator or underestimator to a given rowprep 400 * 401 * \f$ w^h \f$ is special, if all variables are fixed, or it is a constant to estimate, a univariate power to estimate, 402 * or a bivariate power to underestimate. The estimator is multiplied by the multiplier and stored in the rowprep. 403 */ 404 static 405 SCIP_RETCODE estimateSpecialPower( 406 SCIP* scip, /**< SCIP data structure */ 407 SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nonlinear handler expression data */ 408 SCIP_Bool sign, /**< the sign of variables of the power function */ 409 SCIP_Real multiplier, /**< the multiplier of the estimator */ 410 SCIP_Bool overestimate, /**< whether overestimate or underestimator the power function */ 411 SCIP_SOL* sol, /**< solution \f$ w' \f$ to use in estimation */ 412 SCIP_ROWPREP* rowprep, /**< rowprep where to store estimator */ 413 SCIP_Bool* isspecial, /**< buffer to store whether this function is a special case */ 414 SCIP_Bool* success /**< buffer to store whether successful */ 415 ) 416 { 417 int i; 418 SCIP_Bool allfixed = TRUE; 419 int nsignvars; 420 421 assert(scip != NULL); 422 assert(nlhdlrexprdata != NULL); 423 assert(rowprep != NULL); 424 assert(isspecial != NULL); 425 assert(success != NULL); 426 427 *success = FALSE; 428 /* check whether all variables are fixed */ 429 for( i = 0; i < nlhdlrexprdata->nvars; ++i ) 430 { 431 if( nlhdlrexprdata->signs[i] != sign ) 432 continue; 433 434 if( !SCIPisRelEQ(scip, nlhdlrexprdata->intervals[i].inf, nlhdlrexprdata->intervals[i].sup) ) 435 { 436 allfixed = FALSE; 437 break; 438 } 439 } 440 441 /* allfixed is a special case */ 442 if( allfixed ) 443 { 444 /* return a constant */ 445 SCIP_Real funcval = 1.0; 446 SCIP_Real scale; 447 SCIP_Real val; 448 449 for( i = 0; i < nlhdlrexprdata->nvars; ++i ) 450 { 451 if( nlhdlrexprdata->signs[i] != sign ) 452 continue; 453 454 scale = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1.0; 455 val = SCIPgetSolVal(scip, sol, nlhdlrexprdata->vars[i]) / scale; 456 funcval *= pow(val, nlhdlrexprdata->refexponents[i]); 457 } 458 SCIProwprepAddConstant(rowprep, multiplier * funcval); 459 *isspecial = TRUE; 460 461 return SCIP_OKAY; 462 } 463 464 /* number of variables in the power function */ 465 nsignvars = sign ? nlhdlrexprdata->nposvars : nlhdlrexprdata->nnegvars; 466 467 /* if the power function has no more than 2 variables, this a special case */ 468 *isspecial = ( nsignvars <= 1 ) || ( nsignvars == 2 && !overestimate ); 469 if( !*isspecial ) 470 return SCIP_OKAY; 471 472 if( nsignvars == 0 ) 473 { 474 /* constant case */ 475 SCIProwprepAddConstant(rowprep, multiplier); 476 *success = TRUE; 477 } 478 else if( nsignvars == 1 ) 479 { 480 /* univariate case, w^h */ 481 for( i = 0; i < nlhdlrexprdata->nvars; ++i ) 482 { 483 if( nlhdlrexprdata->signs[i] == sign ) 484 { 485 /* the variable w is found */ 486 SCIP_Real scale = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1; 487 SCIP_VAR* var = nlhdlrexprdata->vars[i]; 488 SCIP_Real refexponent = nlhdlrexprdata->refexponents[i]; 489 if( refexponent == 1.0 ) 490 { 491 /* h = 1, a univariate linear function. Only rescale, no need for linearization */ 492 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, var, multiplier / scale) ); 493 } 494 else 495 { 496 /* a univariate power function */ 497 SCIP_Real facetconstant; 498 SCIP_Real facetcoef; 499 SCIP_Real val = SCIPgetSolVal(scip, sol, var) / scale; 500 /* local (using bounds) depends on whether we under- or overestimate */ 501 SCIP_Bool islocal = !overestimate; 502 SCIPestimateRoot(scip, refexponent, overestimate, nlhdlrexprdata->intervals[i].inf, nlhdlrexprdata->intervals[i].sup, 503 val, &facetconstant, &facetcoef, &islocal, success); 504 SCIProwprepAddConstant(rowprep, multiplier * facetconstant); 505 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, var, multiplier * facetcoef / scale) ); 506 } 507 } 508 } 509 *success = TRUE; 510 } 511 else if( nsignvars == 2 && !overestimate ) 512 { 513 /* bivariate case, f(w) = w^h = f_0(w_0) f_1(w_1). The convex under envelope is the maxmimum of 514 * two affine functions, each of which is determined by three extreme points the box bound. 515 * One affine function is supported by three lower left points; the other affine function is 516 * supported by three upper right points. For a point xstar in the box, its corresponding affine function can be determined by 517 * which region (upper right or lower left half space) the point is in. Thus, we can determine the region, and use the 518 * associated three extreme point to interpolate the affine function. 519 */ 520 SCIP_Bool isupperright; 521 SCIP_Real dw0, dw1; 522 SCIP_Real f0w0l, f0w0u, f1w1l, f1w1u; 523 SCIP_Real fw0lw1u, fw0uw1l; 524 SCIP_Real facetconstant; 525 SCIP_Real facetcoefs[2] = {0.0, 0.0}; 526 SCIP_VAR* vars[2] = {NULL, NULL}; 527 SCIP_Real refexponents[2] = {0.0, 0.0}; 528 SCIP_Real xstar[2] = {0.0, 0.0};; 529 SCIP_Real scale[2] = {0.0, 0.0};; 530 SCIP_Real box[4] = {0.0, 0.0, 0.0, 0.0}; 531 int j = 0; 532 533 /* get refexponents, xstar, scale, and box */ 534 for( i = 0; i < nlhdlrexprdata->nvars; ++i ) 535 { 536 if( nlhdlrexprdata->signs[i] != sign ) 537 continue; 538 box[2 * j] = nlhdlrexprdata->intervals[i].inf; 539 box[2 * j + 1] = nlhdlrexprdata->intervals[i].sup; 540 refexponents[j] = nlhdlrexprdata->refexponents[i]; 541 scale[j] = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1; 542 vars[j] = nlhdlrexprdata->vars[i]; 543 xstar[j] = SCIPgetSolVal(scip, sol, nlhdlrexprdata->vars[j]) / scale[j]; 544 j++; 545 } 546 547 /* compute the box length*/ 548 dw0 = box[1] - box[0]; 549 dw1 = box[3] - box[2]; 550 551 /* determine the location (upper right or lower left half sapce) of the xstar. 552 * (dw1, dw0) is the direction vector to the upper right half space. 553 */ 554 isupperright = ( (xstar[0] - box[0]) * dw1 + (xstar[1] - box[3]) * dw0 ) > 0.0; 555 556 /* compute function values of f_0, f_1 at vertices */ 557 f0w0l = pow(box[0], refexponents[0]); 558 f0w0u = pow(box[1], refexponents[0]); 559 f1w1l = pow(box[2], refexponents[1]); 560 f1w1u = pow(box[3], refexponents[1]); 561 fw0lw1u = f0w0l * f1w1u; 562 fw0uw1l = f0w0u * f1w1l; 563 if( isupperright ) 564 { 565 /* underestimator: fw0uw1u + (fw0uw1u - fw0lw1u) / (dw0) * (x0 - w0u) + (fw0uw1u - fw0uw1l) / (dw1) * (x1 - w1u) */ 566 SCIP_Real fw0uw1u = f0w0u * f1w1u; 567 facetcoefs[0] = (fw0uw1u - fw0lw1u) / dw0; 568 facetcoefs[1] = (fw0uw1u - fw0uw1l) / dw1; 569 facetconstant = fw0uw1u - facetcoefs[0] * box[1] - facetcoefs[1] * box[3]; 570 } 571 else 572 { 573 /* underestimator: fw0lw1l + (fw0uw1l - fw0lw1l) / (dw0) * (x0 - w0l) + (fw0lw1u- fw0lw1l) / (dw1) * (x1 - w1l) */ 574 SCIP_Real fw0lw1l = f0w0l * f1w1l; 575 facetcoefs[0] = (fw0uw1l - fw0lw1l) / dw0; 576 facetcoefs[1] = (fw0lw1u - fw0lw1l) / dw1; 577 facetconstant = fw0lw1l - facetcoefs[0] * box[0] - facetcoefs[1] * box[2]; 578 } 579 SCIProwprepAddConstant(rowprep, multiplier * facetconstant); 580 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, vars[0], multiplier * facetcoefs[0] / scale[0]) ); 581 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, vars[1], multiplier * facetcoefs[1] / scale[1]) ); 582 *success = TRUE; 583 } 584 585 return SCIP_OKAY; 586 } 587 588 /** adds an underestimator for a multivariate concave power function \f$ w^h \f$ to a given rowprep 589 * 590 * Calls \ref SCIPcomputeFacetVertexPolyhedralNonlinear() for \f$ w^h \f$ and 591 * box set to local bounds of auxiliary variables. The estimator is multiplied 592 * by the multiplier and stored in the rowprep. 593 */ 594 static 595 SCIP_RETCODE underEstimatePower( 596 SCIP* scip, /**< SCIP data structure */ 597 SCIP_CONSHDLR* conshdlr, /**< nonlinear constraint handler */ 598 SCIP_NLHDLR* nlhdlr, /**< nonlinear handler */ 599 SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nonlinear handler expression data */ 600 SCIP_Bool sign, /**< the sign of variables of the power function */ 601 SCIP_Real multiplier, /**< the multiplier of the estimator */ 602 SCIP_SOL* sol, /**< solution \f$ w' \f$ to use*/ 603 SCIP_Real targetvalue, /**< a target value to achieve; if not reachable, then can give up early */ 604 SCIP_ROWPREP* rowprep, /**< rowprep where to store estimator */ 605 SCIP_Bool* success /**< buffer to store whether successful */ 606 ) 607 { 608 int i; 609 int j; 610 int nsignvars; 611 SCIP_Real facetconstant; 612 SCIP_Real scale; 613 SCIP_Real* box; 614 SCIP_Real* facetcoefs; 615 SCIP_Real* xstar; 616 VERTEXPOLYFUN_EVALDATA evaldata; 617 618 assert(scip != NULL); 619 assert(nlhdlr != NULL); 620 assert(nlhdlrexprdata != NULL); 621 assert(rowprep != NULL); 622 assert(success != NULL); 623 624 *success = FALSE; 625 626 /* number of variables of sign */ 627 nsignvars = sign ? nlhdlrexprdata->nposvars : nlhdlrexprdata->nnegvars; 628 629 /* data structure to evaluate the power function */ 630 evaldata.nlhdlrexprdata = nlhdlrexprdata; 631 evaldata.sign = sign; 632 evaldata.nsignvars = nsignvars; 633 evaldata.scip = scip; 634 635 /* compute box constraints and reference value of variables*/ 636 xstar = nlhdlrexprdata->xstar; 637 box = nlhdlrexprdata->box; 638 j = 0; 639 for( i = 0; i < nlhdlrexprdata->nvars; ++i ) 640 { 641 if( nlhdlrexprdata->signs[i] != sign ) 642 continue; 643 644 box[2 * j] = nlhdlrexprdata->intervals[i].inf; 645 box[2 * j + 1] = nlhdlrexprdata->intervals[i].sup; 646 scale = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1.0; 647 xstar[j] = SCIPgetSolVal(scip, sol, nlhdlrexprdata->vars[i]) / scale; 648 j++; 649 } 650 651 /* find a facet of the underestimator */ 652 facetcoefs = nlhdlrexprdata->facetcoefs; 653 SCIP_CALL( SCIPcomputeFacetVertexPolyhedralNonlinear(scip, conshdlr, FALSE, nlhdlrExprEvalPower, (void*)&evaldata, xstar, box, 654 nsignvars, targetvalue, success, facetcoefs, &facetconstant) ); 655 656 if( !*success ) 657 { 658 SCIPdebugMsg(scip, "failed to compute facet of convex hull\n"); 659 return SCIP_OKAY; 660 } 661 662 /* set coefficients in the rowprep */ 663 SCIProwprepAddConstant(rowprep, multiplier * facetconstant); 664 j = 0; 665 for( i = 0; i < nlhdlrexprdata->nvars; ++i ) 666 { 667 if( nlhdlrexprdata->signs[i] != sign ) 668 continue; 669 scale = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1.0; 670 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, nlhdlrexprdata->vars[i], multiplier * facetcoefs[j] / scale) ); 671 j++; 672 } 673 674 return SCIP_OKAY; 675 } 676 677 /** adds an overestimator for a concave power function \f$ w^h \f$ to a given rowprep 678 * 679 * The estimator is multiplied by the multiplier and stored in the rowprep. 680 */ 681 static 682 SCIP_RETCODE overEstimatePower( 683 SCIP* scip, /**< SCIP data structure */ 684 SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nonlinear handler expression data */ 685 SCIP_Bool sign, /**< the sign of variables of the power function */ 686 SCIP_Real multiplier, /**< the multiplier of the estimator */ 687 SCIP_SOL* sol, /**< solution to use */ 688 SCIP_ROWPREP* rowprep, /**< rowprep where to store estimator */ 689 SCIP_Bool* success /**< buffer to store whether successful */ 690 ) 691 { 692 int i; 693 SCIP_Real facetcoef; 694 SCIP_Real facetconstant; 695 SCIP_Real funcval; 696 SCIP_Real refexponent; 697 SCIP_Real sumrefexponents; 698 SCIP_Real scale; 699 SCIP_Real val; 700 SCIP_VAR* var; 701 assert(scip != NULL); 702 assert(nlhdlrexprdata != NULL); 703 assert(rowprep != NULL); 704 assert(success != NULL); 705 706 *success = FALSE; 707 708 /* compute the value and the sum of reformulated exponents of the power function */ 709 sumrefexponents = 0; 710 funcval = 1; 711 for( i = 0; i < nlhdlrexprdata->nvars; ++i ) 712 { 713 if( nlhdlrexprdata->signs[i] != sign ) 714 continue; 715 scale = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1.0; 716 val = SCIPgetSolVal(scip, sol, nlhdlrexprdata->vars[i]) / scale; 717 val = MAX(val, 0.1); 718 refexponent = nlhdlrexprdata->refexponents[i]; 719 sumrefexponents += refexponent; 720 funcval *= pow(val, refexponent); 721 } 722 723 /* overestimate by gradient cut: w'^h + h w'^(h - 1)(w - w') */ 724 facetconstant = (1 - sumrefexponents) * funcval; 725 SCIProwprepAddConstant(rowprep, multiplier * facetconstant); 726 for( i = 0; i < nlhdlrexprdata->nvars; ++i ) 727 { 728 if( nlhdlrexprdata->signs[i] != sign ) 729 continue; 730 scale = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1.0; 731 var = nlhdlrexprdata->vars[i]; 732 val = SCIPgetSolVal(scip, sol, var) / scale; 733 val = MAX(val, 0.1); 734 facetcoef = nlhdlrexprdata->refexponents[i] * funcval / val; 735 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, var, multiplier * facetcoef / scale) ); 736 } 737 738 *success = TRUE; 739 740 return SCIP_OKAY; 741 } 742 743 /* 744 * Callback methods of nonlinear handler 745 */ 746 747 /** nonlinear handler estimation callback */ 748 static 749 SCIP_DECL_NLHDLRESTIMATE(nlhdlrEstimateSignomial) 750 { /*lint --e{715}*/ 751 SCIP_Bool isboxsignomial; 752 SCIP_Bool isspecial; 753 SCIP_Bool undersign; 754 SCIP_Bool oversign; 755 int i; 756 int nundervars; 757 SCIP_Real undermultiplier; 758 SCIP_Real overmultiplier; 759 SCIP_NLHDLRDATA* nlhdlrdata; 760 SCIP_ROWPREP* rowprep; 761 SCIP_Real targetunder; 762 SCIP_Real scale; 763 764 assert(conshdlr != NULL); 765 assert(nlhdlr != NULL); 766 assert(expr != NULL); 767 assert(nlhdlrexprdata != NULL); 768 assert(rowpreps != NULL); 769 *success = FALSE; 770 *addedbranchscores = FALSE; 771 772 /* store and capture the vars of an expression, if the vars are not stored and captured yet */ 773 if( !nlhdlrexprdata->isstorecapture ) 774 { 775 SCIP_CALL( storeCaptureVars(scip, expr, nlhdlrexprdata) ); 776 } 777 778 /* check whether all variables have finite positive bounds, which is necessary for the expression to be a signomial term */ 779 /* TODO consider allowing 0.0 lower bounds for u variables (and handle gradients at 0.0) */ 780 isboxsignomial = FALSE; 781 checkSignomialBounds(scip, nlhdlrexprdata, &isboxsignomial); 782 783 if( !isboxsignomial ) 784 return SCIP_OKAY; 785 786 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr); 787 788 /* determine the sign of variables for over- and underestimators, and the multiplier for estimators in the rowprep */ 789 if( overestimate ) 790 { 791 /* relax t <= x^a, i.e., 0 <= u^f - v^g, overestimate u^f, underestimate v^g */ 792 oversign = TRUE; 793 undersign = FALSE; 794 overmultiplier = 1.0; 795 undermultiplier = -1.0; 796 nundervars = nlhdlrexprdata->nnegvars; 797 } 798 else 799 { 800 /* relax x^a <= t, i.e., u^f - v^g <= 0, underestimate u^f, overestimate v^g */ 801 oversign = FALSE; 802 undersign = TRUE; 803 overmultiplier = -1.0; 804 undermultiplier = 1.0; 805 nundervars = nlhdlrexprdata->nposvars; 806 } 807 808 /* compute the value of the overestimator, which is a target value for computing the underestimator efficiently */ 809 targetunder = 1.0; 810 for( i = 0; i < nlhdlrexprdata->nvars; i++ ) 811 { 812 if( nlhdlrexprdata->signs[i] == oversign ) 813 { 814 scale = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1.0; 815 targetunder *= pow(SCIPgetSolVal(scip, sol, nlhdlrexprdata->vars[i]) / scale, nlhdlrexprdata->refexponents[i]); 816 } 817 } 818 819 #ifdef SCIP_SIGCUT_DEBUG 820 SCIP_Real targetover = 1.0; 821 822 /* print information on estimators */ 823 SCIP_CALL( printSignomial(scip, expr, nlhdlrexprdata)); 824 SCIPinfoMessage(scip, NULL, " Auxvalue: %f, targetvalue: %f, %sestimate.", auxvalue, targetvalue, overestimate ? "over" : "under"); 825 826 targetunder = 1.0; 827 for( i = 0; i < nlhdlrexprdata->nvars; i++ ) 828 { 829 if( nlhdlrexprdata->signs[i] == undersign ) 830 { 831 scale = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1.0; 832 targetover *= pow(SCIPgetSolVal(scip, sol, nlhdlrexprdata->vars[i]) / scale, nlhdlrexprdata->refexponents[i]); 833 } 834 } 835 SCIPinfoMessage(scip, NULL, " Underestimate: %s, overestimate: %s.", 836 undersign ? "positive" : "negative", oversign ? "positive" : "negative"); 837 SCIPinfoMessage(scip, NULL, " Undervalue (targetover): %f, overvalue (targetunder): %f.", targetover, targetunder); 838 #endif 839 840 /* create a rowprep and allocate memory for it */ 841 SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) ); 842 SCIP_CALL( SCIPensureRowprepSize(scip, rowprep, nlhdlrexprdata->nvars + 1) ); 843 844 /* only underestimate a concave function, if the number of variables is below a given threshold */ 845 if( nundervars <= nlhdlrdata->maxnundervars ) 846 { 847 /* compute underestimator */ 848 isspecial = FALSE; 849 /* check and compute the special case */ 850 SCIP_CALL( estimateSpecialPower(scip, nlhdlrexprdata, undersign, undermultiplier, FALSE, sol, rowprep, &isspecial, success) ); 851 if( !isspecial ) 852 { 853 SCIP_CALL( underEstimatePower(scip, conshdlr, nlhdlr, nlhdlrexprdata, undersign, undermultiplier, sol, targetunder, rowprep, success) ); 854 } 855 } 856 857 if( *success ) 858 { 859 /* compute overestimator */ 860 isspecial = FALSE; 861 /* check and compute the special case */ 862 SCIP_CALL( estimateSpecialPower(scip, nlhdlrexprdata, oversign, overmultiplier, TRUE, sol, rowprep, &isspecial, success) ); 863 if( !isspecial ) 864 { 865 SCIP_CALL( overEstimatePower(scip, nlhdlrexprdata, oversign, overmultiplier, sol, rowprep, success) ); 866 } 867 } 868 869 if( *success ) 870 { 871 reformRowprep(scip, nlhdlrexprdata, rowprep, nlhdlrdata->mincutscale, success); 872 if( *success ) 873 { 874 SCIP_CALL( SCIPsetPtrarrayVal(scip, rowpreps, 0, rowprep) ); 875 } 876 } 877 878 if( !*success ) 879 SCIPfreeRowprep(scip, &rowprep); 880 881 return SCIP_OKAY; 882 } 883 884 /** callback to detect structure in expression tree */ 885 static 886 SCIP_DECL_NLHDLRDETECT(nlhdlrDetectSignomial) 887 { /*lint --e{715}*/ 888 889 assert(expr != NULL); 890 assert(enforcing != NULL); 891 assert(participating != NULL); 892 893 /* for now, we do not get active if separation is already (or does not need to be) provided */ 894 if( (*enforcing & SCIP_NLHDLR_METHOD_SEPABOTH) == SCIP_NLHDLR_METHOD_SEPABOTH ) 895 { 896 return SCIP_OKAY; 897 } 898 899 /* check for product expressions with more than one child */ 900 if( SCIPisExprProduct(scip, expr) && SCIPexprGetNChildren(expr) >= 2 ) 901 { 902 int c; 903 int nf = SCIPexprGetNChildren(expr); 904 int nvars = nf + 1; 905 SCIP_Bool ismultilinear = TRUE; 906 907 /* create expression data for the nonlinear handler */ 908 SCIP_CALL( SCIPallocClearBlockMemory(scip, nlhdlrexprdata) ); 909 (*nlhdlrexprdata)->nfactors = nf; 910 (*nlhdlrexprdata)->nvars = nvars; 911 912 /* allocate memory for expression data */ 913 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlhdlrexprdata)->factors, nf) ); 914 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlhdlrexprdata)->exponents, nf) ); 915 916 /* get monomial information */ 917 SCIP_CALL( SCIPgetExprMonomialData(scip, expr, &((*nlhdlrexprdata)->coef), (*nlhdlrexprdata)->exponents, (*nlhdlrexprdata)->factors) ); 918 919 /* skip multilinear terms, since we wouldn't do better than expr_product */ 920 for( c = 0; c < nf; c++ ) 921 { 922 if( !SCIPisEQ(scip, (*nlhdlrexprdata)->exponents[c], 1.0) ) 923 { 924 ismultilinear = FALSE; 925 break; 926 } 927 } 928 929 if( ismultilinear ) 930 { 931 /* if multilinear, we free memory of the expression data and do nothing */ 932 freeExprDataMem(scip, nlhdlrexprdata, TRUE); 933 return SCIP_OKAY; 934 } 935 else 936 { 937 SCIP_Real normalize; 938 SCIP_Real sumlexponents = 0; 939 SCIP_Real sumrexponents = 1; 940 int nposvars = 0; 941 942 /* allocate more memory for expression data */ 943 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlhdlrexprdata)->signs, nvars) ); 944 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlhdlrexprdata)->refexponents, nvars) ); 945 SCIP_CALL( SCIPallocClearBlockMemoryArray(scip, &(*nlhdlrexprdata)->vars, nvars) ); 946 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlhdlrexprdata)->intervals, nvars) ); 947 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlhdlrexprdata)->xstar, nvars) ); 948 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlhdlrexprdata)->facetcoefs, nvars) ); 949 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlhdlrexprdata)->box, 2 * nvars) ); 950 951 (*nlhdlrexprdata)->isstorecapture = FALSE; 952 953 /* detect more information for the reformulation: we first compute the sum 954 * of positive and negative exponents and determine the sign indicators 955 */ 956 for( c = 0; c < nf; c++ ) 957 { 958 /* capture sub expressions */ 959 SCIPcaptureExpr((*nlhdlrexprdata)->factors[c]); 960 if( (*nlhdlrexprdata)->exponents[c] > 0.0 ) 961 { 962 /* add a positive exponent */ 963 sumlexponents += (*nlhdlrexprdata)->exponents[c]; 964 (*nlhdlrexprdata)->signs[c] = TRUE; 965 nposvars++; 966 } 967 else 968 { 969 /* subtract a negative exponent */ 970 sumrexponents -= (*nlhdlrexprdata)->exponents[c]; 971 (*nlhdlrexprdata)->signs[c] = FALSE; 972 } 973 /* set null to working variables, meaning that they are not stored yet */ 974 } 975 (*nlhdlrexprdata)->signs[nf] = FALSE; 976 (*nlhdlrexprdata)->nposvars = nposvars; 977 (*nlhdlrexprdata)->nnegvars = nf - nposvars + 1; 978 979 /* compute the normalization constant */ 980 normalize = MAX(sumlexponents, sumrexponents); 981 /* normalize positive and negative exponents */ 982 for( c = 0; c < nf; c++ ) 983 { 984 if( (*nlhdlrexprdata)->signs[c] ) 985 (*nlhdlrexprdata)->refexponents[c] = (*nlhdlrexprdata)->exponents[c] / normalize; 986 else 987 (*nlhdlrexprdata)->refexponents[c] = -(*nlhdlrexprdata)->exponents[c] / normalize; 988 } 989 (*nlhdlrexprdata)->refexponents[nf] = 1.0 / normalize; 990 991 /* tell children that we will use their auxvar and use its activity for estimation */ 992 for( c = 0; c < nf; c++ ) 993 { 994 SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, (*nlhdlrexprdata)->factors[c], TRUE, FALSE, TRUE, TRUE) ); 995 } 996 } 997 } 998 999 if( *nlhdlrexprdata != NULL ) 1000 { 1001 *participating = SCIP_NLHDLR_METHOD_SEPABOTH; 1002 } 1003 1004 #ifdef SCIP_SIGCUT_DEBUG 1005 if( *participating ) 1006 { 1007 SCIPdebugMsg(scip, "scip depth: %d, step: %d, expr pointer: %p, expr data pointer: %p, detected expr: total vars (exps) %d ", 1008 SCIPgetSubscipDepth(scip), SCIPgetStage(scip), (void *)expr, (void *)nlhdlrexprdata, (*nlhdlrexprdata)->nfactors); 1009 SCIPprintExpr(scip, expr, NULL); 1010 SCIPinfoMessage(scip, NULL, " participating: %d\n", *participating); 1011 } 1012 #endif 1013 1014 return SCIP_OKAY; 1015 } 1016 1017 /** auxiliary evaluation callback of nonlinear handler */ 1018 static SCIP_DECL_NLHDLREVALAUX(nlhdlrEvalauxSignomial) 1019 { /*lint --e{715}*/ 1020 int c; 1021 SCIP_Real val; 1022 SCIP_VAR* var; 1023 1024 *auxvalue = nlhdlrexprdata->coef; 1025 for( c = 0; c < nlhdlrexprdata->nfactors; ++c ) 1026 { 1027 assert(nlhdlrexprdata->factors[c] != NULL); 1028 1029 var = SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->factors[c]); 1030 assert(var != NULL); 1031 1032 val = SCIPgetSolVal(scip, sol, var); 1033 if( SCIPisPositive(scip, val) ) 1034 { 1035 *auxvalue *= pow(val, nlhdlrexprdata->exponents[c]); 1036 } 1037 else 1038 { 1039 *auxvalue = SCIP_INVALID; 1040 break; 1041 } 1042 } 1043 1044 return SCIP_OKAY; 1045 } 1046 1047 /** nonlinear handler copy callback */ 1048 static 1049 SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrSignomial) 1050 { /*lint --e{715}*/ 1051 assert(targetscip != NULL); 1052 assert(sourcenlhdlr != NULL); 1053 assert(strcmp(SCIPnlhdlrGetName(sourcenlhdlr), NLHDLR_NAME) == 0); 1054 1055 SCIP_CALL( SCIPincludeNlhdlrSignomial(targetscip) ); 1056 1057 return SCIP_OKAY; 1058 } 1059 1060 /** callback to free data of handler */ 1061 static 1062 SCIP_DECL_NLHDLRFREEHDLRDATA(nlhdlrFreehdlrDataSignomial) 1063 { /*lint --e{715}*/ 1064 assert(nlhdlrdata != NULL); 1065 1066 SCIPfreeBlockMemory(scip, nlhdlrdata); 1067 1068 return SCIP_OKAY; 1069 } 1070 1071 /** callback to free expression specific data */ 1072 static 1073 SCIP_DECL_NLHDLRFREEEXPRDATA(nlhdlrFreeExprDataSignomial) 1074 { /*lint --e{715}*/ 1075 int c; 1076 1077 /* release expressions */ 1078 for( c = 0; c < (*nlhdlrexprdata)->nfactors; c++ ) 1079 { 1080 SCIP_CALL( SCIPreleaseExpr(scip, &(*nlhdlrexprdata)->factors[c]) ); 1081 } 1082 1083 /* release variables */ 1084 if( (*nlhdlrexprdata)->isstorecapture ) 1085 { 1086 for( c = 0; c < (*nlhdlrexprdata)->nvars; c++ ) 1087 { 1088 if( (*nlhdlrexprdata)->vars[c] != NULL ) 1089 { 1090 SCIP_CALL( SCIPreleaseVar(scip, &(*nlhdlrexprdata)->vars[c]) ); 1091 } 1092 } 1093 } 1094 1095 /* free memory */ 1096 freeExprDataMem(scip, nlhdlrexprdata, FALSE); 1097 1098 return SCIP_OKAY; 1099 } 1100 1101 1102 /* 1103 * nonlinear handler specific interface methods 1104 */ 1105 1106 /** includes signomial nonlinear handler in nonlinear constraint handler */ 1107 SCIP_RETCODE 1108 SCIPincludeNlhdlrSignomial( 1109 SCIP* scip /**< SCIP data structure */ 1110 ) 1111 { 1112 SCIP_NLHDLRDATA* nlhdlrdata; 1113 SCIP_NLHDLR* nlhdlr; 1114 1115 assert(scip != NULL); 1116 1117 /* create nonlinear handler specific data */ 1118 SCIP_CALL( SCIPallocBlockMemory(scip, &nlhdlrdata)); 1119 BMSclearMemory(nlhdlrdata); 1120 1121 SCIP_CALL( SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, NLHDLR_NAME, NLHDLR_DESC, NLHDLR_DETECTPRIORITY, 1122 NLHDLR_ENFOPRIORITY, nlhdlrDetectSignomial, nlhdlrEvalauxSignomial, nlhdlrdata) ); 1123 assert(nlhdlr != NULL); 1124 1125 SCIPnlhdlrSetCopyHdlr(nlhdlr, nlhdlrCopyhdlrSignomial); 1126 SCIPnlhdlrSetFreeHdlrData(nlhdlr, nlhdlrFreehdlrDataSignomial); 1127 SCIPnlhdlrSetFreeExprData(nlhdlr, nlhdlrFreeExprDataSignomial); 1128 SCIPnlhdlrSetSepa(nlhdlr, NULL, NULL, nlhdlrEstimateSignomial, NULL); 1129 1130 /* parameters */ 1131 SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/maxnundervars", 1132 "maximum number of variables when underestimating a concave power function", 1133 &nlhdlrdata->maxnundervars, TRUE, NLHDLR_MAXNUNDERVARS, 2, 14, NULL, NULL) ); 1134 SCIP_CALL( SCIPaddRealParam(scip, "nlhdlr/" NLHDLR_NAME "/mincutscale", 1135 "minimum scale factor when scaling a cut", 1136 &nlhdlrdata->mincutscale, TRUE, NLHDLR_MINCUTSCALE, 1e-6, 1e6, NULL, NULL) ); 1137 1138 return SCIP_OKAY; 1139 } 1140