1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 2 /* */ 3 /* This file is part of the program and library */ 4 /* SCIP --- Solving Constraint Integer Programs */ 5 /* */ 6 /* Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB) */ 7 /* */ 8 /* Licensed under the Apache License, Version 2.0 (the "License"); */ 9 /* you may not use this file except in compliance with the License. */ 10 /* You may obtain a copy of the License at */ 11 /* */ 12 /* http://www.apache.org/licenses/LICENSE-2.0 */ 13 /* */ 14 /* Unless required by applicable law or agreed to in writing, software */ 15 /* distributed under the License is distributed on an "AS IS" BASIS, */ 16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */ 17 /* See the License for the specific language governing permissions and */ 18 /* limitations under the License. */ 19 /* */ 20 /* You should have received a copy of the Apache-2.0 license */ 21 /* along with SCIP; see the file LICENSE. If not visit scipopt.org. */ 22 /* */ 23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 24 25 /**@file nlhdlr_soc.c 26 * @ingroup DEFPLUGINS_NLHDLR 27 * @brief nonlinear handler for second order cone constraints 28 29 * @author Benjamin Mueller 30 * @author Felipe Serrano 31 * @author Fabian Wegscheider 32 * 33 * This is a nonlinear handler for second order cone constraints of the form 34 * 35 * \f[\sqrt{\sum_{i=1}^{n} (v_i^T x + \beta_i)^2} \leq v_{n+1}^T x + \beta_{n+1}.\f] 36 * 37 * Note that \f$v_i\f$, for \f$i \leq n\f$, could be 0, thus allowing a positive constant term inside the root. 38 * 39 * @todo test if it makes sense to only disaggregate when nterms > some parameter 40 * 41 */ 42 43 #include <string.h> 44 45 #include "scip/nlhdlr_soc.h" 46 #include "scip/cons_nonlinear.h" 47 #include "scip/expr_pow.h" 48 #include "scip/expr_sum.h" 49 #include "scip/expr_var.h" 50 #include "scip/debug.h" 51 #include "scip/pub_nlhdlr.h" 52 #include "scip/lapack_calls.h" 53 54 55 /* fundamental nonlinear handler properties */ 56 #define NLHDLR_NAME "soc" 57 #define NLHDLR_DESC "nonlinear handler for second-order cone structures" 58 #define NLHDLR_DETECTPRIORITY 100 /**< priority of the nonlinear handler for detection */ 59 #define NLHDLR_ENFOPRIORITY 100 /**< priority of the nonlinear handler for enforcement */ 60 #define DEFAULT_MINCUTEFFICACY 1e-5 /**< default value for parameter mincutefficacy */ 61 #define DEFAULT_COMPEIGENVALUES TRUE /**< default value for parameter compeigenvalues */ 62 63 /* 64 * Data structures 65 */ 66 67 /** nonlinear handler expression data. The data is structured in the following way: 68 * 69 * A 'term' is one of the arguments of the quadratic terms, i.e. \f$v_i^T x + beta_i\f$. 70 * The last term is always the one on the right-hand side. This means that nterms is 71 * equal to n+1 in the above description. 72 * 73 * - vars contains a list of all expressions which are treated as variables (no duplicates) 74 * - offsets contains the constants beta_i of each term 75 * - transcoefs contains the non-zero values of the transformation vectors v_i of each term 76 * - transcoefsidx contains for each entry of transcoefs the position of the respective variable in vars 77 * - termbegins contains the index at which the transcoefs of each term start, with a sentinel value 78 * - nterms is the total number of terms appearing on both sides 79 * - nvars is the total number of unique variables appearing (length of vars) 80 * 81 * Note that the numbers of nonzeroes in v_i is termbegins[i+1] - termbegins[i] and that 82 * the total number of entries in transcoefs and transcoefsidx is termbegins[nterms] 83 * 84 * The disaggregation is implicitly stored in the variables disvars and disrow. An SOC as 85 * described above is replaced by n smaller SOCs 86 * 87 * (v_i^T x + beta_i)^2 <= disvar_i * (v_{n+1}^T x + beta_{n+1}) 88 * 89 * and the row sum_i disvar_i <= v_{n+1}^T x + beta_{n+1}. 90 * 91 * The disaggregation only happens if we have more than 3 terms. 92 * 93 * Example: The constraint SQRT(5 + (3x - 4y + 2)^2 + y^2 + 7z^2) <= 5x - y - 1 94 * results in the following nlhdlrexprdata: 95 * 96 * vars = {x, y, z} 97 * offsets = {2, 0, 0, SQRT(5), -1} 98 * transcoefs = {3, -4, 1, SQRT(7), 5, -1} 99 * transcoefsidx = {0, 1, 1, 2, 0, 1} 100 * termbegins = {0, 2, 3, 4, 4, 6} 101 * nvars = 3 102 * nterms = 5 103 * 104 * @note: due to the current implementation, the constant term is the second to last term, except when the SOC was a rotated 105 * SOC, e.g., 1 + x^2 - y*z, i.e., when detected by detectSocQuadraticSimple. In that case, the constant is third to 106 * last term. 107 */ 108 struct SCIP_NlhdlrExprData 109 { 110 SCIP_EXPR** vars; /**< expressions which (aux)variables appear on both sides (x) */ 111 SCIP_Real* offsets; /**< offsets of both sides (beta_i) */ 112 SCIP_Real* transcoefs; /**< non-zeros of linear transformation vectors (v_i) */ 113 int* transcoefsidx; /**< mapping of transformation coefficients to variable indices in vars */ 114 int* termbegins; /**< starting indices of transcoefs for each term */ 115 int nvars; /**< total number of variables appearing */ 116 int nterms; /**< number of summands in the SQRT +1 for RHS (n+1) */ 117 118 /* variables for cone disaggregation */ 119 SCIP_VAR** disvars; /**< disaggregation variables for each term in lhs */ 120 SCIP_ROW* disrow; /**< disaggregation row */ 121 122 /* separation data */ 123 SCIP_Real* varvals; /**< current values for vars */ 124 SCIP_Real* disvarvals; /**< current values for disvars */ 125 }; 126 127 struct SCIP_NlhdlrData 128 { 129 SCIP_Real mincutefficacy; /**< minimum efficacy a cut need to be added */ 130 SCIP_Bool compeigenvalues; /**< whether Eigenvalue computations should be done to detect complex cases */ 131 }; 132 133 /* 134 * Local methods 135 */ 136 137 #ifdef SCIP_DEBUG 138 /** prints the nlhdlr expression data */ 139 static 140 void printNlhdlrExprData( 141 SCIP* scip, /**< SCIP data structure */ 142 SCIP_NLHDLREXPRDATA* nlhdlrexprdata /**< pointer to store nonlinear handler expression data */ 143 ) 144 { 145 int nterms; 146 int i; 147 int j; 148 149 nterms = nlhdlrexprdata->nterms; 150 151 SCIPinfoMessage(scip, NULL, "SQRT( "); 152 153 for( i = 0; i < nterms - 1; ++i ) 154 { 155 int startidx; 156 157 startidx = nlhdlrexprdata->termbegins[i]; 158 159 /* v_i is 0 */ 160 if( startidx == nlhdlrexprdata->termbegins[i + 1] ) 161 { 162 assert(nlhdlrexprdata->offsets[i] != 0.0); 163 164 SCIPinfoMessage(scip, NULL, "%g", SQR(nlhdlrexprdata->offsets[i])); 165 continue; 166 } 167 168 /* v_i is not 0 */ 169 SCIPinfoMessage(scip, NULL, "("); 170 171 for( j = startidx; j < nlhdlrexprdata->termbegins[i + 1]; ++j ) 172 { 173 if( nlhdlrexprdata->transcoefs[j] != 1.0 ) 174 SCIPinfoMessage(scip, NULL, " %+g*", nlhdlrexprdata->transcoefs[j]); 175 else 176 SCIPinfoMessage(scip, NULL, " +"); 177 if( SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]]) != NULL ) 178 { 179 SCIPinfoMessage(scip, NULL, "%s", SCIPvarGetName(SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]]))); 180 SCIPinfoMessage(scip, NULL, "(%p)", (void*)nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]]); 181 } 182 else 183 SCIPinfoMessage(scip, NULL, "%p", (void*)nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]]); 184 } 185 if( nlhdlrexprdata->offsets[i] != 0.0 ) 186 SCIPinfoMessage(scip, NULL, " %+g", nlhdlrexprdata->offsets[i]); 187 188 SCIPinfoMessage(scip, NULL, ")^2"); 189 190 if( i < nterms - 2 ) 191 SCIPinfoMessage(scip, NULL, " + "); 192 } 193 194 SCIPinfoMessage(scip, NULL, " ) <="); 195 196 for( j = nlhdlrexprdata->termbegins[nterms-1]; j < nlhdlrexprdata->termbegins[nterms]; ++j ) 197 { 198 if( nlhdlrexprdata->transcoefs[j] != 1.0 ) 199 SCIPinfoMessage(scip, NULL, " %+g*", nlhdlrexprdata->transcoefs[j]); 200 else 201 SCIPinfoMessage(scip, NULL, " +"); 202 if( SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]]) != NULL ) 203 SCIPinfoMessage(scip, NULL, "%s", SCIPvarGetName(SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]]))); 204 else 205 SCIPinfoMessage(scip, NULL, "%p", (void*)nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]]); 206 } 207 if( nlhdlrexprdata->offsets[nterms-1] != 0.0 ) 208 SCIPinfoMessage(scip, NULL, " %+g", nlhdlrexprdata->offsets[nterms-1]); 209 210 SCIPinfoMessage(scip, NULL, "\n"); 211 } 212 #endif 213 214 /** helper method to create variables for the cone disaggregation */ 215 static 216 SCIP_RETCODE createDisaggrVars( 217 SCIP* scip, /**< SCIP data structure */ 218 SCIP_EXPR* expr, /**< expression */ 219 SCIP_NLHDLREXPRDATA* nlhdlrexprdata /**< nonlinear handler expression data */ 220 ) 221 { 222 char name[SCIP_MAXSTRLEN]; 223 int ndisvars; 224 int i; 225 226 assert(nlhdlrexprdata != NULL); 227 228 ndisvars = nlhdlrexprdata->nterms - 1; 229 230 /* allocate memory */ 231 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &nlhdlrexprdata->disvars, ndisvars) ); 232 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &nlhdlrexprdata->disvarvals, ndisvars) ); 233 234 /* create disaggregation variables representing the epigraph of (v_i^T x + beta_i)^2 / (v_{n+1}^T x + beta_{n+1}) */ 235 for( i = 0; i < ndisvars; ++i ) 236 { 237 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedis_%p_%d", (void*) expr, i); 238 SCIP_CALL( SCIPcreateVarBasic(scip, &nlhdlrexprdata->disvars[i], name, 0.0, SCIPinfinity(scip), 0.0, 239 SCIP_VARTYPE_CONTINUOUS) ); 240 SCIPvarMarkRelaxationOnly(nlhdlrexprdata->disvars[i]); 241 242 SCIP_CALL( SCIPaddVar(scip, nlhdlrexprdata->disvars[i]) ); 243 SCIP_CALL( SCIPaddVarLocksType(scip, nlhdlrexprdata->disvars[i], SCIP_LOCKTYPE_MODEL, 1, 1) ); 244 } 245 246 return SCIP_OKAY; 247 } 248 249 /** helper method to free variables for the cone disaggregation */ 250 static 251 SCIP_RETCODE freeDisaggrVars( 252 SCIP* scip, /**< SCIP data structure */ 253 SCIP_NLHDLREXPRDATA* nlhdlrexprdata /**< nonlinear handler expression data */ 254 ) 255 { 256 int ndisvars; 257 int i; 258 259 assert(nlhdlrexprdata != NULL); 260 261 if( nlhdlrexprdata->disvars == NULL ) 262 return SCIP_OKAY; 263 264 ndisvars = nlhdlrexprdata->nterms - 1; 265 266 /* release variables */ 267 for( i = 0; i < ndisvars; ++i ) 268 { 269 SCIP_CALL( SCIPaddVarLocksType(scip, nlhdlrexprdata->disvars[i], SCIP_LOCKTYPE_MODEL, -1, -1) ); 270 SCIP_CALL( SCIPreleaseVar(scip, &nlhdlrexprdata->disvars[i]) ); 271 } 272 273 /* free memory */ 274 SCIPfreeBlockMemoryArray(scip, &nlhdlrexprdata->disvars, ndisvars); 275 SCIPfreeBlockMemoryArrayNull(scip, &nlhdlrexprdata->disvarvals, ndisvars); 276 277 return SCIP_OKAY; 278 } 279 280 /** helper method to create the disaggregation row \f$\text{disvars}_i \leq v_{n+1}^T x + \beta_{n+1}\f$ */ 281 static 282 SCIP_RETCODE createDisaggrRow( 283 SCIP* scip, /**< SCIP data structure */ 284 SCIP_CONSHDLR* conshdlr, /**< nonlinear constraint handler */ 285 SCIP_EXPR* expr, /**< expression */ 286 SCIP_NLHDLREXPRDATA* nlhdlrexprdata /**< nonlinear handler expression data */ 287 ) 288 { 289 SCIP_Real beta; 290 char name[SCIP_MAXSTRLEN]; 291 int ndisvars; 292 int nterms; 293 int i; 294 295 assert(scip != NULL); 296 assert(expr != NULL); 297 assert(nlhdlrexprdata != NULL); 298 assert(nlhdlrexprdata->disrow == NULL); 299 300 nterms = nlhdlrexprdata->nterms; 301 beta = nlhdlrexprdata->offsets[nterms - 1]; 302 303 ndisvars = nterms - 1; 304 305 /* create row 0 <= beta_{n+1} */ 306 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedis_%p_row", (void*) expr); 307 SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &nlhdlrexprdata->disrow, conshdlr, name, 308 -SCIPinfinity(scip), beta, FALSE, FALSE, TRUE) ); 309 310 /* add disvars to row */ 311 for( i = 0; i < ndisvars; ++i ) 312 { 313 SCIP_CALL( SCIPaddVarToRow(scip, nlhdlrexprdata->disrow, nlhdlrexprdata->disvars[i], 1.0) ); 314 } 315 316 /* add rhs vars to row */ 317 for( i = nlhdlrexprdata->termbegins[nterms - 1]; i < nlhdlrexprdata->termbegins[nterms]; ++i ) 318 { 319 SCIP_VAR* var; 320 SCIP_Real coef; 321 322 var = SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[i]]); 323 assert(var != NULL); 324 325 coef = -nlhdlrexprdata->transcoefs[i]; 326 327 SCIP_CALL( SCIPaddVarToRow(scip, nlhdlrexprdata->disrow, var, coef) ); 328 } 329 330 return SCIP_OKAY; 331 } 332 333 /** helper method to create nonlinear handler expression data */ 334 static 335 SCIP_RETCODE createNlhdlrExprData( 336 SCIP* scip, /**< SCIP data structure */ 337 SCIP_EXPR** vars, /**< expressions which variables appear on both sides (\f$x\f$) */ 338 SCIP_Real* offsets, /**< offsets of bot sides (\f$beta_i\f$) */ 339 SCIP_Real* transcoefs, /**< non-zeroes of linear transformation vectors (\f$v_i\f$) */ 340 int* transcoefsidx, /**< mapping of transformation coefficients to variable indices in vars */ 341 int* termbegins, /**< starting indices of transcoefs for each term */ 342 int nvars, /**< total number of variables appearing */ 343 int nterms, /**< number of summands in the SQRT, +1 for RHS */ 344 SCIP_NLHDLREXPRDATA** nlhdlrexprdata /**< pointer to store nonlinear handler expression data */ 345 ) 346 { 347 int ntranscoefs; 348 349 assert(vars != NULL); 350 assert(offsets != NULL); 351 assert(transcoefs != NULL); 352 assert(transcoefsidx != NULL); 353 assert(termbegins != NULL); 354 assert(nlhdlrexprdata != NULL); 355 356 ntranscoefs = termbegins[nterms]; 357 358 SCIP_CALL( SCIPallocBlockMemory(scip, nlhdlrexprdata) ); 359 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*nlhdlrexprdata)->vars, vars, nvars) ); 360 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*nlhdlrexprdata)->offsets, offsets, nterms) ); 361 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*nlhdlrexprdata)->transcoefs, transcoefs, ntranscoefs) ); 362 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*nlhdlrexprdata)->transcoefsidx, transcoefsidx, ntranscoefs) ); 363 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*nlhdlrexprdata)->termbegins, termbegins, nterms + 1) ); 364 (*nlhdlrexprdata)->nvars = nvars; 365 (*nlhdlrexprdata)->nterms = nterms; 366 367 (*nlhdlrexprdata)->disrow = NULL; 368 (*nlhdlrexprdata)->disvars = NULL; 369 370 (*nlhdlrexprdata)->varvals = NULL; 371 (*nlhdlrexprdata)->disvarvals = NULL; 372 373 #ifdef SCIP_DEBUG 374 SCIPdebugMsg(scip, "created nlhdlr data for the following soc expression:\n"); 375 printNlhdlrExprData(scip, *nlhdlrexprdata); 376 /* SCIPdebugMsg(scip, "x is %p\n", (void *)vars[0]); */ 377 #endif 378 379 return SCIP_OKAY; 380 } 381 382 /** helper method to free nonlinear handler expression data */ 383 static 384 SCIP_RETCODE freeNlhdlrExprData( 385 SCIP* scip, /**< SCIP data structure */ 386 SCIP_NLHDLREXPRDATA** nlhdlrexprdata /**< pointer to free nonlinear handler expression data */ 387 ) 388 { 389 int ntranscoefs; 390 391 assert(nlhdlrexprdata != NULL); 392 assert(*nlhdlrexprdata != NULL); 393 394 /* free variables and row for cone disaggregation */ 395 SCIP_CALL( freeDisaggrVars(scip, *nlhdlrexprdata) ); 396 397 ntranscoefs = (*nlhdlrexprdata)->termbegins[(*nlhdlrexprdata)->nterms]; 398 399 SCIPfreeBlockMemoryArrayNull(scip, &(*nlhdlrexprdata)->varvals, (*nlhdlrexprdata)->nvars); 400 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->termbegins, (*nlhdlrexprdata)->nterms + 1); 401 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->transcoefsidx, ntranscoefs); 402 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->transcoefs, ntranscoefs); 403 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->offsets, (*nlhdlrexprdata)->nterms); 404 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->vars, (*nlhdlrexprdata)->nvars); 405 SCIPfreeBlockMemory(scip, nlhdlrexprdata); 406 407 return SCIP_OKAY; 408 } 409 410 /** set varvalrs in nlhdlrexprdata to values from given SCIP solution */ 411 static 412 void updateVarVals( 413 SCIP* scip, /**< SCIP data structure */ 414 SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nonlinear handler expression data */ 415 SCIP_SOL* sol, /**< SCIP solution */ 416 SCIP_Bool roundtinyfrac /**< whether values close to integers should be rounded */ 417 ) 418 { 419 int i; 420 421 assert(nlhdlrexprdata != NULL); 422 assert(nlhdlrexprdata->varvals != NULL); 423 424 /* update varvals */ 425 for( i = 0; i < nlhdlrexprdata->nvars; ++i ) 426 { 427 nlhdlrexprdata->varvals[i] = SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[i])); 428 if( roundtinyfrac && SCIPisIntegral(scip, nlhdlrexprdata->varvals[i]) ) 429 nlhdlrexprdata->varvals[i] = SCIPround(scip, nlhdlrexprdata->varvals[i]); 430 } 431 432 /* update disvarvals (in unittests, this may be NULL even though nterms > 1 */ 433 if( nlhdlrexprdata->disvarvals != NULL ) 434 for( i = 0; i < nlhdlrexprdata->nterms - 1; ++i ) 435 { 436 nlhdlrexprdata->disvarvals[i] = SCIPgetSolVal(scip, sol, nlhdlrexprdata->disvars[i]); 437 if( roundtinyfrac && SCIPisIntegral(scip, nlhdlrexprdata->disvarvals[i]) ) 438 nlhdlrexprdata->disvarvals[i] = SCIPround(scip, nlhdlrexprdata->disvarvals[i]); 439 } 440 } 441 442 /** evaluate a single term of the form \f$v_i^T x + \beta_i\f$ */ 443 static 444 SCIP_Real evalSingleTerm( 445 SCIP* scip, /**< SCIP data structure */ 446 SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nonlinear handler expression data */ 447 int k /**< term to be evaluated */ 448 ) 449 { 450 SCIP_Real result; 451 int i; 452 453 assert(scip != NULL); 454 assert(nlhdlrexprdata != NULL); 455 assert(0 <= k && k < nlhdlrexprdata->nterms); 456 457 result = nlhdlrexprdata->offsets[k]; 458 459 for( i = nlhdlrexprdata->termbegins[k]; i < nlhdlrexprdata->termbegins[k + 1]; ++i ) 460 result += nlhdlrexprdata->transcoefs[i] * nlhdlrexprdata->varvals[nlhdlrexprdata->transcoefsidx[i]]; 461 462 return result; 463 } 464 465 /** computes gradient cut for a 2D or 3D SOC 466 * 467 * A 3D SOC looks like 468 * \f[ 469 * \sqrt{ (v_1^T x + \beta_1)^2 + (v_2^T x + \beta_2)^2 } \leq v_3^T x + \beta_3 470 * \f] 471 * 472 * Let \f$f(x)\f$ be the left-hand-side. The partial derivatives of \f$f\f$ are given by 473 * \f[ 474 * \frac{\delta f}{\delta x_j} = \frac{(v_1)_j(v_1^T x + \beta_1) + (v_2)_j (v_2^T x + \beta_2)}{f(x)} 475 * \f] 476 * 477 * and the gradient cut is then \f$f(x^*) + \nabla f(x^*)(x - x^*) \leq v_3^T x + \beta_3\f$. 478 * 479 * If \f$\beta_1 = \beta_2 = 0\f$, then the constant on the left-hand-side of the cut becomes zero: 480 * \f[ 481 * f(x^*) - (\frac{(v_1)_j v_1^T x^* + (v_2)_j v_2^T x^*}{f(x^*)})_j^T x^* 482 * = f(x^*) - \frac{1}{f(x^*)} \sum_j ((v_1)_j x_j^* v_1^T x^* + (v_2)_j x_j^* v_2^T x^*) 483 * = f(x^*) - \frac{1}{f(x^*)} ((v_1^T x^*)^2 + (v_2^T x^*)^2) 484 * = f(x^*) - \frac{1}{f(x^*)} f(x^*)^2 = 0 485 * \f] 486 * 487 * A 2D SOC is 488 * \f[ 489 * |v_1^T x + \beta_1| \leq v_2^T x + \beta_2 490 * \f] 491 * but we build the cut using the same procedure as for 3D. 492 */ 493 static 494 SCIP_RETCODE generateCutSolSOC( 495 SCIP* scip, /**< SCIP data structure */ 496 SCIP_ROWPREP** rowprep, /**< buffer to store rowprep with cut data */ 497 SCIP_EXPR* expr, /**< expression */ 498 SCIP_CONS* cons, /**< the constraint that expr is part of */ 499 SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nonlinear handler expression data */ 500 SCIP_Real mincutviolation, /**< minimal required cut violation */ 501 SCIP_Real rhsval /**< value of last term at sol */ 502 ) 503 { 504 SCIP_Real* transcoefs; 505 SCIP_Real cutcoef; 506 SCIP_Real fvalue; 507 SCIP_Real valterms[2] = {0.0, 0.0}; /* for lint */ 508 SCIP_Real cutrhs; 509 SCIP_EXPR** vars; 510 SCIP_VAR* cutvar; 511 SCIP_Bool offsetzero; 512 int* transcoefsidx; 513 int* termbegins; 514 int nterms; 515 int i; 516 int j; 517 518 assert(rowprep != NULL); 519 assert(expr != NULL); 520 assert(cons != NULL); 521 assert(nlhdlrexprdata != NULL); 522 523 vars = nlhdlrexprdata->vars; 524 transcoefs = nlhdlrexprdata->transcoefs; 525 transcoefsidx = nlhdlrexprdata->transcoefsidx; 526 termbegins = nlhdlrexprdata->termbegins; 527 nterms = nlhdlrexprdata->nterms; 528 529 *rowprep = NULL; 530 531 /* evaluate lhs terms and compute f(x*), check whether both beta_1 and beta_2 are zero */ 532 fvalue = 0.0; 533 offsetzero = TRUE; 534 for( i = 0; i < nterms - 1; ++i ) 535 { 536 valterms[i] = evalSingleTerm(scip, nlhdlrexprdata, i); 537 fvalue += SQR( valterms[i] ); 538 if( nlhdlrexprdata->offsets[i] != 0.0 ) 539 offsetzero = FALSE; 540 } 541 fvalue = SQRT( fvalue ); 542 543 /* don't generate cut if we are not violated @todo: remove this once core detects better when a nlhdlr's cons is 544 * violated 545 */ 546 if( fvalue - rhsval <= mincutviolation ) 547 { 548 SCIPdebugMsg(scip, "do not generate cut: rhsval %g, fvalue %g violation is %g\n", rhsval, fvalue, fvalue - rhsval); 549 return SCIP_OKAY; 550 } 551 552 /* if f(x*) = 0 then we are at top of cone, where we cannot generate cut */ 553 if( SCIPisZero(scip, fvalue) ) 554 { 555 SCIPdebugMsg(scip, "do not generate cut for lhs=%g, cannot linearize at top of cone\n", fvalue); 556 return SCIP_OKAY; 557 } 558 559 /* create cut */ 560 SCIP_CALL( SCIPcreateRowprep(scip, rowprep, SCIP_SIDETYPE_RIGHT, FALSE) ); 561 SCIP_CALL( SCIPensureRowprepSize(scip, *rowprep, termbegins[nterms]) ); 562 563 /* cut is f(x*) + \nabla f(x*)^T (x - x*) \leq v_n^T x + \beta_n, i.e., 564 * \nabla f(x*)^T x - v_n^T x \leq \beta_n + \nabla f(x*)^T x* - f(x*) 565 * thus cutrhs is \beta_n - f(x*) + \nabla f(x*)^T x* 566 * if offsetzero, then we make sure that cutrhs is exactly \beta_n 567 */ 568 cutrhs = nlhdlrexprdata->offsets[nterms - 1]; 569 if( !offsetzero ) 570 cutrhs -= fvalue; 571 572 /* add cut coefficients from lhs terms and compute cut's rhs */ 573 for( j = 0; j < nterms - 1; ++j ) 574 { 575 for( i = termbegins[j]; i < termbegins[j + 1]; ++i ) 576 { 577 cutvar = SCIPgetExprAuxVarNonlinear(vars[transcoefsidx[i]]); 578 579 /* cutcoef is (the first part of) the partial derivative w.r.t cutvar */ 580 cutcoef = transcoefs[i] * valterms[j] / fvalue; 581 582 SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, cutvar, cutcoef) ); 583 584 if( !offsetzero ) 585 cutrhs += cutcoef * nlhdlrexprdata->varvals[transcoefsidx[i]]; 586 } 587 } 588 589 /* add terms for v_n */ 590 for( i = termbegins[nterms - 1]; i < termbegins[nterms]; ++i ) 591 { 592 cutvar = SCIPgetExprAuxVarNonlinear(vars[transcoefsidx[i]]); 593 SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, cutvar, -transcoefs[i]) ); 594 } 595 596 /* add side */ 597 SCIProwprepAddSide(*rowprep, cutrhs); 598 599 /* set name */ 600 (void) SCIPsnprintf(SCIProwprepGetName(*rowprep), SCIP_MAXSTRLEN, "soc%d_%p_%d", nterms, (void*) expr, SCIPgetNLPs(scip)); 601 602 return SCIP_OKAY; 603 } 604 605 /** helper method to compute and add a gradient cut for the k-th cone disaggregation 606 * 607 * After the SOC constraint \f$\sqrt{\sum_{i = 0}^{n-1} (v_i^T x + \beta_i)^2} \leq v_n^T x + \beta_n\f$ 608 * has been disaggregated into the row \f$\sum_{i = 0}^{n-1} y_i \leq v_n^T x + \beta_n\f$ and the smaller SOC constraints 609 * \f[ 610 * (v_i^T x + \beta_i)^2 \leq (v_n^T x + \beta_n) y_i \text{ for } i \in \{0, \ldots, n -1\}, 611 * \f] 612 * we want to separate one of the small rotated cones. 613 * We first transform it into standard form: 614 * \f[ 615 * \sqrt{4(v_i^T x + \beta_i)^2 + (v_n^T x + \beta_n - y_i)^2} - v_n^T x - \beta_n - y_i \leq 0. 616 * \f] 617 * Let \f$f(x,y)\f$ be the left-hand-side. We now compute the gradient by 618 * \f{align*}{ 619 * \frac{\delta f}{\delta x_j} &= \frac{(v_i)_j(4v_i^T x + 4\beta_i) + (v_n)_j(v_n^T x + \beta_n - y_i)}{\sqrt{4(v_i^T x + \beta_i)^2 + (v_n^T x + \beta_n - y_i)^2}} - (v_n)_j \\ 620 * \frac{\delta f}{\delta y_i} &= \frac{y_i - v_n^T x -\beta_n}{\sqrt{4(v_i^T x + \beta_i)^2 + (v_n^T x + \beta_n - y_i)^2}} - 1 621 * \f} 622 * and the gradient cut is then \f$f(x^*, y^*) + \nabla f(x^*,y^*)((x,y) - (x^*, y^*)) \leq 0\f$. 623 * 624 * As in \ref generateCutSolSOC(), the cut constant is zero if \f$\beta_i = \beta_n = 0\f$. 625 */ 626 static 627 SCIP_RETCODE generateCutSolDisagg( 628 SCIP* scip, /**< SCIP data structure */ 629 SCIP_ROWPREP** rowprep, /**< buffer to store rowprep with cut data */ 630 SCIP_EXPR* expr, /**< expression */ 631 SCIP_CONS* cons, /**< the constraint that expr is part of */ 632 SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nonlinear handler expression data */ 633 int disaggidx, /**< index of disaggregation to separate */ 634 SCIP_Real mincutviolation, /**< minimal required cut violation */ 635 SCIP_Real rhsval /**< value of the rhs term */ 636 ) 637 { 638 SCIP_EXPR** vars; 639 SCIP_VAR** disvars; 640 SCIP_Real* transcoefs; 641 int* transcoefsidx; 642 int* termbegins; 643 SCIP_VAR* cutvar; 644 SCIP_Real cutcoef; 645 SCIP_Real fvalue; 646 SCIP_Real disvarval; 647 SCIP_Real lhsval; 648 SCIP_Real constant; 649 SCIP_Real denominator; 650 SCIP_Bool offsetzero; 651 int ncutvars; 652 int nterms; 653 int i; 654 655 assert(rowprep != NULL); 656 assert(expr != NULL); 657 assert(cons != NULL); 658 assert(nlhdlrexprdata != NULL); 659 assert(disaggidx < nlhdlrexprdata->nterms-1); 660 661 vars = nlhdlrexprdata->vars; 662 disvars = nlhdlrexprdata->disvars; 663 transcoefs = nlhdlrexprdata->transcoefs; 664 transcoefsidx = nlhdlrexprdata->transcoefsidx; 665 termbegins = nlhdlrexprdata->termbegins; 666 nterms = nlhdlrexprdata->nterms; 667 668 /* nterms is equal to n in the description and disaggidx is in {0, ..., n - 1} */ 669 670 *rowprep = NULL; 671 672 disvarval = nlhdlrexprdata->disvarvals[disaggidx]; 673 674 lhsval = evalSingleTerm(scip, nlhdlrexprdata, disaggidx); 675 676 denominator = SQRT(4.0 * SQR(lhsval) + SQR(rhsval - disvarval)); 677 678 /* compute value of function to be separated (f(x*,y*)) */ 679 fvalue = denominator - rhsval - disvarval; 680 681 /* if the disagg soc is not violated don't compute cut */ 682 if( fvalue <= mincutviolation ) 683 { 684 SCIPdebugMsg(scip, "skip cut on disaggregation index %d as violation=%g below minviolation %g\n", disaggidx, 685 fvalue, mincutviolation); 686 return SCIP_OKAY; 687 } 688 689 /* if the denominator is 0 -> the constraint can't be violated, and the gradient is infinite */ 690 if( SCIPisZero(scip, denominator) ) 691 { 692 SCIPdebugMsg(scip, "skip cut on disaggregation index %d as we are on top of cone (denom=%g)\n", disaggidx, denominator); 693 return SCIP_OKAY; 694 } 695 696 /* compute upper bound on the number of variables in cut: vars in rhs + vars in term + disagg var */ 697 ncutvars = (termbegins[nterms] - termbegins[nterms-1]) + (termbegins[disaggidx + 1] - termbegins[disaggidx]) + 1; 698 699 /* create cut */ 700 SCIP_CALL( SCIPcreateRowprep(scip, rowprep, SCIP_SIDETYPE_RIGHT, FALSE) ); 701 SCIP_CALL( SCIPensureRowprepSize(scip, *rowprep, ncutvars) ); 702 703 /* check whether offsets (beta) are zero, so we can know cut constant will be zero */ 704 offsetzero = nlhdlrexprdata->offsets[disaggidx] == 0.0 && nlhdlrexprdata->offsets[nterms-1] == 0.0; 705 706 /* constant will be grad_f(x*,y*)^T (x*, y*) */ 707 constant = 0.0; 708 709 /* a variable could appear on the lhs and rhs, but we add the coefficients separately */ 710 711 /* add terms for v_disaggidx */ 712 for( i = termbegins[disaggidx]; i < termbegins[disaggidx + 1]; ++i ) 713 { 714 cutvar = SCIPgetExprAuxVarNonlinear(vars[transcoefsidx[i]]); 715 assert(cutvar != NULL); 716 717 /* cutcoef is (the first part of) the partial derivative w.r.t cutvar */ 718 cutcoef = 4.0 * lhsval * transcoefs[i] / denominator; 719 720 SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, cutvar, cutcoef) ); 721 722 if( !offsetzero ) 723 constant += cutcoef * nlhdlrexprdata->varvals[transcoefsidx[i]]; 724 } 725 726 /* add terms for v_n */ 727 for( i = termbegins[nterms - 1]; i < termbegins[nterms]; ++i ) 728 { 729 cutvar = SCIPgetExprAuxVarNonlinear(vars[transcoefsidx[i]]); 730 assert(cutvar != NULL); 731 732 /* cutcoef is the (second part of) the partial derivative w.r.t cutvar */ 733 cutcoef = (rhsval - disvarval) * transcoefs[i] / denominator - transcoefs[i]; 734 735 SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, cutvar, cutcoef) ); 736 737 if( !offsetzero ) 738 constant += cutcoef * nlhdlrexprdata->varvals[transcoefsidx[i]]; 739 } 740 741 /* add term for disvar: cutcoef is the the partial derivative w.r.t. the disaggregation variable */ 742 cutcoef = (disvarval - rhsval) / denominator - 1.0; 743 cutvar = disvars[disaggidx]; 744 745 SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, cutvar, cutcoef) ); 746 747 if( !offsetzero ) 748 { 749 constant += cutcoef * nlhdlrexprdata->disvarvals[disaggidx]; 750 751 /* add side */ 752 SCIProwprepAddSide(*rowprep, constant - fvalue); 753 } 754 755 /* set name */ 756 (void) SCIPsnprintf(SCIProwprepGetName(*rowprep), SCIP_MAXSTRLEN, "soc_%p_%d_%d", (void*) expr, disaggidx, SCIPgetNLPs(scip)); 757 758 return SCIP_OKAY; 759 } 760 761 /** given a rowprep, does a number of cleanup and checks and, if successful, generate a cut to be added to the sepastorage */ 762 static 763 SCIP_RETCODE addCut( 764 SCIP* scip, /**< SCIP data structure */ 765 SCIP_NLHDLRDATA* nlhdlrdata, /**< nonlinear handler data */ 766 SCIP_ROWPREP* rowprep, /**< rowprep from which to generate row and add as cut */ 767 SCIP_SOL* sol, /**< solution to be separated */ 768 SCIP_CONS* cons, /**< constraint for which cut is generated, or NULL */ 769 SCIP_Bool allowweakcuts, /**< whether weak cuts are allowed */ 770 SCIP_RESULT* result /**< result pointer to update (set to SCIP_CUTOFF or SCIP_SEPARATED if cut is added) */ 771 ) 772 { 773 SCIP_ROW* cut; 774 SCIP_Real cutefficacy; 775 SCIP_Bool success; 776 777 assert(scip != NULL); 778 assert(nlhdlrdata != NULL); 779 assert(rowprep != NULL); 780 assert(result != NULL); 781 782 SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, sol, SCIPgetHugeValue(scip), &success) ); 783 784 if( !success ) 785 { 786 SCIPdebugMsg(scip, "rowprep cleanup failed, skip cut\n"); 787 return SCIP_OKAY; 788 } 789 790 if( SCIPgetRowprepViolation(scip, rowprep, sol, NULL) <= SCIPgetLPFeastol(scip) ) 791 { 792 SCIPdebugMsg(scip, "rowprep violation %g below LP feastol %g, skip cut\n", 793 SCIPgetRowprepViolation(scip, rowprep, sol, NULL), SCIPgetLPFeastol(scip)); 794 return SCIP_OKAY; 795 } 796 797 SCIP_CALL( SCIPgetRowprepRowCons(scip, &cut, rowprep, cons) ); 798 799 cutefficacy = SCIPgetCutEfficacy(scip, sol, cut); 800 801 SCIPdebugMsg(scip, "generated row for SOC, efficacy=%g, minefficacy=%g, allowweakcuts=%u\n", 802 cutefficacy, nlhdlrdata->mincutefficacy, allowweakcuts); 803 804 /* check whether cut is applicable */ 805 if( SCIPisCutApplicable(scip, cut) && (allowweakcuts || cutefficacy >= nlhdlrdata->mincutefficacy) ) 806 { 807 SCIP_Bool infeasible; 808 809 SCIP_CALL( SCIPaddRow(scip, cut, FALSE, &infeasible) ); 810 811 #ifdef SCIP_CONSNONLINEAR_ROWNOTREMOVABLE 812 /* mark row as not removable from LP for current node, if in enforcement (==addbranchscores) (this can prevent some cycling) */ 813 if( addbranchscores ) 814 SCIPmarkRowNotRemovableLocal(scip, row); 815 #endif 816 817 if( infeasible ) 818 *result = SCIP_CUTOFF; 819 else 820 *result = SCIP_SEPARATED; 821 } 822 823 /* release row */ 824 SCIP_CALL( SCIPreleaseRow(scip, &cut) ); 825 826 return SCIP_OKAY; 827 } 828 829 /** given a rowprep, does a number of cleanup and checks and, if successful, generate a cut to be added to the cutpool */ 830 static 831 SCIP_RETCODE addCutPool( 832 SCIP* scip, /**< SCIP data structure */ 833 SCIP_NLHDLRDATA* nlhdlrdata, /**< nonlinear handler data */ 834 SCIP_ROWPREP* rowprep, /**< rowprep from which to generate row and add as cut */ 835 SCIP_SOL* sol, /**< solution to be separated */ 836 SCIP_CONS* cons /**< constraint for which cut is generated, or NULL */ 837 ) 838 { 839 SCIP_ROW* cut; 840 SCIP_Bool success; 841 842 assert(scip != NULL); 843 assert(nlhdlrdata != NULL); 844 assert(rowprep != NULL); 845 846 assert(!SCIProwprepIsLocal(rowprep)); 847 848 SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, sol, SCIPgetHugeValue(scip), &success) ); 849 /* if failed or cut is only locally valid now, then skip */ 850 if( !success || SCIProwprepIsLocal(rowprep) ) 851 return SCIP_OKAY; 852 853 /* if row after cleanup is just a boundchange, then skip */ 854 if( SCIProwprepGetNVars(rowprep) <= 1 ) 855 return SCIP_OKAY; 856 857 /* generate row and add to cutpool */ 858 SCIP_CALL( SCIPgetRowprepRowCons(scip, &cut, rowprep, cons) ); 859 860 SCIP_CALL( SCIPaddPoolCut(scip, cut) ); 861 862 SCIP_CALL( SCIPreleaseRow(scip, &cut) ); 863 864 return SCIP_OKAY; 865 } 866 867 /** checks if an expression is quadratic and collects all occurring expressions 868 * 869 * @pre `expr2idx` and `occurringexprs` need to be initialized with capacity 2 * nchildren 870 * 871 * @note We assume that a linear term always appears before its corresponding 872 * quadratic term in quadexpr; this should be ensured by canonicalize 873 */ 874 static 875 SCIP_RETCODE checkAndCollectQuadratic( 876 SCIP* scip, /**< SCIP data structure */ 877 SCIP_EXPR* quadexpr, /**< candidate for a quadratic expression */ 878 SCIP_HASHMAP* expr2idx, /**< hashmap to store expressions */ 879 SCIP_EXPR** occurringexprs, /**< array to store expressions */ 880 int* nexprs, /**< buffer to store number of expressions */ 881 SCIP_Bool* success /**< buffer to store whether the check was successful */ 882 ) 883 { 884 SCIP_EXPR** children; 885 int nchildren; 886 int i; 887 888 assert(scip != NULL); 889 assert(quadexpr != NULL); 890 assert(expr2idx != NULL); 891 assert(occurringexprs != NULL); 892 assert(nexprs != NULL); 893 assert(success != NULL); 894 895 *nexprs = 0; 896 *success = FALSE; 897 children = SCIPexprGetChildren(quadexpr); 898 nchildren = SCIPexprGetNChildren(quadexpr); 899 900 /* iterate in reverse order to ensure that quadratic terms are found before linear terms */ 901 for( i = nchildren - 1; i >= 0; --i ) 902 { 903 SCIP_EXPR* child; 904 905 child = children[i]; 906 if( SCIPisExprPower(scip, child) ) 907 { 908 SCIP_EXPR* childarg; 909 910 if( SCIPgetExponentExprPow(child) != 2.0 ) 911 return SCIP_OKAY; 912 913 childarg = SCIPexprGetChildren(child)[0]; 914 915 if( !SCIPhashmapExists(expr2idx, (void*) childarg) ) 916 { 917 SCIP_CALL( SCIPhashmapInsertInt(expr2idx, (void*) childarg, *nexprs) ); 918 919 /* store the expression so we know it later */ 920 assert(*nexprs < 2 * nchildren); 921 occurringexprs[*nexprs] = childarg; 922 923 ++(*nexprs); 924 } 925 } 926 else if( SCIPisExprVar(scip, child) && SCIPvarIsBinary(SCIPgetVarExprVar(child)) ) 927 { 928 if( !SCIPhashmapExists(expr2idx, (void*) child) ) 929 { 930 SCIP_CALL( SCIPhashmapInsertInt(expr2idx, (void*) child, *nexprs) ); 931 932 /* store the expression so we know it later */ 933 assert(*nexprs < 2 * nchildren); 934 occurringexprs[*nexprs] = child; 935 936 ++(*nexprs); 937 } 938 } 939 else if( SCIPisExprProduct(scip, child) ) 940 { 941 SCIP_EXPR* childarg1; 942 SCIP_EXPR* childarg2; 943 944 if( SCIPexprGetNChildren(child) != 2 ) 945 return SCIP_OKAY; 946 947 childarg1 = SCIPexprGetChildren(child)[0]; 948 childarg2 = SCIPexprGetChildren(child)[1]; 949 950 if( !SCIPhashmapExists(expr2idx, (void*) childarg1) ) 951 { 952 SCIP_CALL( SCIPhashmapInsertInt(expr2idx, (void*) childarg1, *nexprs) ); 953 954 /* store the expression so we know it later */ 955 assert(*nexprs < 2 * nchildren); 956 occurringexprs[*nexprs] = childarg1; 957 958 ++(*nexprs); 959 } 960 961 if( !SCIPhashmapExists(expr2idx, (void*) childarg2) ) 962 { 963 SCIP_CALL( SCIPhashmapInsertInt(expr2idx, (void*) childarg2, *nexprs) ); 964 965 /* store the expression so we know it later */ 966 assert(*nexprs < 2 * nchildren); 967 occurringexprs[*nexprs] = childarg2; 968 969 ++(*nexprs); 970 } 971 } 972 else 973 { 974 /* if there is a linear term without corresponding quadratic term, it is not a SOC */ 975 if( !SCIPhashmapExists(expr2idx, (void*) child) ) 976 return SCIP_OKAY; 977 } 978 } 979 980 *success = TRUE; 981 982 return SCIP_OKAY; 983 } 984 985 /* builds the constraint defining matrix and vector of a quadratic expression 986 * 987 * @pre `quadmatrix` and `linvector` need to be initialized with size `nexprs`^2 and `nexprs`, resp. 988 */ 989 static 990 void buildQuadExprMatrix( 991 SCIP* scip, /**< SCIP data structure */ 992 SCIP_EXPR* quadexpr, /**< the quadratic expression */ 993 SCIP_HASHMAP* expr2idx, /**< hashmap mapping the occurring expressions to their index */ 994 int nexprs, /**< number of occurring expressions */ 995 SCIP_Real* quadmatrix, /**< pointer to store (the lower-left triangle of) the quadratic matrix */ 996 SCIP_Real* linvector /**< pointer to store the linear vector */ 997 ) 998 { 999 SCIP_EXPR** children; 1000 SCIP_Real* childcoefs; 1001 int nchildren; 1002 int i; 1003 1004 assert(scip != NULL); 1005 assert(quadexpr != NULL); 1006 assert(expr2idx != NULL); 1007 assert(quadmatrix != NULL); 1008 assert(linvector != NULL); 1009 1010 children = SCIPexprGetChildren(quadexpr); 1011 nchildren = SCIPexprGetNChildren(quadexpr); 1012 childcoefs = SCIPgetCoefsExprSum(quadexpr); 1013 1014 /* iterate over children to build the constraint defining matrix and vector */ 1015 for( i = 0; i < nchildren; ++i ) 1016 { 1017 int varpos; 1018 1019 if( SCIPisExprPower(scip, children[i]) ) 1020 { 1021 assert(SCIPgetExponentExprPow(children[i]) == 2.0); 1022 assert(SCIPhashmapExists(expr2idx, (void*) SCIPexprGetChildren(children[i])[0])); 1023 1024 varpos = SCIPhashmapGetImageInt(expr2idx, (void*) SCIPexprGetChildren(children[i])[0]); 1025 assert(0 <= varpos && varpos < nexprs); 1026 1027 quadmatrix[varpos * nexprs + varpos] = childcoefs[i]; 1028 } 1029 else if( SCIPisExprVar(scip, children[i]) && SCIPvarIsBinary(SCIPgetVarExprVar(children[i])) ) 1030 { 1031 assert(SCIPhashmapExists(expr2idx, (void*) children[i])); 1032 1033 varpos = SCIPhashmapGetImageInt(expr2idx, (void*) children[i]); 1034 assert(0 <= varpos && varpos < nexprs); 1035 1036 quadmatrix[varpos * nexprs + varpos] = childcoefs[i]; 1037 } 1038 else if( SCIPisExprProduct(scip, children[i]) ) 1039 { 1040 int varpos2; 1041 1042 assert(SCIPexprGetNChildren(children[i]) == 2); 1043 assert(SCIPhashmapExists(expr2idx, (void*) SCIPexprGetChildren(children[i])[0])); 1044 assert(SCIPhashmapExists(expr2idx, (void*) SCIPexprGetChildren(children[i])[1])); 1045 1046 varpos = SCIPhashmapGetImageInt(expr2idx, (void*) SCIPexprGetChildren(children[i])[0]); 1047 assert(0 <= varpos && varpos < nexprs); 1048 1049 varpos2 = SCIPhashmapGetImageInt(expr2idx, (void*) SCIPexprGetChildren(children[i])[1]); 1050 assert(0 <= varpos2 && varpos2 < nexprs); 1051 assert(varpos != varpos2); 1052 1053 /* Lapack uses only the lower left triangle of the symmetric matrix */ 1054 quadmatrix[MIN(varpos, varpos2) * nexprs + MAX(varpos, varpos2)] = childcoefs[i] / 2.0; 1055 } 1056 else 1057 { 1058 varpos = SCIPhashmapGetImageInt(expr2idx, (void*) children[i]); 1059 assert(0 <= varpos && varpos < nexprs); 1060 1061 linvector[varpos] = childcoefs[i]; 1062 } 1063 } 1064 } 1065 1066 /** tries to fill the nlhdlrexprdata for a potential quadratic SOC expression 1067 * 1068 * We say "try" because the expression might still turn out not to be a SOC at this point. 1069 */ 1070 static 1071 SCIP_RETCODE tryFillNlhdlrExprDataQuad( 1072 SCIP* scip, /**< SCIP data structure */ 1073 SCIP_EXPR** occurringexprs, /**< array of all occurring expressions (nvars many) */ 1074 SCIP_Real* eigvecmatrix, /**< array containing the Eigenvectors */ 1075 SCIP_Real* eigvals, /**< array containing the Eigenvalues */ 1076 SCIP_Real* bp, /**< product of linear vector b * P (eigvecmatrix^t) */ 1077 int nvars, /**< number of variables */ 1078 int* termbegins, /**< pointer to store the termbegins */ 1079 SCIP_Real* transcoefs, /**< pointer to store the transcoefs */ 1080 int* transcoefsidx, /**< pointer to store the transcoefsidx */ 1081 SCIP_Real* offsets, /**< pointer to store the offsets */ 1082 SCIP_Real* lhsconstant, /**< pointer to store the lhsconstant */ 1083 int* nterms, /**< pointer to store the total number of terms */ 1084 SCIP_Bool* success /**< whether the expression is indeed a SOC */ 1085 ) 1086 { 1087 SCIP_Real sqrteigval; 1088 int nextterm = 0; 1089 int nexttranscoef = 0; 1090 int specialtermidx; 1091 int i; 1092 int j; 1093 1094 assert(scip != NULL); 1095 assert(occurringexprs != NULL); 1096 assert(eigvecmatrix != NULL); 1097 assert(eigvals != NULL); 1098 assert(bp != NULL); 1099 assert(termbegins != NULL); 1100 assert(transcoefs != NULL); 1101 assert(transcoefsidx != NULL); 1102 assert(offsets != NULL); 1103 assert(lhsconstant != NULL); 1104 assert(success != NULL); 1105 1106 *success = FALSE; 1107 *nterms = 0; 1108 1109 /* we have lhsconstant + x^t A x + b x <= 0 and A has a single negative eigenvalue; try to build soc; 1110 * we now store all the v_i^T x + beta_i on the lhs, and compute the constant 1111 */ 1112 specialtermidx = -1; 1113 for( i = 0; i < nvars; ++i ) 1114 { 1115 if( SCIPisZero(scip, eigvals[i]) ) 1116 continue; 1117 1118 if( eigvals[i] < 0.0 ) 1119 { 1120 assert(specialtermidx == -1); /* there should only be one negative eigenvalue */ 1121 1122 specialtermidx = i; 1123 1124 *lhsconstant -= bp[i] * bp[i] / (4.0 * eigvals[i]); 1125 1126 continue; 1127 } 1128 1129 assert(eigvals[i] > 0.0); 1130 sqrteigval = SQRT(eigvals[i]); 1131 1132 termbegins[nextterm] = nexttranscoef; 1133 offsets[nextterm] = bp[i] / (2.0 * sqrteigval); 1134 *lhsconstant -= bp[i] * bp[i] / (4.0 * eigvals[i]); 1135 1136 /* set transcoefs */ 1137 for( j = 0; j < nvars; ++j ) 1138 { 1139 if( !SCIPisZero(scip, eigvecmatrix[i * nvars + j]) ) 1140 { 1141 transcoefs[nexttranscoef] = sqrteigval * eigvecmatrix[i * nvars + j]; 1142 transcoefsidx[nexttranscoef] = j; 1143 1144 ++nexttranscoef; 1145 } 1146 } 1147 ++nextterm; 1148 } 1149 assert(specialtermidx > -1); 1150 1151 /* process constant; if constant is negative -> no soc */ 1152 if( SCIPisNegative(scip, *lhsconstant) ) 1153 return SCIP_OKAY; 1154 1155 /* we need lhsconstant to be >= 0 */ 1156 if( *lhsconstant < 0.0 ) 1157 *lhsconstant = 0.0; 1158 1159 /* store constant term */ 1160 if( *lhsconstant > 0.0 ) 1161 { 1162 termbegins[nextterm] = nexttranscoef; 1163 offsets[nextterm] = SQRT( *lhsconstant ); 1164 ++nextterm; 1165 } 1166 1167 /* now process rhs */ 1168 { 1169 SCIP_Real rhstermlb; 1170 SCIP_Real rhstermub; 1171 SCIP_Real signfactor; 1172 1173 assert(-eigvals[specialtermidx] > 0.0); 1174 sqrteigval = SQRT(-eigvals[specialtermidx]); 1175 1176 termbegins[nextterm] = nexttranscoef; 1177 offsets[nextterm] = -bp[specialtermidx] / (2.0 * sqrteigval); 1178 1179 /* the expression can only be an soc if the resulting rhs term does not change sign; 1180 * the rhs term is a linear combination of variables, so estimate its bounds 1181 */ 1182 rhstermlb = offsets[nextterm]; 1183 for( j = 0; j < nvars; ++j ) 1184 { 1185 SCIP_INTERVAL activity; 1186 SCIP_Real aux; 1187 1188 if( SCIPisZero(scip, eigvecmatrix[specialtermidx * nvars + j]) ) 1189 continue; 1190 1191 SCIP_CALL( SCIPevalExprActivity(scip, occurringexprs[j]) ); 1192 activity = SCIPexprGetActivity(occurringexprs[j]); 1193 1194 if( eigvecmatrix[specialtermidx * nvars + j] > 0.0 ) 1195 { 1196 aux = activity.inf; 1197 assert(!SCIPisInfinity(scip, aux)); 1198 } 1199 else 1200 { 1201 aux = activity.sup; 1202 assert(!SCIPisInfinity(scip, -aux)); 1203 } 1204 1205 if( SCIPisInfinity(scip, aux) || SCIPisInfinity(scip, -aux) ) 1206 { 1207 rhstermlb = -SCIPinfinity(scip); 1208 break; 1209 } 1210 else 1211 rhstermlb += sqrteigval * eigvecmatrix[specialtermidx * nvars + j] * aux; 1212 } 1213 1214 rhstermub = offsets[nextterm]; 1215 for( j = 0; j < nvars; ++j ) 1216 { 1217 SCIP_INTERVAL activity; 1218 SCIP_Real aux; 1219 1220 if( SCIPisZero(scip, eigvecmatrix[specialtermidx * nvars + j]) ) 1221 continue; 1222 1223 SCIP_CALL( SCIPevalExprActivity(scip, occurringexprs[j]) ); 1224 activity = SCIPexprGetActivity(occurringexprs[j]); 1225 1226 if( eigvecmatrix[specialtermidx * nvars + j] > 0.0 ) 1227 { 1228 aux = activity.sup; 1229 assert(!SCIPisInfinity(scip, -aux)); 1230 } 1231 else 1232 { 1233 aux = activity.inf; 1234 assert(!SCIPisInfinity(scip, aux)); 1235 } 1236 1237 if( SCIPisInfinity(scip, aux) || SCIPisInfinity(scip, -aux) ) 1238 { 1239 rhstermub = SCIPinfinity(scip); 1240 break; 1241 } 1242 else 1243 rhstermub += sqrteigval * eigvecmatrix[specialtermidx * nvars + j] * aux; 1244 } 1245 1246 /* since we are just interested in obtaining an interval that contains the real bounds 1247 * and is tight enough so that we can identify that the rhsvar does not change sign, 1248 * we swap the bounds in case of numerical troubles 1249 */ 1250 if( rhstermub < rhstermlb ) 1251 { 1252 assert(SCIPisEQ(scip, rhstermub, rhstermlb)); 1253 SCIPswapReals(&rhstermub, &rhstermlb); 1254 } 1255 1256 /* if rhs changes sign -> not a SOC */ 1257 if( SCIPisLT(scip, rhstermlb, 0.0) && SCIPisGT(scip, rhstermub, 0.0) ) 1258 return SCIP_OKAY; 1259 1260 signfactor = SCIPisLE(scip, rhstermub, 0.0) ? -1.0 : 1.0; 1261 1262 offsets[nextterm] *= signfactor; 1263 1264 /* set transcoefs for rhs term */ 1265 for( j = 0; j < nvars; ++j ) 1266 { 1267 if( SCIPisZero(scip, eigvecmatrix[specialtermidx * nvars + j]) ) 1268 continue; 1269 1270 transcoefs[nexttranscoef] = signfactor * sqrteigval * eigvecmatrix[specialtermidx * nvars + j]; 1271 transcoefsidx[nexttranscoef] = j; 1272 1273 ++nexttranscoef; 1274 } 1275 1276 /* if rhs is a constant this method shouldn't have been called */ 1277 assert(nexttranscoef > termbegins[nextterm]); 1278 1279 /* finish processing term */ 1280 ++nextterm; 1281 } 1282 1283 *nterms = nextterm; 1284 1285 /* sentinel value */ 1286 termbegins[nextterm] = nexttranscoef; 1287 1288 *success = TRUE; 1289 1290 return SCIP_OKAY; 1291 } 1292 1293 /** detects if expr ≤ auxvar is of the form SQRT(sum_i coef_i (expr_i + shift_i)^2 + const) ≤ auxvar 1294 * 1295 * @note if a user inputs the above expression with `const` = -epsilon, then `const` is going to be set to 0. 1296 */ 1297 static 1298 SCIP_RETCODE detectSocNorm( 1299 SCIP* scip, /**< SCIP data structure */ 1300 SCIP_EXPR* expr, /**< expression */ 1301 SCIP_NLHDLREXPRDATA** nlhdlrexprdata, /**< pointer to store nonlinear handler expression data */ 1302 SCIP_Bool* success /**< pointer to store whether SOC structure has been detected */ 1303 ) 1304 { 1305 SCIP_EXPR** children; 1306 SCIP_EXPR* child; 1307 SCIP_EXPR** vars; 1308 SCIP_HASHMAP* expr2idx; 1309 SCIP_HASHSET* linexprs; 1310 SCIP_Real* childcoefs; 1311 SCIP_Real* offsets; 1312 SCIP_Real* transcoefs; 1313 SCIP_Real constant; 1314 SCIP_Bool issoc; 1315 int* transcoefsidx; 1316 int* termbegins; 1317 int nchildren; 1318 int nterms; 1319 int nvars; 1320 int nextentry; 1321 int i; 1322 1323 assert(expr != NULL); 1324 assert(success != NULL); 1325 1326 *success = FALSE; 1327 issoc = TRUE; 1328 1329 /* relation is not "<=" -> skip */ 1330 if( SCIPgetExprNLocksPosNonlinear(expr) == 0 ) 1331 return SCIP_OKAY; 1332 1333 assert(SCIPexprGetNChildren(expr) > 0); 1334 1335 child = SCIPexprGetChildren(expr)[0]; 1336 assert(child != NULL); 1337 1338 /* check whether expression is a SQRT and has a sum as child with at least 2 children and a non-negative constant */ 1339 if( ! SCIPisExprPower(scip, expr) 1340 || SCIPgetExponentExprPow(expr) != 0.5 1341 || !SCIPisExprSum(scip, child) 1342 || SCIPexprGetNChildren(child) < 2 1343 || SCIPgetConstantExprSum(child) < 0.0) 1344 { 1345 return SCIP_OKAY; 1346 } 1347 1348 /* assert(SCIPvarGetLbLocal(auxvar) >= 0.0); */ 1349 1350 /* get children of the sum */ 1351 children = SCIPexprGetChildren(child); 1352 nchildren = SCIPexprGetNChildren(child); 1353 childcoefs = SCIPgetCoefsExprSum(child); 1354 1355 /* TODO: should we initialize the hashmap with size SCIPgetNVars() so that it never has to be resized? */ 1356 SCIP_CALL( SCIPhashmapCreate(&expr2idx, SCIPblkmem(scip), nchildren) ); 1357 SCIP_CALL( SCIPhashsetCreate(&linexprs, SCIPblkmem(scip), nchildren) ); 1358 1359 /* we create coefs array here already, since we have to fill it in first loop in case of success 1360 * +1 for auxvar 1361 */ 1362 SCIP_CALL( SCIPallocBufferArray(scip, &transcoefs, nchildren+1) ); 1363 1364 nterms = 0; 1365 1366 /* check if all children are squares or linear terms with matching square term: 1367 * if the i-th child is (pow, expr, 2) we store the association <|expr -> i|> in expr2idx and if expr was in 1368 * linexprs, we remove it from there. 1369 * if the i-th child is expr' (different from (pow, expr, 2)) and expr' is not a key of expr2idx, we add it 1370 * to linexprs. 1371 * if at the end there is any expr in linexpr -> we do not have a separable quadratic function. 1372 */ 1373 for( i = 0; i < nchildren; ++i ) 1374 { 1375 /* handle quadratic expressions children */ 1376 if( SCIPisExprPower(scip, children[i]) && SCIPgetExponentExprPow(children[i]) == 2.0 ) 1377 { 1378 SCIP_EXPR* squarearg = SCIPexprGetChildren(children[i])[0]; 1379 1380 if( !SCIPhashmapExists(expr2idx, (void*) squarearg) ) 1381 { 1382 SCIP_CALL( SCIPhashmapInsertInt(expr2idx, (void *) squarearg, nterms) ); 1383 } 1384 1385 if( childcoefs[i] < 0.0 ) 1386 { 1387 issoc = FALSE; 1388 break; 1389 } 1390 transcoefs[nterms] = SQRT(childcoefs[i]); 1391 1392 SCIP_CALL( SCIPhashsetRemove(linexprs, (void*) squarearg) ); 1393 ++nterms; 1394 } 1395 /* handle binary variable children */ 1396 else if( SCIPisExprVar(scip, children[i]) && SCIPvarIsBinary(SCIPgetVarExprVar(children[i])) ) 1397 { 1398 assert(!SCIPhashmapExists(expr2idx, (void*) children[i])); 1399 assert(!SCIPhashsetExists(linexprs, (void*) children[i])); 1400 1401 SCIP_CALL( SCIPhashmapInsertInt(expr2idx, (void *) children[i], nterms) ); 1402 1403 if( childcoefs[i] < 0.0 ) 1404 { 1405 issoc = FALSE; 1406 break; 1407 } 1408 transcoefs[nterms] = SQRT(childcoefs[i]); 1409 1410 ++nterms; 1411 } 1412 else 1413 { 1414 if( !SCIPhashmapExists(expr2idx, (void*) children[i]) ) 1415 { 1416 SCIP_CALL( SCIPhashsetInsert(linexprs, SCIPblkmem(scip), (void*) children[i]) ); 1417 } 1418 } 1419 } 1420 1421 /* there are linear terms without corresponding quadratic terms or it was detected not to be soc */ 1422 if( SCIPhashsetGetNElements(linexprs) > 0 || ! issoc ) 1423 { 1424 SCIPfreeBufferArray(scip, &transcoefs); 1425 SCIPhashsetFree(&linexprs, SCIPblkmem(scip) ); 1426 SCIPhashmapFree(&expr2idx); 1427 return SCIP_OKAY; 1428 } 1429 1430 /* add one to terms counter for auxvar */ 1431 ++nterms; 1432 1433 constant = SCIPgetConstantExprSum(child); 1434 1435 /* compute constant of possible soc expression to check its sign */ 1436 for( i = 0; i < nchildren; ++i ) 1437 { 1438 if( ! SCIPisExprPower(scip, children[i]) || SCIPgetExponentExprPow(children[i]) != 2.0 ) 1439 { 1440 int auxvarpos; 1441 1442 assert(SCIPhashmapExists(expr2idx, (void*) children[i]) ); 1443 auxvarpos = SCIPhashmapGetImageInt(expr2idx, (void*) children[i]); 1444 1445 constant -= SQR(0.5 * childcoefs[i] / transcoefs[auxvarpos]); 1446 } 1447 } 1448 1449 /* if the constant is negative -> no SOC */ 1450 if( SCIPisNegative(scip, constant) ) 1451 { 1452 SCIPfreeBufferArray(scip, &transcoefs); 1453 SCIPhashsetFree(&linexprs, SCIPblkmem(scip) ); 1454 SCIPhashmapFree(&expr2idx); 1455 return SCIP_OKAY; 1456 } 1457 else if( SCIPisZero(scip, constant) ) 1458 constant = 0.0; 1459 assert(constant >= 0.0); 1460 1461 /* at this point, we have found an SOC structure */ 1462 *success = TRUE; 1463 1464 nvars = nterms; 1465 1466 /* add one to terms counter for constant term */ 1467 if( constant > 0.0 ) 1468 ++nterms; 1469 1470 /* allocate temporary memory to collect data */ 1471 SCIP_CALL( SCIPallocBufferArray(scip, &vars, nvars) ); 1472 SCIP_CALL( SCIPallocBufferArray(scip, &offsets, nterms) ); 1473 SCIP_CALL( SCIPallocBufferArray(scip, &transcoefsidx, nvars) ); 1474 SCIP_CALL( SCIPallocBufferArray(scip, &termbegins, nterms + 1) ); 1475 1476 /* fill in data for non constant terms of lhs; initialize their offsets */ 1477 for( i = 0; i < nvars - 1; ++i ) 1478 { 1479 transcoefsidx[i] = i; 1480 termbegins[i] = i; 1481 offsets[i] = 0.0; 1482 } 1483 1484 /* add constant term and rhs */ 1485 vars[nvars - 1] = expr; 1486 if( constant > 0.0 ) 1487 { 1488 /* constant term */ 1489 termbegins[nterms - 2] = nterms - 2; 1490 offsets[nterms - 2] = SQRT(constant); 1491 1492 /* rhs */ 1493 termbegins[nterms - 1] = nterms - 2; 1494 offsets[nterms - 1] = 0.0; 1495 transcoefsidx[nterms - 2] = nvars - 1; 1496 transcoefs[nterms - 2] = 1.0; 1497 1498 /* sentinel value */ 1499 termbegins[nterms] = nterms - 1; 1500 } 1501 else 1502 { 1503 /* rhs */ 1504 termbegins[nterms - 1] = nterms - 1; 1505 offsets[nterms - 1] = 0.0; 1506 transcoefsidx[nterms - 1] = nvars - 1; 1507 transcoefs[nterms - 1] = 1.0; 1508 1509 /* sentinel value */ 1510 termbegins[nterms] = nterms; 1511 } 1512 1513 /* request required auxiliary variables and fill vars and offsets array */ 1514 nextentry = 0; 1515 for( i = 0; i < nchildren; ++i ) 1516 { 1517 if( SCIPisExprPower(scip, children[i]) && SCIPgetExponentExprPow(children[i]) == 2.0 ) 1518 { 1519 SCIP_EXPR* squarearg; 1520 1521 squarearg = SCIPexprGetChildren(children[i])[0]; 1522 assert(SCIPhashmapGetImageInt(expr2idx, (void*) squarearg) == nextentry); 1523 1524 SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, squarearg, TRUE, FALSE, FALSE, FALSE) ); 1525 1526 vars[nextentry] = squarearg; 1527 ++nextentry; 1528 } 1529 else if( SCIPisExprVar(scip, children[i]) && SCIPvarIsBinary(SCIPgetVarExprVar(children[i])) ) 1530 { 1531 /* handle binary variable children: no need to request auxvar */ 1532 assert(SCIPhashmapGetImageInt(expr2idx, (void*) children[i]) == nextentry); 1533 vars[nextentry] = children[i]; 1534 ++nextentry; 1535 } 1536 else 1537 { 1538 int auxvarpos; 1539 1540 assert(SCIPhashmapExists(expr2idx, (void*) children[i])); 1541 auxvarpos = SCIPhashmapGetImageInt(expr2idx, (void*) children[i]); 1542 1543 SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, children[i], TRUE, FALSE, FALSE, FALSE) ); 1544 1545 offsets[auxvarpos] = 0.5 * childcoefs[i] / transcoefs[auxvarpos]; 1546 } 1547 } 1548 assert(nextentry == nvars - 1); 1549 1550 #ifdef SCIP_DEBUG 1551 SCIPdebugMsg(scip, "found SOC structure for expression %p\n", (void*)expr); 1552 SCIPprintExpr(scip, expr, NULL); 1553 SCIPinfoMessage(scip, NULL, " <= auxvar\n"); 1554 #endif 1555 1556 /* create and store nonlinear handler expression data */ 1557 SCIP_CALL( createNlhdlrExprData(scip, vars, offsets, transcoefs, transcoefsidx, termbegins, 1558 nvars, nterms, nlhdlrexprdata) ); 1559 assert(*nlhdlrexprdata != NULL); 1560 1561 /* free memory */ 1562 SCIPhashsetFree(&linexprs, SCIPblkmem(scip) ); 1563 SCIPhashmapFree(&expr2idx); 1564 SCIPfreeBufferArray(scip, &termbegins); 1565 SCIPfreeBufferArray(scip, &transcoefsidx); 1566 SCIPfreeBufferArray(scip, &offsets); 1567 SCIPfreeBufferArray(scip, &vars); 1568 SCIPfreeBufferArray(scip, &transcoefs); 1569 1570 return SCIP_OKAY; 1571 } 1572 1573 /** helper method to detect c + sum_i coef_i expr_i^2 - coef_k expr_k^2 ≤ 0 1574 * and c + sum_i coef_i expr_i^2 - coef_k expr_k expr_l ≤ 0 1575 * 1576 * binary linear variables are interpreted as quadratic terms 1577 * 1578 * @todo: extend this function to detect c + sum_i coef_i (expr_i + const_i)^2 - ... 1579 * this would probably share a lot of code with detectSocNorm 1580 */ 1581 static 1582 SCIP_RETCODE detectSocQuadraticSimple( 1583 SCIP* scip, /**< SCIP data structure */ 1584 SCIP_EXPR* expr, /**< expression */ 1585 SCIP_Real conslhs, /**< lhs of the constraint that the expression defines (or SCIP_INVALID) */ 1586 SCIP_Real consrhs, /**< rhs of the constraint that the expression defines (or SCIP_INVALID) */ 1587 SCIP_NLHDLREXPRDATA** nlhdlrexprdata, /**< pointer to store nonlinear handler expression data */ 1588 SCIP_Bool* enforcebelow, /**< pointer to store whether we enforce <= (TRUE) or >= (FALSE); only valid when success is TRUE */ 1589 SCIP_Bool* success /**< pointer to store whether SOC structure has been detected */ 1590 ) 1591 { 1592 SCIP_EXPR** children; 1593 SCIP_EXPR** vars = NULL; 1594 SCIP_Real* childcoefs; 1595 SCIP_Real* offsets = NULL; 1596 SCIP_Real* transcoefs = NULL; 1597 int* transcoefsidx = NULL; 1598 int* termbegins = NULL; 1599 SCIP_Real constant; 1600 SCIP_Real lhsconstant; 1601 SCIP_Real lhs; 1602 SCIP_Real rhs; 1603 SCIP_Real rhssign; 1604 SCIP_INTERVAL expractivity; 1605 int ntranscoefs; 1606 int nposquadterms; 1607 int nnegquadterms; 1608 int nposbilinterms; 1609 int nnegbilinterms; 1610 int rhsidx; 1611 int lhsidx; 1612 int specialtermidx; 1613 int nchildren; 1614 int nnzinterms; 1615 int nterms; 1616 int nvars; 1617 int nextentry; 1618 int i; 1619 SCIP_Bool ishyperbolic; 1620 1621 assert(expr != NULL); 1622 assert(success != NULL); 1623 1624 *success = FALSE; 1625 1626 /* check whether expression is a sum with at least 2 children */ 1627 if( ! SCIPisExprSum(scip, expr) || SCIPexprGetNChildren(expr) < 2 ) 1628 return SCIP_OKAY; 1629 1630 /* get children of the sum */ 1631 children = SCIPexprGetChildren(expr); 1632 nchildren = SCIPexprGetNChildren(expr); 1633 constant = SCIPgetConstantExprSum(expr); 1634 1635 /* we duplicate the child coefficients since we have to manipulate them */ 1636 SCIP_CALL( SCIPduplicateBufferArray(scip, &childcoefs, SCIPgetCoefsExprSum(expr), nchildren) ); /*lint !e666*/ 1637 1638 /* initialize data */ 1639 lhsidx = -1; 1640 rhsidx = -1; 1641 nposquadterms = 0; 1642 nnegquadterms = 0; 1643 nposbilinterms = 0; 1644 nnegbilinterms = 0; 1645 1646 /* check if all children are quadratic or binary linear and count number of positive and negative terms */ 1647 for( i = 0; i < nchildren; ++i ) 1648 { 1649 if( SCIPisExprPower(scip, children[i]) && SCIPgetExponentExprPow(children[i]) == 2.0 ) 1650 { 1651 if( childcoefs[i] > 0.0 ) 1652 { 1653 ++nposquadterms; 1654 lhsidx = i; 1655 } 1656 else 1657 { 1658 ++nnegquadterms; 1659 rhsidx = i; 1660 } 1661 } 1662 else if( SCIPisExprVar(scip, children[i]) && SCIPvarIsBinary(SCIPgetVarExprVar(children[i])) ) 1663 { 1664 if( childcoefs[i] > 0.0 ) 1665 { 1666 ++nposquadterms; 1667 lhsidx = i; 1668 } 1669 else 1670 { 1671 ++nnegquadterms; 1672 rhsidx = i; 1673 } 1674 } 1675 else if( SCIPisExprProduct(scip, children[i]) && SCIPexprGetNChildren(children[i]) == 2 ) 1676 { 1677 if( childcoefs[i] > 0.0 ) 1678 { 1679 ++nposbilinterms; 1680 lhsidx = i; 1681 } 1682 else 1683 { 1684 ++nnegbilinterms; 1685 rhsidx = i; 1686 } 1687 } 1688 else 1689 { 1690 goto CLEANUP; 1691 } 1692 1693 /* more than one positive eigenvalue and more than one negative eigenvalue -> can't be convex */ 1694 if( nposquadterms > 1 && nnegquadterms > 1 ) 1695 goto CLEANUP; 1696 1697 /* more than one bilinear term -> can't be handled by this method */ 1698 if( nposbilinterms + nnegbilinterms > 1 ) 1699 goto CLEANUP; 1700 1701 /* one positive bilinear term and also at least one positive quadratic term -> not a simple SOC */ 1702 if( nposbilinterms > 0 && nposquadterms > 0 ) 1703 goto CLEANUP; 1704 1705 /* one negative bilinear term and also at least one negative quadratic term -> not a simple SOC */ 1706 if( nnegbilinterms > 0 && nnegquadterms > 0 ) 1707 goto CLEANUP; 1708 } 1709 1710 if( nposquadterms == nchildren || nnegquadterms == nchildren ) 1711 goto CLEANUP; 1712 1713 assert(nposquadterms <= 1 || nnegquadterms <= 1); 1714 assert(nposbilinterms + nnegbilinterms <= 1); 1715 assert(nposbilinterms == 0 || nposquadterms == 0); 1716 assert(nnegbilinterms == 0 || nnegquadterms == 0); 1717 1718 /* if a bilinear term is involved, it is a hyperbolic expression */ 1719 ishyperbolic = (nposbilinterms + nnegbilinterms > 0); 1720 1721 if( conslhs == SCIP_INVALID || consrhs == SCIP_INVALID ) /*lint !e777*/ 1722 { 1723 SCIP_CALL( SCIPevalExprActivity(scip, expr) ); 1724 expractivity = SCIPexprGetActivity(expr); 1725 1726 lhs = (conslhs == SCIP_INVALID ? expractivity.inf : conslhs); /*lint !e777*/ 1727 rhs = (consrhs == SCIP_INVALID ? expractivity.sup : consrhs); /*lint !e777*/ 1728 } 1729 else 1730 { 1731 lhs = conslhs; 1732 rhs = consrhs; 1733 } 1734 1735 /* detect case and store lhs/rhs information */ 1736 if( (ishyperbolic && nnegbilinterms > 0) || (!ishyperbolic && nnegquadterms < 2) ) 1737 { 1738 /* we have -x*y + z^2 ... -> we want to write z^2 ... <= x*y; 1739 * or we have -x^2 + y^2 ... -> we want to write y^2 ... <= x^2; 1740 * in any case, we need a finite rhs 1741 */ 1742 assert(nnegbilinterms == 1 || nnegquadterms == 1); 1743 assert(rhsidx != -1); 1744 1745 /* if rhs is infinity, it can't be soc 1746 * TODO: if it can't be soc, then we should enforce the caller so that we do not try the more complex quadratic 1747 * method 1748 */ 1749 if( SCIPisInfinity(scip, rhs) ) 1750 goto CLEANUP; 1751 1752 specialtermidx = rhsidx; 1753 lhsconstant = constant - rhs; 1754 *enforcebelow = TRUE; /* enforce expr <= rhs */ 1755 } 1756 else 1757 { 1758 assert(lhsidx != -1); 1759 1760 /* if lhs is infinity, it can't be soc */ 1761 if( SCIPisInfinity(scip, -lhs) ) 1762 goto CLEANUP; 1763 1764 specialtermidx = lhsidx; 1765 lhsconstant = lhs - constant; 1766 1767 /* negate all coefficients */ 1768 for( i = 0; i < nchildren; ++i ) 1769 childcoefs[i] = -childcoefs[i]; 1770 *enforcebelow = FALSE; /* enforce lhs <= expr */ 1771 } 1772 assert(childcoefs[specialtermidx] != 0.0); 1773 1774 if( ishyperbolic ) 1775 { 1776 SCIP_INTERVAL yactivity; 1777 SCIP_INTERVAL zactivity; 1778 1779 assert(SCIPexprGetNChildren(children[specialtermidx]) == 2); 1780 1781 SCIP_CALL( SCIPevalExprActivity(scip, SCIPexprGetChildren(children[specialtermidx])[0]) ); 1782 yactivity = SCIPexprGetActivity(SCIPexprGetChildren(children[specialtermidx])[0]); 1783 1784 SCIP_CALL( SCIPevalExprActivity(scip, SCIPexprGetChildren(children[specialtermidx])[1]) ); 1785 zactivity = SCIPexprGetActivity(SCIPexprGetChildren(children[specialtermidx])[1]); 1786 1787 if( SCIPisNegative(scip, yactivity.inf + zactivity.inf) ) 1788 { 1789 /* the sum of the expressions in the bilinear term changes sign -> no SOC */ 1790 if( SCIPisPositive(scip, yactivity.sup + zactivity.sup) ) 1791 goto CLEANUP; 1792 1793 rhssign = -1.0; 1794 } 1795 else 1796 rhssign = 1.0; 1797 1798 lhsconstant *= 4.0 / -childcoefs[specialtermidx]; 1799 } 1800 else if( SCIPisExprVar(scip, children[specialtermidx]) ) 1801 { 1802 /* children[specialtermidx] can be a variable, in which case we treat it as if it is squared */ 1803 rhssign = 1.0; 1804 } 1805 else 1806 { 1807 SCIP_INTERVAL rhsactivity; 1808 1809 assert(SCIPexprGetNChildren(children[specialtermidx]) == 1); 1810 SCIP_CALL( SCIPevalExprActivity(scip, SCIPexprGetChildren(children[specialtermidx])[0]) ); 1811 rhsactivity = SCIPexprGetActivity(SCIPexprGetChildren(children[specialtermidx])[0]); 1812 1813 if( rhsactivity.inf < 0.0 ) 1814 { 1815 /* rhs variable changes sign -> no SOC */ 1816 if( rhsactivity.sup > 0.0 ) 1817 goto CLEANUP; 1818 1819 rhssign = -1.0; 1820 } 1821 else 1822 rhssign = 1.0; 1823 } 1824 1825 if( SCIPisNegative(scip, lhsconstant) ) 1826 goto CLEANUP; 1827 1828 if( SCIPisZero(scip, lhsconstant) ) 1829 lhsconstant = 0.0; 1830 1831 /* 1832 * we have found an SOC-representable expression. Now build the nlhdlrexprdata 1833 * 1834 * in the non-hyperbolic case, c + sum_i coef_i expr_i^2 - coef_k expr_k^2 <= 0 is transformed to 1835 * SQRT( c + sum_i coef_i expr_i^2 ) <= coef_k expr_k 1836 * so there are nchildren many vars, nchildren (+ 1 if c != 0) many terms, nchildren many coefficients in the vs 1837 * in SOC representation 1838 * 1839 * in the hyperbolic case, c + sum_i coef_i expr_i^2 - coef_k expr_k expr_l <= 0 is transformed to 1840 * SQRT( 4(c + sum_i coef_i expr_i^2) + (expr_k - expr_l)^2 ) <= expr_k + expr_l 1841 * so there are nchildren + 1many vars, nchildren + 1(+ 1 if c != 0) many terms, nchildren + 3 many coefficients in 1842 * the vs in SOC representation 1843 */ 1844 1845 ntranscoefs = ishyperbolic ? nchildren + 3 : nchildren; 1846 nvars = ishyperbolic ? nchildren + 1 : nchildren; 1847 nterms = nvars; 1848 1849 /* constant term */ 1850 if( lhsconstant > 0.0 ) 1851 nterms++; 1852 1853 /* SOC was detected, allocate temporary memory for data to collect */ 1854 SCIP_CALL( SCIPallocBufferArray(scip, &vars, nvars) ); 1855 SCIP_CALL( SCIPallocBufferArray(scip, &offsets, nterms) ); 1856 SCIP_CALL( SCIPallocBufferArray(scip, &transcoefs, ntranscoefs) ); 1857 SCIP_CALL( SCIPallocBufferArray(scip, &transcoefsidx, ntranscoefs) ); 1858 SCIP_CALL( SCIPallocBufferArray(scip, &termbegins, nterms + 1) ); 1859 1860 *success = TRUE; 1861 nextentry = 0; 1862 1863 /* collect all the v_i and beta_i */ 1864 nnzinterms = 0; 1865 for( i = 0; i < nchildren; ++i ) 1866 { 1867 /* variable and coef for rhs have to be set to the last entry */ 1868 if( i == specialtermidx ) 1869 continue; 1870 1871 /* extract (unique) variable appearing in term */ 1872 if( SCIPisExprVar(scip, children[i]) ) 1873 { 1874 vars[nextentry] = children[i]; 1875 1876 assert(SCIPvarIsBinary(SCIPgetVarExprVar(vars[nextentry]))); 1877 } 1878 else 1879 { 1880 assert(SCIPisExprPower(scip, children[i])); 1881 1882 /* notify that we will require auxiliary variable */ 1883 SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, SCIPexprGetChildren(children[i])[0], TRUE, FALSE, FALSE, FALSE) ); 1884 vars[nextentry] = SCIPexprGetChildren(children[i])[0]; 1885 } 1886 assert(vars[nextentry] != NULL); 1887 1888 /* store v_i and beta_i */ 1889 termbegins[nextentry] = nnzinterms; 1890 offsets[nextentry] = 0.0; 1891 1892 transcoefsidx[nnzinterms] = nextentry; 1893 if( ishyperbolic ) 1894 { 1895 /* we eliminate the coefficient of the bilinear term to arrive at standard form */ 1896 assert(4.0 * childcoefs[i] / -childcoefs[specialtermidx] > 0.0); 1897 transcoefs[nnzinterms] = SQRT(4.0 * childcoefs[i] / -childcoefs[specialtermidx]); 1898 } 1899 else 1900 { 1901 assert(childcoefs[i] > 0.0); 1902 transcoefs[nnzinterms] = SQRT(childcoefs[i]); 1903 } 1904 1905 /* finish adding nonzeros */ 1906 ++nnzinterms; 1907 1908 /* finish processing term */ 1909 ++nextentry; 1910 } 1911 assert(nextentry == nchildren - 1); 1912 1913 /* store term for constant (v_i = 0) */ 1914 if( lhsconstant > 0.0 ) 1915 { 1916 termbegins[nextentry] = nnzinterms; 1917 offsets[nextentry] = SQRT(lhsconstant); 1918 1919 /* finish processing term; this term has 0 nonzero thus we do not increase nnzinterms */ 1920 ++nextentry; 1921 } 1922 1923 if( !ishyperbolic ) 1924 { 1925 /* store rhs term */ 1926 if( SCIPisExprVar(scip, children[specialtermidx]) ) 1927 { 1928 /* this should be the "children[specialtermidx] can be a variable, in which case we treat it as if it is squared" case */ 1929 SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, children[specialtermidx], TRUE, FALSE, FALSE, FALSE) ); 1930 vars[nvars - 1] = children[specialtermidx]; 1931 } 1932 else 1933 { 1934 assert(SCIPisExprPower(scip, children[specialtermidx])); 1935 assert(SCIPexprGetChildren(children[specialtermidx]) != NULL); 1936 SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, SCIPexprGetChildren(children[specialtermidx])[0], TRUE, FALSE, FALSE, FALSE) ); 1937 vars[nvars - 1] = SCIPexprGetChildren(children[specialtermidx])[0]; 1938 } 1939 1940 assert(childcoefs[specialtermidx] < 0.0); 1941 1942 termbegins[nextentry] = nnzinterms; 1943 offsets[nextentry] = 0.0; 1944 transcoefs[nnzinterms] = rhssign * SQRT(-childcoefs[specialtermidx]); 1945 transcoefsidx[nnzinterms] = nvars - 1; 1946 1947 /* finish adding nonzeros */ 1948 ++nnzinterms; 1949 1950 /* finish processing term */ 1951 ++nextentry; 1952 } 1953 else 1954 { 1955 /* store last lhs term and rhs term coming from the bilinear term */ 1956 SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, SCIPexprGetChildren(children[specialtermidx])[0], TRUE, FALSE, FALSE, FALSE) ); 1957 vars[nvars - 2] = SCIPexprGetChildren(children[specialtermidx])[0]; 1958 1959 SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, SCIPexprGetChildren(children[specialtermidx])[1], TRUE, FALSE, FALSE, FALSE) ); 1960 vars[nvars - 1] = SCIPexprGetChildren(children[specialtermidx])[1]; 1961 1962 /* at this point, vars[nvars - 2] = expr_k and vars[nvars - 1] = expr_l; 1963 * on the lhs we have the term (expr_k - expr_l)^2 1964 */ 1965 termbegins[nextentry] = nnzinterms; 1966 offsets[nextentry] = 0.0; 1967 1968 /* expr_k */ 1969 transcoefsidx[nnzinterms] = nvars - 2; 1970 transcoefs[nnzinterms] = 1.0; 1971 ++nnzinterms; 1972 1973 /* - expr_l */ 1974 transcoefsidx[nnzinterms] = nvars - 1; 1975 transcoefs[nnzinterms] = -1.0; 1976 ++nnzinterms; 1977 1978 /* finish processing term */ 1979 ++nextentry; 1980 1981 /* on rhs we have +/-(expr_k + expr_l) */ 1982 termbegins[nextentry] = nnzinterms; 1983 offsets[nextentry] = 0.0; 1984 1985 /* rhssing * expr_k */ 1986 transcoefsidx[nnzinterms] = nvars - 2; 1987 transcoefs[nnzinterms] = rhssign; 1988 ++nnzinterms; 1989 1990 /* rhssing * expr_l */ 1991 transcoefsidx[nnzinterms] = nvars - 1; 1992 transcoefs[nnzinterms] = rhssign; 1993 ++nnzinterms; 1994 1995 /* finish processing term */ 1996 ++nextentry; 1997 } 1998 assert(nextentry == nterms); 1999 assert(nnzinterms == ntranscoefs); 2000 2001 /* sentinel value */ 2002 termbegins[nextentry] = nnzinterms; 2003 2004 #ifdef SCIP_DEBUG 2005 SCIPdebugMsg(scip, "found SOC structure for expression %p\n %g <= ", (void*)expr, lhs); 2006 SCIPprintExpr(scip, expr, NULL); 2007 SCIPinfoMessage(scip, NULL, " <= %g\n", rhs); 2008 #endif 2009 2010 /* create and store nonlinear handler expression data */ 2011 SCIP_CALL( createNlhdlrExprData(scip, vars, offsets, transcoefs, transcoefsidx, termbegins, nvars, nterms, 2012 nlhdlrexprdata) ); 2013 assert(*nlhdlrexprdata != NULL); 2014 2015 CLEANUP: 2016 SCIPfreeBufferArrayNull(scip, &termbegins); 2017 SCIPfreeBufferArrayNull(scip, &transcoefsidx); 2018 SCIPfreeBufferArrayNull(scip, &transcoefs); 2019 SCIPfreeBufferArrayNull(scip, &offsets); 2020 SCIPfreeBufferArrayNull(scip, &vars); 2021 SCIPfreeBufferArrayNull(scip, &childcoefs); 2022 2023 return SCIP_OKAY; 2024 } 2025 2026 /** detects complex quadratic expressions that can be represented as SOC constraints 2027 * 2028 * These are quadratic expressions with either exactly one positive or exactly one negative eigenvalue, 2029 * in addition to some extra conditions. One needs to write the quadratic as 2030 * sum eigval_i (eigvec_i . x)^2 + c ≤ -eigval_k (eigvec_k . x)^2, where eigval_k is the negative eigenvalue, 2031 * and c must be positive and (eigvec_k . x) must not change sign. 2032 * This is described in more details in 2033 * Mahajan, Ashutosh & Munson, Todd, Exploiting Second-Order Cone Structure for Global Optimization, 2010. 2034 * 2035 * The eigen-decomposition is computed using Lapack. 2036 * Binary linear variables are interpreted as quadratic terms. 2037 * 2038 * @todo: In the case -b <= a + x^2 - y^2 <= b, it is possible to represent both sides by SOC. Currently, the 2039 * datastructure can only handle one SOC. If this should appear more often, it could be worth to extend it, 2040 * such that both sides can be handled (see e.g. instance chp_partload). 2041 * FS: this shouldn't be possible. For a <= b + x^2 - y^2 <= c to be SOC representable on both sides, we would need 2042 * that a - b >= 0 and b -c >= 0, but this implies that a >= c and assuming the constraint is not trivially infeasible, 2043 * a <= b. Thus, a = b = c and the constraint is x^2 == y^2. 2044 * 2045 * @todo: Since cons_nonlinear multiplies as many terms out as possible during presolving, some SOC-representable 2046 * structures cannot be detected, (see e.g. instances bearing or wager). There is currently no obvious way 2047 * to handle this. 2048 */ 2049 static 2050 SCIP_RETCODE detectSocQuadraticComplex( 2051 SCIP* scip, /**< SCIP data structure */ 2052 SCIP_EXPR* expr, /**< expression */ 2053 SCIP_Real conslhs, /**< lhs of the constraint that the expression defines (or SCIP_INVALID) */ 2054 SCIP_Real consrhs, /**< rhs of the constraint that the expression defines (or SCIP_INVALID) */ 2055 SCIP_NLHDLREXPRDATA** nlhdlrexprdata, /**< pointer to store nonlinear handler expression data */ 2056 SCIP_Bool* enforcebelow, /**< pointer to store whether we enforce <= (TRUE) or >= (FALSE); only 2057 * valid when success is TRUE */ 2058 SCIP_Bool* success /**< pointer to store whether SOC structure has been detected */ 2059 ) 2060 { 2061 SCIP_EXPR** occurringexprs; 2062 SCIP_HASHMAP* expr2idx; 2063 SCIP_Real* offsets; 2064 SCIP_Real* transcoefs; 2065 SCIP_Real* eigvecmatrix; 2066 SCIP_Real* eigvals; 2067 SCIP_Real* lincoefs; 2068 SCIP_Real* bp; 2069 int* transcoefsidx; 2070 int* termbegins; 2071 SCIP_Real constant; 2072 SCIP_Real lhsconstant; 2073 SCIP_Real lhs; 2074 SCIP_Real rhs; 2075 SCIP_INTERVAL expractivity; 2076 int nvars; 2077 int nterms; 2078 int nchildren; 2079 int npos; 2080 int nneg; 2081 int ntranscoefs; 2082 int i; 2083 int j; 2084 SCIP_Bool rhsissoc; 2085 SCIP_Bool lhsissoc; 2086 SCIP_Bool isquadratic; 2087 2088 assert(expr != NULL); 2089 assert(success != NULL); 2090 2091 *success = FALSE; 2092 2093 /* check whether expression is a sum with at least 2 children */ 2094 if( ! SCIPisExprSum(scip, expr) || SCIPexprGetNChildren(expr) < 2 ) 2095 { 2096 return SCIP_OKAY; 2097 } 2098 2099 /* we need Lapack to compute eigenvalues/vectors below */ 2100 if( ! SCIPlapackIsAvailable() ) 2101 return SCIP_OKAY; 2102 2103 /* get children of the sum */ 2104 nchildren = SCIPexprGetNChildren(expr); 2105 constant = SCIPgetConstantExprSum(expr); 2106 2107 /* initialize data */ 2108 offsets = NULL; 2109 transcoefs = NULL; 2110 transcoefsidx = NULL; 2111 termbegins = NULL; 2112 bp = NULL; 2113 2114 SCIP_CALL( SCIPhashmapCreate(&expr2idx, SCIPblkmem(scip), 2 * nchildren) ); 2115 SCIP_CALL( SCIPallocBufferArray(scip, &occurringexprs, 2 * nchildren) ); 2116 2117 /* check if the expression is quadratic and collect all occurring expressions */ 2118 SCIP_CALL( checkAndCollectQuadratic(scip, expr, expr2idx, occurringexprs, &nvars, &isquadratic) ); 2119 2120 if( !isquadratic ) 2121 { 2122 SCIPfreeBufferArray(scip, &occurringexprs); 2123 SCIPhashmapFree(&expr2idx); 2124 return SCIP_OKAY; 2125 } 2126 2127 /* check that nvars*nvars doesn't get too large, see also SCIPcomputeExprQuadraticCurvature() */ 2128 if( nvars > 7000 ) 2129 { 2130 SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, "nlhdlr_soc - number of quadratic variables is too large (%d) to check the curvature\n", nvars); 2131 SCIPfreeBufferArray(scip, &occurringexprs); 2132 SCIPhashmapFree(&expr2idx); 2133 return SCIP_OKAY; 2134 } 2135 2136 assert(SCIPhashmapGetNElements(expr2idx) == nvars); 2137 2138 /* create datastructures for constaint defining matrix and vector */ 2139 SCIP_CALL( SCIPallocClearBufferArray(scip, &eigvecmatrix, nvars * nvars) ); /*lint !e647*/ 2140 SCIP_CALL( SCIPallocClearBufferArray(scip, &lincoefs, nvars) ); 2141 2142 /* build constraint defining matrix (stored in eigvecmatrix) and vector (stored in lincoefs) */ 2143 buildQuadExprMatrix(scip, expr, expr2idx, nvars, eigvecmatrix, lincoefs); 2144 2145 SCIP_CALL( SCIPallocBufferArray(scip, &eigvals, nvars) ); 2146 2147 /* compute eigenvalues and vectors, A = PDP^t 2148 * note: eigvecmatrix stores P^t, i.e., P^t_{i,j} = eigvecmatrix[i*nvars+j] 2149 */ 2150 if( SCIPlapackComputeEigenvalues(SCIPbuffer(scip), TRUE, nvars, eigvecmatrix, eigvals) != SCIP_OKAY ) 2151 { 2152 SCIPdebugMsg(scip, "Failed to compute eigenvalues and eigenvectors for expression:\n"); 2153 2154 #ifdef SCIP_DEBUG 2155 SCIPdismantleExpr(scip, NULL, expr); 2156 #endif 2157 2158 goto CLEANUP; 2159 } 2160 2161 SCIP_CALL( SCIPallocClearBufferArray(scip, &bp, nvars) ); 2162 2163 nneg = 0; 2164 npos = 0; 2165 ntranscoefs = 0; 2166 2167 /* set small eigenvalues to 0 and compute b*P */ 2168 for( i = 0; i < nvars; ++i ) 2169 { 2170 for( j = 0; j < nvars; ++j ) 2171 { 2172 bp[i] += lincoefs[j] * eigvecmatrix[i * nvars + j]; 2173 2174 /* count the number of transcoefs to be used later */ 2175 if( !SCIPisZero(scip, eigvals[i]) && !SCIPisZero(scip, eigvecmatrix[i * nvars + j]) ) 2176 ++ntranscoefs; 2177 } 2178 2179 if( SCIPisZero(scip, eigvals[i]) ) 2180 { 2181 /* if there is a purely linear variable, the constraint can't be written as a SOC */ 2182 if( !SCIPisZero(scip, bp[i]) ) 2183 goto CLEANUP; 2184 2185 bp[i] = 0.0; 2186 eigvals[i] = 0.0; 2187 } 2188 else if( eigvals[i] > 0.0 ) 2189 npos++; 2190 else 2191 nneg++; 2192 } 2193 2194 /* a proper SOC constraint needs at least 2 variables */ 2195 if( npos + nneg < 2 ) 2196 goto CLEANUP; 2197 2198 /* determine whether rhs or lhs of cons is potentially SOC, if any */ 2199 rhsissoc = (nneg == 1 && SCIPgetExprNLocksPosNonlinear(expr) > 0); 2200 lhsissoc = (npos == 1 && SCIPgetExprNLocksNegNonlinear(expr) > 0); 2201 2202 if( rhsissoc || lhsissoc ) 2203 { 2204 if( conslhs == SCIP_INVALID || consrhs == SCIP_INVALID ) /*lint !e777*/ 2205 { 2206 SCIP_CALL( SCIPevalExprActivity(scip, expr) ); 2207 expractivity = SCIPexprGetActivity(expr); 2208 lhs = (conslhs == SCIP_INVALID ? expractivity.inf : conslhs); /*lint !e777*/ 2209 rhs = (consrhs == SCIP_INVALID ? expractivity.sup : consrhs); /*lint !e777*/ 2210 } 2211 else 2212 { 2213 lhs = conslhs; 2214 rhs = consrhs; 2215 } 2216 } 2217 else 2218 { 2219 /* if none of the sides is potentially SOC, stop */ 2220 goto CLEANUP; 2221 } 2222 2223 /* @TODO: what do we do if both sides are possible? */ 2224 if( !rhsissoc ) 2225 { 2226 assert(lhsissoc); 2227 2228 /* lhs is potentially SOC, change signs */ 2229 lhsconstant = lhs - constant; /*lint !e644*/ 2230 2231 for( i = 0; i < nvars; ++i ) 2232 { 2233 eigvals[i] = -eigvals[i]; 2234 bp[i] = -bp[i]; 2235 } 2236 *enforcebelow = FALSE; /* enforce lhs <= expr */ 2237 } 2238 else 2239 { 2240 lhsconstant = constant - rhs; /*lint !e644*/ 2241 *enforcebelow = TRUE; /* enforce expr <= rhs */ 2242 } 2243 2244 /* initialize remaining datastructures for nonlinear handler */ 2245 SCIP_CALL( SCIPallocBufferArray(scip, &offsets, npos + nneg + 1) ); 2246 SCIP_CALL( SCIPallocBufferArray(scip, &transcoefs, ntranscoefs) ); 2247 SCIP_CALL( SCIPallocBufferArray(scip, &transcoefsidx, ntranscoefs) ); 2248 SCIP_CALL( SCIPallocBufferArray(scip, &termbegins, npos + nneg + 2) ); 2249 2250 /* try to fill the nlhdlrexprdata (at this point, it can still fail) */ 2251 SCIP_CALL( tryFillNlhdlrExprDataQuad(scip, occurringexprs, eigvecmatrix, eigvals, bp, nvars, termbegins, transcoefs, 2252 transcoefsidx, offsets, &lhsconstant, &nterms, success) ); 2253 2254 if( !(*success) ) 2255 goto CLEANUP; 2256 2257 assert(0 < nterms && nterms <= npos + nneg + 1); 2258 assert(ntranscoefs == termbegins[nterms]); 2259 2260 /* 2261 * at this point, the expression passed all checks and is SOC-representable 2262 */ 2263 2264 /* register all requests for auxiliary variables */ 2265 for( i = 0; i < nvars; ++i ) 2266 { 2267 SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, occurringexprs[i], TRUE, FALSE, FALSE, FALSE) ); 2268 } 2269 2270 #ifdef SCIP_DEBUG 2271 SCIPdebugMsg(scip, "found SOC structure for expression %p\n%f <= ", (void*)expr, lhs); 2272 SCIPprintExpr(scip, expr, NULL); 2273 SCIPinfoMessage(scip, NULL, "<= %f\n", rhs); 2274 #endif 2275 2276 /* finally, create and store nonlinear handler expression data */ 2277 SCIP_CALL( createNlhdlrExprData(scip, occurringexprs, offsets, transcoefs, transcoefsidx, termbegins, nvars, nterms, 2278 nlhdlrexprdata) ); 2279 assert(*nlhdlrexprdata != NULL); 2280 2281 CLEANUP: 2282 SCIPfreeBufferArrayNull(scip, &termbegins); 2283 SCIPfreeBufferArrayNull(scip, &transcoefsidx); 2284 SCIPfreeBufferArrayNull(scip, &transcoefs); 2285 SCIPfreeBufferArrayNull(scip, &offsets); 2286 SCIPfreeBufferArrayNull(scip, &bp); 2287 SCIPfreeBufferArray(scip, &eigvals); 2288 SCIPfreeBufferArray(scip, &lincoefs); 2289 SCIPfreeBufferArray(scip, &eigvecmatrix); 2290 SCIPfreeBufferArray(scip, &occurringexprs); 2291 SCIPhashmapFree(&expr2idx); 2292 2293 return SCIP_OKAY; 2294 } 2295 2296 /** helper method to detect SOC structures 2297 * 2298 * The detection runs in 3 steps: 2299 * 1. check if expression is a norm of the form \f$\sqrt{\sum_i (\text{sqrcoef}_i\, \text{expr}_i^2 + \text{lincoef}_i\, \text{expr}_i) + c}\f$ 2300 * which can be transformed to the form \f$\sqrt{\sum_i (\text{coef}_i \text{expr}_i + \text{const}_i)^2 + c^*}\f$ with \f$c^* \geq 0\f$.\n 2301 * -> this results in the SOC expr ≤ auxvar(expr) 2302 * 2303 * TODO we should generalize and check for sqrt(positive-semidefinite-quadratic) 2304 * 2305 * 2. check if expression represents a quadratic function of one of the following forms (all coefs > 0) 2306 * 1. \f$(\sum_i \text{coef}_i \text{expr}_i^2) - \text{coef}_k \text{expr}_k^2 \leq \text{RHS}\f$ or 2307 * 2. \f$(\sum_i - \text{coef}_i \text{expr}_i^2) + \text{coef}_k \text{expr}_k^2 \geq \text{LHS}\f$ or 2308 * 3. \f$(\sum_i \text{coef}_i \text{expr}_i^2) - \text{coef}_k \text{expr}_k \text{expr}_l \leq \text{RHS}\f$ or 2309 * 4. \f$(\sum_i - \text{coef}_i \text{expr}_i^2) + \text{coef}_k \text{expr}_k \text{expr}_l \geq \text{LHS}\f$, 2310 * 2311 * where RHS ≥ 0 or LHS ≤ 0, respectively. For LHS and RHS we use the constraint sides if it is a root expr 2312 * and the bounds of the auxiliary variable otherwise. 2313 * The last two cases are called hyperbolic or rotated second order cone.\n 2314 * -> this results in the SOC \f$\sqrt{(\sum_i \text{coef}_i \text{expr}_i^2) - \text{RHS}} \leq \sqrt{\text{coef}_k} \text{expr}_k\f$ 2315 * or \f$\sqrt{4(\sum_i \text{coef}_i \text{expr}_i^2) - 4\text{RHS} + (\text{expr}_k - \text{expr}_l)^2)} \leq \text{expr}_k + \text{expr}_l\f$. 2316 * (analogously for the LHS cases) 2317 * 2318 * 3. check if expression represents a quadratic inequality of the form \f$f(x) = x^TAx + b^Tx + c \leq 0\f$ such that \f$f(x)\f$ 2319 * has exactly one negative eigenvalue plus some extra conditions, see detectSocQuadraticComplex(). 2320 * 2321 * Note that step 3 is only performed if parameter `compeigenvalues` is set to TRUE. 2322 */ 2323 static 2324 SCIP_RETCODE detectSOC( 2325 SCIP* scip, /**< SCIP data structure */ 2326 SCIP_NLHDLRDATA* nlhdlrdata, /**< nonlinear handler data */ 2327 SCIP_EXPR* expr, /**< expression */ 2328 SCIP_Real conslhs, /**< lhs of the constraint that the expression defines (or SCIP_INVALID) */ 2329 SCIP_Real consrhs, /**< rhs of the constraint that the expression defines (or SCIP_INVALID) */ 2330 SCIP_NLHDLREXPRDATA** nlhdlrexprdata, /**< pointer to store nonlinear handler expression data */ 2331 SCIP_Bool* enforcebelow, /**< pointer to store whether we enforce <= (TRUE) or >= (FALSE); only 2332 * valid when success is TRUE */ 2333 SCIP_Bool* success /**< pointer to store whether SOC structure has been detected */ 2334 ) 2335 { 2336 assert(expr != NULL); 2337 assert(nlhdlrdata != NULL); 2338 assert(nlhdlrexprdata != NULL); 2339 assert(success != NULL); 2340 2341 *success = FALSE; 2342 2343 /* check whether expression is given as norm as described in case 1 above: if we have a constraint 2344 * sqrt(sum x_i^2) <= constant, then it might be better not to handle this here; thus, we only call detectSocNorm 2345 * when the expr is _not_ the root of a constraint 2346 */ 2347 if( conslhs == SCIP_INVALID && consrhs == SCIP_INVALID ) /*lint !e777*/ 2348 { 2349 SCIP_CALL( detectSocNorm(scip, expr, nlhdlrexprdata, success) ); 2350 *enforcebelow = *success; 2351 } 2352 2353 if( !(*success) ) 2354 { 2355 /* check whether expression is a simple soc-respresentable quadratic expression as described in case 2 above */ 2356 SCIP_CALL( detectSocQuadraticSimple(scip, expr, conslhs, consrhs, nlhdlrexprdata, enforcebelow, success) ); 2357 } 2358 2359 if( !(*success) && nlhdlrdata->compeigenvalues ) 2360 { 2361 /* check whether expression is a more complex soc-respresentable quadratic expression as described in case 3 */ 2362 SCIP_CALL( detectSocQuadraticComplex(scip, expr, conslhs, consrhs, nlhdlrexprdata, enforcebelow, success) ); 2363 } 2364 2365 return SCIP_OKAY; 2366 } 2367 2368 /* 2369 * Callback methods of nonlinear handler 2370 */ 2371 2372 /** nonlinear handler copy callback */ 2373 static 2374 SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrSoc) 2375 { /*lint --e{715}*/ 2376 assert(targetscip != NULL); 2377 assert(sourcenlhdlr != NULL); 2378 assert(strcmp(SCIPnlhdlrGetName(sourcenlhdlr), NLHDLR_NAME) == 0); 2379 2380 SCIP_CALL( SCIPincludeNlhdlrSoc(targetscip) ); 2381 2382 return SCIP_OKAY; 2383 } 2384 2385 /** callback to free data of handler */ 2386 static 2387 SCIP_DECL_NLHDLRFREEHDLRDATA(nlhdlrFreehdlrdataSoc) 2388 { /*lint --e{715}*/ 2389 assert(nlhdlrdata != NULL); 2390 2391 SCIPfreeBlockMemory(scip, nlhdlrdata); 2392 2393 return SCIP_OKAY; 2394 } 2395 2396 /** callback to free expression specific data */ 2397 static 2398 SCIP_DECL_NLHDLRFREEEXPRDATA(nlhdlrFreeExprDataSoc) 2399 { /*lint --e{715}*/ 2400 assert(*nlhdlrexprdata != NULL); 2401 2402 SCIP_CALL( freeNlhdlrExprData(scip, nlhdlrexprdata) ); 2403 2404 return SCIP_OKAY; 2405 } 2406 2407 2408 /** callback to be called in initialization */ 2409 #if 0 2410 static 2411 SCIP_DECL_NLHDLRINIT(nlhdlrInitSoc) 2412 { /*lint --e{715}*/ 2413 SCIPerrorMessage("method of soc nonlinear handler not implemented yet\n"); 2414 SCIPABORT(); /*lint --e{527}*/ 2415 2416 return SCIP_OKAY; 2417 } 2418 #else 2419 #define nlhdlrInitSoc NULL 2420 #endif 2421 2422 2423 /** callback to be called in deinitialization */ 2424 #if 0 2425 static 2426 SCIP_DECL_NLHDLREXIT(nlhdlrExitSoc) 2427 { /*lint --e{715}*/ 2428 SCIPerrorMessage("method of soc nonlinear handler not implemented yet\n"); 2429 SCIPABORT(); /*lint --e{527}*/ 2430 2431 return SCIP_OKAY; 2432 } 2433 #else 2434 #define nlhdlrExitSoc NULL 2435 #endif 2436 2437 2438 /** callback to detect structure in expression tree */ 2439 static 2440 SCIP_DECL_NLHDLRDETECT(nlhdlrDetectSoc) 2441 { /*lint --e{715}*/ 2442 SCIP_Real conslhs; 2443 SCIP_Real consrhs; 2444 SCIP_Bool enforcebelow; 2445 SCIP_Bool success; 2446 SCIP_NLHDLRDATA* nlhdlrdata; 2447 2448 assert(expr != NULL); 2449 2450 /* don't try if no sepa is required 2451 * TODO implement some bound strengthening 2452 */ 2453 if( (*enforcing & SCIP_NLHDLR_METHOD_SEPABOTH) == SCIP_NLHDLR_METHOD_SEPABOTH ) 2454 return SCIP_OKAY; 2455 2456 assert(SCIPgetExprNAuxvarUsesNonlinear(expr) > 0); /* since some sepa is required, there should have been demand for it */ 2457 2458 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr); 2459 assert(nlhdlrdata != NULL); 2460 2461 conslhs = (cons == NULL ? SCIP_INVALID : SCIPgetLhsNonlinear(cons)); 2462 consrhs = (cons == NULL ? SCIP_INVALID : SCIPgetRhsNonlinear(cons)); 2463 2464 SCIP_CALL( detectSOC(scip, nlhdlrdata, expr, conslhs, consrhs, nlhdlrexprdata, &enforcebelow, &success) ); 2465 2466 if( !success ) 2467 return SCIP_OKAY; 2468 2469 /* inform what we can do */ 2470 *participating = enforcebelow ? SCIP_NLHDLR_METHOD_SEPABELOW : SCIP_NLHDLR_METHOD_SEPAABOVE; 2471 2472 /* if we have been successful on sqrt(...) <= auxvar, then we enforce 2473 * otherwise, expr is quadratic and we separate for expr <= ub(auxvar) only 2474 * in that case, we enforce only if expr is the root of a constraint, since then replacing auxvar by up(auxvar) does not relax anything (auxvar <= ub(auxvar) is the only constraint on auxvar) 2475 */ 2476 if( (SCIPisExprPower(scip, expr) && SCIPgetExponentExprPow(expr) == 0.5) || (cons != NULL) ) 2477 *enforcing |= *participating; 2478 2479 return SCIP_OKAY; 2480 } 2481 2482 2483 /** auxiliary evaluation callback of nonlinear handler 2484 * @todo: remember if we are in the original variables and avoid reevaluating 2485 */ 2486 static 2487 SCIP_DECL_NLHDLREVALAUX(nlhdlrEvalauxSoc) 2488 { /*lint --e{715}*/ 2489 int i; 2490 2491 assert(nlhdlrexprdata != NULL); 2492 assert(nlhdlrexprdata->vars != NULL); 2493 assert(nlhdlrexprdata->transcoefs != NULL); 2494 assert(nlhdlrexprdata->transcoefsidx != NULL); 2495 assert(nlhdlrexprdata->nterms > 1); 2496 2497 /* if the original expression is a norm, evaluate w.r.t. the auxiliary variables */ 2498 if( SCIPisExprPower(scip, expr) ) 2499 { 2500 assert(SCIPgetExponentExprPow(expr) == 0.5); 2501 2502 updateVarVals(scip, nlhdlrexprdata, sol, FALSE); 2503 2504 /* compute sum_i coef_i expr_i^2 */ 2505 *auxvalue = 0.0; 2506 for( i = 0; i < nlhdlrexprdata->nterms - 1; ++i ) 2507 { 2508 SCIP_Real termval; 2509 2510 termval = evalSingleTerm(scip, nlhdlrexprdata, i); 2511 *auxvalue += SQR(termval); 2512 } 2513 2514 assert(*auxvalue >= 0.0); 2515 2516 /* compute SQRT(sum_i coef_i expr_i^2) */ 2517 *auxvalue = SQRT(*auxvalue); 2518 } 2519 /* otherwise, evaluate the original quadratic expression w.r.t. the created auxvars of the children */ 2520 else 2521 { 2522 SCIP_EXPR** children; 2523 SCIP_Real* childcoefs; 2524 int nchildren; 2525 2526 assert(SCIPisExprSum(scip, expr)); 2527 2528 children = SCIPexprGetChildren(expr); 2529 childcoefs = SCIPgetCoefsExprSum(expr); 2530 nchildren = SCIPexprGetNChildren(expr); 2531 2532 *auxvalue = SCIPgetConstantExprSum(expr); 2533 2534 for( i = 0; i < nchildren; ++i ) 2535 { 2536 if( SCIPisExprPower(scip, children[i]) ) 2537 { 2538 SCIP_VAR* argauxvar; 2539 SCIP_Real solval; 2540 2541 assert(SCIPgetExponentExprPow(children[i]) == 2.0); 2542 2543 argauxvar = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(children[i])[0]); 2544 assert(argauxvar != NULL); 2545 2546 solval = SCIPgetSolVal(scip, sol, argauxvar); 2547 *auxvalue += childcoefs[i] * SQR( solval ); 2548 } 2549 else if( SCIPisExprProduct(scip, children[i]) ) 2550 { 2551 SCIP_VAR* argauxvar1; 2552 SCIP_VAR* argauxvar2; 2553 2554 assert(SCIPexprGetNChildren(children[i]) == 2); 2555 2556 argauxvar1 = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(children[i])[0]); 2557 argauxvar2 = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(children[i])[1]); 2558 assert(argauxvar1 != NULL); 2559 assert(argauxvar2 != NULL); 2560 2561 *auxvalue += childcoefs[i] * SCIPgetSolVal(scip, sol, argauxvar1) * SCIPgetSolVal(scip, sol, argauxvar2); 2562 } 2563 else 2564 { 2565 SCIP_VAR* argauxvar; 2566 2567 argauxvar = SCIPgetExprAuxVarNonlinear(children[i]); 2568 assert(argauxvar != NULL); 2569 2570 *auxvalue += childcoefs[i] * SCIPgetSolVal(scip, sol, argauxvar); 2571 } 2572 } 2573 } 2574 2575 return SCIP_OKAY; 2576 } 2577 2578 2579 /** separation deinitialization method of a nonlinear handler (called during CONSINITLP) */ 2580 static 2581 SCIP_DECL_NLHDLRINITSEPA(nlhdlrInitSepaSoc) 2582 { /*lint --e{715}*/ 2583 SCIP_ROWPREP* rowprep; 2584 2585 assert(conshdlr != NULL); 2586 assert(expr != NULL); 2587 assert(nlhdlrexprdata != NULL); 2588 2589 /* already needed for debug solution */ 2590 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &nlhdlrexprdata->varvals, nlhdlrexprdata->nvars) ); 2591 2592 /* if we have 3 or more terms in lhs create variable and row for disaggregation */ 2593 if( nlhdlrexprdata->nterms > 3 ) 2594 { 2595 /* create variables for cone disaggregation */ 2596 SCIP_CALL( createDisaggrVars(scip, expr, nlhdlrexprdata) ); 2597 2598 #ifdef WITH_DEBUG_SOLUTION 2599 if( SCIPdebugIsMainscip(scip) ) 2600 { 2601 SCIP_Real lhsval; 2602 SCIP_Real rhsval; 2603 SCIP_Real disvarval; 2604 int ndisvars; 2605 int nterms; 2606 int i; 2607 2608 for( i = 0; i < nlhdlrexprdata->nvars; ++i ) 2609 { 2610 SCIP_CALL( SCIPdebugGetSolVal(scip, SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[i]), &nlhdlrexprdata->varvals[i]) ); 2611 } 2612 2613 /* the debug solution value of the disaggregation variables is set to 2614 * (v_i^T x + beta_i)^2 / (v_{n+1}^T x + beta_{n+1}) 2615 * if (v_{n+1}^T x + beta_{n+1}) is different from 0. 2616 * Otherwise, the debug solution value is set to 0. 2617 */ 2618 2619 nterms = nlhdlrexprdata->nterms; 2620 2621 /* find value of rhs */ 2622 rhsval = evalSingleTerm(scip, nlhdlrexprdata, nterms - 1); 2623 2624 /* set value of disaggregation vars */ 2625 ndisvars = nlhdlrexprdata->nterms - 1; 2626 2627 if( SCIPisZero(scip, rhsval) ) 2628 { 2629 for( i = 0; i < ndisvars; ++i ) 2630 { 2631 SCIP_CALL( SCIPdebugAddSolVal(scip, nlhdlrexprdata->disvars[i], 0.0) ); 2632 } 2633 } 2634 else 2635 { 2636 /* set value for each disaggregation variable corresponding to quadratic term */ 2637 for( i = 0; i < ndisvars; ++i ) 2638 { 2639 lhsval = evalSingleTerm(scip, nlhdlrexprdata, i); 2640 2641 disvarval = SQR(lhsval) / rhsval; 2642 2643 SCIP_CALL( SCIPdebugAddSolVal(scip, nlhdlrexprdata->disvars[i], disvarval) ); 2644 } 2645 } 2646 } 2647 #endif 2648 2649 /* create the disaggregation row and store it in nlhdlrexprdata */ 2650 SCIP_CALL( createDisaggrRow(scip, conshdlr, expr, nlhdlrexprdata) ); 2651 } 2652 2653 #ifdef SCIP_DEBUG 2654 SCIPdebugMsg(scip, "initlp for \n"); 2655 printNlhdlrExprData(scip, nlhdlrexprdata); 2656 #endif 2657 2658 /* add some initial cuts on well-selected coordinates */ 2659 if( nlhdlrexprdata->nterms == 2 ) 2660 { 2661 /* we have |v_1^T x + \beta_1| \leq v_2^T x + \beta_2 2662 * 2663 * we should linearize at points where the first term is -1 or 1, so we can take 2664 * 2665 * x = v_1 / ||v_1||^2 (+/-1 - beta_1) 2666 */ 2667 SCIP_Real plusminus1; 2668 SCIP_Real norm; 2669 int i; 2670 2671 /* calculate ||v_1||^2 */ 2672 norm = 0.0; 2673 for( i = nlhdlrexprdata->termbegins[0]; i < nlhdlrexprdata->termbegins[1]; ++i ) 2674 norm += SQR(nlhdlrexprdata->transcoefs[i]); 2675 assert(norm > 0.0); 2676 2677 BMSclearMemoryArray(nlhdlrexprdata->varvals, nlhdlrexprdata->nvars); 2678 2679 for( plusminus1 = -1.0; plusminus1 <= 1.0; plusminus1 += 2.0 ) 2680 { 2681 /* set x = v_1 / ||v_1||^2 (plusminus1 - beta_1) */ 2682 for( i = nlhdlrexprdata->termbegins[0]; i < nlhdlrexprdata->termbegins[1]; ++i ) 2683 nlhdlrexprdata->varvals[nlhdlrexprdata->transcoefsidx[i]] = nlhdlrexprdata->transcoefs[i] / norm * (plusminus1 - nlhdlrexprdata->offsets[0]); 2684 assert(SCIPisEQ(scip, evalSingleTerm(scip, nlhdlrexprdata, 0), plusminus1)); 2685 2686 /* compute gradient cut */ 2687 SCIP_CALL( generateCutSolSOC(scip, &rowprep, expr, cons, nlhdlrexprdata, -SCIPinfinity(scip), evalSingleTerm(scip, nlhdlrexprdata, 1)) ); 2688 2689 if( rowprep != NULL ) 2690 { 2691 SCIP_Bool success = FALSE; 2692 2693 SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, NULL, SCIPgetHugeValue(scip), &success) ); 2694 if( success ) 2695 { 2696 SCIP_ROW* cut; 2697 SCIP_CALL( SCIPgetRowprepRowCons(scip, &cut, rowprep, cons) ); 2698 SCIP_CALL( SCIPaddRow(scip, cut, FALSE, infeasible) ); 2699 SCIP_CALL( SCIPreleaseRow(scip, &cut) ); 2700 } 2701 2702 SCIPfreeRowprep(scip, &rowprep); 2703 } 2704 2705 if( *infeasible ) 2706 return SCIP_OKAY; 2707 } 2708 } 2709 else if( nlhdlrexprdata->nterms == 3 && nlhdlrexprdata->termbegins[0] != nlhdlrexprdata->termbegins[1] ) 2710 { 2711 /* we have \sqrt{ (v_1^T x + \beta_1)^2 + (v_2^T x + \beta_2)^2 } \leq v_3^T x + \beta_3 2712 * with v_1 != 0 2713 * 2714 * we should linearize at points where the first and second term are (-1,0), (1,1), (1,-1), i.e., 2715 * 2716 * v_1^T x + \beta_1 = -1 1 1 2717 * v_2^T x + \beta_2 = 0 1 -1 2718 * 2719 * Let i be the index of the first nonzero element in v_1. 2720 * Let j != i be an index of v_2, to be determined. 2721 * Assume all other entries of x will be set to 0. 2722 * Then we have 2723 * (v_1)_i x_i + (v_1)_j x_j = c with c = -1 - beta_1 2724 * (v_2)_i x_i + (v_2)_j x_j = d with d = 0 - beta_2 2725 * 2726 * Since (v_1)_i != 0, this gives 2727 * x_i = 1/(v_1)_i (c - (v_1)_j x_j) 2728 * Substituting in the 2nd equation, we get 2729 * (v_2)_i/(v_1)_i (c - (v_1)_j x_j) + (v_2)_j x_j = d 2730 * -> ((v_2)_j - (v_2)_i (v_1)_j / (v_1)_i) x_j = d - (v_2)_i/(v_1)_i c 2731 * Now find j such that (v_2)_j - (v_2)_i (v_1)_j / (v_1)_i != 0. 2732 * 2733 * If v_2 = 0, then linearize only for first term being -1 or 1 and don't care about value of second term. 2734 * We then set j arbitrary, x_i = 1/(v_1)_i c, other coordinates of x = zero. 2735 */ 2736 static const SCIP_Real refpoints[3][2] = { {-1.0, 0.0}, {1.0, 1.0}, {1.0, -1.0} }; 2737 SCIP_Real v1i, v1j = 0.0; 2738 SCIP_Real v2i, v2j = 0.0; 2739 SCIP_Bool v2zero; 2740 int i; 2741 int j = -1; 2742 int k; 2743 int pos; 2744 2745 i = nlhdlrexprdata->transcoefsidx[nlhdlrexprdata->termbegins[0]]; 2746 v1i = nlhdlrexprdata->transcoefs[nlhdlrexprdata->termbegins[0]]; 2747 assert(v1i != 0.0); 2748 2749 v2zero = nlhdlrexprdata->termbegins[1] == nlhdlrexprdata->termbegins[2]; 2750 2751 /* get (v_2)_i; as it is a sparse vector, we need to search for i in transcoefsidx */ 2752 v2i = 0.0; 2753 if( !v2zero && SCIPsortedvecFindInt(nlhdlrexprdata->transcoefsidx + nlhdlrexprdata->termbegins[1], i, nlhdlrexprdata->termbegins[2] - nlhdlrexprdata->termbegins[1], &pos) ) 2754 { 2755 assert(nlhdlrexprdata->transcoefsidx[nlhdlrexprdata->termbegins[1] + pos] == i); 2756 v2i = nlhdlrexprdata->transcoefs[nlhdlrexprdata->termbegins[1] + pos]; 2757 } 2758 2759 /* find a j, for now look only into indices with (v_2)_j != 0 */ 2760 for( k = nlhdlrexprdata->termbegins[1]; k < nlhdlrexprdata->termbegins[2]; ++k ) 2761 { 2762 /* check whether transcoefsidx[k] could be a good j */ 2763 2764 if( nlhdlrexprdata->transcoefsidx[k] == i ) /* i == j is not good */ 2765 continue; 2766 2767 /* get (v_1)_j; as it is a sparse vector, we need to search for j in transcoefsidx */ 2768 v1j = 0.0; 2769 if( SCIPsortedvecFindInt(nlhdlrexprdata->transcoefsidx + nlhdlrexprdata->termbegins[0], nlhdlrexprdata->transcoefsidx[k], nlhdlrexprdata->termbegins[1] - nlhdlrexprdata->termbegins[0], &pos) ) 2770 { 2771 assert(nlhdlrexprdata->transcoefsidx[nlhdlrexprdata->termbegins[0] + pos] == nlhdlrexprdata->transcoefsidx[k]); 2772 v1j = nlhdlrexprdata->transcoefs[nlhdlrexprdata->termbegins[0] + pos]; 2773 } 2774 2775 v2j = nlhdlrexprdata->transcoefs[k]; 2776 2777 if( SCIPisZero(scip, v2j - v2i * v1j / v1i) ) /* (v_2)_j - (v_2)_i (v_1)_j / (v_1)_i = 0 is also not good */ 2778 continue; 2779 2780 j = nlhdlrexprdata->transcoefsidx[k]; 2781 break; 2782 } 2783 2784 if( v2zero ) 2785 { 2786 j = 0; 2787 v1j = 0.0; 2788 v2j = 0.0; 2789 } 2790 2791 if( j != -1 ) 2792 { 2793 SCIP_Real c, d; 2794 int point; 2795 2796 BMSclearMemoryArray(nlhdlrexprdata->varvals, nlhdlrexprdata->nvars); 2797 2798 for( point = 0; point < (v2zero ? 2 : 3); ++point ) 2799 { 2800 c = refpoints[point][0] - nlhdlrexprdata->offsets[0]; 2801 2802 if( !v2zero ) 2803 { 2804 /* set x_j and x_i */ 2805 d = refpoints[point][1] - nlhdlrexprdata->offsets[1]; 2806 nlhdlrexprdata->varvals[j] = (d - v2i/v1i*c) / (v2j - v2i * v1j / v1i); 2807 nlhdlrexprdata->varvals[i] = (c - v1j * nlhdlrexprdata->varvals[j]) / v1i; 2808 2809 SCIPdebugMsg(scip, "<%s>(%d) = %g, <%s>(%d) = %g\n", 2810 SCIPvarGetName(SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[i])), i, nlhdlrexprdata->varvals[i], 2811 SCIPvarGetName(SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[j])), j, nlhdlrexprdata->varvals[j]); 2812 } 2813 else 2814 { 2815 /* set x_i */ 2816 nlhdlrexprdata->varvals[i] = c / v1i; 2817 2818 SCIPdebugMsg(scip, "<%s>(%d) = %g\n", 2819 SCIPvarGetName(SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[i])), i, nlhdlrexprdata->varvals[i]); 2820 } 2821 2822 assert(SCIPisEQ(scip, evalSingleTerm(scip, nlhdlrexprdata, 0), refpoints[point][0])); 2823 assert(v2zero || SCIPisEQ(scip, evalSingleTerm(scip, nlhdlrexprdata, 1), refpoints[point][1])); 2824 2825 /* compute gradient cut */ 2826 SCIP_CALL( generateCutSolSOC(scip, &rowprep, expr, cons, nlhdlrexprdata, -SCIPinfinity(scip), evalSingleTerm(scip, nlhdlrexprdata, 2)) ); 2827 2828 if( rowprep != NULL ) 2829 { 2830 SCIP_Bool success = FALSE; 2831 2832 SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, NULL, SCIPgetHugeValue(scip), &success) ); 2833 if( success ) 2834 { 2835 SCIP_ROW* cut; 2836 SCIP_CALL( SCIPgetRowprepRowCons(scip, &cut, rowprep, cons) ); 2837 SCIP_CALL( SCIPaddRow(scip, cut, FALSE, infeasible) ); 2838 SCIP_CALL( SCIPreleaseRow(scip, &cut) ); 2839 } 2840 2841 SCIPfreeRowprep(scip, &rowprep); 2842 } 2843 2844 if( *infeasible ) 2845 return SCIP_OKAY; 2846 } 2847 } 2848 } 2849 else if( nlhdlrexprdata->nterms == 3 ) 2850 { 2851 /* we have \sqrt{ \beta_1^2 + (v_2^T x + \beta_2)^2 } \leq v_3^T x + \beta_3 2852 * with v_2 != 0 2853 * 2854 * we should linearize at points where the second term is -1 or 1 2855 * 2856 * set x = v_2 / ||v_2||^2 (+/-1 - beta_2) 2857 */ 2858 SCIP_Real plusminus1; 2859 SCIP_Real norm; 2860 int i; 2861 2862 /* calculate ||v_2||^2 */ 2863 norm = 0.0; 2864 for( i = nlhdlrexprdata->termbegins[1]; i < nlhdlrexprdata->termbegins[2]; ++i ) 2865 norm += SQR(nlhdlrexprdata->transcoefs[i]); 2866 assert(norm > 0.0); 2867 2868 BMSclearMemoryArray(nlhdlrexprdata->varvals, nlhdlrexprdata->nvars); 2869 2870 for( plusminus1 = -1.0; plusminus1 <= 1.0; plusminus1 += 2.0 ) 2871 { 2872 /* set x = v_2 / ||v_2||^2 (plusminus1 - beta_2) */ 2873 for( i = nlhdlrexprdata->termbegins[1]; i < nlhdlrexprdata->termbegins[2]; ++i ) 2874 nlhdlrexprdata->varvals[nlhdlrexprdata->transcoefsidx[i]] = nlhdlrexprdata->transcoefs[i] / norm * (plusminus1 - nlhdlrexprdata->offsets[1]); 2875 assert(SCIPisEQ(scip, evalSingleTerm(scip, nlhdlrexprdata, 1), plusminus1)); 2876 2877 /* compute gradient cut */ 2878 SCIP_CALL( generateCutSolSOC(scip, &rowprep, expr, cons, nlhdlrexprdata, -SCIPinfinity(scip), evalSingleTerm(scip, nlhdlrexprdata, 2)) ); 2879 2880 if( rowprep != NULL ) 2881 { 2882 SCIP_Bool success = FALSE; 2883 2884 SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, NULL, SCIPgetHugeValue(scip), &success) ); 2885 if( success ) 2886 { 2887 SCIP_ROW* cut; 2888 SCIP_CALL( SCIPgetRowprepRowCons(scip, &cut, rowprep, cons) ); 2889 SCIP_CALL( SCIPaddRow(scip, cut, FALSE, infeasible) ); 2890 SCIP_CALL( SCIPreleaseRow(scip, &cut) ); 2891 } 2892 2893 SCIPfreeRowprep(scip, &rowprep); 2894 } 2895 2896 if( *infeasible ) 2897 return SCIP_OKAY; 2898 } 2899 } 2900 else 2901 { 2902 /* generate gradient cuts for the small rotated cones 2903 * \f[ 2904 * \sqrt{4(v_k^T x + \beta_k)^2 + (v_n^T x + \beta_n - y_k)^2} - v_n^T x - \beta_n - y_k \leq 0. 2905 * \f] 2906 * 2907 * we should linearize again at points where the first and second term (inside sqr) are (-1/2,0), (1/2,1), (1/2,-1). 2908 * Since we have y_k, we can achieve this more easily here via 2909 * x = v_k/||v_k||^2 (+/-0.5 - beta_k) 2910 * y_k = v_n^T x + beta_n + 0/1/-1 2911 * 2912 * If v_k = 0, then we use x = 0 and linearize for second term being 1 and -1 only 2913 */ 2914 static const SCIP_Real refpoints[3][2] = { {-0.5, 0.0}, {0.5, 1.0}, {0.5, -1.0} }; 2915 SCIP_Real rhsval; 2916 SCIP_Real norm; 2917 int point; 2918 int i; 2919 int k; 2920 SCIP_Bool vkzero; 2921 2922 /* add disaggregation row to LP */ 2923 SCIP_CALL( SCIPaddRow(scip, nlhdlrexprdata->disrow, FALSE, infeasible) ); 2924 2925 if( *infeasible ) 2926 return SCIP_OKAY; 2927 2928 for( k = 0; k < nlhdlrexprdata->nterms - 1; ++k ) 2929 { 2930 vkzero = nlhdlrexprdata->termbegins[k+1] == nlhdlrexprdata->termbegins[k]; 2931 assert(!vkzero || nlhdlrexprdata->offsets[k] != 0.0); 2932 2933 /* calculate ||v_k||^2 */ 2934 norm = 0.0; 2935 for( i = nlhdlrexprdata->termbegins[k]; i < nlhdlrexprdata->termbegins[k+1]; ++i ) 2936 norm += SQR(nlhdlrexprdata->transcoefs[i]); 2937 assert(vkzero || norm > 0.0); 2938 2939 BMSclearMemoryArray(nlhdlrexprdata->varvals, nlhdlrexprdata->nvars); 2940 2941 for( point = vkzero ? 1 : 0; point < 3; ++point ) 2942 { 2943 /* set x = v_k / ||v_k||^2 (refpoints[point][0] - beta_k) / 2 */ 2944 for( i = nlhdlrexprdata->termbegins[k]; i < nlhdlrexprdata->termbegins[k+1]; ++i ) 2945 nlhdlrexprdata->varvals[nlhdlrexprdata->transcoefsidx[i]] = nlhdlrexprdata->transcoefs[i] / norm * (refpoints[point][0] - nlhdlrexprdata->offsets[k]); /*lint !e795*/ 2946 assert(vkzero || SCIPisEQ(scip, evalSingleTerm(scip, nlhdlrexprdata, k), refpoints[point][0])); 2947 2948 /* set y_k = v_n^T x + beta_n + 0/1/-1 */ 2949 rhsval = evalSingleTerm(scip, nlhdlrexprdata, nlhdlrexprdata->nterms - 1); 2950 nlhdlrexprdata->disvarvals[k] = rhsval + refpoints[point][1]; 2951 2952 /* compute gradient cut */ 2953 SCIP_CALL( generateCutSolDisagg(scip, &rowprep, expr, cons, nlhdlrexprdata, k, -SCIPinfinity(scip), rhsval) ); 2954 2955 if( rowprep != NULL ) 2956 { 2957 SCIP_Bool success = FALSE; 2958 2959 SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, NULL, SCIPgetHugeValue(scip), &success) ); 2960 if( success ) 2961 { 2962 SCIP_ROW* cut; 2963 SCIP_CALL( SCIPgetRowprepRowCons(scip, &cut, rowprep, cons) ); 2964 SCIP_CALL( SCIPaddRow(scip, cut, FALSE, infeasible) ); 2965 SCIP_CALL( SCIPreleaseRow(scip, &cut) ); 2966 } 2967 2968 SCIPfreeRowprep(scip, &rowprep); 2969 } 2970 2971 if( *infeasible ) 2972 return SCIP_OKAY; 2973 } 2974 } 2975 } 2976 2977 return SCIP_OKAY; 2978 } 2979 2980 2981 /** separation deinitialization method of a nonlinear handler (called during CONSEXITSOL) */ 2982 static 2983 SCIP_DECL_NLHDLREXITSEPA(nlhdlrExitSepaSoc) 2984 { /*lint --e{715}*/ 2985 assert(nlhdlrexprdata != NULL); 2986 2987 /* free disaggreagation row */ 2988 if( nlhdlrexprdata->disrow != NULL ) 2989 { 2990 SCIP_CALL( SCIPreleaseRow(scip, &nlhdlrexprdata->disrow) ); 2991 } 2992 2993 SCIPfreeBlockMemoryArray(scip, &nlhdlrexprdata->varvals, nlhdlrexprdata->nvars); 2994 2995 return SCIP_OKAY; 2996 } 2997 2998 2999 /** nonlinear handler separation callback */ 3000 static 3001 SCIP_DECL_NLHDLRENFO(nlhdlrEnfoSoc) 3002 { /*lint --e{715}*/ 3003 SCIP_NLHDLRDATA* nlhdlrdata; 3004 SCIP_Real rhsval; 3005 int ndisaggrs; 3006 int k; 3007 SCIP_Bool infeasible; 3008 3009 assert(nlhdlrexprdata != NULL); 3010 assert(nlhdlrexprdata->nterms < 4 || nlhdlrexprdata->disrow != NULL); 3011 assert(nlhdlrexprdata->nterms > 1); 3012 3013 *result = SCIP_DIDNOTFIND; 3014 3015 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr); 3016 assert(nlhdlrdata != NULL); 3017 3018 /* update varvals 3019 * set variables close to integer to integer, in particular when close to zero 3020 * for simple soc's (no large v_i, no offsets), variables close to zero would give coefficients close to zero in the cut, 3021 * which the cut cleanup may have problems to relax (and we end up with local or much relaxed cuts) 3022 * also when close to other integers, rounding now may prevent some relaxation in cut cleanup 3023 */ 3024 updateVarVals(scip, nlhdlrexprdata, sol, TRUE); 3025 3026 rhsval = evalSingleTerm(scip, nlhdlrexprdata, nlhdlrexprdata->nterms - 1); 3027 3028 /* if there are three or two terms just compute gradient cut */ 3029 if( nlhdlrexprdata->nterms < 4 ) 3030 { 3031 SCIP_ROWPREP* rowprep; 3032 3033 /* compute gradient cut */ 3034 SCIP_CALL( generateCutSolSOC(scip, &rowprep, expr, cons, nlhdlrexprdata, SCIPgetLPFeastol(scip), rhsval) ); 3035 3036 if( rowprep != NULL ) 3037 { 3038 SCIP_CALL( addCut(scip, nlhdlrdata, rowprep, sol, cons, allowweakcuts, result) ); 3039 3040 SCIPfreeRowprep(scip, &rowprep); 3041 } 3042 else 3043 { 3044 SCIPdebugMsg(scip, "failed to generate cut for SOC\n"); 3045 } 3046 3047 return SCIP_OKAY; 3048 } 3049 3050 ndisaggrs = nlhdlrexprdata->nterms - 1; 3051 3052 /* check whether the aggregation row is in the LP */ 3053 if( !SCIProwIsInLP(nlhdlrexprdata->disrow) && -SCIPgetRowSolFeasibility(scip, nlhdlrexprdata->disrow, sol) > SCIPgetLPFeastol(scip) ) 3054 { 3055 SCIP_CALL( SCIPaddRow(scip, nlhdlrexprdata->disrow, TRUE, &infeasible) ); 3056 SCIPdebugMsg(scip, "added disaggregation row to LP, cutoff=%u\n", infeasible); 3057 3058 if( infeasible ) 3059 { 3060 *result = SCIP_CUTOFF; 3061 return SCIP_OKAY; 3062 } 3063 3064 *result = SCIP_SEPARATED; 3065 } 3066 3067 for( k = 0; k < ndisaggrs && *result != SCIP_CUTOFF; ++k ) 3068 { 3069 SCIP_ROWPREP* rowprep; 3070 3071 /* compute gradient cut */ 3072 SCIP_CALL( generateCutSolDisagg(scip, &rowprep, expr, cons, nlhdlrexprdata, k, SCIPgetLPFeastol(scip), rhsval) ); 3073 3074 if( rowprep != NULL ) 3075 { 3076 SCIP_CALL( addCut(scip, nlhdlrdata, rowprep, sol, cons, allowweakcuts, result) ); 3077 3078 SCIPfreeRowprep(scip, &rowprep); 3079 } 3080 } 3081 3082 return SCIP_OKAY; 3083 } 3084 3085 static 3086 SCIP_DECL_NLHDLRSOLLINEARIZE(nlhdlrSollinearizeSoc) 3087 { /*lint --e{715}*/ 3088 SCIP_NLHDLRDATA* nlhdlrdata; 3089 SCIP_Real rhsval; 3090 int k; 3091 3092 assert(sol != NULL); 3093 assert(nlhdlrexprdata != NULL); 3094 assert(nlhdlrexprdata->nterms < 4 || nlhdlrexprdata->disrow != NULL); 3095 assert(nlhdlrexprdata->nterms > 1); 3096 3097 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr); 3098 assert(nlhdlrdata != NULL); 3099 3100 /* update varvals */ 3101 updateVarVals(scip, nlhdlrexprdata, sol, TRUE); 3102 3103 rhsval = evalSingleTerm(scip, nlhdlrexprdata, nlhdlrexprdata->nterms - 1); 3104 3105 /* if there are three or two terms just compute gradient cut */ 3106 if( nlhdlrexprdata->nterms < 4 ) 3107 { 3108 SCIP_ROWPREP* rowprep; 3109 3110 /* compute gradient cut */ 3111 SCIP_CALL( generateCutSolSOC(scip, &rowprep, expr, cons, nlhdlrexprdata, -SCIPinfinity(scip), rhsval) ); 3112 3113 if( rowprep != NULL ) 3114 { 3115 SCIP_CALL( addCutPool(scip, nlhdlrdata, rowprep, sol, cons) ); 3116 3117 SCIPfreeRowprep(scip, &rowprep); 3118 } 3119 3120 return SCIP_OKAY; 3121 } 3122 3123 for( k = 0; k < nlhdlrexprdata->nterms - 1; ++k ) 3124 { 3125 SCIP_ROWPREP* rowprep; 3126 3127 /* compute gradient cut */ 3128 SCIP_CALL( generateCutSolDisagg(scip, &rowprep, expr, cons, nlhdlrexprdata, k, -SCIPinfinity(scip), rhsval) ); 3129 3130 if( rowprep != NULL ) 3131 { 3132 SCIP_CALL( addCutPool(scip, nlhdlrdata, rowprep, sol, cons) ); 3133 3134 SCIPfreeRowprep(scip, &rowprep); 3135 } 3136 } 3137 3138 return SCIP_OKAY; 3139 } 3140 3141 /* 3142 * nonlinear handler specific interface methods 3143 */ 3144 3145 /** includes SOC nonlinear handler in nonlinear constraint handler */ 3146 SCIP_RETCODE SCIPincludeNlhdlrSoc( 3147 SCIP* scip /**< SCIP data structure */ 3148 ) 3149 { 3150 SCIP_NLHDLRDATA* nlhdlrdata; 3151 SCIP_NLHDLR* nlhdlr; 3152 3153 assert(scip != NULL); 3154 3155 /* create nonlinear handler data */ 3156 SCIP_CALL( SCIPallocClearBlockMemory(scip, &nlhdlrdata) ); 3157 3158 SCIP_CALL( SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, NLHDLR_NAME, NLHDLR_DESC, NLHDLR_DETECTPRIORITY, NLHDLR_ENFOPRIORITY, nlhdlrDetectSoc, nlhdlrEvalauxSoc, nlhdlrdata) ); 3159 assert(nlhdlr != NULL); 3160 3161 SCIPnlhdlrSetCopyHdlr(nlhdlr, nlhdlrCopyhdlrSoc); 3162 SCIPnlhdlrSetFreeHdlrData(nlhdlr, nlhdlrFreehdlrdataSoc); 3163 SCIPnlhdlrSetFreeExprData(nlhdlr, nlhdlrFreeExprDataSoc); 3164 SCIPnlhdlrSetInitExit(nlhdlr, nlhdlrInitSoc, nlhdlrExitSoc); 3165 SCIPnlhdlrSetSepa(nlhdlr, nlhdlrInitSepaSoc, nlhdlrEnfoSoc, NULL, nlhdlrExitSepaSoc); 3166 SCIPnlhdlrSetSollinearize(nlhdlr, nlhdlrSollinearizeSoc); 3167 3168 /* add soc nlhdlr parameters */ 3169 /* TODO should we get rid of this and use separating/mineffiacy(root) instead, which is 1e-4? */ 3170 SCIP_CALL( SCIPaddRealParam(scip, "nlhdlr/" NLHDLR_NAME "/mincutefficacy", 3171 "Minimum efficacy which a cut needs in order to be added.", 3172 &nlhdlrdata->mincutefficacy, FALSE, DEFAULT_MINCUTEFFICACY, 0.0, SCIPinfinity(scip), NULL, NULL) ); 3173 3174 SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/compeigenvalues", 3175 "Should Eigenvalue computations be done to detect complex cases in quadratic constraints?", 3176 &nlhdlrdata->compeigenvalues, FALSE, DEFAULT_COMPEIGENVALUES, NULL, NULL) ); 3177 3178 return SCIP_OKAY; 3179 } 3180 3181 /** checks whether constraint is SOC representable in original variables and returns the SOC representation 3182 * 3183 * The SOC representation has the form: 3184 * \f$\sqrt{\sum_{i=1}^{n} (v_i^T x + \beta_i)^2} - v_{n+1}^T x - \beta_{n+1} \lessgtr 0\f$, 3185 * where \f$n+1 = \text{nterms}\f$ and the inequality type is given by sidetype (`SCIP_SIDETYPE_RIGHT` if inequality 3186 * is \f$\leq\f$, `SCIP_SIDETYPE_LEFT` if \f$\geq\f$). 3187 * 3188 * For each term (i.e. for each \f$i\f$ in the above notation as well as \f$n+1\f$), the constant \f$\beta_i\f$ is given by the 3189 * corresponding element `offsets[i-1]` and `termbegins[i-1]` is the starting position of the term in arrays 3190 * `transcoefs` and `transcoefsidx`. The overall number of nonzeros is `termbegins[nterms]`. 3191 * 3192 * Arrays `transcoefs` and `transcoefsidx` have size `termbegins[nterms]` and define the linear expressions \f$v_i^T x\f$ 3193 * for each term. For a term \f$i\f$ in the above notation, the nonzeroes are given by elements 3194 * `termbegins[i-1]...termbegins[i]` of `transcoefs` and `transcoefsidx`. There may be no nonzeroes for some term (i.e., 3195 * constant terms are possible). `transcoefs` contains the coefficients \f$v_i\f$ and `transcoefsidx` contains positions of 3196 * variables in the `vars` array. 3197 * 3198 * The `vars` array has size `nvars` and contains \f$x\f$ variables; each variable is included at most once. 3199 * 3200 * The arrays should be freed by calling SCIPfreeSOCArraysNonlinear(). 3201 * 3202 * This function uses the methods that are used in the detection algorithm of the SOC nonlinear handler. 3203 */ 3204 SCIP_RETCODE SCIPisSOCNonlinear( 3205 SCIP* scip, /**< SCIP data structure */ 3206 SCIP_CONS* cons, /**< nonlinear constraint */ 3207 SCIP_Bool compeigenvalues, /**< whether eigenvalues should be computed to detect complex cases */ 3208 SCIP_Bool* success, /**< pointer to store whether SOC structure has been detected */ 3209 SCIP_SIDETYPE* sidetype, /**< pointer to store which side of cons is SOC representable; only 3210 * valid when success is TRUE */ 3211 SCIP_VAR*** vars, /**< variables (x) that appear on both sides; no duplicates are allowed */ 3212 SCIP_Real** offsets, /**< offsets of both sides (beta_i) */ 3213 SCIP_Real** transcoefs, /**< non-zeros of linear transformation vectors (v_i) */ 3214 int** transcoefsidx, /**< mapping of transformation coefficients to variable indices in vars */ 3215 int** termbegins, /**< starting indices of transcoefs for each term */ 3216 int* nvars, /**< total number of variables appearing (i.e. size of vars) */ 3217 int* nterms /**< number of summands in the SQRT +1 for RHS (n+1) */ 3218 ) 3219 { 3220 SCIP_NLHDLRDATA nlhdlrdata; 3221 SCIP_NLHDLREXPRDATA *nlhdlrexprdata; 3222 SCIP_Real conslhs; 3223 SCIP_Real consrhs; 3224 SCIP_EXPR* expr; 3225 SCIP_Bool enforcebelow; 3226 int i; 3227 3228 assert(cons != NULL); 3229 3230 expr = SCIPgetExprNonlinear(cons); 3231 assert(expr != NULL); 3232 3233 nlhdlrdata.mincutefficacy = 0.0; 3234 nlhdlrdata.compeigenvalues = compeigenvalues; 3235 3236 conslhs = SCIPgetLhsNonlinear(cons); 3237 consrhs = SCIPgetRhsNonlinear(cons); 3238 3239 SCIP_CALL( detectSOC(scip, &nlhdlrdata, expr, conslhs, consrhs, &nlhdlrexprdata, &enforcebelow, success) ); 3240 3241 /* the constraint must be SOC representable in original variables */ 3242 if( *success ) 3243 { 3244 assert(nlhdlrexprdata != NULL); 3245 3246 for( i = 0; i < nlhdlrexprdata->nvars; ++i ) 3247 { 3248 if( !SCIPisExprVar(scip, nlhdlrexprdata->vars[i]) ) 3249 { 3250 *success = FALSE; 3251 break; 3252 } 3253 } 3254 } 3255 3256 if( *success ) 3257 { 3258 *sidetype = enforcebelow ? SCIP_SIDETYPE_RIGHT : SCIP_SIDETYPE_LEFT; 3259 SCIP_CALL( SCIPallocBlockMemoryArray(scip, vars, nlhdlrexprdata->nvars) ); 3260 3261 for( i = 0; i < nlhdlrexprdata->nvars; ++i ) 3262 { 3263 (*vars)[i] = SCIPgetVarExprVar(nlhdlrexprdata->vars[i]); 3264 assert((*vars)[i] != NULL); 3265 } 3266 SCIPfreeBlockMemoryArray(scip, &nlhdlrexprdata->vars, nlhdlrexprdata->nvars); 3267 *offsets = nlhdlrexprdata->offsets; 3268 *transcoefs = nlhdlrexprdata->transcoefs; 3269 *transcoefsidx = nlhdlrexprdata->transcoefsidx; 3270 *termbegins = nlhdlrexprdata->termbegins; 3271 *nvars = nlhdlrexprdata->nvars; 3272 *nterms = nlhdlrexprdata->nterms; 3273 SCIPfreeBlockMemory(scip, &nlhdlrexprdata); 3274 } 3275 else 3276 { 3277 if( nlhdlrexprdata != NULL ) 3278 { 3279 SCIP_CALL( freeNlhdlrExprData(scip, &nlhdlrexprdata) ); 3280 } 3281 *vars = NULL; 3282 *offsets = NULL; 3283 *transcoefs = NULL; 3284 *transcoefsidx = NULL; 3285 *termbegins = NULL; 3286 *nvars = 0; 3287 *nterms = 0; 3288 } 3289 3290 return SCIP_OKAY; 3291 } 3292 3293 /** frees arrays created by SCIPisSOCNonlinear() */ 3294 void SCIPfreeSOCArraysNonlinear( 3295 SCIP* scip, /**< SCIP data structure */ 3296 SCIP_VAR*** vars, /**< variables that appear on both sides (x) */ 3297 SCIP_Real** offsets, /**< offsets of both sides (beta_i) */ 3298 SCIP_Real** transcoefs, /**< non-zeros of linear transformation vectors (v_i) */ 3299 int** transcoefsidx, /**< mapping of transformation coefficients to variable indices in vars */ 3300 int** termbegins, /**< starting indices of transcoefs for each term */ 3301 int nvars, /**< total number of variables appearing */ 3302 int nterms /**< number of summands in the SQRT +1 for RHS (n+1) */ 3303 ) 3304 { 3305 int ntranscoefs; 3306 3307 if( nvars == 0 ) 3308 return; 3309 3310 assert(vars != NULL); 3311 assert(offsets != NULL); 3312 assert(transcoefs != NULL); 3313 assert(transcoefsidx != NULL); 3314 assert(termbegins != NULL); 3315 3316 ntranscoefs = (*termbegins)[nterms]; 3317 3318 SCIPfreeBlockMemoryArray(scip, termbegins, nterms + 1); 3319 SCIPfreeBlockMemoryArray(scip, transcoefsidx, ntranscoefs); 3320 SCIPfreeBlockMemoryArray(scip, transcoefs, ntranscoefs); 3321 SCIPfreeBlockMemoryArray(scip, offsets, nterms); 3322 SCIPfreeBlockMemoryArray(scip, vars, nvars); 3323 } 3324