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_convexproj.c 26 * @ingroup DEFPLUGINS_SEPA 27 * @brief convexproj separator 28 * @author Felipe Serrano 29 * 30 * @todo should separator only be run when SCIPallColsInLP is true? 31 * @todo check if it makes sense to implement the copy callback 32 * @todo add SCIPisStopped(scip) to the condition of time consuming loops 33 */ 34 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 35 36 #include <assert.h> 37 #include <string.h> 38 39 #include "blockmemshell/memory.h" 40 #include "scip/scip_expr.h" 41 #include "scip/scip_nlpi.h" 42 #include "scip/expr_varidx.h" 43 #include "scip/expr_pow.h" 44 #include "scip/expr_sum.h" 45 #include "scip/pub_message.h" 46 #include "scip/pub_misc.h" 47 #include "scip/pub_nlp.h" 48 #include "scip/pub_sepa.h" 49 #include "scip/pub_var.h" 50 #include "scip/scip_cut.h" 51 #include "scip/scip_general.h" 52 #include "scip/scip_lp.h" 53 #include "scip/scip_mem.h" 54 #include "scip/scip_message.h" 55 #include "scip/scip_nlp.h" 56 #include "scip/scip_numerics.h" 57 #include "scip/scip_param.h" 58 #include "scip/scip_prob.h" 59 #include "scip/scip_sepa.h" 60 #include "scip/scip_sol.h" 61 #include "scip/scip_solvingstats.h" 62 #include "scip/scip_timing.h" 63 #include "scip/scip_tree.h" 64 #include "scip/sepa_convexproj.h" 65 #include <string.h> 66 67 68 #define SEPA_NAME "convexproj" 69 #define SEPA_DESC "separate at projection of point onto convex region" 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 TRUE /**< should separation method be delayed, if other separators found cuts? */ 75 76 #define DEFAULT_MAXDEPTH -1 /**< maximum depth at which the separator is applied; -1 means no limit */ 77 #define DEFAULT_NLPITERLIM 250 /**< default NLP iteration limit */ 78 79 #define VIOLATIONFAC 100 /**< points regarded violated if max violation > VIOLATIONFAC*SCIPfeastol() */ 80 81 /* 82 * Data structures 83 */ 84 85 /** side that makes an nlrow convex */ 86 enum ConvexSide 87 { 88 LHS = 0, /**< left hand side */ 89 RHS = 1 /**< right hand side */ 90 }; 91 typedef enum ConvexSide CONVEXSIDE; 92 93 /** separator data 94 * it keeps the nlpi which represents the projection problem (see sepa_convexproj.h); it also keeps the convex nlrows 95 * and the side which actually makes them convex; when separating, we use the nlpi to compute the projection and then 96 * the convex nlrows to compute the actual gradient cuts */ 97 struct SCIP_SepaData 98 { 99 SCIP_NLPI* nlpi; /**< nlpi used to create the nlpi problem */ 100 SCIP_NLPIPROBLEM* nlpiprob; /**< nlpi problem representing the convex NLP relaxation */ 101 SCIP_VAR** nlpivars; /**< array containing all variables of the nlpi */ 102 SCIP_HASHMAP* var2nlpiidx; /**< mapping between variables and nlpi indices */ 103 int nlpinvars; /**< total number of nlpi variables */ 104 105 SCIP_Bool skipsepa; /**< should separator be skipped? */ 106 107 SCIP_NLROW** nlrows; /**< convex nlrows */ 108 CONVEXSIDE* convexsides; /**< which sides make the nlrows convex */ 109 SCIP_Real* constraintviolation;/**< array storing the violation of constraint by current solution; 0.0 if it is not violated */ 110 int nnlrows; /**< total number of nlrows */ 111 int nlrowssize; /**< memory allocated for nlrows, convexsides and nlrowsidx */ 112 113 /* parameter */ 114 int nlpiterlimit; /**< iteration limit of NLP solver; 0 for no limit */ 115 int maxdepth; /**< maximal depth at which the separator is applied */ 116 117 int ncuts; /**< number of cuts generated */ 118 }; 119 120 121 /* 122 * Local methods 123 */ 124 125 /** clears the sepadata data */ 126 static 127 SCIP_RETCODE sepadataClear( 128 SCIP* scip, /**< SCIP data structure */ 129 SCIP_SEPADATA* sepadata /**< separator data */ 130 ) 131 { 132 assert(sepadata != NULL); 133 134 /* nlrowssize gets allocated first and then its decided whether to create the nlpiprob */ 135 if( sepadata->nlrowssize > 0 ) 136 { 137 SCIPfreeBlockMemoryArray(scip, &sepadata->constraintviolation, sepadata->nlrowssize); 138 SCIPfreeBlockMemoryArray(scip, &sepadata->convexsides, sepadata->nlrowssize); 139 SCIPfreeBlockMemoryArray(scip, &sepadata->nlrows, sepadata->nlrowssize); 140 sepadata->nlrowssize = 0; 141 } 142 143 if( sepadata->nlpiprob != NULL ) 144 { 145 assert(sepadata->nlpi != NULL); 146 147 SCIPfreeBlockMemoryArray(scip, &sepadata->nlpivars, sepadata->nlpinvars); 148 149 SCIPhashmapFree(&sepadata->var2nlpiidx); 150 SCIP_CALL( SCIPfreeNlpiProblem(scip, sepadata->nlpi, &sepadata->nlpiprob) ); 151 152 sepadata->nlpinvars = 0; 153 sepadata->nnlrows = 0; 154 } 155 assert(sepadata->nlpinvars == 0); 156 assert(sepadata->nnlrows == 0); 157 assert(sepadata->nlrowssize == 0); 158 159 sepadata->skipsepa = FALSE; 160 161 return SCIP_OKAY; 162 } 163 164 /** computes gradient cut (linearization) of nlrow at projection */ 165 static 166 SCIP_RETCODE generateCut( 167 SCIP* scip, /**< SCIP data structure */ 168 SCIP_SEPA* sepa, /**< the cut separator itself */ 169 SCIP_SOL* projection, /**< point where we compute gradient cut */ 170 SCIP_NLROW* nlrow, /**< constraint for which we generate gradient cut */ 171 CONVEXSIDE convexside, /**< which side makes the nlrow convex */ 172 SCIP_Real activity, /**< activity of constraint at projection */ 173 SCIP_EXPRITER* exprit, /**< expression iterator that can be used */ 174 SCIP_ROW** row /**< storage for cut */ 175 ) 176 { 177 char rowname[SCIP_MAXSTRLEN]; 178 SCIP_SEPADATA* sepadata; 179 SCIP_Real gradx0; /* <grad f(x_0), x_0> */ 180 SCIP_EXPR* expr; 181 int i; 182 183 assert(scip != NULL); 184 assert(sepa != NULL); 185 assert(nlrow != NULL); 186 assert(row != NULL); 187 188 sepadata = SCIPsepaGetData(sepa); 189 190 assert(sepadata != NULL); 191 192 gradx0 = 0.0; 193 194 /* an nlrow has a linear part and expression; ideally one would just build the gradient but we 195 * do not know if the different parts share variables or not, so we can't just build the gradient; for this reason 196 * we create the row right away and compute the gradients of each part independently and add them to the row; the 197 * row takes care to add coeffs corresponding to the same variable when they appear in different parts of the nlrow 198 * NOTE: a gradient cut is globally valid whenever the constraint from which it is deduced is globally valid; since 199 * we build the convex relaxation using only globally valid constraints, the cuts are globally valid 200 */ 201 (void) SCIPsnprintf(rowname, SCIP_MAXSTRLEN, "proj_cut_%s_%u", SCIPnlrowGetName(nlrow), ++(sepadata->ncuts)); 202 SCIP_CALL( SCIPcreateEmptyRowSepa(scip, row, sepa, rowname, -SCIPinfinity(scip), SCIPinfinity(scip), TRUE, FALSE , 203 TRUE) ); 204 205 SCIP_CALL( SCIPcacheRowExtensions(scip, *row) ); 206 207 /* linear part */ 208 for( i = 0; i < SCIPnlrowGetNLinearVars(nlrow); i++ ) 209 { 210 gradx0 += SCIPgetSolVal(scip, projection, SCIPnlrowGetLinearVars(nlrow)[i]) * SCIPnlrowGetLinearCoefs(nlrow)[i]; 211 SCIP_CALL( SCIPaddVarToRow(scip, *row, SCIPnlrowGetLinearVars(nlrow)[i], SCIPnlrowGetLinearCoefs(nlrow)[i]) ); 212 } 213 214 expr = SCIPnlrowGetExpr(nlrow); 215 assert(expr != NULL); 216 217 SCIP_CALL( SCIPevalExprGradient(scip, expr, projection, 0L) ); 218 219 SCIP_CALL( SCIPexpriterInit(exprit, expr, SCIP_EXPRITER_DFS, FALSE) ); 220 for( ; !SCIPexpriterIsEnd(exprit); expr = SCIPexpriterGetNext(exprit) ) /*lint !e441*/ /*lint !e440*/ 221 { 222 SCIP_Real grad; 223 SCIP_VAR* var; 224 225 if( !SCIPisExprVar(scip, expr) ) 226 continue; 227 228 grad = SCIPexprGetDerivative(expr); 229 var = SCIPgetVarExprVar(expr); 230 assert(var != NULL); 231 232 gradx0 += grad * SCIPgetSolVal(scip, projection, var); 233 SCIP_CALL( SCIPaddVarToRow(scip, *row, var, grad) ); 234 } 235 236 SCIP_CALL( SCIPflushRowExtensions(scip, *row) ); 237 238 SCIPdebugPrintf("gradient: "); 239 SCIPdebug( SCIP_CALL( SCIPprintRow(scip, *row, NULL) ) ); 240 SCIPdebugPrintf("gradient dot x_0: %g\n", gradx0); 241 242 /* gradient cut is f(x_0) - <grad f(x_0), x_0> + <grad f(x_0), x> <= rhs or >= lhs */ 243 if( convexside == RHS ) 244 { 245 assert(!SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrow))); 246 SCIP_CALL( SCIPchgRowRhs(scip, *row, SCIPnlrowGetRhs(nlrow) - activity + gradx0) ); 247 } 248 else 249 { 250 assert(convexside == LHS); 251 assert(!SCIPisInfinity(scip, -SCIPnlrowGetLhs(nlrow))); 252 SCIP_CALL( SCIPchgRowLhs(scip, *row, SCIPnlrowGetLhs(nlrow) - activity + gradx0) ); 253 } 254 255 SCIPdebugPrintf("gradient cut: "); 256 SCIPdebug( SCIP_CALL( SCIPprintRow(scip, *row, NULL) ) ); 257 258 return SCIP_OKAY; 259 } 260 261 /** set quadratic part of objective function: \f$ \sum_i x_i^2 \f$ 262 * 263 * the objective function is \f$ ||x - x_0||^2 \f$, 264 * where \f$ x_0 \f$ is the point to separate; the only part that changes is the term \f$ -2 \langle x_0, x \rangle \f$ 265 * which is linear and is set every time we want to separate a point, see separateCuts() 266 */ 267 static 268 SCIP_RETCODE setQuadraticObj( 269 SCIP* scip, /**< SCIP data structure */ 270 SCIP_SEPADATA* sepadata /**< the cut separator data */ 271 ) 272 { 273 SCIP_EXPR* exprsum; 274 SCIP_EXPR** exprspow; 275 int i; 276 277 assert(scip != NULL); 278 assert(sepadata != NULL); 279 assert(sepadata->nlpi != NULL); 280 assert(sepadata->nlpiprob != NULL); 281 assert(sepadata->var2nlpiidx != NULL); 282 assert(sepadata->nlpinvars > 0); 283 284 SCIP_CALL( SCIPallocBufferArray(scip, &exprspow, sepadata->nlpinvars) ); 285 for( i = 0; i < sepadata->nlpinvars; i++ ) 286 { 287 SCIP_VAR* var; 288 SCIP_EXPR* varexpr; 289 290 var = sepadata->nlpivars[i]; 291 assert(SCIPhashmapExists(sepadata->var2nlpiidx, (void*)var) ); 292 293 SCIP_CALL( SCIPcreateExprVaridx(scip, &varexpr, SCIPhashmapGetImageInt(sepadata->var2nlpiidx, (void*)var), NULL, NULL) ); 294 SCIP_CALL( SCIPcreateExprPow(scip, &exprspow[i], varexpr, 2.0, NULL, NULL) ); 295 SCIP_CALL( SCIPreleaseExpr(scip, &varexpr) ); 296 } 297 298 SCIP_CALL( SCIPcreateExprSum(scip, &exprsum, sepadata->nlpinvars, exprspow, NULL, 0.0, NULL, NULL) ); 299 300 /* set quadratic part of objective function */ 301 SCIP_CALL( SCIPsetNlpiObjective(scip, sepadata->nlpi, sepadata->nlpiprob, 0, NULL, NULL, exprsum, 0.0) ); 302 303 /* free memory */ 304 SCIP_CALL( SCIPreleaseExpr(scip, &exprsum) ); 305 for( i = sepadata->nlpinvars-1; i >= 0; --i ) 306 { 307 SCIP_CALL( SCIPreleaseExpr(scip, &exprspow[i]) ); 308 } 309 SCIPfreeBufferArray(scip, &exprspow); 310 311 return SCIP_OKAY; 312 } 313 314 /** projects sol onto convex relaxation (stored in sepadata) and tries to generate gradient cuts at the projection 315 * 316 * it generates cuts only for the constraints that were violated by the LP solution and are now active or still 317 * violated (in case we don't solve to optimality). 318 * @todo: store a feasible solution if one is found to use as warmstart 319 */ 320 static 321 SCIP_RETCODE separateCuts( 322 SCIP* scip, /**< SCIP data structure */ 323 SCIP_SEPA* sepa, /**< the cut separator itself */ 324 SCIP_SOL* sol, /**< solution that should be separated */ 325 SCIP_RESULT* result /**< pointer to store the result of the separation call */ 326 ) 327 { 328 SCIP_SEPADATA* sepadata; 329 SCIP_SOL* projection; 330 SCIP_Real* linvals; 331 SCIP_Real* nlpisol; 332 int nlpinvars; 333 int i; 334 int* lininds; 335 SCIP_Bool nlpunstable; 336 SCIP_EXPRITER* exprit; 337 338 nlpunstable = FALSE; 339 340 assert(sepa != NULL); 341 342 sepadata = SCIPsepaGetData(sepa); 343 344 assert(result != NULL); 345 assert(sepadata != NULL); 346 assert(sepadata->nnlrows > 0); 347 assert(sepadata->nlpi != NULL); 348 assert(sepadata->nlpinvars > 0); 349 assert(sepadata->nlrows != NULL); 350 assert(sepadata->nlpiprob != NULL); 351 assert(sepadata->var2nlpiidx != NULL); 352 assert(sepadata->convexsides != NULL); 353 assert(sepadata->constraintviolation != NULL); 354 355 nlpinvars = sepadata->nlpinvars; 356 /* set linear part of objective function: \norm(x - x^0)^2 = \norm(x)^2 - \sum 2 * x_i * x^0_i + const 357 * we ignore the constant; x0 is `sol` 358 */ 359 SCIP_CALL( SCIPallocBufferArray(scip, &linvals, nlpinvars) ); 360 SCIP_CALL( SCIPallocBufferArray(scip, &lininds, nlpinvars) ); 361 for( i = 0; i < nlpinvars; i++ ) 362 { 363 SCIP_VAR* var; 364 365 var = sepadata->nlpivars[i]; 366 assert(SCIPhashmapExists(sepadata->var2nlpiidx, (void*)var) ); 367 368 lininds[i] = SCIPhashmapGetImageInt(sepadata->var2nlpiidx, (void*)var); 369 linvals[i] = - 2.0 * SCIPgetSolVal(scip, sol, var); 370 371 /* if coefficient is too large, don't separate */ 372 if( SCIPisHugeValue(scip, REALABS(linvals[i])) ) 373 { 374 SCIPdebugMsg(scip, "Don't separate points too close to infinity\n"); 375 goto CLEANUP; 376 } 377 } 378 379 /* set linear part of objective function */ 380 SCIP_CALL( SCIPchgNlpiLinearCoefs(scip, sepadata->nlpi, sepadata->nlpiprob, -1, nlpinvars, lininds, linvals) ); 381 382 /* compute the projection onto the convex NLP relaxation */ 383 SCIP_CALL( SCIPsolveNlpi(scip, sepadata->nlpi, sepadata->nlpiprob, 384 .iterlimit = sepadata->nlpiterlimit > 0 ? sepadata->nlpiterlimit : INT_MAX, 385 .feastol = SCIPfeastol(scip) / 10.0, /* use tighter tolerances for the NLP solver */ 386 .opttol = MAX(SCIPfeastol(scip), SCIPdualfeastol(scip))) ); /*lint !e666*/ 387 SCIPdebugMsg(scip, "NLP solstat = %d\n", SCIPgetNlpiSolstat(scip, sepadata->nlpi, sepadata->nlpiprob)); 388 389 /* if solution is feasible, add cuts */ 390 switch( SCIPgetNlpiSolstat(scip, sepadata->nlpi, sepadata->nlpiprob) ) 391 { 392 case SCIP_NLPSOLSTAT_GLOBOPT: 393 case SCIP_NLPSOLSTAT_LOCOPT: 394 /* @todo: if solution is optimal, we might as well add the cut <x - P(x_0), x_0 - P(x_0)> <= 0 395 * even though this cut is implied by all the gradient cuts of the rows active at the projection, 396 * we do not add them all (only the gradient cuts of constraints that violated the LP solution */ 397 case SCIP_NLPSOLSTAT_FEASIBLE: 398 /* generate cuts for violated constraints (at sol) that are active or still violated at the projection, since 399 * a suboptimal solution or numerical issues could give a solution of the projection problem where constraints 400 * are not active; if the solution of the projection problem is in the interior of the region, we do nothing 401 */ 402 403 /* get solution: build SCIP_SOL out of nlpi sol */ 404 SCIP_CALL( SCIPgetNlpiSolution(scip, sepadata->nlpi, sepadata->nlpiprob, &nlpisol, NULL, NULL, NULL, NULL) ); 405 assert(nlpisol != NULL); 406 407 SCIP_CALL( SCIPcreateSol(scip, &projection, NULL) ); 408 for( i = 0; i < nlpinvars; i++ ) 409 { 410 SCIP_VAR* var; 411 412 var = sepadata->nlpivars[i]; 413 assert(SCIPhashmapExists(sepadata->var2nlpiidx, (void*)var) ); 414 415 SCIP_CALL( SCIPsetSolVal(scip, projection, var, 416 nlpisol[SCIPhashmapGetImageInt(sepadata->var2nlpiidx, (void *)var)]) ); 417 } 418 SCIPdebug( SCIPprintSol(scip, projection, NULL, TRUE) ); 419 420 /** @todo this could just be created inside generateCut and the extra argument removed */ 421 SCIP_CALL( SCIPcreateExpriter(scip, &exprit) ); 422 423 /* check for active or violated constraints */ 424 for( i = 0; i < sepadata->nnlrows; ++i ) 425 { 426 SCIP_NLROW* nlrow; 427 CONVEXSIDE convexside; 428 SCIP_Real activity; 429 430 /* ignore constraints that are not violated by `sol` */ 431 if( SCIPisFeasZero(scip, sepadata->constraintviolation[i]) ) 432 continue; 433 434 convexside = sepadata->convexsides[i]; 435 nlrow = sepadata->nlrows[i]; 436 assert(nlrow != NULL); 437 438 /* check for currently active constraints at projected point */ 439 SCIP_CALL( SCIPgetNlRowSolActivity(scip, nlrow, projection, &activity) ); 440 441 SCIPdebugMsg(scip, "NlRow activity at nlpi solution: %g <= %g <= %g\n", SCIPnlrowGetLhs(nlrow), activity, 442 SCIPnlrowGetRhs(nlrow) ); 443 444 /* if nlrow is active or violates the projection, build gradient cut at projection */ 445 if( (convexside == RHS && SCIPisFeasGE(scip, activity, SCIPnlrowGetRhs(nlrow))) 446 || (convexside == LHS && SCIPisFeasLE(scip, activity, SCIPnlrowGetLhs(nlrow))) ) 447 { 448 SCIP_ROW* row; 449 450 SCIP_CALL( generateCut(scip, sepa, projection, nlrow, convexside, activity, exprit, 451 &row) ); 452 453 SCIPdebugMsg(scip, "active or violated nlrow: (sols vio: %e)\n", sepadata->constraintviolation[i]); 454 SCIPdebug( SCIP_CALL( SCIPprintNlRow(scip, nlrow, NULL) ) ); 455 SCIPdebugMsg(scip, "cut with efficacy %g generated\n", SCIPgetCutEfficacy(scip, sol, row)); 456 SCIPdebug( SCIPprintRow(scip, row, NULL) ); 457 458 /* add cut if it is efficacious for the point we want to separate (sol) */ 459 if( SCIPisCutEfficacious(scip, sol, row) ) 460 { 461 SCIP_Bool infeasible; 462 463 SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) ); 464 465 if( infeasible ) 466 { 467 *result = SCIP_CUTOFF; 468 SCIP_CALL( SCIPreleaseRow(scip, &row) ); 469 break; 470 } 471 else 472 { 473 *result = SCIP_SEPARATED; 474 } 475 } 476 477 /* release the row */ 478 SCIP_CALL( SCIPreleaseRow(scip, &row) ); 479 } 480 } 481 482 SCIPfreeExpriter(&exprit); 483 484 #ifdef SCIP_DEBUG 485 { 486 SCIP_Real distance; 487 488 /* compute distance between LP sol and its projection (only makes sense when it is optimal) */ 489 distance = 0.0; 490 for( i = 0; i < SCIPgetNNLPVars(scip); ++i ) 491 { 492 SCIP_VAR* var; 493 494 var = SCIPgetNLPVars(scip)[i]; 495 assert(var != NULL); 496 497 /* assert NLP solution is within the bounds of the variable (only make sense when sol is optimal) */ 498 if( !SCIPisInfinity(scip, -SCIPvarGetLbLocal(var)) ) 499 assert(SCIPisFeasLE(scip, SCIPvarGetLbLocal(var), SCIPvarGetNLPSol(var))); 500 if( !SCIPisInfinity(scip, SCIPvarGetUbLocal(var)) ) 501 assert(SCIPisFeasLE(scip, SCIPvarGetNLPSol(var), SCIPvarGetUbLocal(var))); 502 503 /*SCIPdebugMsg(scip, "NLP sol (LP sol): %s = %f (%g)\n", SCIPvarGetName(var), 504 * SCIPvarGetNLPSol(var), SCIPgetSolVal(scip, sol, var)); 505 */ 506 507 distance += SQR( SCIPvarGetNLPSol(var) - SCIPgetSolVal(scip, sol, var) ); 508 } 509 510 SCIPdebugMsg(scip, "NLP objval: %e, distance: %e\n", SCIPgetNLPObjval(scip), distance); 511 } 512 #endif 513 514 /* free solution */ 515 SCIP_CALL( SCIPfreeSol(scip, &projection) ); 516 break; 517 518 case SCIP_NLPSOLSTAT_GLOBINFEASIBLE: 519 case SCIP_NLPSOLSTAT_LOCINFEASIBLE: 520 /* fallthrough; 521 * @todo: write what it means to be locinfeasible and why it can't be used to cutoff the node */ 522 case SCIP_NLPSOLSTAT_UNKNOWN: 523 /* unknown... assume numerical issues */ 524 nlpunstable = TRUE; 525 break; 526 527 case SCIP_NLPSOLSTAT_UNBOUNDED: 528 default: 529 SCIPerrorMessage("Projection NLP is not unbounded by construction, should not get here!\n"); 530 SCIPABORT(); 531 nlpunstable = TRUE; 532 } 533 534 /* if nlp is detected to be unstable, don't try to separate again */ 535 if( nlpunstable ) 536 { 537 /* @todo: maybe change objective function to \sum [(x_i - x_i^*)/max(|x_i^*|, 1)]^2 538 * or some other scaling when unstable and try again. 539 * maybe free it here */ 540 sepadata->skipsepa = TRUE; 541 } 542 543 /* reset objective */ 544 BMSclearMemoryArray(linvals, nlpinvars); 545 SCIP_CALL( SCIPchgNlpiLinearCoefs(scip, sepadata->nlpi, sepadata->nlpiprob, -1, nlpinvars, lininds, linvals) ); 546 547 CLEANUP: 548 /* free memory */ 549 SCIPfreeBufferArray(scip, &lininds); 550 SCIPfreeBufferArray(scip, &linvals); 551 552 return SCIP_OKAY; 553 } 554 555 /** computes the violation and maximum violation of the convex nlrows stored in sepadata wrt sol */ 556 static 557 SCIP_RETCODE computeMaxViolation( 558 SCIP* scip, /**< SCIP data structure */ 559 SCIP_SEPADATA* sepadata, /**< separator data */ 560 SCIP_SOL* sol, /**< solution that should be separated */ 561 SCIP_Real* maxviolation /**< buffer to store maximum violation */ 562 ) 563 { 564 SCIP_NLROW* nlrow; 565 int i; 566 567 assert(sepadata != NULL); 568 assert(sepadata->nnlrows > 0); 569 assert(sepadata->nlrows != NULL); 570 assert(sepadata->convexsides != NULL); 571 assert(sepadata->constraintviolation != NULL); 572 573 *maxviolation = 0.0; 574 for( i = 0; i < sepadata->nnlrows; i++ ) 575 { 576 SCIP_Real activity; 577 SCIP_Real violation; 578 579 nlrow = sepadata->nlrows[i]; 580 581 /* get activity of nlrow */ 582 SCIP_CALL( SCIPgetNlRowSolActivity(scip, nlrow, sol, &activity) ); 583 584 /* violation = max{activity - rhs, 0.0} when convex and max{lhs - activity, 0.0} when concave */ 585 if( sepadata->convexsides[i] == RHS ) 586 { 587 assert(SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_CONVEX); 588 assert(!SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrow))); 589 590 violation = activity - SCIPnlrowGetRhs(nlrow); 591 sepadata->constraintviolation[i] = MAX(violation, 0.0); 592 } 593 if( sepadata->convexsides[i] == LHS ) 594 { 595 assert(SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_CONCAVE); 596 assert(!SCIPisInfinity(scip, -SCIPnlrowGetLhs(nlrow))); 597 598 violation = SCIPnlrowGetLhs(nlrow) - activity; 599 sepadata->constraintviolation[i] = MAX(violation, 0.0); 600 } 601 602 /* compute maximum */ 603 if( *maxviolation < sepadata->constraintviolation[i] ) 604 *maxviolation = sepadata->constraintviolation[i]; 605 } 606 607 SCIPdebugMsg(scip, "Maximum violation %g\n", *maxviolation); 608 609 return SCIP_OKAY; 610 } 611 612 613 /** stores, from the constraints represented by nlrows, the nonlinear convex ones in sepadata */ 614 static 615 SCIP_RETCODE storeNonlinearConvexNlrows( 616 SCIP* scip, /**< SCIP data structure */ 617 SCIP_SEPADATA* sepadata, /**< separator data */ 618 SCIP_NLROW** nlrows, /**< nlrows from which to store convex ones */ 619 int nnlrows /**< number of nlrows */ 620 ) 621 { 622 int i; 623 624 assert(scip != NULL); 625 assert(sepadata != NULL); 626 627 SCIPdebugMsg(scip, "storing convex nlrows\n"); 628 629 sepadata->nlrowssize = nnlrows; 630 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(sepadata->nlrows), nnlrows) ); 631 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(sepadata->convexsides), nnlrows) ); 632 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(sepadata->constraintviolation), nnlrows) ); 633 634 /* count the number of nonlinear convex rows and store them */ 635 sepadata->nnlrows = 0; 636 for( i = 0; i < nnlrows; ++i ) 637 { 638 SCIP_NLROW* nlrow; 639 640 nlrow = nlrows[i]; 641 assert(nlrow != NULL); 642 643 /* linear case */ 644 if( SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_LINEAR || SCIPnlrowGetExpr(nlrow) == NULL ) 645 continue; 646 647 /* nonlinear case */ 648 if( !SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrow)) && SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_CONVEX ) 649 { 650 sepadata->convexsides[sepadata->nnlrows] = RHS; 651 sepadata->nlrows[sepadata->nnlrows] = nlrow; 652 ++(sepadata->nnlrows); 653 } 654 else if( !SCIPisInfinity(scip, -SCIPnlrowGetLhs(nlrow)) && SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_CONCAVE ) 655 { 656 sepadata->convexsides[sepadata->nnlrows] = LHS; 657 sepadata->nlrows[sepadata->nnlrows] = nlrow; 658 ++(sepadata->nnlrows); 659 } 660 } 661 662 return SCIP_OKAY; 663 } 664 665 666 /* 667 * Callback methods of separator 668 */ 669 670 671 /** destructor of separator to free user data (called when SCIP is exiting) */ 672 static 673 SCIP_DECL_SEPAFREE(sepaFreeConvexproj) 674 { /*lint --e{715}*/ 675 SCIP_SEPADATA* sepadata; 676 677 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0); 678 679 /* free separator data */ 680 sepadata = SCIPsepaGetData(sepa); 681 assert(sepadata != NULL); 682 683 SCIP_CALL( sepadataClear(scip, sepadata) ); 684 685 SCIPfreeBlockMemory(scip, &sepadata); 686 687 SCIPsepaSetData(sepa, NULL); 688 689 return SCIP_OKAY; 690 } 691 692 /** solving process deinitialization method of separator (called before branch and bound process data is freed) */ 693 static 694 SCIP_DECL_SEPAEXITSOL(sepaExitsolConvexproj) 695 { /*lint --e{715}*/ 696 SCIP_SEPADATA* sepadata; 697 698 assert(sepa != NULL); 699 700 sepadata = SCIPsepaGetData(sepa); 701 702 assert(sepadata != NULL); 703 704 SCIP_CALL( sepadataClear(scip, sepadata) ); 705 706 return SCIP_OKAY; 707 } 708 709 710 /** LP solution separation method of separator */ 711 static 712 SCIP_DECL_SEPAEXECLP(sepaExeclpConvexproj) 713 { /*lint --e{715}*/ 714 SCIP_Real maxviolation; 715 SCIP_SOL* lpsol; 716 SCIP_SEPADATA* sepadata; 717 718 *result = SCIP_DIDNOTRUN; 719 720 sepadata = SCIPsepaGetData(sepa); 721 assert(sepadata != NULL); 722 723 /* do not run if there is no interesting convex relaxation (with at least one nonlinear convex constraint), 724 * or if we have found it to be numerically unstable 725 * @todo: should it be with at least 2 nonlinear convex constraints? 726 */ 727 if( sepadata->skipsepa ) 728 { 729 SCIPdebugMsg(scip, "not running because convex relaxation is uninteresting or numerically unstable\n"); 730 return SCIP_OKAY; 731 } 732 733 /* the separator needs an NLP solver */ 734 if( SCIPgetNNlpis(scip) == 0 ) 735 return SCIP_OKAY; 736 737 /* only call separator up to a maximum depth */ 738 if( sepadata->maxdepth >= 0 && depth > sepadata->maxdepth ) 739 return SCIP_OKAY; 740 741 /* only call separator, if we are not close to terminating */ 742 if( SCIPisStopped(scip) ) 743 return SCIP_OKAY; 744 745 /* do not run if SCIP does not have constructed an NLP */ 746 if( !SCIPisNLPConstructed(scip) ) 747 { 748 SCIPdebugMsg(scip, "NLP not constructed, skipping convex projection separator\n"); 749 return SCIP_OKAY; 750 } 751 752 /* recompute convex NLP relaxation if the variable set changed and we are still at the root node */ 753 if( sepadata->nlpiprob != NULL && SCIPgetNVars(scip) != sepadata->nlpinvars && SCIPgetDepth(scip) == 0 ) 754 { 755 SCIP_CALL( sepadataClear(scip, sepadata) ); 756 assert(sepadata->nlpiprob == NULL); 757 } 758 759 /* create or update convex NLP relaxation */ 760 if( sepadata->nlpiprob == NULL ) 761 { 762 /* store convex nonlinear constraints */ 763 SCIP_CALL( storeNonlinearConvexNlrows(scip, sepadata, SCIPgetNLPNlRows(scip), SCIPgetNNLPNlRows(scip)) ); 764 765 /* check that convex NLP relaxation is interesting (more than one nonlinear constraint) */ 766 if( sepadata->nnlrows < 1 ) 767 { 768 SCIPdebugMsg(scip, "convex relaxation uninteresting, don't run\n"); 769 sepadata->skipsepa = TRUE; 770 return SCIP_OKAY; 771 } 772 773 sepadata->nlpinvars = SCIPgetNVars(scip); 774 sepadata->nlpi = SCIPgetNlpis(scip)[0]; 775 assert(sepadata->nlpi != NULL); 776 777 SCIP_CALL( SCIPhashmapCreate(&sepadata->var2nlpiidx, SCIPblkmem(scip), sepadata->nlpinvars) ); 778 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &sepadata->nlpivars, SCIPgetVars(scip), sepadata->nlpinvars) ); /*lint !e666*/ 779 780 SCIP_CALL( SCIPcreateNlpiProblemFromNlRows(scip, sepadata->nlpi, &sepadata->nlpiprob, "convexproj-nlp", SCIPgetNLPNlRows(scip), SCIPgetNNLPNlRows(scip), 781 sepadata->var2nlpiidx, NULL, NULL, SCIPgetCutoffbound(scip), FALSE, TRUE) ); 782 783 /* add rows of the LP 784 * we do not sue the depth argument of the callback because we want to build a globally valid initia lrelaxation 785 */ 786 if( SCIPgetDepth(scip) == 0 ) 787 { 788 SCIP_CALL( SCIPaddNlpiProblemRows(scip, sepadata->nlpi, sepadata->nlpiprob, sepadata->var2nlpiidx, 789 SCIPgetLPRows(scip), SCIPgetNLPRows(scip)) ); 790 } 791 792 /* set quadratic part of objective function */ 793 SCIP_CALL( setQuadraticObj(scip, sepadata) ); 794 } 795 else 796 { 797 SCIP_CALL( SCIPupdateNlpiProblem(scip, sepadata->nlpi, sepadata->nlpiprob, sepadata->var2nlpiidx, 798 sepadata->nlpivars, sepadata->nlpinvars, SCIPgetCutoffbound(scip)) ); 799 } 800 801 /* assert that the lp solution satisfies the cutoff bound; if this fails then we shouldn't have a cutoff bound in the 802 * nlpi, since then the projection could be in the interior of the actual convex relaxation */ 803 assert(SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL || 804 SCIPisLE(scip, SCIPgetLPObjval(scip), SCIPgetCutoffbound(scip))); 805 806 /* get current sol: LP or pseudo solution if LP sol is not available */ 807 SCIP_CALL( SCIPcreateCurrentSol(scip, &lpsol, NULL) ); 808 809 /* do not run if current solution's violation is small */ 810 SCIP_CALL( computeMaxViolation(scip, sepadata, lpsol, &maxviolation) ); 811 if( maxviolation < VIOLATIONFAC * SCIPfeastol(scip) ) 812 { 813 SCIPdebugMsg(scip, "solution doesn't violate constraints enough, do not separate\n"); 814 SCIP_CALL( SCIPfreeSol(scip, &lpsol) ); 815 return SCIP_OKAY; 816 } 817 818 /* run the separator */ 819 *result = SCIP_DIDNOTFIND; 820 821 /* separateCuts computes the projection and then gradient cuts on each constraint that was originally violated */ 822 SCIP_CALL( separateCuts(scip, sepa, lpsol, result) ); 823 824 /* free memory */ 825 SCIP_CALL( SCIPfreeSol(scip, &lpsol) ); 826 827 return SCIP_OKAY; 828 } 829 830 831 /* 832 * separator specific interface methods 833 */ 834 835 /** creates the convexproj separator and includes it in SCIP */ 836 SCIP_RETCODE SCIPincludeSepaConvexproj( 837 SCIP* scip /**< SCIP data structure */ 838 ) 839 { 840 SCIP_SEPADATA* sepadata; 841 SCIP_SEPA* sepa; 842 843 /* create convexproj separator data */ 844 SCIP_CALL( SCIPallocBlockMemory(scip, &sepadata) ); 845 846 /* this sets all data in sepadata to 0 */ 847 BMSclearMemory(sepadata); 848 849 /* include separator */ 850 SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST, 851 SEPA_USESSUBSCIP, SEPA_DELAY, 852 sepaExeclpConvexproj, NULL, 853 sepadata) ); 854 assert(sepa != NULL); 855 856 /* set non fundamental callbacks via setter functions */ 857 SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeConvexproj) ); 858 SCIP_CALL( SCIPsetSepaExitsol(scip, sepa, sepaExitsolConvexproj) ); 859 860 /* add convexproj separator parameters */ 861 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxdepth", 862 "maximal depth at which the separator is applied (-1: unlimited)", 863 &sepadata->maxdepth, FALSE, DEFAULT_MAXDEPTH, -1, INT_MAX, NULL, NULL) ); 864 865 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nlpiterlimit", 866 "iteration limit of NLP solver; 0 for no limit", 867 &sepadata->nlpiterlimit, TRUE, DEFAULT_NLPITERLIM, 0, INT_MAX, NULL, NULL) ); 868 869 return SCIP_OKAY; 870 } 871