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 sepa_gauge.c 26 * @ingroup DEFPLUGINS_SEPA 27 * @brief gauge separator 28 * @author Felipe Serrano 29 * 30 * @todo should separator only be run when SCIPallColsInLP is true? 31 * @todo add SCIPisStopped(scip) to the condition of time consuming loops 32 * @todo check if it makes sense to implement the copy callback 33 */ 34 35 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 36 37 #include <assert.h> 38 #include <string.h> 39 40 #include "blockmemshell/memory.h" 41 #include "scip/scip_nlpi.h" 42 #include "scip/nlpi_ipopt.h" 43 #include "scip/nlpioracle.h" 44 #include "scip/scip_expr.h" 45 #include "scip/pub_expr.h" 46 #include "scip/pub_lp.h" 47 #include "scip/pub_message.h" 48 #include "scip/pub_misc.h" 49 #include "scip/pub_nlp.h" 50 #include "scip/pub_sepa.h" 51 #include "scip/pub_var.h" 52 #include "scip/scip_cut.h" 53 #include "scip/scip_lp.h" 54 #include "scip/scip_mem.h" 55 #include "scip/scip_message.h" 56 #include "scip/scip_nlp.h" 57 #include "scip/scip_numerics.h" 58 #include "scip/scip_param.h" 59 #include "scip/scip_prob.h" 60 #include "scip/scip_sepa.h" 61 #include "scip/scip_sol.h" 62 #include "scip/scip_solvingstats.h" 63 #include "scip/scip_timing.h" 64 #include "scip/sepa_gauge.h" 65 #include <string.h> 66 67 68 #define SEPA_NAME "gauge" 69 #define SEPA_DESC "gauge separator" 70 #define SEPA_PRIORITY 0 71 #define SEPA_FREQ -1 72 #define SEPA_MAXBOUNDDIST 1.0 73 #define SEPA_USESSUBSCIP FALSE /**< does the separator use a secondary SCIP instance? */ 74 #define SEPA_DELAY FALSE /**< should separation method be delayed, if other separators found cuts? */ 75 76 #define VIOLATIONFAC 100 /**< constraints regarded as violated when violation > VIOLATIONFAC*SCIPfeastol */ 77 #define MAX_ITER 75 /**< maximum number of iterations for the line search */ 78 79 #define DEFAULT_NLPITERLIM 1000 /**< default NLP iteration limit */ 80 81 #define NLPFEASFAC 1e-1/**< NLP feasibility tolerance = NLPFEASFAC * SCIP's feasibility tolerance */ 82 83 #define INTERIOROBJVARLB -100 /**< lower bound of the objective variable when computing interior point */ 84 85 /* 86 * Data structures 87 */ 88 89 /** side that makes a nlrow convex */ 90 enum ConvexSide 91 { 92 LHS = 0, /**< left hand side */ 93 RHS = 1 /**< right hand side */ 94 }; 95 typedef enum ConvexSide CONVEXSIDE; 96 97 /** position of a point */ 98 enum Position 99 { 100 INTERIOR = 0, /**< point is in the interior of the region */ 101 BOUNDARY = 1, /**< point is in the boundary of the region */ 102 EXTERIOR = 2 /**< point is in the exterior of the region */ 103 }; 104 typedef enum Position POSITION; 105 106 /** separator data */ 107 struct SCIP_SepaData 108 { 109 SCIP_NLROW** nlrows; /**< stores convex nlrows */ 110 CONVEXSIDE* convexsides; /**< which sides make the nlrows convex */ 111 int* nlrowsidx; /**< indices of nlrows that violate the current lp solution */ 112 int nnlrowsidx; /**< total number of convex nonlinear nlrows that violate the current lp solution */ 113 int nnlrows; /**< total number of convex nonlinear nlrows */ 114 int nlrowssize; /**< memory allocated for nlrows, convexsides and nlrowsidx */ 115 116 SCIP_Bool isintsolavailable; /**< do we have an interior point available? */ 117 SCIP_Bool skipsepa; /**< whether separator should be skipped */ 118 SCIP_SOL* intsol; /**< stores interior point */ 119 120 int ncuts; /**< number of cuts generated */ 121 122 /* parameters */ 123 int nlpiterlimit; /**< iteration limit of NLP solver; 0 for no limit */ 124 }; 125 126 /* 127 * Local methods 128 */ 129 130 /** stores, from the constraints represented by nlrows, the nonlinear convex ones in sepadata */ 131 static 132 SCIP_RETCODE storeNonlinearConvexNlrows( 133 SCIP* scip, /**< SCIP data structure */ 134 SCIP_SEPADATA* sepadata, /**< separator data */ 135 SCIP_NLROW** nlrows, /**< nlrows from which to store convex ones */ 136 int nnlrows /**< number of nlrows */ 137 ) 138 { 139 int i; 140 141 assert(scip != NULL); 142 assert(sepadata != NULL); 143 assert(nlrows != NULL); 144 assert(nnlrows > 0); 145 146 SCIPdebugMsg(scip, "storing convex nlrows\n"); 147 148 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(sepadata->nlrows), nnlrows) ); 149 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(sepadata->convexsides), nnlrows) ); 150 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(sepadata->nlrowsidx), nnlrows) ); 151 sepadata->nlrowssize = nnlrows; 152 153 sepadata->nnlrows = 0; 154 for( i = 0; i < nnlrows; ++i ) 155 { 156 SCIP_NLROW* nlrow; 157 158 nlrow = nlrows[i]; 159 assert(nlrow != NULL); 160 161 /* linear case */ 162 if( SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_LINEAR || SCIPnlrowGetExpr(nlrow) == NULL ) 163 continue; 164 165 /* nonlinear case */ 166 if( !SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrow)) && SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_CONVEX ) 167 { 168 sepadata->convexsides[sepadata->nnlrows] = RHS; 169 sepadata->nlrows[sepadata->nnlrows] = nlrow; 170 ++(sepadata->nnlrows); 171 } 172 else if( !SCIPisInfinity(scip, -SCIPnlrowGetLhs(nlrow)) && SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_CONCAVE ) 173 { 174 sepadata->convexsides[sepadata->nnlrows] = LHS; 175 sepadata->nlrows[sepadata->nnlrows] = nlrow; 176 ++(sepadata->nnlrows); 177 } 178 } 179 180 return SCIP_OKAY; 181 } 182 183 /** computes an interior point of a convex NLP relaxation 184 * 185 * builds the convex relaxation, modifies it to find an interior 186 * point, solves it and frees it; more details in @ref sepa_gauge.h 187 * 188 * @note the method also counts the number of nonlinear convex constraints and if there are < 2, then the convex 189 * relaxation is not interesting and the separator will not run again 190 */ 191 static 192 SCIP_RETCODE computeInteriorPoint( 193 SCIP* scip, /**< SCIP data structure */ 194 SCIP_SEPADATA* sepadata /**< separator data */ 195 ) 196 { 197 SCIP_NLPIORACLE* nlpioracle; 198 SCIP_NLPIPROBLEM* nlpiprob; 199 SCIP_NLPI* nlpi; 200 SCIP_HASHMAP* var2nlpiidx; 201 SCIP_Real objvarlb; 202 SCIP_Real minusone; 203 SCIP_Real one; 204 int nconvexnlrows; 205 int objvaridx; 206 int nconss; 207 int nvars; 208 int i; 209 210 assert(scip != NULL); 211 assert(sepadata != NULL); 212 assert(!sepadata->skipsepa); 213 214 SCIPdebugMsg(scip, "Computing interior point\n"); 215 216 /* create convex relaxation NLP */ 217 assert(SCIPgetNNlpis(scip) > 0); 218 219 nlpi = SCIPgetNlpis(scip)[0]; 220 assert(nlpi != NULL); 221 222 nvars = SCIPgetNVars(scip); 223 SCIP_CALL( SCIPhashmapCreate(&var2nlpiidx, SCIPblkmem(scip), nvars) ); 224 SCIP_CALL( SCIPcreateNlpiProblemFromNlRows(scip, nlpi, &nlpiprob, "gauge-interiorpoint-nlp", SCIPgetNLPNlRows(scip), SCIPgetNNLPNlRows(scip), var2nlpiidx, 225 NULL, NULL, SCIPgetCutoffbound(scip), FALSE, TRUE) ); 226 227 /* add objective variable; the problem is \min t, s.t. g(x) <= t, l(x) <= 0, where g are nonlinear and l linear */ 228 objvaridx = nvars; 229 objvarlb = INTERIOROBJVARLB; 230 one = 1.0; 231 SCIP_CALL( SCIPaddNlpiVars(scip, nlpi, nlpiprob, 1, &objvarlb, NULL, NULL) ); 232 SCIP_CALL( SCIPsetNlpiObjective(scip, nlpi, nlpiprob, 1, &objvaridx, &one, NULL, 0.0) ); 233 234 /* add objective variables to constraints; for this we need to get nlpi oracle to have access to number of 235 * constraints and which constraints are nonlinear 236 */ 237 /* @todo: this code is only valid when using IPOPT and needs to be changed when new NLP solvers get interfaced */ 238 assert(strcmp(SCIPnlpiGetName(nlpi), "ipopt") == 0); 239 nlpioracle = (SCIP_NLPIORACLE *)SCIPgetNlpiOracleIpopt(nlpiprob); 240 assert(nlpioracle != NULL); 241 assert(SCIPnlpiOracleGetNVars(nlpioracle) == objvaridx + 1); 242 243 minusone = -1.0; 244 nconvexnlrows = 0; 245 nconss = SCIPnlpiOracleGetNConstraints(nlpioracle); 246 for( i = 0; i < nconss; i++ ) 247 { 248 if( SCIPnlpiOracleIsConstraintNonlinear(nlpioracle, i) ) 249 { 250 SCIP_CALL( SCIPchgNlpiLinearCoefs(scip, nlpi, nlpiprob, i, 1, &objvaridx, &minusone) ); 251 ++nconvexnlrows; 252 } 253 } 254 SCIPdebug( SCIP_CALL( SCIPnlpiOraclePrintProblem(scip, nlpioracle, NULL) ) ); 255 256 /* check if convex relaxation is interesting */ 257 if( nconvexnlrows < 2 ) 258 { 259 SCIPdebugMsg(scip, "convex relaxation is not interesting, only %d nonlinear convex rows; abort\n", nconvexnlrows); 260 sepadata->skipsepa = TRUE; 261 goto CLEANUP; 262 } 263 264 /* add linear rows */ 265 SCIP_CALL( SCIPaddNlpiProblemRows(scip, nlpi, nlpiprob, var2nlpiidx, SCIPgetLPRows(scip), SCIPgetNLPRows(scip)) ); 266 267 /* compute interior point */ 268 SCIPdebugMsg(scip, "starting interior point computation\n"); 269 SCIP_CALL( SCIPsolveNlpi(scip, nlpi, nlpiprob, 270 .iterlimit = sepadata->nlpiterlimit > 0 ? sepadata->nlpiterlimit : INT_MAX, 271 .feastol = NLPFEASFAC * SCIPfeastol(scip), 272 .opttol = MAX(SCIPfeastol(scip), SCIPdualfeastol(scip))) ); /*lint !e666*/ 273 SCIPdebugMsg(scip, "finish interior point computation\n"); 274 275 #ifdef SCIP_DEBUG 276 { 277 SCIP_NLPSTATISTICS nlpstatistics; 278 279 /* get statistics */ 280 SCIP_CALL( SCIPgetNlpiStatistics(scip, nlpi, nlpiprob, &nlpstatistics) ); 281 282 SCIPdebugMsg(scip, "nlpi took iters %d, time %g searching for an find interior point: solstat %d\n", 283 nlpstatistics.niterations, nlpstatistics.totaltime, 284 SCIPgetNlpiSolstat(scip, nlpi, nlpiprob)); 285 } 286 #endif 287 288 if( SCIPgetNlpiSolstat(scip, nlpi, nlpiprob) <= SCIP_NLPSOLSTAT_FEASIBLE ) 289 { 290 SCIP_Real* nlpisol; 291 292 SCIP_CALL( SCIPgetNlpiSolution(scip, nlpi, nlpiprob, &nlpisol, NULL, NULL, NULL, NULL) ); 293 294 assert(nlpisol != NULL); 295 SCIPdebugMsg(scip, "NLP solved: sol found has objvalue = %g\n", nlpisol[objvaridx]); 296 297 /* if we found an interior point store it */ 298 if( SCIPisFeasNegative(scip, nlpisol[objvaridx]) ) 299 { 300 SCIPdebugMsg(scip, "Interior point found!, storing it\n"); 301 SCIP_CALL( SCIPcreateSol(scip, &sepadata->intsol, NULL) ); 302 for( i = 0; i < nvars; i ++ ) 303 { 304 SCIP_VAR* var; 305 306 var = SCIPgetVars(scip)[i]; 307 assert(SCIPhashmapExists(var2nlpiidx, (void*)var) ); 308 309 /* @todo: filter zero? */ 310 SCIP_CALL( SCIPsetSolVal(scip, sepadata->intsol, var, 311 nlpisol[SCIPhashmapGetImageInt(var2nlpiidx, (void *)var)]) ); 312 } 313 314 sepadata->isintsolavailable = TRUE; 315 } 316 else 317 { 318 SCIPdebugMsg(scip, "We got a feasible point but not interior (objval: %g)\n", nlpisol[objvaridx]); 319 sepadata->skipsepa = TRUE; 320 } 321 } 322 else 323 { 324 SCIPdebugMsg(scip, "We couldn't get an interior point (stat: %d)\n", SCIPgetNlpiSolstat(scip, nlpi, nlpiprob)); 325 sepadata->skipsepa = TRUE; 326 } 327 328 CLEANUP: 329 /* free memory */ 330 SCIPhashmapFree(&var2nlpiidx); 331 SCIP_CALL( SCIPfreeNlpiProblem(scip, nlpi, &nlpiprob) ); 332 333 return SCIP_OKAY; 334 } 335 336 337 /** find whether point is in the interior, at the boundary, or in the exterior of the region described by the 338 * intersection of `nlrows[i]` ≤ rhs if `convexsides[i]` = RHS or lhs ≤ `nlrows[i]` if `convexsides[i]` = LHS 339 * 340 * @note point corresponds to a convex combination between the LP solution and the interior point 341 */ 342 static 343 SCIP_RETCODE findPointPosition( 344 SCIP* scip, /**< SCIP data structure */ 345 SCIP_NLROW** nlrows, /**< nlrows defining the region */ 346 int* nlrowsidx, /**< indices of nlrows defining the region */ 347 int nnlrowsidx, /**< number of nlrows indices */ 348 CONVEXSIDE* convexsides, /**< sides of the nlrows involved in the region */ 349 SCIP_SOL* point, /**< point for which we want to know its position */ 350 POSITION* position /**< buffer to store position of sol */ 351 ) 352 { 353 int i; 354 355 assert(scip != NULL); 356 assert(nlrows != NULL); 357 assert(convexsides != NULL); 358 assert(nnlrowsidx > 0); 359 assert(point != NULL); 360 assert(position != NULL); 361 362 *position = INTERIOR; 363 for( i = 0; i < nnlrowsidx; i++ ) 364 { 365 SCIP_NLROW* nlrow; 366 SCIP_Real activity; 367 CONVEXSIDE convexside; 368 369 nlrow = nlrows[nlrowsidx[i]]; 370 convexside = convexsides[nlrowsidx[i]]; 371 372 /* compute activity of nlrow at point */ 373 SCIP_CALL( SCIPgetNlRowSolActivity(scip, nlrow, point, &activity) ); 374 375 if( convexside == RHS ) 376 { 377 assert(!SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrow))); 378 379 /* if nlrow <= rhs is violated, then we are in the exterior */ 380 if( SCIPisFeasGT(scip, activity, SCIPnlrowGetRhs(nlrow)) ) 381 { 382 *position = EXTERIOR; 383 SCIPdebugMsg(scip, "exterior because cons <%s> has activity %g. rhs: %g\n", SCIPnlrowGetName(nlrow), 384 activity, SCIPnlrowGetRhs(nlrow)); 385 SCIPdebug( SCIPprintNlRow(scip, nlrow, NULL) ); 386 387 return SCIP_OKAY; 388 } 389 390 /* if nlrow(point) == rhs, then we are currently at the boundary */ 391 if( SCIPisFeasEQ(scip, activity, SCIPnlrowGetRhs(nlrow)) ) 392 *position = BOUNDARY; 393 } 394 else 395 { 396 assert(!SCIPisInfinity(scip, -SCIPnlrowGetLhs(nlrow))); 397 assert(convexside == LHS); 398 399 /* if lhs <= nlrow is violated, then we are in the exterior */ 400 if( SCIPisFeasLT(scip, activity, SCIPnlrowGetLhs(nlrow)) ) 401 { 402 *position = EXTERIOR; 403 return SCIP_OKAY; 404 } 405 406 /* if lhs == nlrow(point), then we are currently at the boundary */ 407 if( SCIPisFeasEQ(scip, activity, SCIPnlrowGetLhs(nlrow)) ) 408 *position = BOUNDARY; 409 } 410 } 411 412 return SCIP_OKAY; 413 } 414 415 416 /** returns, in convexcomb, the convex combination 417 * \f$ \lambda\, \text{endpoint} + (1 - \lambda) \text{startpoint} = \text{startpoint} + \lambda (\text{endpoint} - \text{startpoint})\f$ 418 */ 419 static 420 SCIP_RETCODE buildConvexCombination( 421 SCIP* scip, /**< SCIP data structure */ 422 SCIP_Real lambda, /**< convex combination multiplier */ 423 SCIP_SOL* startpoint, /**< point corresponding to \f$ \lambda = 0 \f$ */ 424 SCIP_SOL* endpoint, /**< point corresponding to \f$ \lambda = 1 \f$ */ 425 SCIP_SOL* convexcomb /**< solution to store convex combination of intsol and tosepasol */ 426 ) 427 { 428 SCIP_VAR** vars; 429 int nvars; 430 int i; 431 432 assert(scip != NULL); 433 assert(startpoint != NULL); 434 assert(endpoint != NULL); 435 assert(convexcomb != NULL); 436 437 vars = SCIPgetVars(scip); 438 nvars = SCIPgetNVars(scip); 439 440 for( i = 0; i < nvars; i++ ) 441 { 442 SCIP_Real val; 443 SCIP_VAR* var; 444 445 var = vars[i]; 446 val = lambda * SCIPgetSolVal(scip, endpoint, var) + (1.0 - lambda) * SCIPgetSolVal(scip, startpoint, var); 447 448 if( !SCIPisZero(scip, val) ) 449 { 450 SCIP_CALL( SCIPsetSolVal(scip, convexcomb, var, val) ); 451 } 452 else 453 { 454 SCIP_CALL( SCIPsetSolVal(scip, convexcomb, var, 0.0) ); 455 } 456 } 457 458 return SCIP_OKAY; 459 } 460 461 462 /** performs binary search to find the point belonging to the segment [`intsol`, `tosepasol`] that intersects the boundary 463 * of the region described by the intersection of `nlrows[i]` ≤ rhs if `convexsides[i] = RHS` or lhs ≤ `nlrows[i]` if not, 464 * for i in `nlrowsidx` 465 */ 466 static 467 SCIP_RETCODE findBoundaryPoint( 468 SCIP* scip, /**< SCIP data structure */ 469 SCIP_NLROW** nlrows, /**< nlrows defining the region */ 470 int* nlrowsidx, /**< indices of nlrows defining the region */ 471 int nnlrowsidx, /**< number of nlrows indices */ 472 CONVEXSIDE* convexsides, /**< sides of the nlrows involved in the region */ 473 SCIP_SOL* intsol, /**< point acting as 'interior point' */ 474 SCIP_SOL* tosepasol, /**< solution that should be separated */ 475 SCIP_SOL* sol, /**< convex combination of intsol and lpsol */ 476 POSITION* position /**< buffer to store position of sol */ 477 ) 478 { 479 SCIP_Real lb; 480 SCIP_Real ub; 481 int i; 482 483 assert(scip != NULL); 484 assert(nlrows != NULL); 485 assert(nlrowsidx != NULL); 486 assert(convexsides != NULL); 487 assert(intsol != NULL); 488 assert(tosepasol != NULL); 489 assert(sol != NULL); 490 assert(position != NULL); 491 492 SCIPdebugMsg(scip, "starting binary search\n"); 493 lb = 0.0; /* corresponds to intsol */ 494 ub = 1.0; /* corresponds to tosepasol */ 495 for( i = 0; i < MAX_ITER; i++ ) 496 { 497 /* sol = (ub+lb)/2 * lpsol + (1 - (ub+lb)/2) * intsol */ 498 SCIP_CALL( buildConvexCombination(scip, (ub + lb)/2.0, intsol, tosepasol, sol) ); 499 500 /* find poisition of point: boundary, interior, exterior */ 501 SCIP_CALL( findPointPosition(scip, nlrows, nlrowsidx, nnlrowsidx, convexsides, sol, position) ); 502 SCIPdebugMsg(scip, "Position: %d, lambda: %g\n", *position, (ub + lb)/2.0); 503 504 switch( *position ) 505 { 506 case BOUNDARY: 507 SCIPdebugMsg(scip, "Done\n"); 508 return SCIP_OKAY; 509 510 case INTERIOR: 511 /* want to be closer to tosepasol */ 512 lb = (ub + lb)/2.0; 513 break; 514 515 case EXTERIOR: 516 /* want to be closer to intsol */ 517 ub = (ub + lb)/2.0; 518 break; 519 } 520 } 521 SCIPdebugMsg(scip, "Done\n"); 522 return SCIP_OKAY; 523 } 524 525 526 /** computes gradient cut (linearization) of nlrow at sol */ 527 static 528 SCIP_RETCODE generateCut( 529 SCIP* scip, /**< SCIP data structure */ 530 SCIP_SOL* sol, /**< point used to construct gradient cut (x_0) */ 531 SCIP_NLROW* nlrow, /**< constraint */ 532 CONVEXSIDE convexside, /**< whether we use rhs or lhs of nlrow */ 533 SCIP_EXPRITER* exprit, /**< expression iterator that can be used */ 534 SCIP_ROW* row, /**< storage for cut */ 535 SCIP_Bool* success /**< buffer to store whether the gradient was finite */ 536 ) 537 { 538 SCIP_EXPR* expr; 539 SCIP_Real exprval; 540 SCIP_Real gradx0; /* <grad f(x_0), x_0> */ 541 int i; 542 543 assert(scip != NULL); 544 assert(nlrow != NULL); 545 assert(row != NULL); 546 547 gradx0 = 0.0; 548 *success = TRUE; 549 550 SCIP_CALL( SCIPcacheRowExtensions(scip, row) ); 551 552 #ifdef CUT_DEBUG 553 SCIPdebug( SCIP_CALL( SCIPprintNlRow(scip, nlrow, NULL) ) ); 554 #endif 555 556 /* linear part */ 557 for( i = 0; i < SCIPnlrowGetNLinearVars(nlrow); i++ ) 558 { 559 SCIP_CALL( SCIPaddVarToRow(scip, row, SCIPnlrowGetLinearVars(nlrow)[i], SCIPnlrowGetLinearCoefs(nlrow)[i]) ); 560 } 561 562 expr = SCIPnlrowGetExpr(nlrow); 563 assert(expr != NULL); 564 565 SCIP_CALL( SCIPevalExprGradient(scip, expr, sol, 0L) ); 566 567 SCIP_CALL( SCIPexpriterInit(exprit, expr, SCIP_EXPRITER_DFS, FALSE) ); 568 for( ; !SCIPexpriterIsEnd(exprit); expr = SCIPexpriterGetNext(exprit) ) /*lint !e441*/ /*lint !e440*/ 569 { 570 SCIP_Real grad; 571 SCIP_VAR* var; 572 573 if( !SCIPisExprVar(scip, expr) ) 574 continue; 575 576 grad = SCIPexprGetDerivative(expr); 577 var = SCIPgetVarExprVar(expr); 578 assert(var != NULL); 579 580 /* check gradient entries: function might not be differentiable */ 581 if( !SCIPisFinite(grad) || grad == SCIP_INVALID ) /*lint !e777*/ 582 { 583 *success = FALSE; 584 break; 585 } 586 /* SCIPdebugMsg(scip, "grad w.r.t. <%s> (%g) = %g, gradx0 += %g\n", SCIPvarGetName(var), SCIPgetSolVal(scip, sol, var), grad, grad * SCIPgetSolVal(scip, sol, var)); */ 587 588 gradx0 += grad * SCIPgetSolVal(scip, sol, var); 589 SCIP_CALL( SCIPaddVarToRow(scip, row, var, grad) ); 590 } 591 592 SCIP_CALL( SCIPflushRowExtensions(scip, row) ); 593 594 /* if there was a problem computing the cut -> return */ 595 if( ! *success ) 596 return SCIP_OKAY; 597 598 #ifdef CUT_DEBUG 599 SCIPdebugMsg(scip, "gradient: "); 600 SCIPdebug( SCIP_CALL( SCIPprintRow(scip, row, NULL) ) ); 601 SCIPdebugMsg(scip, "gradient dot x_0: %g\n", gradx0); 602 #endif 603 604 /* gradient cut is linear part + f(x_0) - <grad f(x_0), x_0> + <grad f(x_0), x> <= rhs or >= lhs */ 605 exprval = SCIPexprGetEvalValue(SCIPnlrowGetExpr(nlrow)); 606 assert(exprval != SCIP_INVALID); /* we should have noticed a domain error above */ /*lint !e777*/ 607 if( convexside == RHS ) 608 { 609 assert(!SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrow))); 610 SCIP_CALL( SCIPchgRowRhs(scip, row, SCIPnlrowGetRhs(nlrow) - SCIPnlrowGetConstant(nlrow) - exprval + gradx0) ); 611 } 612 else 613 { 614 assert(convexside == LHS); 615 assert(!SCIPisInfinity(scip, -SCIPnlrowGetLhs(nlrow))); 616 SCIP_CALL( SCIPchgRowLhs(scip, row, SCIPnlrowGetLhs(nlrow) - SCIPnlrowGetConstant(nlrow) - exprval + gradx0) ); 617 } 618 619 #ifdef CUT_DEBUG 620 SCIPdebugMsg(scip, "gradient cut: "); 621 SCIPdebug( SCIP_CALL( SCIPprintRow(scip, row, NULL) ) ); 622 #endif 623 624 return SCIP_OKAY; 625 } 626 627 /** tries to generate gradient cuts at the point on the segment [`intsol`, `tosepasol`] that intersecs the boundary of the 628 * convex relaxation 629 * 630 * -# checks that the relative interior of the segment actually intersects the boundary 631 * (this check is needed since `intsol` is not necessarily an interior point) 632 * -# finds point on the boundary 633 * -# generates gradient cut at point on the boundary 634 */ 635 static 636 SCIP_RETCODE separateCuts( 637 SCIP* scip, /**< SCIP data structure */ 638 SCIP_SEPA* sepa, /**< the cut separator itself */ 639 SCIP_SOL* tosepasol, /**< solution that should be separated */ 640 SCIP_RESULT* result /**< pointer to store the result of the separation call */ 641 ) 642 { 643 SCIP_SEPADATA* sepadata; 644 SCIP_NLROW** nlrows; 645 CONVEXSIDE* convexsides; 646 SCIP_SOL* sol; 647 SCIP_SOL* intsol; 648 POSITION position; 649 int* nlrowsidx; 650 int nnlrowsidx; 651 int i; 652 SCIP_EXPRITER* exprit; 653 654 assert(sepa != NULL); 655 656 sepadata = SCIPsepaGetData(sepa); 657 assert(sepadata != NULL); 658 659 intsol = sepadata->intsol; 660 nlrows = sepadata->nlrows; 661 nlrowsidx = sepadata->nlrowsidx; 662 nnlrowsidx = sepadata->nnlrowsidx; 663 convexsides = sepadata->convexsides; 664 665 assert(intsol != NULL); 666 assert(nlrows != NULL); 667 assert(nlrowsidx != NULL); 668 assert(nnlrowsidx > 0); 669 assert(convexsides != NULL); 670 671 /* to evaluate the nlrow one needs a solution */ 672 SCIP_CALL( SCIPcreateSol(scip, &sol, NULL) ); 673 674 /* don't separate if, under SCIP tolerances, only a slight perturbation of the interior point in the direction of 675 * tosepasol gives a point that is in the exterior */ 676 SCIP_CALL( buildConvexCombination(scip, VIOLATIONFAC * SCIPfeastol(scip), intsol, tosepasol, sol) ); 677 SCIP_CALL( findPointPosition(scip, nlrows, nlrowsidx, nnlrowsidx, convexsides, sol, &position) ); 678 679 if( position == EXTERIOR ) 680 { 681 #ifdef SCIP_DEBUG 682 SCIPdebugMsg(scip, "segment joining intsol and tosepasol seems to be contained in the exterior of the region, can't separate\n"); 683 /* move from intsol in the direction of -tosepasol to check if we are really tangent to the region */ 684 SCIP_CALL( buildConvexCombination(scip, -1e-3, intsol, tosepasol, sol) ); 685 SCIP_CALL( findPointPosition(scip, nlrows, nlrowsidx, nnlrowsidx, convexsides, sol, &position) ); 686 if( position == EXTERIOR ) 687 { 688 SCIPdebugMsg(scip, "line through intsol and tosepasol is tangent to region; can't separate\n"); 689 } 690 SCIP_CALL( findPointPosition(scip, nlrows, nlrowsidx, nnlrowsidx, convexsides, intsol, &position) ); 691 printf("Position of intsol is %s\n", 692 position == EXTERIOR ? "exterior" : position == INTERIOR ? "interior": "boundary"); 693 SCIP_CALL( findPointPosition(scip, nlrows, nlrowsidx, nnlrowsidx, convexsides, tosepasol, &position) ); 694 printf("Position of tosepasol is %s\n", 695 position == EXTERIOR ? "exterior" : position == INTERIOR ? "interior": "boundary"); 696 697 /* slightly move from intsol in the direction of +-tosepasol */ 698 SCIP_CALL( buildConvexCombination(scip, 1e-5, intsol, tosepasol, sol) ); 699 SCIP_CALL( findPointPosition(scip, nlrows, nlrowsidx, nnlrowsidx, convexsides, sol, &position) ); 700 printf("Position of intsol + 0.00001(tosepasol - inisol) is %s\n", 701 position == EXTERIOR ? "exterior" : position == INTERIOR ? "interior": "boundary"); 702 SCIPdebug( SCIPprintSol(scip, sol, NULL, FALSE) ); 703 704 SCIP_CALL( buildConvexCombination(scip, -1e-5, intsol, tosepasol, sol) ); 705 SCIP_CALL( findPointPosition(scip, nlrows, nlrowsidx, nnlrowsidx, convexsides, sol, &position) ); 706 printf("Position of intsol - 0.00001(tosepasol - inisol) is %s\n", 707 position == EXTERIOR ? "exterior" : position == INTERIOR ? "interior": "boundary"); 708 SCIPdebug( SCIPprintSol(scip, sol, NULL, FALSE) ); 709 #endif 710 *result = SCIP_DIDNOTFIND; 711 goto CLEANUP; 712 } 713 714 /* find point on boundary */ 715 if( position != BOUNDARY ) 716 { 717 SCIP_CALL( findBoundaryPoint(scip, nlrows, nlrowsidx, nnlrowsidx, convexsides, intsol, tosepasol, sol, 718 &position) ); 719 720 /* if MAX_ITER weren't enough to find a point in the boundary we don't separate */ 721 if( position != BOUNDARY ) 722 { 723 SCIPdebugMsg(scip, "couldn't find boundary point, don't separate\n"); 724 goto CLEANUP; 725 } 726 } 727 728 /** @todo: could probably be moved inside generateCut */ 729 SCIP_CALL( SCIPcreateExpriter(scip, &exprit) ); 730 731 /* generate cuts at sol */ 732 for( i = 0; i < nnlrowsidx; i++ ) 733 { 734 SCIP_NLROW* nlrow; 735 SCIP_ROW* row; 736 SCIP_Real activity; 737 CONVEXSIDE convexside; 738 SCIP_Bool success; 739 char rowname[SCIP_MAXSTRLEN]; 740 741 nlrow = nlrows[nlrowsidx[i]]; 742 convexside = convexsides[nlrowsidx[i]]; 743 744 (void) SCIPsnprintf(rowname, SCIP_MAXSTRLEN, "%s_%u", SCIPnlrowGetName(nlrow), ++(sepadata->ncuts)); 745 746 /* only separate nlrows that are tight at the boundary point */ 747 SCIP_CALL( SCIPgetNlRowSolActivity(scip, nlrow, sol, &activity) ); 748 SCIPdebugMsg(scip, "cons <%s> at boundary point has activity: %g\n", SCIPnlrowGetName(nlrow), activity); 749 750 if( (convexside == RHS && !SCIPisFeasEQ(scip, activity, SCIPnlrowGetRhs(nlrow))) 751 || (convexside == LHS && !SCIPisFeasEQ(scip, activity, SCIPnlrowGetLhs(nlrow))) ) 752 continue; 753 754 /* cut is globally valid, since we work on nlrows from the NLP built at the root node, which are globally valid */ 755 /* @todo: when local nlrows get supported in SCIP, one can think of recomputing the interior point */ 756 SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &row, sepa, rowname, -SCIPinfinity(scip), SCIPinfinity(scip), 757 FALSE, FALSE , TRUE) ); 758 SCIP_CALL( generateCut(scip, sol, nlrow, convexside, exprit, row, &success) ); 759 760 /* add cut */ 761 SCIPdebugMsg(scip, "cut <%s> has efficacy %g\n", SCIProwGetName(row), SCIPgetCutEfficacy(scip, NULL, row)); 762 if( success && SCIPisCutEfficacious(scip, NULL, row) ) 763 { 764 SCIP_Bool infeasible; 765 766 SCIPdebugMsg(scip, "adding cut\n"); 767 SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) ); 768 769 if( infeasible ) 770 { 771 *result = SCIP_CUTOFF; 772 SCIP_CALL( SCIPreleaseRow(scip, &row) ); 773 break; 774 } 775 else 776 { 777 *result = SCIP_SEPARATED; 778 } 779 } 780 781 /* release the row */ 782 SCIP_CALL( SCIPreleaseRow(scip, &row) ); 783 } 784 785 SCIPfreeExpriter(&exprit); 786 787 CLEANUP: 788 SCIP_CALL( SCIPfreeSol(scip, &sol) ); 789 790 return SCIP_OKAY; 791 } 792 793 /* 794 * Callback methods of separator 795 */ 796 797 798 /** destructor of separator to free user data (called when SCIP is exiting) */ 799 static 800 SCIP_DECL_SEPAFREE(sepaFreeGauge) 801 { /*lint --e{715}*/ 802 SCIP_SEPADATA* sepadata; 803 804 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0); 805 806 /* free separator data */ 807 sepadata = SCIPsepaGetData(sepa); 808 assert(sepadata != NULL); 809 810 SCIPfreeBlockMemory(scip, &sepadata); 811 812 SCIPsepaSetData(sepa, NULL); 813 814 return SCIP_OKAY; 815 } 816 817 818 /** solving process deinitialization method of separator (called before branch and bound process data is freed) */ 819 static 820 SCIP_DECL_SEPAEXITSOL(sepaExitsolGauge) 821 { /*lint --e{715}*/ 822 SCIP_SEPADATA* sepadata; 823 824 assert(sepa != NULL); 825 826 sepadata = SCIPsepaGetData(sepa); 827 828 assert(sepadata != NULL); 829 830 /* free memory and reset data */ 831 if( sepadata->isintsolavailable ) 832 { 833 SCIPfreeBlockMemoryArray(scip, &sepadata->nlrowsidx, sepadata->nlrowssize); 834 SCIPfreeBlockMemoryArray(scip, &sepadata->convexsides, sepadata->nlrowssize); 835 SCIPfreeBlockMemoryArray(scip, &sepadata->nlrows, sepadata->nlrowssize); 836 SCIP_CALL( SCIPfreeSol(scip, &sepadata->intsol) ); 837 838 sepadata->nnlrows = 0; 839 sepadata->nnlrowsidx = 0; 840 sepadata->nlrowssize = 0; 841 sepadata->isintsolavailable = FALSE; 842 } 843 assert(sepadata->nnlrows == 0); 844 assert(sepadata->nnlrowsidx == 0); 845 assert(sepadata->nlrowssize == 0); 846 assert(sepadata->isintsolavailable == FALSE); 847 848 sepadata->skipsepa = FALSE; 849 850 return SCIP_OKAY; 851 } 852 853 854 /** LP solution separation method of separator */ 855 static 856 SCIP_DECL_SEPAEXECLP(sepaExeclpGauge) 857 { /*lint --e{715}*/ 858 SCIP_SEPADATA* sepadata; 859 SCIP_SOL* lpsol; 860 int i; 861 862 assert(scip != NULL); 863 assert(sepa != NULL); 864 865 sepadata = SCIPsepaGetData(sepa); 866 867 assert(sepadata != NULL); 868 869 *result = SCIP_DIDNOTRUN; 870 871 /* do not run if there is no interesting convex relaxation (with at least two nonlinear convex constraint) */ 872 if( sepadata->skipsepa ) 873 { 874 SCIPdebugMsg(scip, "not running because convex relaxation is uninteresting\n"); 875 return SCIP_OKAY; 876 } 877 878 /* do not run if SCIP has not constructed an NLP */ 879 if( !SCIPisNLPConstructed(scip) ) 880 { 881 SCIPdebugMsg(scip, "NLP not constructed, skipping gauge separator\n"); 882 return SCIP_OKAY; 883 } 884 885 /* do not run if SCIP has no way of solving nonlinear problems */ 886 if( SCIPgetNNlpis(scip) == 0 ) 887 { 888 SCIPdebugMsg(scip, "Skip gauge separator: no nlpi and SCIP can't solve nonlinear problems without a nlpi\n"); 889 return SCIP_OKAY; 890 } 891 892 /* if we don't have an interior point compute one; if we fail to compute one, then separator will not be run again; 893 * otherwise, we also store the convex nlrows in sepadata 894 */ 895 if( !sepadata->isintsolavailable ) 896 { 897 /* @todo: one could store the convex nonlinear rows inside computeInteriorPoint */ 898 SCIP_CALL( computeInteriorPoint(scip, sepadata) ); 899 assert(sepadata->skipsepa || sepadata->isintsolavailable); 900 901 if( sepadata->skipsepa ) 902 return SCIP_OKAY; 903 904 SCIP_CALL( storeNonlinearConvexNlrows(scip, sepadata, SCIPgetNLPNlRows(scip), SCIPgetNNLPNlRows(scip)) ); 905 } 906 907 #ifdef SCIP_DISABLED_CODE 908 /* get interior point: try to compute an interior point, otherwise use primal solution, otherwise use NLP solution */ 909 /* @todo: - decide order: 910 * - we can also use convex combination of solutions; there is a function SCIPvarGetAvgSol! 911 * - can add an event handler to only update when a new solution has been found 912 */ 913 if( !sepadata->isintsolavailable ) 914 { 915 if( SCIPgetNSols(scip) > 0 ) 916 { 917 SCIPdebugMsg(scip, "Using current primal solution as interior point!\n"); 918 SCIP_CALL( SCIPcreateSolCopy(scip, &sepadata->intsol, SCIPgetBestSol(scip)) ); 919 sepadata->isintsolavailable = TRUE; 920 } 921 else if( SCIPnlpGetSolstat(scip) <= SCIP_NLPSOLSTAT_FEASIBLE ) 922 { 923 SCIPdebugMsg(scip, "Using NLP solution as interior point!\n"); 924 SCIP_CALL( SCIPcreateNLPSol(scip, &sepadata->intsol, NULL) ); 925 sepadata->isintsolavailable = TRUE; 926 } 927 else 928 { 929 SCIPdebugMsg(scip, "We couldn't find an interior point, don't have a feasible nor an NLP solution; skip separator\n"); 930 return SCIP_OKAY; 931 } 932 } 933 #endif 934 935 /* store lp sol (or pseudo sol when lp is not solved) to be able to use it to compute nlrows' activities */ 936 SCIP_CALL( SCIPcreateCurrentSol(scip, &lpsol, NULL) ); 937 938 /* store indices of relevant constraints, ie, the ones that violate the lp sol */ 939 sepadata->nnlrowsidx = 0; 940 for( i = 0; i < sepadata->nnlrows; i++ ) 941 { 942 SCIP_NLROW* nlrow; 943 SCIP_Real activity; 944 945 nlrow = sepadata->nlrows[i]; 946 947 SCIP_CALL( SCIPgetNlRowSolActivity(scip, nlrow, lpsol, &activity) ); 948 949 if( sepadata->convexsides[i] == RHS ) 950 { 951 assert(!SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrow))); 952 953 if( activity - SCIPnlrowGetRhs(nlrow) < VIOLATIONFAC * SCIPfeastol(scip) ) 954 continue; 955 } 956 else 957 { 958 assert(sepadata->convexsides[i] == LHS); 959 assert(!SCIPisInfinity(scip, -SCIPnlrowGetLhs(nlrow))); 960 961 if( SCIPnlrowGetLhs(nlrow) - activity < VIOLATIONFAC * SCIPfeastol(scip) ) 962 continue; 963 } 964 965 sepadata->nlrowsidx[sepadata->nnlrowsidx] = i; 966 ++(sepadata->nnlrowsidx); 967 } 968 969 /* separate only if there are violated nlrows */ 970 SCIPdebugMsg(scip, "there are %d violated nlrows\n", sepadata->nnlrowsidx); 971 if( sepadata->nnlrowsidx > 0 ) 972 { 973 SCIP_CALL( separateCuts(scip, sepa, lpsol, result) ); 974 } 975 976 /* free lpsol */ 977 SCIP_CALL( SCIPfreeSol(scip, &lpsol) ); 978 979 return SCIP_OKAY; 980 } 981 982 983 /* 984 * separator specific interface methods 985 */ 986 987 /** creates the gauge separator and includes it in SCIP */ 988 SCIP_RETCODE SCIPincludeSepaGauge( 989 SCIP* scip /**< SCIP data structure */ 990 ) 991 { 992 SCIP_SEPADATA* sepadata; 993 SCIP_SEPA* sepa; 994 995 /* create gauge separator data */ 996 SCIP_CALL( SCIPallocBlockMemory(scip, &sepadata) ); 997 998 /* this sets all data in sepadata to 0 */ 999 BMSclearMemory(sepadata); 1000 1001 /* include separator */ 1002 SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST, 1003 SEPA_USESSUBSCIP, SEPA_DELAY, 1004 sepaExeclpGauge, NULL, 1005 sepadata) ); 1006 1007 assert(sepa != NULL); 1008 1009 /* set non fundamental callbacks via setter functions */ 1010 SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeGauge) ); 1011 SCIP_CALL( SCIPsetSepaExitsol(scip, sepa, sepaExitsolGauge) ); 1012 1013 /* add gauge separator parameters */ 1014 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nlpiterlimit", 1015 "iteration limit of NLP solver; 0 for no limit", 1016 &sepadata->nlpiterlimit, TRUE, DEFAULT_NLPITERLIM, 0, INT_MAX, NULL, NULL) ); 1017 1018 return SCIP_OKAY; 1019 } 1020