SCIP Doxygen Documentation
 
Loading...
Searching...
No Matches
nlhdlr_soc.c
Go to the documentation of this file.
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/nlpi_ipopt.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 */
108struct 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
123struct SCIP_NlhdlrData
124{
125 SCIP_Real mincutefficacy; /**< minimum efficacy a cut need to be added */
126 SCIP_Bool compeigenvalues; /**< whether Eigenvalue computations should be done to detect complex cases */
127};
128
129/*
130 * Local methods
131 */
132
133#ifdef SCIP_DEBUG
134/** prints the nlhdlr expression data */
135static
137 SCIP* scip, /**< SCIP data structure */
138 SCIP_NLHDLREXPRDATA* nlhdlrexprdata /**< pointer to store nonlinear handler expression data */
139 )
140{
141 int nterms;
142 int i;
143 int j;
144
145 nterms = nlhdlrexprdata->nterms;
146
147 SCIPinfoMessage(scip, NULL, "SQRT( ");
148
149 for( i = 0; i < nterms - 1; ++i )
150 {
151 int startidx;
152
153 startidx = nlhdlrexprdata->termbegins[i];
154
155 /* v_i is 0 */
156 if( startidx == nlhdlrexprdata->termbegins[i + 1] )
157 {
158 assert(nlhdlrexprdata->offsets[i] != 0.0);
159
160 SCIPinfoMessage(scip, NULL, "%f", SQR(nlhdlrexprdata->offsets[i]));
161 continue;
162 }
163
164 /* v_i is not 0 */
166
167 for( j = startidx; j < nlhdlrexprdata->termbegins[i + 1]; ++j )
168 {
169 if( nlhdlrexprdata->transcoefs[j] != 1.0 )
170 SCIPinfoMessage(scip, NULL, "%f*", nlhdlrexprdata->transcoefs[j]);
171 if( SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]]) != NULL )
172 {
173 SCIPinfoMessage(scip, NULL, "%s", SCIPvarGetName(SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]])));
174 SCIPinfoMessage(scip, NULL, "(%p)", (void*)nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]]);
175 }
176 else
177 SCIPinfoMessage(scip, NULL, "%p", (void*)nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]]);
178
179 if( j < nlhdlrexprdata->termbegins[i + 1] - 1 )
180 SCIPinfoMessage(scip, NULL, " + ");
181 else if( nlhdlrexprdata->offsets[i] != 0.0 )
182 SCIPinfoMessage(scip, NULL, " + %f", nlhdlrexprdata->offsets[i]);
183 }
184
185 SCIPinfoMessage(scip, NULL, ")^2");
186
187 if( i < nterms - 2 )
188 SCIPinfoMessage(scip, NULL, " + ");
189 }
190
191 SCIPinfoMessage(scip, NULL, " ) <= ");
192
193 for( j = nlhdlrexprdata->termbegins[nterms-1]; j < nlhdlrexprdata->termbegins[nterms]; ++j )
194 {
195 if( nlhdlrexprdata->transcoefs[j] != 1.0 )
196 SCIPinfoMessage(scip, NULL, "%f*", nlhdlrexprdata->transcoefs[j]);
197 if( SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]]) != NULL )
198 SCIPinfoMessage(scip, NULL, "%s", SCIPvarGetName(SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]])));
199 else
200 SCIPinfoMessage(scip, NULL, "%p", (void*)nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[j]]);
201
202 if( j < nlhdlrexprdata->termbegins[nterms] - 1 )
203 SCIPinfoMessage(scip, NULL, " + ");
204 else if( nlhdlrexprdata->offsets[nterms-1] != 0.0 )
205 SCIPinfoMessage(scip, NULL, " + %f", nlhdlrexprdata->offsets[nterms-1]);
206 }
207
208 SCIPinfoMessage(scip, NULL, "\n");
209}
210#endif
211
212/** helper method to create variables for the cone disaggregation */
213static
215 SCIP* scip, /**< SCIP data structure */
216 SCIP_EXPR* expr, /**< expression */
217 SCIP_NLHDLREXPRDATA* nlhdlrexprdata /**< nonlinear handler expression data */
218 )
219{
220 char name[SCIP_MAXSTRLEN];
221 int ndisvars;
222 int i;
223
224 assert(nlhdlrexprdata != NULL);
225
226 ndisvars = nlhdlrexprdata->nterms - 1;
227
228 /* allocate memory */
229 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &nlhdlrexprdata->disvars, ndisvars) );
230
231 /* create disaggregation variables representing the epigraph of (v_i^T x + beta_i)^2 / (v_{n+1}^T x + beta_{n+1}) */
232 for( i = 0; i < ndisvars; ++i )
233 {
234 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedis_%p_%d", (void*) expr, i);
235 SCIP_CALL( SCIPcreateVarBasic(scip, &nlhdlrexprdata->disvars[i], name, 0.0, SCIPinfinity(scip), 0.0,
237 SCIPvarMarkRelaxationOnly(nlhdlrexprdata->disvars[i]);
238
239 SCIP_CALL( SCIPaddVar(scip, nlhdlrexprdata->disvars[i]) );
240 SCIP_CALL( SCIPaddVarLocksType(scip, nlhdlrexprdata->disvars[i], SCIP_LOCKTYPE_MODEL, 1, 1) );
241 }
242
243 return SCIP_OKAY;
244}
245
246/** helper method to free variables for the cone disaggregation */
247static
249 SCIP* scip, /**< SCIP data structure */
250 SCIP_NLHDLREXPRDATA* nlhdlrexprdata /**< nonlinear handler expression data */
251 )
252{
253 int ndisvars;
254 int i;
255
256 assert(nlhdlrexprdata != NULL);
257
258 if( nlhdlrexprdata->disvars == NULL )
259 return SCIP_OKAY;
260
261 ndisvars = nlhdlrexprdata->nterms - 1;
262
263 /* release variables */
264 for( i = 0; i < ndisvars; ++i )
265 {
266 SCIP_CALL( SCIPaddVarLocksType(scip, nlhdlrexprdata->disvars[i], SCIP_LOCKTYPE_MODEL, -1, -1) );
267 SCIP_CALL( SCIPreleaseVar(scip, &nlhdlrexprdata->disvars[i]) );
268 }
269
270 /* free memory */
271 SCIPfreeBlockMemoryArrayNull(scip, &nlhdlrexprdata->disvars, ndisvars);
272
273 return SCIP_OKAY;
274}
275
276/** helper method to create the disaggregation row \f$\text{disvars}_i \leq v_{n+1}^T x + \beta_{n+1}\f$ */
277static
279 SCIP* scip, /**< SCIP data structure */
280 SCIP_CONSHDLR* conshdlr, /**< nonlinear constraint handler */
281 SCIP_EXPR* expr, /**< expression */
282 SCIP_NLHDLREXPRDATA* nlhdlrexprdata /**< nonlinear handler expression data */
283 )
284{
285 SCIP_Real beta;
286 char name[SCIP_MAXSTRLEN];
287 int ndisvars;
288 int nterms;
289 int i;
290
291 assert(scip != NULL);
292 assert(expr != NULL);
293 assert(nlhdlrexprdata != NULL);
294 assert(nlhdlrexprdata->disrow == NULL);
295
296 nterms = nlhdlrexprdata->nterms;
297 beta = nlhdlrexprdata->offsets[nterms - 1];
298
299 ndisvars = nterms - 1;
300
301 /* create row 0 <= beta_{n+1} */
302 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedis_%p_row", (void*) expr);
303 SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &nlhdlrexprdata->disrow, conshdlr, name,
304 -SCIPinfinity(scip), beta, FALSE, FALSE, TRUE) );
305
306 /* add disvars to row */
307 for( i = 0; i < ndisvars; ++i )
308 {
309 SCIP_CALL( SCIPaddVarToRow(scip, nlhdlrexprdata->disrow, nlhdlrexprdata->disvars[i], 1.0) );
310 }
311
312 /* add rhs vars to row */
313 for( i = nlhdlrexprdata->termbegins[nterms - 1]; i < nlhdlrexprdata->termbegins[nterms]; ++i )
314 {
315 SCIP_VAR* var;
316 SCIP_Real coef;
317
318 var = SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[i]]);
319 assert(var != NULL);
320
321 coef = -nlhdlrexprdata->transcoefs[i];
322
323 SCIP_CALL( SCIPaddVarToRow(scip, nlhdlrexprdata->disrow, var, coef) );
324 }
325
326 return SCIP_OKAY;
327}
328
329/** helper method to create nonlinear handler expression data */
330static
332 SCIP* scip, /**< SCIP data structure */
333 SCIP_EXPR** vars, /**< expressions which variables appear on both sides (\f$x\f$) */
334 SCIP_Real* offsets, /**< offsets of bot sides (\f$beta_i\f$) */
335 SCIP_Real* transcoefs, /**< non-zeroes of linear transformation vectors (\f$v_i\f$) */
336 int* transcoefsidx, /**< mapping of transformation coefficients to variable indices in vars */
337 int* termbegins, /**< starting indices of transcoefs for each term */
338 int nvars, /**< total number of variables appearing */
339 int nterms, /**< number of summands in the SQRT, +1 for RHS */
340 SCIP_NLHDLREXPRDATA** nlhdlrexprdata /**< pointer to store nonlinear handler expression data */
341 )
342{
343 int ntranscoefs;
344
345 assert(vars != NULL);
346 assert(offsets != NULL);
347 assert(transcoefs != NULL);
348 assert(transcoefsidx != NULL);
349 assert(termbegins != NULL);
350 assert(nlhdlrexprdata != NULL);
351
352 ntranscoefs = termbegins[nterms];
353
354 SCIP_CALL( SCIPallocBlockMemory(scip, nlhdlrexprdata) );
355 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*nlhdlrexprdata)->vars, vars, nvars) );
356 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*nlhdlrexprdata)->offsets, offsets, nterms) );
357 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*nlhdlrexprdata)->transcoefs, transcoefs, ntranscoefs) );
358 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*nlhdlrexprdata)->transcoefsidx, transcoefsidx, ntranscoefs) );
359 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*nlhdlrexprdata)->termbegins, termbegins, nterms + 1) );
360 (*nlhdlrexprdata)->nvars = nvars;
361 (*nlhdlrexprdata)->nterms = nterms;
362
363 (*nlhdlrexprdata)->disrow = NULL;
364 (*nlhdlrexprdata)->disvars = NULL;
365
366#ifdef SCIP_DEBUG
367 SCIPdebugMsg(scip, "created nlhdlr data for the following soc expression:\n");
368 printNlhdlrExprData(scip, *nlhdlrexprdata);
369 printf("x is %p\n", (void *)vars[0]);
370#endif
371
372 return SCIP_OKAY;
373}
374
375/** helper method to free nonlinear handler expression data */
376static
378 SCIP* scip, /**< SCIP data structure */
379 SCIP_NLHDLREXPRDATA** nlhdlrexprdata /**< pointer to free nonlinear handler expression data */
380 )
381{
382 int ntranscoefs;
383
384 assert(nlhdlrexprdata != NULL);
385 assert(*nlhdlrexprdata != NULL);
386
387 /* free variables and row for cone disaggregation */
388 SCIP_CALL( freeDisaggrVars(scip, *nlhdlrexprdata) );
389
390 ntranscoefs = (*nlhdlrexprdata)->termbegins[(*nlhdlrexprdata)->nterms];
391
392 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->termbegins, (*nlhdlrexprdata)->nterms + 1);
393 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->transcoefsidx, ntranscoefs);
394 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->transcoefs, ntranscoefs);
395 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->offsets, (*nlhdlrexprdata)->nterms);
396 SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->vars, (*nlhdlrexprdata)->nvars);
397 SCIPfreeBlockMemory(scip, nlhdlrexprdata);
398
399 return SCIP_OKAY;
400}
401
402/** evaluate a single term of the form \f$v_i^T x + \beta_i\f$ */
403static
405 SCIP* scip, /**< SCIP data structure */
406 SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nonlinear handler expression data */
407 SCIP_SOL* sol, /**< solution */
408 int k /**< term to be evaluated */
409 )
410{
411 SCIP_Real result;
412 int i;
413
414 assert(scip != NULL);
415 assert(nlhdlrexprdata != NULL);
417
418 result = nlhdlrexprdata->offsets[k];
419
420 for( i = nlhdlrexprdata->termbegins[k]; i < nlhdlrexprdata->termbegins[k + 1]; ++i )
421 {
422 SCIP_Real varval = SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[i]]));
423 result += nlhdlrexprdata->transcoefs[i] * varval;
424 }
425
426 return result;
427}
428
429/** computes gradient cut for a 2D or 3D SOC
430 *
431 * A 3D SOC looks like
432 * \f[
433 * \sqrt{ (v_1^T x + \beta_1)^2 + (v_2^T x + \beta_2)^2 } \leq v_3^T x + \beta_3
434 * \f]
435 *
436 * Let \f$f(x)\f$ be the left-hand-side. The partial derivatives of \f$f\f$ are given by
437 * \f[
438 * \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)}
439 * \f]
440 *
441 * and the gradient cut is then \f$f(x^*) + \nabla f(x^*)(x - x^*) \leq v_3^T x + \beta_3\f$.
442 *
443 * A 2D SOC is
444 * \f[
445 * |v_1^T x + \beta_1| \leq v_2^T x + \beta_2
446 * \f]
447 * but we build the cut using the same procedure as for 3D.
448 */
449static
451 SCIP* scip, /**< SCIP data structure */
452 SCIP_EXPR* expr, /**< expression */
453 SCIP_CONS* cons, /**< the constraint that expr is part of */
454 SCIP_SOL* sol, /**< solution to separate or NULL for the LP solution */
455 SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nonlinear handler expression data */
456 SCIP_Real mincutviolation, /**< minimal required cut violation */
457 SCIP_Real rhsval, /**< value of last term at sol */
458 SCIP_ROW** cut /**< pointer to store a cut */
459 )
460{
462 SCIP_Real* transcoefs;
463 SCIP_Real cutcoef;
464 SCIP_Real fvalue;
465 SCIP_Real valterms[2] = {0.0, 0.0}; /* for lint */
466 SCIP_Real cutrhs;
467 SCIP_EXPR** vars;
469 int* transcoefsidx;
470 int* termbegins;
471 int nterms;
472 int i;
473 int j;
474
475 assert(expr != NULL);
476 assert(cons != NULL);
477 assert(nlhdlrexprdata != NULL);
478 assert(mincutviolation >= 0.0);
479 assert(cut != NULL);
480
481 vars = nlhdlrexprdata->vars;
482 transcoefs = nlhdlrexprdata->transcoefs;
483 transcoefsidx = nlhdlrexprdata->transcoefsidx;
484 termbegins = nlhdlrexprdata->termbegins;
485 nterms = nlhdlrexprdata->nterms;
486
487 *cut = NULL;
488
489 /* evaluate lhs terms and compute f(x*) */
490 fvalue = 0.0;
491 for( i = 0; i < nterms - 1; ++i )
492 {
493 valterms[i] = evalSingleTerm(scip, nlhdlrexprdata, sol, i);
494 fvalue += SQR( valterms[i] );
495 }
496 fvalue = SQRT( fvalue );
497
498 /* don't generate cut if we are not violated @todo: remove this once core detects better when a nlhdlr's cons is
499 * violated
500 */
501 if( fvalue - rhsval <= mincutviolation )
502 {
503 SCIPdebugMsg(scip, "do not generate cut: rhsval %g, fvalue %g violation is %g\n", rhsval, fvalue, fvalue - rhsval);
504 return SCIP_OKAY;
505 }
506
507 /* if f(x*) = 0 then SOC can't be violated and we shouldn't be here */
508 assert(fvalue > 0.0);
509
510 /* create cut */
513
514 /* cut is f(x*) + \nabla f(x*)^T (x - x*) \leq v_n^T x + \beta_n, i.e.,
515 * \nabla f(x*)^T x - v_n^T x \leq \beta_n + \nabla f(x*)^T x* - f(x*)
516 * thus cutrhs is \beta_n - f(x*) + \nabla f(x*)^T x*
517 */
518 cutrhs = nlhdlrexprdata->offsets[nterms - 1] - fvalue;
519
520 /* add cut coefficients from lhs terms and compute cut's rhs */
521 for( j = 0; j < nterms - 1; ++j )
522 {
523 for( i = termbegins[j]; i < termbegins[j + 1]; ++i )
524 {
525 cutvar = SCIPgetExprAuxVarNonlinear(vars[transcoefsidx[i]]);
526
527 /* cutcoef is (the first part of) the partial derivative w.r.t cutvar */
528 cutcoef = transcoefs[i] * valterms[j] / fvalue;
529
531
533 }
534 }
535
536 /* add terms for v_n */
537 for( i = termbegins[nterms - 1]; i < termbegins[nterms]; ++i )
538 {
539 cutvar = SCIPgetExprAuxVarNonlinear(vars[transcoefsidx[i]]);
540 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, cutvar, -transcoefs[i]) );
541 }
542
543 /* add side */
545
547
548 if( SCIPgetRowprepViolation(scip, rowprep, sol, NULL) >= mincutviolation )
549 {
552 }
553 else
554 {
555 SCIPdebugMsg(scip, "%d-SOC rowprep violation %g below mincutviolation %g\n", nterms, SCIPgetRowprepViolation(scip,
556 rowprep, sol, NULL), mincutviolation);
557 /* SCIPprintRowprep(scip, rowprep, NULL); */
558 }
559
560 /* free memory */
562
563 return SCIP_OKAY;
564}
565
566/** helper method to compute and add a gradient cut for the k-th cone disaggregation
567 *
568 * 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$
569 * 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
570 * \f[
571 * (v_i^T x + \beta_i)^2 \leq (v_n^T x + \beta_n) y_i \text{ for } i \in \{0, \ldots, n -1\},
572 * \f]
573 * we want to separate one of the small rotated cones.
574 * We first transform it into standard form:
575 * \f[
576 * \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.
577 * \f]
578 * Let \f$f(x,y)\f$ be the left-hand-side. We now compute the gradient by
579 * \f{align*}{
580 * \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 \\
581 * \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
582 * \f}
583 * and the gradient cut is then \f$f(x^*, y^*) + \nabla f(x^*,y^*)((x,y) - (x^*, y^*)) \leq 0\f$.
584 */
585static
587 SCIP* scip, /**< SCIP data structure */
588 SCIP_EXPR* expr, /**< expression */
589 SCIP_CONS* cons, /**< the constraint that expr is part of */
590 SCIP_SOL* sol, /**< solution to separate or NULL for the LP solution */
591 SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nonlinear handler expression data */
592 int disaggidx, /**< index of disaggregation to separate */
593 SCIP_Real mincutviolation, /**< minimal required cut violation */
594 SCIP_Real rhsval, /**< value of the rhs term */
595 SCIP_ROW** cut /**< pointer to store a cut */
596 )
597{
598 SCIP_EXPR** vars;
599 SCIP_VAR** disvars;
600 SCIP_Real* transcoefs;
601 int* transcoefsidx;
602 int* termbegins;
605 SCIP_Real cutcoef;
606 SCIP_Real fvalue;
607 SCIP_Real disvarval;
608 SCIP_Real lhsval;
609 SCIP_Real constant;
610 SCIP_Real denominator;
611 int ncutvars;
612 int nterms;
613 int i;
614
615 assert(expr != NULL);
616 assert(cons != NULL);
617 assert(nlhdlrexprdata != NULL);
619 assert(mincutviolation >= 0.0);
620 assert(cut != NULL);
621
622 vars = nlhdlrexprdata->vars;
623 disvars = nlhdlrexprdata->disvars;
624 transcoefs = nlhdlrexprdata->transcoefs;
625 transcoefsidx = nlhdlrexprdata->transcoefsidx;
626 termbegins = nlhdlrexprdata->termbegins;
627 nterms = nlhdlrexprdata->nterms;
628
629 /* nterms is equal to n in the description and disaggidx is in {0, ..., n - 1} */
630
631 *cut = NULL;
632
634
635 lhsval = evalSingleTerm(scip, nlhdlrexprdata, sol, disaggidx);
636
638
639 /* compute value of function to be separated (f(x*,y*)) */
641
642 /* if the disagg soc is not violated don't compute cut */
643 if( fvalue <= mincutviolation )
644 {
645 SCIPdebugMsg(scip, "skip cut on disaggregation index %d as violation=%g below minviolation %g\n", disaggidx,
646 fvalue, mincutviolation);
647 return SCIP_OKAY;
648 }
649
650 /* if the denominator is 0 -> the constraint can't be violated */
652 /* if v_disaggidx^T x + beta_disaggidx is 0 -> the constraint can't be violated */
654
655 /* compute upper bound on the number of variables in cut: vars in rhs + vars in term + disagg var */
656 ncutvars = (termbegins[nterms] - termbegins[nterms-1]) + (termbegins[disaggidx + 1] - termbegins[disaggidx]) + 1;
657
658 /* create cut */
661
662 /* constant will be grad_f(x*,y*)^T (x*, y*) */
663 constant = 0.0;
664
665 /* a variable could appear on the lhs and rhs, but we add the coefficients separately */
666
667 /* add terms for v_disaggidx */
668 for( i = termbegins[disaggidx]; i < termbegins[disaggidx + 1]; ++i )
669 {
670 cutvar = SCIPgetExprAuxVarNonlinear(vars[transcoefsidx[i]]);
671 assert(cutvar != NULL);
672
673 /* cutcoef is (the first part of) the partial derivative w.r.t cutvar */
674 cutcoef = 4.0 * lhsval * transcoefs[i] / denominator;
675
677
678 constant += cutcoef * SCIPgetSolVal(scip, sol, cutvar);
679 }
680
681 /* add terms for v_n */
682 for( i = termbegins[nterms - 1]; i < termbegins[nterms]; ++i )
683 {
684 cutvar = SCIPgetExprAuxVarNonlinear(vars[transcoefsidx[i]]);
685 assert(cutvar != NULL);
686
687 /* cutcoef is the (second part of) the partial derivative w.r.t cutvar */
688 cutcoef = (rhsval - disvarval) * transcoefs[i] / denominator - transcoefs[i];
689
691
692 constant += cutcoef * SCIPgetSolVal(scip, sol, cutvar);
693 }
694
695 /* add term for disvar: cutcoef is the the partial derivative w.r.t. the disaggregation variable */
696 cutcoef = (disvarval - rhsval) / denominator - 1.0;
697 cutvar = disvars[disaggidx];
698
700
701 constant += cutcoef * SCIPgetSolVal(scip, sol, cutvar);
702
703 /* add side */
704 SCIProwprepAddSide(rowprep, constant - fvalue);
705
707
708 if( SCIPgetRowprepViolation(scip, rowprep, sol, NULL) >= mincutviolation )
709 {
712 }
713 else
714 {
715 SCIPdebugMsg(scip, "rowprep violation %g below mincutviolation %g\n", SCIPgetRowprepViolation(scip, rowprep, sol,
716 NULL), mincutviolation);
717 /* SCIPprintRowprep(scip, rowprep, NULL); */
718 }
719
720 /* free memory */
722
723 return SCIP_OKAY;
724}
725
726/** checks if an expression is quadratic and collects all occurring expressions
727 *
728 * @pre `expr2idx` and `occurringexprs` need to be initialized with capacity 2 * nchildren
729 *
730 * @note We assume that a linear term always appears before its corresponding
731 * quadratic term in quadexpr; this should be ensured by canonicalize
732 */
733static
735 SCIP* scip, /**< SCIP data structure */
736 SCIP_EXPR* quadexpr, /**< candidate for a quadratic expression */
737 SCIP_HASHMAP* expr2idx, /**< hashmap to store expressions */
738 SCIP_EXPR** occurringexprs, /**< array to store expressions */
739 int* nexprs, /**< buffer to store number of expressions */
740 SCIP_Bool* success /**< buffer to store whether the check was successful */
741 )
742{
743 SCIP_EXPR** children;
744 int nchildren;
745 int i;
746
747 assert(scip != NULL);
748 assert(quadexpr != NULL);
749 assert(expr2idx != NULL);
751 assert(nexprs != NULL);
752 assert(success != NULL);
753
754 *nexprs = 0;
755 *success = FALSE;
756 children = SCIPexprGetChildren(quadexpr);
757 nchildren = SCIPexprGetNChildren(quadexpr);
758
759 /* iterate in reverse order to ensure that quadratic terms are found before linear terms */
760 for( i = nchildren - 1; i >= 0; --i )
761 {
762 SCIP_EXPR* child;
763
764 child = children[i];
765 if( SCIPisExprPower(scip, child) )
766 {
768
769 if( SCIPgetExponentExprPow(child) != 2.0 )
770 return SCIP_OKAY;
771
772 childarg = SCIPexprGetChildren(child)[0];
773
774 if( !SCIPhashmapExists(expr2idx, (void*) childarg) )
775 {
776 SCIP_CALL( SCIPhashmapInsertInt(expr2idx, (void*) childarg, *nexprs) );
777
778 /* store the expression so we know it later */
779 assert(*nexprs < 2 * nchildren);
780 occurringexprs[*nexprs] = childarg;
781
782 ++(*nexprs);
783 }
784 }
785 else if( SCIPisExprVar(scip, child) && SCIPvarIsBinary(SCIPgetVarExprVar(child)) )
786 {
787 if( !SCIPhashmapExists(expr2idx, (void*) child) )
788 {
789 SCIP_CALL( SCIPhashmapInsertInt(expr2idx, (void*) child, *nexprs) );
790
791 /* store the expression so we know it later */
792 assert(*nexprs < 2 * nchildren);
793 occurringexprs[*nexprs] = child;
794
795 ++(*nexprs);
796 }
797 }
798 else if( SCIPisExprProduct(scip, child) )
799 {
802
803 if( SCIPexprGetNChildren(child) != 2 )
804 return SCIP_OKAY;
805
806 childarg1 = SCIPexprGetChildren(child)[0];
807 childarg2 = SCIPexprGetChildren(child)[1];
808
809 if( !SCIPhashmapExists(expr2idx, (void*) childarg1) )
810 {
811 SCIP_CALL( SCIPhashmapInsertInt(expr2idx, (void*) childarg1, *nexprs) );
812
813 /* store the expression so we know it later */
814 assert(*nexprs < 2 * nchildren);
815 occurringexprs[*nexprs] = childarg1;
816
817 ++(*nexprs);
818 }
819
820 if( !SCIPhashmapExists(expr2idx, (void*) childarg2) )
821 {
822 SCIP_CALL( SCIPhashmapInsertInt(expr2idx, (void*) childarg2, *nexprs) );
823
824 /* store the expression so we know it later */
825 assert(*nexprs < 2 * nchildren);
826 occurringexprs[*nexprs] = childarg2;
827
828 ++(*nexprs);
829 }
830 }
831 else
832 {
833 /* if there is a linear term without corresponding quadratic term, it is not a SOC */
834 if( !SCIPhashmapExists(expr2idx, (void*) child) )
835 return SCIP_OKAY;
836 }
837 }
838
839 *success = TRUE;
840
841 return SCIP_OKAY;
842}
843
844/* builds the constraint defining matrix and vector of a quadratic expression
845 *
846 * @pre `quadmatrix` and `linvector` need to be initialized with size `nexprs`^2 and `nexprs`, resp.
847 */
848static
850 SCIP* scip, /**< SCIP data structure */
851 SCIP_EXPR* quadexpr, /**< the quadratic expression */
852 SCIP_HASHMAP* expr2idx, /**< hashmap mapping the occurring expressions to their index */
853 int nexprs, /**< number of occurring expressions */
854 SCIP_Real* quadmatrix, /**< pointer to store (the lower-left triangle of) the quadratic matrix */
855 SCIP_Real* linvector /**< pointer to store the linear vector */
856 )
857{
858 SCIP_EXPR** children;
859 SCIP_Real* childcoefs;
860 int nchildren;
861 int i;
862
863 assert(scip != NULL);
864 assert(quadexpr != NULL);
865 assert(expr2idx != NULL);
868
869 children = SCIPexprGetChildren(quadexpr);
870 nchildren = SCIPexprGetNChildren(quadexpr);
872
873 /* iterate over children to build the constraint defining matrix and vector */
874 for( i = 0; i < nchildren; ++i )
875 {
876 int varpos;
877
878 if( SCIPisExprPower(scip, children[i]) )
879 {
880 assert(SCIPgetExponentExprPow(children[i]) == 2.0);
881 assert(SCIPhashmapExists(expr2idx, (void*) SCIPexprGetChildren(children[i])[0]));
882
883 varpos = SCIPhashmapGetImageInt(expr2idx, (void*) SCIPexprGetChildren(children[i])[0]);
884 assert(0 <= varpos && varpos < nexprs);
885
886 quadmatrix[varpos * nexprs + varpos] = childcoefs[i];
887 }
888 else if( SCIPisExprVar(scip, children[i]) && SCIPvarIsBinary(SCIPgetVarExprVar(children[i])) )
889 {
890 assert(SCIPhashmapExists(expr2idx, (void*) children[i]));
891
892 varpos = SCIPhashmapGetImageInt(expr2idx, (void*) children[i]);
893 assert(0 <= varpos && varpos < nexprs);
894
895 quadmatrix[varpos * nexprs + varpos] = childcoefs[i];
896 }
897 else if( SCIPisExprProduct(scip, children[i]) )
898 {
899 int varpos2;
900
901 assert(SCIPexprGetNChildren(children[i]) == 2);
902 assert(SCIPhashmapExists(expr2idx, (void*) SCIPexprGetChildren(children[i])[0]));
903 assert(SCIPhashmapExists(expr2idx, (void*) SCIPexprGetChildren(children[i])[1]));
904
905 varpos = SCIPhashmapGetImageInt(expr2idx, (void*) SCIPexprGetChildren(children[i])[0]);
906 assert(0 <= varpos && varpos < nexprs);
907
909 assert(0 <= varpos2 && varpos2 < nexprs);
910 assert(varpos != varpos2);
911
912 /* Lapack uses only the lower left triangle of the symmetric matrix */
913 quadmatrix[MIN(varpos, varpos2) * nexprs + MAX(varpos, varpos2)] = childcoefs[i] / 2.0;
914 }
915 else
916 {
917 varpos = SCIPhashmapGetImageInt(expr2idx, (void*) children[i]);
918 assert(0 <= varpos && varpos < nexprs);
919
920 linvector[varpos] = childcoefs[i];
921 }
922 }
923}
924
925/** tries to fill the nlhdlrexprdata for a potential quadratic SOC expression
926 *
927 * We say "try" because the expression might still turn out not to be a SOC at this point.
928 */
929static
931 SCIP* scip, /**< SCIP data structure */
932 SCIP_EXPR** occurringexprs, /**< array of all occurring expressions (nvars many) */
933 SCIP_Real* eigvecmatrix, /**< array containing the Eigenvectors */
934 SCIP_Real* eigvals, /**< array containing the Eigenvalues */
935 SCIP_Real* bp, /**< product of linear vector b * P (eigvecmatrix^t) */
936 int nvars, /**< number of variables */
937 int* termbegins, /**< pointer to store the termbegins */
938 SCIP_Real* transcoefs, /**< pointer to store the transcoefs */
939 int* transcoefsidx, /**< pointer to store the transcoefsidx */
940 SCIP_Real* offsets, /**< pointer to store the offsets */
941 SCIP_Real* lhsconstant, /**< pointer to store the lhsconstant */
942 int* nterms, /**< pointer to store the total number of terms */
943 SCIP_Bool* success /**< whether the expression is indeed a SOC */
944 )
945{
946 SCIP_Real sqrteigval;
947 int nextterm = 0;
948 int nexttranscoef = 0;
949 int specialtermidx;
950 int i;
951 int j;
952
953 assert(scip != NULL);
956 assert(eigvals != NULL);
957 assert(bp != NULL);
958 assert(termbegins != NULL);
959 assert(transcoefs != NULL);
960 assert(transcoefsidx != NULL);
961 assert(offsets != NULL);
963 assert(success != NULL);
964
965 *success = FALSE;
966 *nterms = 0;
967
968 /* we have lhsconstant + x^t A x + b x <= 0 and A has a single negative eigenvalue; try to build soc;
969 * we now store all the v_i^T x + beta_i on the lhs, and compute the constant
970 */
971 specialtermidx = -1;
972 for( i = 0; i < nvars; ++i )
973 {
974 if( SCIPisZero(scip, eigvals[i]) )
975 continue;
976
977 if( eigvals[i] < 0.0 )
978 {
979 assert(specialtermidx == -1); /* there should only be one negative eigenvalue */
980
982
983 *lhsconstant -= bp[i] * bp[i] / (4.0 * eigvals[i]);
984
985 continue;
986 }
987
988 assert(eigvals[i] > 0.0);
990
991 termbegins[nextterm] = nexttranscoef;
992 offsets[nextterm] = bp[i] / (2.0 * sqrteigval);
993 *lhsconstant -= bp[i] * bp[i] / (4.0 * eigvals[i]);
994
995 /* set transcoefs */
996 for( j = 0; j < nvars; ++j )
997 {
998 if( !SCIPisZero(scip, eigvecmatrix[i * nvars + j]) )
999 {
1000 transcoefs[nexttranscoef] = sqrteigval * eigvecmatrix[i * nvars + j];
1001 transcoefsidx[nexttranscoef] = j;
1002
1003 ++nexttranscoef;
1004 }
1005 }
1006 ++nextterm;
1007 }
1008 assert(specialtermidx > -1);
1009
1010 /* process constant; if constant is negative -> no soc */
1012 return SCIP_OKAY;
1013
1014 /* we need lhsconstant to be >= 0 */
1015 if( *lhsconstant < 0.0 )
1016 *lhsconstant = 0.0;
1017
1018 /* store constant term */
1019 if( *lhsconstant > 0.0 )
1020 {
1021 termbegins[nextterm] = nexttranscoef;
1022 offsets[nextterm] = SQRT( *lhsconstant );
1023 ++nextterm;
1024 }
1025
1026 /* now process rhs */
1027 {
1028 SCIP_Real rhstermlb;
1029 SCIP_Real rhstermub;
1030 SCIP_Real signfactor;
1031
1032 assert(-eigvals[specialtermidx] > 0.0);
1034
1035 termbegins[nextterm] = nexttranscoef;
1036 offsets[nextterm] = -bp[specialtermidx] / (2.0 * sqrteigval);
1037
1038 /* the expression can only be an soc if the resulting rhs term does not change sign;
1039 * the rhs term is a linear combination of variables, so estimate its bounds
1040 */
1041 rhstermlb = offsets[nextterm];
1042 for( j = 0; j < nvars; ++j )
1043 {
1044 SCIP_INTERVAL activity;
1045 SCIP_Real aux;
1046
1048 continue;
1049
1052
1053 if( eigvecmatrix[specialtermidx * nvars + j] > 0.0 )
1054 {
1055 aux = activity.inf;
1056 assert(!SCIPisInfinity(scip, aux));
1057 }
1058 else
1059 {
1060 aux = activity.sup;
1061 assert(!SCIPisInfinity(scip, -aux));
1062 }
1063
1064 if( SCIPisInfinity(scip, aux) || SCIPisInfinity(scip, -aux) )
1065 {
1067 break;
1068 }
1069 else
1071 }
1072
1073 rhstermub = offsets[nextterm];
1074 for( j = 0; j < nvars; ++j )
1075 {
1076 SCIP_INTERVAL activity;
1077 SCIP_Real aux;
1078
1080 continue;
1081
1084
1085 if( eigvecmatrix[specialtermidx * nvars + j] > 0.0 )
1086 {
1087 aux = activity.sup;
1088 assert(!SCIPisInfinity(scip, -aux));
1089 }
1090 else
1091 {
1092 aux = activity.inf;
1093 assert(!SCIPisInfinity(scip, aux));
1094 }
1095
1096 if( SCIPisInfinity(scip, aux) || SCIPisInfinity(scip, -aux) )
1097 {
1099 break;
1100 }
1101 else
1103 }
1104
1105 /* since we are just interested in obtaining an interval that contains the real bounds
1106 * and is tight enough so that we can identify that the rhsvar does not change sign,
1107 * we swap the bounds in case of numerical troubles
1108 */
1109 if( rhstermub < rhstermlb )
1110 {
1113 }
1114
1115 /* if rhs changes sign -> not a SOC */
1116 if( SCIPisLT(scip, rhstermlb, 0.0) && SCIPisGT(scip, rhstermub, 0.0) )
1117 return SCIP_OKAY;
1118
1119 signfactor = SCIPisLE(scip, rhstermub, 0.0) ? -1.0 : 1.0;
1120
1121 offsets[nextterm] *= signfactor;
1122
1123 /* set transcoefs for rhs term */
1124 for( j = 0; j < nvars; ++j )
1125 {
1127 continue;
1128
1130 transcoefsidx[nexttranscoef] = j;
1131
1132 ++nexttranscoef;
1133 }
1134
1135 /* if rhs is a constant this method shouldn't have been called */
1136 assert(nexttranscoef > termbegins[nextterm]);
1137
1138 /* finish processing term */
1139 ++nextterm;
1140 }
1141
1142 *nterms = nextterm;
1143
1144 /* sentinel value */
1145 termbegins[nextterm] = nexttranscoef;
1146
1147 *success = TRUE;
1148
1149 return SCIP_OKAY;
1150}
1151
1152/** detects if expr &le; auxvar is of the form SQRT(sum_i coef_i (expr_i + shift_i)^2 + const) &le; auxvar
1153 *
1154 * @note if a user inputs the above expression with `const` = -epsilon, then `const` is going to be set to 0.
1155 */
1156static
1158 SCIP* scip, /**< SCIP data structure */
1159 SCIP_EXPR* expr, /**< expression */
1160 SCIP_NLHDLREXPRDATA** nlhdlrexprdata, /**< pointer to store nonlinear handler expression data */
1161 SCIP_Bool* success /**< pointer to store whether SOC structure has been detected */
1162 )
1163{
1164 SCIP_EXPR** children;
1165 SCIP_EXPR* child;
1166 SCIP_EXPR** vars;
1168 SCIP_HASHSET* linexprs;
1169 SCIP_Real* childcoefs;
1170 SCIP_Real* offsets;
1171 SCIP_Real* transcoefs;
1172 SCIP_Real constant;
1173 SCIP_Bool issoc;
1174 int* transcoefsidx;
1175 int* termbegins;
1176 int nchildren;
1177 int nterms;
1178 int nvars;
1179 int nextentry;
1180 int i;
1181
1182 assert(expr != NULL);
1183 assert(success != NULL);
1184
1185 *success = FALSE;
1186 issoc = TRUE;
1187
1188 /* relation is not "<=" -> skip */
1189 if( SCIPgetExprNLocksPosNonlinear(expr) == 0 )
1190 return SCIP_OKAY;
1191
1192 assert(SCIPexprGetNChildren(expr) > 0);
1193
1194 child = SCIPexprGetChildren(expr)[0];
1195 assert(child != NULL);
1196
1197 /* check whether expression is a SQRT and has a sum as child with at least 2 children and a non-negative constant */
1198 if( ! SCIPisExprPower(scip, expr)
1199 || SCIPgetExponentExprPow(expr) != 0.5
1200 || !SCIPisExprSum(scip, child)
1201 || SCIPexprGetNChildren(child) < 2
1202 || SCIPgetConstantExprSum(child) < 0.0)
1203 {
1204 return SCIP_OKAY;
1205 }
1206
1207 /* assert(SCIPvarGetLbLocal(auxvar) >= 0.0); */
1208
1209 /* get children of the sum */
1210 children = SCIPexprGetChildren(child);
1211 nchildren = SCIPexprGetNChildren(child);
1213
1214 /* TODO: should we initialize the hashmap with size SCIPgetNVars() so that it never has to be resized? */
1216 SCIP_CALL( SCIPhashsetCreate(&linexprs, SCIPblkmem(scip), nchildren) );
1217
1218 /* we create coefs array here already, since we have to fill it in first loop in case of success
1219 * +1 for auxvar
1220 */
1221 SCIP_CALL( SCIPallocBufferArray(scip, &transcoefs, nchildren+1) );
1222
1223 nterms = 0;
1224
1225 /* check if all children are squares or linear terms with matching square term:
1226 * if the i-th child is (pow, expr, 2) we store the association <|expr -> i|> in expr2idx and if expr was in
1227 * linexprs, we remove it from there.
1228 * if the i-th child is expr' (different from (pow, expr, 2)) and expr' is not a key of expr2idx, we add it
1229 * to linexprs.
1230 * if at the end there is any expr in linexpr -> we do not have a separable quadratic function.
1231 */
1232 for( i = 0; i < nchildren; ++i )
1233 {
1234 /* handle quadratic expressions children */
1235 if( SCIPisExprPower(scip, children[i]) && SCIPgetExponentExprPow(children[i]) == 2.0 )
1236 {
1237 SCIP_EXPR* squarearg = SCIPexprGetChildren(children[i])[0];
1238
1239 if( !SCIPhashmapExists(expr2idx, (void*) squarearg) )
1240 {
1242 }
1243
1244 if( childcoefs[i] < 0.0 )
1245 {
1246 issoc = FALSE;
1247 break;
1248 }
1249 transcoefs[nterms] = SQRT(childcoefs[i]);
1250
1251 SCIP_CALL( SCIPhashsetRemove(linexprs, (void*) squarearg) );
1252 ++nterms;
1253 }
1254 /* handle binary variable children */
1255 else if( SCIPisExprVar(scip, children[i]) && SCIPvarIsBinary(SCIPgetVarExprVar(children[i])) )
1256 {
1257 assert(!SCIPhashmapExists(expr2idx, (void*) children[i]));
1258 assert(!SCIPhashsetExists(linexprs, (void*) children[i]));
1259
1260 SCIP_CALL( SCIPhashmapInsertInt(expr2idx, (void *) children[i], nterms) );
1261
1262 if( childcoefs[i] < 0.0 )
1263 {
1264 issoc = FALSE;
1265 break;
1266 }
1267 transcoefs[nterms] = SQRT(childcoefs[i]);
1268
1269 ++nterms;
1270 }
1271 else
1272 {
1273 if( !SCIPhashmapExists(expr2idx, (void*) children[i]) )
1274 {
1275 SCIP_CALL( SCIPhashsetInsert(linexprs, SCIPblkmem(scip), (void*) children[i]) );
1276 }
1277 }
1278 }
1279
1280 /* there are linear terms without corresponding quadratic terms or it was detected not to be soc */
1281 if( SCIPhashsetGetNElements(linexprs) > 0 || ! issoc )
1282 {
1283 SCIPfreeBufferArray(scip, &transcoefs);
1284 SCIPhashsetFree(&linexprs, SCIPblkmem(scip) );
1286 return SCIP_OKAY;
1287 }
1288
1289 /* add one to terms counter for auxvar */
1290 ++nterms;
1291
1292 constant = SCIPgetConstantExprSum(child);
1293
1294 /* compute constant of possible soc expression to check its sign */
1295 for( i = 0; i < nchildren; ++i )
1296 {
1297 if( ! SCIPisExprPower(scip, children[i]) || SCIPgetExponentExprPow(children[i]) != 2.0 )
1298 {
1299 int auxvarpos;
1300
1301 assert(SCIPhashmapExists(expr2idx, (void*) children[i]) );
1302 auxvarpos = SCIPhashmapGetImageInt(expr2idx, (void*) children[i]);
1303
1304 constant -= SQR(0.5 * childcoefs[i] / transcoefs[auxvarpos]);
1305 }
1306 }
1307
1308 /* if the constant is negative -> no SOC */
1309 if( SCIPisNegative(scip, constant) )
1310 {
1311 SCIPfreeBufferArray(scip, &transcoefs);
1312 SCIPhashsetFree(&linexprs, SCIPblkmem(scip) );
1314 return SCIP_OKAY;
1315 }
1316 else if( SCIPisZero(scip, constant) )
1317 constant = 0.0;
1318 assert(constant >= 0.0);
1319
1320 /* at this point, we have found an SOC structure */
1321 *success = TRUE;
1322
1323 nvars = nterms;
1324
1325 /* add one to terms counter for constant term */
1326 if( constant > 0.0 )
1327 ++nterms;
1328
1329 /* allocate temporary memory to collect data */
1332 SCIP_CALL( SCIPallocBufferArray(scip, &transcoefsidx, nvars) );
1333 SCIP_CALL( SCIPallocBufferArray(scip, &termbegins, nterms + 1) );
1334
1335 /* fill in data for non constant terms of lhs; initialize their offsets */
1336 for( i = 0; i < nvars - 1; ++i )
1337 {
1338 transcoefsidx[i] = i;
1339 termbegins[i] = i;
1340 offsets[i] = 0.0;
1341 }
1342
1343 /* add constant term and rhs */
1344 vars[nvars - 1] = expr;
1345 if( constant > 0.0 )
1346 {
1347 /* constant term */
1348 termbegins[nterms - 2] = nterms - 2;
1349 offsets[nterms - 2] = SQRT(constant);
1350
1351 /* rhs */
1352 termbegins[nterms - 1] = nterms - 2;
1353 offsets[nterms - 1] = 0.0;
1354 transcoefsidx[nterms - 2] = nvars - 1;
1355 transcoefs[nterms - 2] = 1.0;
1356
1357 /* sentinel value */
1358 termbegins[nterms] = nterms - 1;
1359 }
1360 else
1361 {
1362 /* rhs */
1363 termbegins[nterms - 1] = nterms - 1;
1364 offsets[nterms - 1] = 0.0;
1365 transcoefsidx[nterms - 1] = nvars - 1;
1366 transcoefs[nterms - 1] = 1.0;
1367
1368 /* sentinel value */
1369 termbegins[nterms] = nterms;
1370 }
1371
1372 /* request required auxiliary variables and fill vars and offsets array */
1373 nextentry = 0;
1374 for( i = 0; i < nchildren; ++i )
1375 {
1376 if( SCIPisExprPower(scip, children[i]) && SCIPgetExponentExprPow(children[i]) == 2.0 )
1377 {
1379
1380 squarearg = SCIPexprGetChildren(children[i])[0];
1382
1384
1386 ++nextentry;
1387 }
1388 else if( SCIPisExprVar(scip, children[i]) && SCIPvarIsBinary(SCIPgetVarExprVar(children[i])) )
1389 {
1390 /* handle binary variable children: no need to request auxvar */
1391 assert(SCIPhashmapGetImageInt(expr2idx, (void*) children[i]) == nextentry);
1392 vars[nextentry] = children[i];
1393 ++nextentry;
1394 }
1395 else
1396 {
1397 int auxvarpos;
1398
1399 assert(SCIPhashmapExists(expr2idx, (void*) children[i]));
1400 auxvarpos = SCIPhashmapGetImageInt(expr2idx, (void*) children[i]);
1401
1403
1404 offsets[auxvarpos] = 0.5 * childcoefs[i] / transcoefs[auxvarpos];
1405 }
1406 }
1407 assert(nextentry == nvars - 1);
1408
1409#ifdef SCIP_DEBUG
1410 SCIPdebugMsg(scip, "found SOC structure for expression %p\n", (void*)expr);
1411 SCIPprintExpr(scip, expr, NULL);
1412 SCIPinfoMessage(scip, NULL, " <= auxvar\n");
1413#endif
1414
1415 /* create and store nonlinear handler expression data */
1416 SCIP_CALL( createNlhdlrExprData(scip, vars, offsets, transcoefs, transcoefsidx, termbegins,
1417 nvars, nterms, nlhdlrexprdata) );
1418 assert(*nlhdlrexprdata != NULL);
1419
1420 /* free memory */
1421 SCIPhashsetFree(&linexprs, SCIPblkmem(scip) );
1423 SCIPfreeBufferArray(scip, &termbegins);
1424 SCIPfreeBufferArray(scip, &transcoefsidx);
1425 SCIPfreeBufferArray(scip, &offsets);
1427 SCIPfreeBufferArray(scip, &transcoefs);
1428
1429 return SCIP_OKAY;
1430}
1431
1432/** helper method to detect c + sum_i coef_i expr_i^2 - coef_k expr_k^2 &le; 0
1433 * and c + sum_i coef_i expr_i^2 - coef_k expr_k expr_l &le; 0
1434 *
1435 * binary linear variables are interpreted as quadratic terms
1436 *
1437 * @todo: extend this function to detect c + sum_i coef_i (expr_i + const_i)^2 - ...
1438 * this would probably share a lot of code with detectSocNorm
1439 */
1440static
1442 SCIP* scip, /**< SCIP data structure */
1443 SCIP_EXPR* expr, /**< expression */
1444 SCIP_Real conslhs, /**< lhs of the constraint that the expression defines (or SCIP_INVALID) */
1445 SCIP_Real consrhs, /**< rhs of the constraint that the expression defines (or SCIP_INVALID) */
1446 SCIP_NLHDLREXPRDATA** nlhdlrexprdata, /**< pointer to store nonlinear handler expression data */
1447 SCIP_Bool* enforcebelow, /**< pointer to store whether we enforce <= (TRUE) or >= (FALSE); only valid when success is TRUE */
1448 SCIP_Bool* success /**< pointer to store whether SOC structure has been detected */
1449 )
1450{
1451 SCIP_EXPR** children;
1452 SCIP_EXPR** vars = NULL;
1453 SCIP_Real* childcoefs;
1454 SCIP_Real* offsets = NULL;
1455 SCIP_Real* transcoefs = NULL;
1456 int* transcoefsidx = NULL;
1457 int* termbegins = NULL;
1458 SCIP_Real constant;
1459 SCIP_Real lhsconstant;
1460 SCIP_Real lhs;
1461 SCIP_Real rhs;
1462 SCIP_Real rhssign;
1464 int ntranscoefs;
1465 int nposquadterms;
1466 int nnegquadterms;
1467 int nposbilinterms;
1468 int nnegbilinterms;
1469 int rhsidx;
1470 int lhsidx;
1471 int specialtermidx;
1472 int nchildren;
1473 int nnzinterms;
1474 int nterms;
1475 int nvars;
1476 int nextentry;
1477 int i;
1478 SCIP_Bool ishyperbolic;
1479
1480 assert(expr != NULL);
1481 assert(success != NULL);
1482
1483 *success = FALSE;
1484
1485 /* check whether expression is a sum with at least 2 children */
1486 if( ! SCIPisExprSum(scip, expr) || SCIPexprGetNChildren(expr) < 2 )
1487 return SCIP_OKAY;
1488
1489 /* get children of the sum */
1490 children = SCIPexprGetChildren(expr);
1491 nchildren = SCIPexprGetNChildren(expr);
1492 constant = SCIPgetConstantExprSum(expr);
1493
1494 /* we duplicate the child coefficients since we have to manipulate them */
1495 SCIP_CALL( SCIPduplicateBufferArray(scip, &childcoefs, SCIPgetCoefsExprSum(expr), nchildren) ); /*lint !e666*/
1496
1497 /* initialize data */
1498 lhsidx = -1;
1499 rhsidx = -1;
1500 nposquadterms = 0;
1501 nnegquadterms = 0;
1502 nposbilinterms = 0;
1503 nnegbilinterms = 0;
1504
1505 /* check if all children are quadratic or binary linear and count number of positive and negative terms */
1506 for( i = 0; i < nchildren; ++i )
1507 {
1508 if( SCIPisExprPower(scip, children[i]) && SCIPgetExponentExprPow(children[i]) == 2.0 )
1509 {
1510 if( childcoefs[i] > 0.0 )
1511 {
1512 ++nposquadterms;
1513 lhsidx = i;
1514 }
1515 else
1516 {
1517 ++nnegquadterms;
1518 rhsidx = i;
1519 }
1520 }
1521 else if( SCIPisExprVar(scip, children[i]) && SCIPvarIsBinary(SCIPgetVarExprVar(children[i])) )
1522 {
1523 if( childcoefs[i] > 0.0 )
1524 {
1525 ++nposquadterms;
1526 lhsidx = i;
1527 }
1528 else
1529 {
1530 ++nnegquadterms;
1531 rhsidx = i;
1532 }
1533 }
1534 else if( SCIPisExprProduct(scip, children[i]) && SCIPexprGetNChildren(children[i]) == 2 )
1535 {
1536 if( childcoefs[i] > 0.0 )
1537 {
1539 lhsidx = i;
1540 }
1541 else
1542 {
1544 rhsidx = i;
1545 }
1546 }
1547 else
1548 {
1549 goto CLEANUP;
1550 }
1551
1552 /* more than one positive eigenvalue and more than one negative eigenvalue -> can't be convex */
1553 if( nposquadterms > 1 && nnegquadterms > 1 )
1554 goto CLEANUP;
1555
1556 /* more than one bilinear term -> can't be handled by this method */
1557 if( nposbilinterms + nnegbilinterms > 1 )
1558 goto CLEANUP;
1559
1560 /* one positive bilinear term and also at least one positive quadratic term -> not a simple SOC */
1561 if( nposbilinterms > 0 && nposquadterms > 0 )
1562 goto CLEANUP;
1563
1564 /* one negative bilinear term and also at least one negative quadratic term -> not a simple SOC */
1565 if( nnegbilinterms > 0 && nnegquadterms > 0 )
1566 goto CLEANUP;
1567 }
1568
1569 if( nposquadterms == nchildren || nnegquadterms == nchildren )
1570 goto CLEANUP;
1571
1572 assert(nposquadterms <= 1 || nnegquadterms <= 1);
1574 assert(nposbilinterms == 0 || nposquadterms == 0);
1575 assert(nnegbilinterms == 0 || nnegquadterms == 0);
1576
1577 /* if a bilinear term is involved, it is a hyperbolic expression */
1579
1580 if( conslhs == SCIP_INVALID || consrhs == SCIP_INVALID ) /*lint !e777*/
1581 {
1584
1585 lhs = (conslhs == SCIP_INVALID ? expractivity.inf : conslhs); /*lint !e777*/
1586 rhs = (consrhs == SCIP_INVALID ? expractivity.sup : consrhs); /*lint !e777*/
1587 }
1588 else
1589 {
1590 lhs = conslhs;
1591 rhs = consrhs;
1592 }
1593
1594 /* detect case and store lhs/rhs information */
1595 if( (ishyperbolic && nnegbilinterms > 0) || (!ishyperbolic && nnegquadterms < 2) )
1596 {
1597 /* we have -x*y + z^2 ... -> we want to write z^2 ... <= x*y;
1598 * or we have -x^2 + y^2 ... -> we want to write y^2 ... <= x^2;
1599 * in any case, we need a finite rhs
1600 */
1601 assert(nnegbilinterms == 1 || nnegquadterms == 1);
1602 assert(rhsidx != -1);
1603
1604 /* if rhs is infinity, it can't be soc
1605 * TODO: if it can't be soc, then we should enforce the caller so that we do not try the more complex quadratic
1606 * method
1607 */
1608 if( SCIPisInfinity(scip, rhs) )
1609 goto CLEANUP;
1610
1611 specialtermidx = rhsidx;
1612 lhsconstant = constant - rhs;
1613 *enforcebelow = TRUE; /* enforce expr <= rhs */
1614 }
1615 else
1616 {
1617 assert(lhsidx != -1);
1618
1619 /* if lhs is infinity, it can't be soc */
1620 if( SCIPisInfinity(scip, -lhs) )
1621 goto CLEANUP;
1622
1624 lhsconstant = lhs - constant;
1625
1626 /* negate all coefficients */
1627 for( i = 0; i < nchildren; ++i )
1628 childcoefs[i] = -childcoefs[i];
1629 *enforcebelow = FALSE; /* enforce lhs <= expr */
1630 }
1632
1633 if( ishyperbolic )
1634 {
1637
1639
1642
1645
1646 if( SCIPisNegative(scip, yactivity.inf + zactivity.inf) )
1647 {
1648 /* the sum of the expressions in the bilinear term changes sign -> no SOC */
1649 if( SCIPisPositive(scip, yactivity.sup + zactivity.sup) )
1650 goto CLEANUP;
1651
1652 rhssign = -1.0;
1653 }
1654 else
1655 rhssign = 1.0;
1656
1658 }
1659 else if( SCIPisExprVar(scip, children[specialtermidx]) )
1660 {
1661 /* children[specialtermidx] can be a variable, in which case we treat it as if it is squared */
1662 rhssign = 1.0;
1663 }
1664 else
1665 {
1667
1671
1672 if( rhsactivity.inf < 0.0 )
1673 {
1674 /* rhs variable changes sign -> no SOC */
1675 if( rhsactivity.sup > 0.0 )
1676 goto CLEANUP;
1677
1678 rhssign = -1.0;
1679 }
1680 else
1681 rhssign = 1.0;
1682 }
1683
1685 goto CLEANUP;
1686
1688 lhsconstant = 0.0;
1689
1690 /*
1691 * we have found an SOC-representable expression. Now build the nlhdlrexprdata
1692 *
1693 * in the non-hyperbolic case, c + sum_i coef_i expr_i^2 - coef_k expr_k^2 <= 0 is transformed to
1694 * SQRT( c + sum_i coef_i expr_i^2 ) <= coef_k expr_k
1695 * so there are nchildren many vars, nchildren (+ 1 if c != 0) many terms, nchildren many coefficients in the vs
1696 * in SOC representation
1697 *
1698 * in the hyperbolic case, c + sum_i coef_i expr_i^2 - coef_k expr_k expr_l <= 0 is transformed to
1699 * SQRT( 4(c + sum_i coef_i expr_i^2) + (expr_k - expr_l)^2 ) <= expr_k + expr_l
1700 * so there are nchildren + 1many vars, nchildren + 1(+ 1 if c != 0) many terms, nchildren + 3 many coefficients in
1701 * the vs in SOC representation
1702 */
1703
1704 ntranscoefs = ishyperbolic ? nchildren + 3 : nchildren;
1705 nvars = ishyperbolic ? nchildren + 1 : nchildren;
1706 nterms = nvars;
1707
1708 /* constant term */
1709 if( lhsconstant > 0.0 )
1710 nterms++;
1711
1712 /* SOC was detected, allocate temporary memory for data to collect */
1716 SCIP_CALL( SCIPallocBufferArray(scip, &transcoefsidx, ntranscoefs) );
1717 SCIP_CALL( SCIPallocBufferArray(scip, &termbegins, nterms + 1) );
1718
1719 *success = TRUE;
1720 nextentry = 0;
1721
1722 /* collect all the v_i and beta_i */
1723 nnzinterms = 0;
1724 for( i = 0; i < nchildren; ++i )
1725 {
1726 /* variable and coef for rhs have to be set to the last entry */
1727 if( i == specialtermidx )
1728 continue;
1729
1730 /* extract (unique) variable appearing in term */
1731 if( SCIPisExprVar(scip, children[i]) )
1732 {
1733 vars[nextentry] = children[i];
1734
1736 }
1737 else
1738 {
1739 assert(SCIPisExprPower(scip, children[i]));
1740
1741 /* notify that we will require auxiliary variable */
1743 vars[nextentry] = SCIPexprGetChildren(children[i])[0];
1744 }
1746
1747 /* store v_i and beta_i */
1748 termbegins[nextentry] = nnzinterms;
1749 offsets[nextentry] = 0.0;
1750
1751 transcoefsidx[nnzinterms] = nextentry;
1752 if( ishyperbolic )
1753 {
1754 /* we eliminate the coefficient of the bilinear term to arrive at standard form */
1755 assert(4.0 * childcoefs[i] / -childcoefs[specialtermidx] > 0.0);
1756 transcoefs[nnzinterms] = SQRT(4.0 * childcoefs[i] / -childcoefs[specialtermidx]);
1757 }
1758 else
1759 {
1760 assert(childcoefs[i] > 0.0);
1761 transcoefs[nnzinterms] = SQRT(childcoefs[i]);
1762 }
1763
1764 /* finish adding nonzeros */
1765 ++nnzinterms;
1766
1767 /* finish processing term */
1768 ++nextentry;
1769 }
1770 assert(nextentry == nchildren - 1);
1771
1772 /* store term for constant (v_i = 0) */
1773 if( lhsconstant > 0.0 )
1774 {
1775 termbegins[nextentry] = nnzinterms;
1776 offsets[nextentry] = SQRT(lhsconstant);
1777
1778 /* finish processing term; this term has 0 nonzero thus we do not increase nnzinterms */
1779 ++nextentry;
1780 }
1781
1782 if( !ishyperbolic )
1783 {
1784 /* store rhs term */
1785 if( SCIPisExprVar(scip, children[specialtermidx]) )
1786 {
1787 /* this should be the "children[specialtermidx] can be a variable, in which case we treat it as if it is squared" case */
1789 vars[nvars - 1] = children[specialtermidx];
1790 }
1791 else
1792 {
1796 vars[nvars - 1] = SCIPexprGetChildren(children[specialtermidx])[0];
1797 }
1798
1800
1801 termbegins[nextentry] = nnzinterms;
1802 offsets[nextentry] = 0.0;
1803 transcoefs[nnzinterms] = rhssign * SQRT(-childcoefs[specialtermidx]);
1804 transcoefsidx[nnzinterms] = nvars - 1;
1805
1806 /* finish adding nonzeros */
1807 ++nnzinterms;
1808
1809 /* finish processing term */
1810 ++nextentry;
1811 }
1812 else
1813 {
1814 /* store last lhs term and rhs term coming from the bilinear term */
1816 vars[nvars - 2] = SCIPexprGetChildren(children[specialtermidx])[0];
1817
1819 vars[nvars - 1] = SCIPexprGetChildren(children[specialtermidx])[1];
1820
1821 /* at this point, vars[nvars - 2] = expr_k and vars[nvars - 1] = expr_l;
1822 * on the lhs we have the term (expr_k - expr_l)^2
1823 */
1824 termbegins[nextentry] = nnzinterms;
1825 offsets[nextentry] = 0.0;
1826
1827 /* expr_k */
1828 transcoefsidx[nnzinterms] = nvars - 2;
1829 transcoefs[nnzinterms] = 1.0;
1830 ++nnzinterms;
1831
1832 /* - expr_l */
1833 transcoefsidx[nnzinterms] = nvars - 1;
1834 transcoefs[nnzinterms] = -1.0;
1835 ++nnzinterms;
1836
1837 /* finish processing term */
1838 ++nextentry;
1839
1840 /* on rhs we have +/-(expr_k + expr_l) */
1841 termbegins[nextentry] = nnzinterms;
1842 offsets[nextentry] = 0.0;
1843
1844 /* rhssing * expr_k */
1845 transcoefsidx[nnzinterms] = nvars - 2;
1846 transcoefs[nnzinterms] = rhssign;
1847 ++nnzinterms;
1848
1849 /* rhssing * expr_l */
1850 transcoefsidx[nnzinterms] = nvars - 1;
1851 transcoefs[nnzinterms] = rhssign;
1852 ++nnzinterms;
1853
1854 /* finish processing term */
1855 ++nextentry;
1856 }
1859
1860 /* sentinel value */
1861 termbegins[nextentry] = nnzinterms;
1862
1863#ifdef SCIP_DEBUG
1864 SCIPdebugMsg(scip, "found SOC structure for expression %p\n%f <= ", (void*)expr, lhs);
1865 SCIPprintExpr(scip, expr, NULL);
1866 SCIPinfoMessage(scip, NULL, "<= %f\n", rhs);
1867#endif
1868
1869 /* create and store nonlinear handler expression data */
1870 SCIP_CALL( createNlhdlrExprData(scip, vars, offsets, transcoefs, transcoefsidx, termbegins, nvars, nterms,
1871 nlhdlrexprdata) );
1872 assert(*nlhdlrexprdata != NULL);
1873
1874CLEANUP:
1875 SCIPfreeBufferArrayNull(scip, &termbegins);
1876 SCIPfreeBufferArrayNull(scip, &transcoefsidx);
1877 SCIPfreeBufferArrayNull(scip, &transcoefs);
1878 SCIPfreeBufferArrayNull(scip, &offsets);
1881
1882 return SCIP_OKAY;
1883}
1884
1885/** detects complex quadratic expressions that can be represented as SOC constraints
1886 *
1887 * These are quadratic expressions with either exactly one positive or exactly one negative eigenvalue,
1888 * in addition to some extra conditions. One needs to write the quadratic as
1889 * sum eigval_i (eigvec_i . x)^2 + c &le; -eigval_k (eigvec_k . x)^2, where eigval_k is the negative eigenvalue,
1890 * and c must be positive and (eigvec_k . x) must not change sign.
1891 * This is described in more details in
1892 * Mahajan, Ashutosh & Munson, Todd, Exploiting Second-Order Cone Structure for Global Optimization, 2010.
1893 *
1894 * The eigen-decomposition is computed using Lapack.
1895 * Binary linear variables are interpreted as quadratic terms.
1896 *
1897 * @todo: In the case -b <= a + x^2 - y^2 <= b, it is possible to represent both sides by SOC. Currently, the
1898 * datastructure can only handle one SOC. If this should appear more often, it could be worth to extend it,
1899 * such that both sides can be handled (see e.g. instance chp_partload).
1900 * FS: this shouldn't be possible. For a <= b + x^2 - y^2 <= c to be SOC representable on both sides, we would need
1901 * that a - b >= 0 and b -c >= 0, but this implies that a >= c and assuming the constraint is not trivially infeasible,
1902 * a <= b. Thus, a = b = c and the constraint is x^2 == y^2.
1903 *
1904 * @todo: Since cons_nonlinear multiplies as many terms out as possible during presolving, some SOC-representable
1905 * structures cannot be detected, (see e.g. instances bearing or wager). There is currently no obvious way
1906 * to handle this.
1907 */
1908static
1910 SCIP* scip, /**< SCIP data structure */
1911 SCIP_EXPR* expr, /**< expression */
1912 SCIP_Real conslhs, /**< lhs of the constraint that the expression defines (or SCIP_INVALID) */
1913 SCIP_Real consrhs, /**< rhs of the constraint that the expression defines (or SCIP_INVALID) */
1914 SCIP_NLHDLREXPRDATA** nlhdlrexprdata, /**< pointer to store nonlinear handler expression data */
1915 SCIP_Bool* enforcebelow, /**< pointer to store whether we enforce <= (TRUE) or >= (FALSE); only
1916 * valid when success is TRUE */
1917 SCIP_Bool* success /**< pointer to store whether SOC structure has been detected */
1918 )
1919{
1922 SCIP_Real* offsets;
1923 SCIP_Real* transcoefs;
1924 SCIP_Real* eigvecmatrix;
1925 SCIP_Real* eigvals;
1926 SCIP_Real* lincoefs;
1927 SCIP_Real* bp;
1928 int* transcoefsidx;
1929 int* termbegins;
1930 SCIP_Real constant;
1931 SCIP_Real lhsconstant;
1932 SCIP_Real lhs;
1933 SCIP_Real rhs;
1935 int nvars;
1936 int nterms;
1937 int nchildren;
1938 int npos;
1939 int nneg;
1940 int ntranscoefs;
1941 int i;
1942 int j;
1943 SCIP_Bool rhsissoc;
1944 SCIP_Bool lhsissoc;
1945 SCIP_Bool isquadratic;
1946
1947 assert(expr != NULL);
1948 assert(success != NULL);
1949
1950 *success = FALSE;
1951
1952 /* check whether expression is a sum with at least 2 children */
1953 if( ! SCIPisExprSum(scip, expr) || SCIPexprGetNChildren(expr) < 2 )
1954 {
1955 return SCIP_OKAY;
1956 }
1957
1958 /* we need Lapack to compute eigenvalues/vectors below */
1960 return SCIP_OKAY;
1961
1962 /* get children of the sum */
1963 nchildren = SCIPexprGetNChildren(expr);
1964 constant = SCIPgetConstantExprSum(expr);
1965
1966 /* initialize data */
1967 offsets = NULL;
1968 transcoefs = NULL;
1969 transcoefsidx = NULL;
1970 termbegins = NULL;
1971 bp = NULL;
1972
1973 SCIP_CALL( SCIPhashmapCreate(&expr2idx, SCIPblkmem(scip), 2 * nchildren) );
1975
1976 /* check if the expression is quadratic and collect all occurring expressions */
1978
1979 if( !isquadratic )
1980 {
1983 return SCIP_OKAY;
1984 }
1985
1986 /* check that nvars*nvars doesn't get too large, see also SCIPcomputeExprQuadraticCurvature() */
1987 if( nvars > 7000 )
1988 {
1989 SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, "nlhdlr_soc - number of quadratic variables is too large (%d) to check the curvature\n", nvars);
1992 return SCIP_OKAY;
1993 }
1994
1996
1997 /* create datastructures for constaint defining matrix and vector */
2000
2001 /* build constraint defining matrix (stored in eigvecmatrix) and vector (stored in lincoefs) */
2003
2005
2006 /* compute eigenvalues and vectors, A = PDP^t
2007 * note: eigvecmatrix stores P^t, i.e., P^t_{i,j} = eigvecmatrix[i*nvars+j]
2008 */
2010 {
2011 SCIPdebugMsg(scip, "Failed to compute eigenvalues and eigenvectors for expression:\n");
2012
2013#ifdef SCIP_DEBUG
2014 SCIPdismantleExpr(scip, NULL, expr);
2015#endif
2016
2017 goto CLEANUP;
2018 }
2019
2021
2022 nneg = 0;
2023 npos = 0;
2024 ntranscoefs = 0;
2025
2026 /* set small eigenvalues to 0 and compute b*P */
2027 for( i = 0; i < nvars; ++i )
2028 {
2029 for( j = 0; j < nvars; ++j )
2030 {
2031 bp[i] += lincoefs[j] * eigvecmatrix[i * nvars + j];
2032
2033 /* count the number of transcoefs to be used later */
2035 ++ntranscoefs;
2036 }
2037
2038 if( SCIPisZero(scip, eigvals[i]) )
2039 {
2040 /* if there is a purely linear variable, the constraint can't be written as a SOC */
2041 if( !SCIPisZero(scip, bp[i]) )
2042 goto CLEANUP;
2043
2044 bp[i] = 0.0;
2045 eigvals[i] = 0.0;
2046 }
2047 else if( eigvals[i] > 0.0 )
2048 npos++;
2049 else
2050 nneg++;
2051 }
2052
2053 /* a proper SOC constraint needs at least 2 variables */
2054 if( npos + nneg < 2 )
2055 goto CLEANUP;
2056
2057 /* determine whether rhs or lhs of cons is potentially SOC, if any */
2058 rhsissoc = (nneg == 1 && SCIPgetExprNLocksPosNonlinear(expr) > 0);
2059 lhsissoc = (npos == 1 && SCIPgetExprNLocksNegNonlinear(expr) > 0);
2060
2061 if( rhsissoc || lhsissoc )
2062 {
2063 if( conslhs == SCIP_INVALID || consrhs == SCIP_INVALID ) /*lint !e777*/
2064 {
2067 lhs = (conslhs == SCIP_INVALID ? expractivity.inf : conslhs); /*lint !e777*/
2068 rhs = (consrhs == SCIP_INVALID ? expractivity.sup : consrhs); /*lint !e777*/
2069 }
2070 else
2071 {
2072 lhs = conslhs;
2073 rhs = consrhs;
2074 }
2075 }
2076 else
2077 {
2078 /* if none of the sides is potentially SOC, stop */
2079 goto CLEANUP;
2080 }
2081
2082 /* @TODO: what do we do if both sides are possible? */
2083 if( !rhsissoc )
2084 {
2086
2087 /* lhs is potentially SOC, change signs */
2088 lhsconstant = lhs - constant; /*lint !e644*/
2089
2090 for( i = 0; i < nvars; ++i )
2091 {
2092 eigvals[i] = -eigvals[i];
2093 bp[i] = -bp[i];
2094 }
2095 *enforcebelow = FALSE; /* enforce lhs <= expr */
2096 }
2097 else
2098 {
2099 lhsconstant = constant - rhs; /*lint !e644*/
2100 *enforcebelow = TRUE; /* enforce expr <= rhs */
2101 }
2102
2103 /* initialize remaining datastructures for nonlinear handler */
2104 SCIP_CALL( SCIPallocBufferArray(scip, &offsets, npos + nneg + 1) );
2106 SCIP_CALL( SCIPallocBufferArray(scip, &transcoefsidx, ntranscoefs) );
2107 SCIP_CALL( SCIPallocBufferArray(scip, &termbegins, npos + nneg + 2) );
2108
2109 /* try to fill the nlhdlrexprdata (at this point, it can still fail) */
2111 transcoefsidx, offsets, &lhsconstant, &nterms, success) );
2112
2113 if( !(*success) )
2114 goto CLEANUP;
2115
2116 assert(0 < nterms && nterms <= npos + nneg + 1);
2117 assert(ntranscoefs == termbegins[nterms]);
2118
2119 /*
2120 * at this point, the expression passed all checks and is SOC-representable
2121 */
2122
2123 /* register all requests for auxiliary variables */
2124 for( i = 0; i < nvars; ++i )
2125 {
2127 }
2128
2129#ifdef SCIP_DEBUG
2130 SCIPdebugMsg(scip, "found SOC structure for expression %p\n%f <= ", (void*)expr, lhs);
2131 SCIPprintExpr(scip, expr, NULL);
2132 SCIPinfoMessage(scip, NULL, "<= %f\n", rhs);
2133#endif
2134
2135 /* finally, create and store nonlinear handler expression data */
2136 SCIP_CALL( createNlhdlrExprData(scip, occurringexprs, offsets, transcoefs, transcoefsidx, termbegins, nvars, nterms,
2137 nlhdlrexprdata) );
2138 assert(*nlhdlrexprdata != NULL);
2139
2140CLEANUP:
2141 SCIPfreeBufferArrayNull(scip, &termbegins);
2142 SCIPfreeBufferArrayNull(scip, &transcoefsidx);
2143 SCIPfreeBufferArrayNull(scip, &transcoefs);
2144 SCIPfreeBufferArrayNull(scip, &offsets);
2147 SCIPfreeBufferArray(scip, &lincoefs);
2151
2152 return SCIP_OKAY;
2153}
2154
2155/** helper method to detect SOC structures
2156 *
2157 * The detection runs in 3 steps:
2158 * 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$
2159 * 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
2160 * -> this results in the SOC expr &le; auxvar(expr)
2161 *
2162 * TODO we should generalize and check for sqrt(positive-semidefinite-quadratic)
2163 *
2164 * 2. check if expression represents a quadratic function of one of the following forms (all coefs > 0)
2165 * 1. \f$(\sum_i \text{coef}_i \text{expr}_i^2) - \text{coef}_k \text{expr}_k^2 \leq \text{RHS}\f$ or
2166 * 2. \f$(\sum_i - \text{coef}_i \text{expr}_i^2) + \text{coef}_k \text{expr}_k^2 \geq \text{LHS}\f$ or
2167 * 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
2168 * 4. \f$(\sum_i - \text{coef}_i \text{expr}_i^2) + \text{coef}_k \text{expr}_k \text{expr}_l \geq \text{LHS}\f$,
2169 *
2170 * where RHS &ge; 0 or LHS &le; 0, respectively. For LHS and RHS we use the constraint sides if it is a root expr
2171 * and the bounds of the auxiliary variable otherwise.
2172 * The last two cases are called hyperbolic or rotated second order cone.\n
2173 * -> 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$
2174 * 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$.
2175 * (analogously for the LHS cases)
2176 *
2177 * 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$
2178 * has exactly one negative eigenvalue plus some extra conditions, see detectSocQuadraticComplex().
2179 *
2180 * Note that step 3 is only performed if parameter `compeigenvalues` is set to TRUE.
2181 */
2182static
2184 SCIP* scip, /**< SCIP data structure */
2185 SCIP_NLHDLRDATA* nlhdlrdata, /**< nonlinear handler data */
2186 SCIP_EXPR* expr, /**< expression */
2187 SCIP_Real conslhs, /**< lhs of the constraint that the expression defines (or SCIP_INVALID) */
2188 SCIP_Real consrhs, /**< rhs of the constraint that the expression defines (or SCIP_INVALID) */
2189 SCIP_NLHDLREXPRDATA** nlhdlrexprdata, /**< pointer to store nonlinear handler expression data */
2190 SCIP_Bool* enforcebelow, /**< pointer to store whether we enforce <= (TRUE) or >= (FALSE); only
2191 * valid when success is TRUE */
2192 SCIP_Bool* success /**< pointer to store whether SOC structure has been detected */
2193 )
2194{
2195 assert(expr != NULL);
2197 assert(nlhdlrexprdata != NULL);
2198 assert(success != NULL);
2199
2200 *success = FALSE;
2201
2202 /* check whether expression is given as norm as described in case 1 above: if we have a constraint
2203 * sqrt(sum x_i^2) <= constant, then it might be better not to handle this here; thus, we only call detectSocNorm
2204 * when the expr is _not_ the root of a constraint
2205 */
2206 if( conslhs == SCIP_INVALID && consrhs == SCIP_INVALID ) /*lint !e777*/
2207 {
2208 SCIP_CALL( detectSocNorm(scip, expr, nlhdlrexprdata, success) );
2210 }
2211
2212 if( !(*success) )
2213 {
2214 /* check whether expression is a simple soc-respresentable quadratic expression as described in case 2 above */
2216 }
2217
2218 if( !(*success) && nlhdlrdata->compeigenvalues )
2219 {
2220 /* check whether expression is a more complex soc-respresentable quadratic expression as described in case 3 */
2222 }
2223
2224 return SCIP_OKAY;
2225}
2226
2227/*
2228 * Callback methods of nonlinear handler
2229 */
2230
2231/** nonlinear handler copy callback */
2232static
2243
2244/** callback to free data of handler */
2245static
2247{ /*lint --e{715}*/
2249
2251
2252 return SCIP_OKAY;
2253}
2254
2255/** callback to free expression specific data */
2256static
2258{ /*lint --e{715}*/
2259 assert(*nlhdlrexprdata != NULL);
2260
2261 SCIP_CALL( freeNlhdlrExprData(scip, nlhdlrexprdata) );
2262
2263 return SCIP_OKAY;
2264}
2265
2266
2267/** callback to be called in initialization */
2268#if 0
2269static
2271{ /*lint --e{715}*/
2272 SCIPerrorMessage("method of soc nonlinear handler not implemented yet\n");
2273 SCIPABORT(); /*lint --e{527}*/
2274
2275 return SCIP_OKAY;
2276}
2277#else
2278#define nlhdlrInitSoc NULL
2279#endif
2280
2281
2282/** callback to be called in deinitialization */
2283#if 0
2284static
2286{ /*lint --e{715}*/
2287 SCIPerrorMessage("method of soc nonlinear handler not implemented yet\n");
2288 SCIPABORT(); /*lint --e{527}*/
2289
2290 return SCIP_OKAY;
2291}
2292#else
2293#define nlhdlrExitSoc NULL
2294#endif
2295
2296
2297/** callback to detect structure in expression tree */
2298static
2300{ /*lint --e{715}*/
2301 SCIP_Real conslhs;
2302 SCIP_Real consrhs;
2303 SCIP_Bool enforcebelow;
2304 SCIP_Bool success;
2306
2307 assert(expr != NULL);
2308
2309 /* don't try if no sepa is required
2310 * TODO implement some bound strengthening
2311 */
2313 return SCIP_OKAY;
2314
2315 assert(SCIPgetExprNAuxvarUsesNonlinear(expr) > 0); /* since some sepa is required, there should have been demand for it */
2316
2317 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
2319
2320 conslhs = (cons == NULL ? SCIP_INVALID : SCIPgetLhsNonlinear(cons));
2321 consrhs = (cons == NULL ? SCIP_INVALID : SCIPgetRhsNonlinear(cons));
2322
2323 SCIP_CALL( detectSOC(scip, nlhdlrdata, expr, conslhs, consrhs, nlhdlrexprdata, &enforcebelow, &success) );
2324
2325 if( !success )
2326 return SCIP_OKAY;
2327
2328 /* inform what we can do */
2330
2331 /* if we have been successful on sqrt(...) <= auxvar, then we enforce
2332 * otherwise, expr is quadratic and we separate for expr <= ub(auxvar) only
2333 * 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)
2334 */
2335 if( (SCIPisExprPower(scip, expr) && SCIPgetExponentExprPow(expr) == 0.5) || (cons != NULL) )
2337
2338 return SCIP_OKAY;
2339}
2340
2341
2342/** auxiliary evaluation callback of nonlinear handler
2343 * @todo: remember if we are in the original variables and avoid reevaluating
2344 */
2345static
2347{ /*lint --e{715}*/
2348 int i;
2349
2350 assert(nlhdlrexprdata != NULL);
2351 assert(nlhdlrexprdata->vars != NULL);
2352 assert(nlhdlrexprdata->transcoefs != NULL);
2353 assert(nlhdlrexprdata->transcoefsidx != NULL);
2354 assert(nlhdlrexprdata->nterms > 1);
2355
2356 /* if the original expression is a norm, evaluate w.r.t. the auxiliary variables */
2357 if( SCIPisExprPower(scip, expr) )
2358 {
2359 assert(SCIPgetExponentExprPow(expr) == 0.5);
2360
2361 /* compute sum_i coef_i expr_i^2 */
2362 *auxvalue = 0.0;
2363 for( i = 0; i < nlhdlrexprdata->nterms - 1; ++i )
2364 {
2365 SCIP_Real termval;
2366
2367 termval = evalSingleTerm(scip, nlhdlrexprdata, sol, i);
2368 *auxvalue += SQR(termval);
2369 }
2370
2371 assert(*auxvalue >= 0.0);
2372
2373 /* compute SQRT(sum_i coef_i expr_i^2) */
2374 *auxvalue = SQRT(*auxvalue);
2375 }
2376 /* otherwise, evaluate the original quadratic expression w.r.t. the created auxvars of the children */
2377 else
2378 {
2379 SCIP_EXPR** children;
2380 SCIP_Real* childcoefs;
2381 int nchildren;
2382
2383 assert(SCIPisExprSum(scip, expr));
2384
2385 children = SCIPexprGetChildren(expr);
2387 nchildren = SCIPexprGetNChildren(expr);
2388
2389 *auxvalue = SCIPgetConstantExprSum(expr);
2390
2391 for( i = 0; i < nchildren; ++i )
2392 {
2393 if( SCIPisExprPower(scip, children[i]) )
2394 {
2396 SCIP_Real solval;
2397
2398 assert(SCIPgetExponentExprPow(children[i]) == 2.0);
2399
2401 assert(argauxvar != NULL);
2402
2403 solval = SCIPgetSolVal(scip, sol, argauxvar);
2404 *auxvalue += childcoefs[i] * SQR( solval );
2405 }
2406 else if( SCIPisExprProduct(scip, children[i]) )
2407 {
2410
2411 assert(SCIPexprGetNChildren(children[i]) == 2);
2412
2417
2419 }
2420 else
2421 {
2423
2425 assert(argauxvar != NULL);
2426
2427 *auxvalue += childcoefs[i] * SCIPgetSolVal(scip, sol, argauxvar);
2428 }
2429 }
2430 }
2431
2432 return SCIP_OKAY;
2433}
2434
2435
2436/** separation deinitialization method of a nonlinear handler (called during CONSINITLP) */
2437static
2439{ /*lint --e{715}*/
2440 assert(conshdlr != NULL);
2441 assert(expr != NULL);
2442 assert(nlhdlrexprdata != NULL);
2443
2444 /* if we have 3 or more terms in lhs create variable and row for disaggregation */
2445 if( nlhdlrexprdata->nterms > 3 )
2446 {
2447 /* create variables for cone disaggregation */
2448 SCIP_CALL( createDisaggrVars(scip, expr, nlhdlrexprdata) );
2449
2450#ifdef WITH_DEBUG_SOLUTION
2452 {
2453 SCIP_Real lhsval;
2454 SCIP_Real rhsval;
2455 SCIP_Real disvarval;
2456 int ndisvars;
2457 int nterms;
2458 int i;
2459 int k;
2460
2461 /* the debug solution value of the disaggregation variables is set to
2462 * (v_i^T x + beta_i)^2 / (v_{n+1}^T x + beta_{n+1})
2463 * if (v_{n+1}^T x + beta_{n+1}) is different from 0.
2464 * Otherwise, the debug solution value is set to 0.
2465 */
2466
2467 nterms = nlhdlrexprdata->nterms;
2468
2469 /* find value of rhs */
2470 rhsval = nlhdlrexprdata->offsets[nterms - 1];
2471 for( i = nlhdlrexprdata->termbegins[nterms - 1]; i < nlhdlrexprdata->termbegins[nterms]; ++i )
2472 {
2473 SCIP_VAR* var;
2474 SCIP_Real varval;
2475
2476 var = SCIPgetVarExprVar(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[i]]);
2477
2479 rhsval += nlhdlrexprdata->transcoefs[i] * varval;
2480 }
2481
2482 /* set value of disaggregation vars */
2483 ndisvars = nlhdlrexprdata->nterms - 1;
2484
2485 if( SCIPisZero(scip, rhsval) )
2486 {
2487 for( i = 0; i < ndisvars; ++i )
2488 {
2489 SCIP_CALL( SCIPdebugAddSolVal(scip, nlhdlrexprdata->disvars[i], 0.0) );
2490 }
2491 }
2492 else
2493 {
2494 /* set value for each disaggregation variable corresponding to quadratic term */
2495 for( k = 0; k < ndisvars; ++k )
2496 {
2497 lhsval = nlhdlrexprdata->offsets[k];
2498
2499 for( i = nlhdlrexprdata->termbegins[k]; i < nlhdlrexprdata->termbegins[k + 1]; ++i )
2500 {
2501 SCIP_VAR* var;
2502 SCIP_Real varval;
2503
2504 var = SCIPgetVarExprVar(nlhdlrexprdata->vars[nlhdlrexprdata->transcoefsidx[i]]);
2505
2507 lhsval += nlhdlrexprdata->transcoefs[i] * varval;
2508 }
2509
2511
2512 SCIP_CALL( SCIPdebugAddSolVal(scip, nlhdlrexprdata->disvars[k], disvarval) );
2513 }
2514 }
2515 }
2516#endif
2517
2518 /* create the disaggregation row and store it in nlhdlrexprdata */
2519 SCIP_CALL( createDisaggrRow(scip, conshdlr, expr, nlhdlrexprdata) );
2520 }
2521
2522 /* TODO add something to the LP as well, at least the disaggregation row */
2523
2524 return SCIP_OKAY;
2525}
2526
2527
2528/** separation deinitialization method of a nonlinear handler (called during CONSEXITSOL) */
2529static
2531{ /*lint --e{715}*/
2532 assert(nlhdlrexprdata != NULL);
2533
2534 /* free disaggreagation row */
2535 if( nlhdlrexprdata->disrow != NULL )
2536 {
2537 SCIP_CALL( SCIPreleaseRow(scip, &nlhdlrexprdata->disrow) );
2538 }
2539
2540 return SCIP_OKAY;
2541}
2542
2543
2544/** nonlinear handler separation callback */
2545static
2547{ /*lint --e{715}*/
2549 SCIP_Real rhsval;
2550 int ndisaggrs;
2551 int k;
2552 SCIP_Bool infeasible;
2553
2554 assert(nlhdlrexprdata != NULL);
2555 assert(nlhdlrexprdata->nterms < 4 || nlhdlrexprdata->disrow != NULL);
2556 assert(nlhdlrexprdata->nterms > 1);
2557
2559
2560 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
2562
2563 rhsval = evalSingleTerm(scip, nlhdlrexprdata, sol, nlhdlrexprdata->nterms - 1);
2564
2565 /* if there are three or two terms just compute gradient cut */
2566 if( nlhdlrexprdata->nterms < 4 )
2567 {
2568 SCIP_ROW* row;
2569
2570 /* compute gradient cut */
2571 SCIP_CALL( generateCutSolSOC(scip, expr, cons, sol, nlhdlrexprdata, SCIPgetLPFeastol(scip), rhsval, &row) );
2572
2573 /* TODO this code repeats below, factorize out */
2574 if( row != NULL )
2575 {
2576 SCIP_Real cutefficacy;
2577
2579
2580 SCIPdebugMsg(scip, "generated row for %d-SOC, efficacy=%g, minefficacy=%g, allowweakcuts=%u\n",
2581 nlhdlrexprdata->nterms, cutefficacy, nlhdlrdata->mincutefficacy, allowweakcuts);
2582
2583 /* check whether cut is applicable */
2584 if( SCIPisCutApplicable(scip, row) && (allowweakcuts || cutefficacy >= nlhdlrdata->mincutefficacy) )
2585 {
2586 SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
2587
2588#ifdef SCIP_CONSNONLINEAR_ROWNOTREMOVABLE
2589 /* mark row as not removable from LP for current node, if in enforcement (==addbranchscores) (this can prevent some cycling) */
2590 if( addbranchscores )
2592#endif
2593
2594 if( infeasible )
2596 else
2598 }
2599
2600 /* release row */
2601 SCIP_CALL( SCIPreleaseRow(scip, &row) );
2602 }
2603#ifdef SCIP_DEBUG
2604 else
2605 {
2606 SCIPdebugMsg(scip, "failed to generate SOC\n");
2607 }
2608#endif
2609
2610 return SCIP_OKAY;
2611 }
2612
2613 ndisaggrs = nlhdlrexprdata->nterms - 1;
2614
2615 /* check whether the aggregation row is in the LP */
2616 if( !SCIProwIsInLP(nlhdlrexprdata->disrow) && -SCIPgetRowSolFeasibility(scip, nlhdlrexprdata->disrow, sol) > SCIPgetLPFeastol(scip) )
2617 {
2618 SCIP_CALL( SCIPaddRow(scip, nlhdlrexprdata->disrow, TRUE, &infeasible) );
2619 SCIPdebugMsg(scip, "added disaggregation row to LP, cutoff=%u\n", infeasible);
2620
2621 if( infeasible )
2622 {
2624 return SCIP_OKAY;
2625 }
2626
2628 }
2629
2630 for( k = 0; k < ndisaggrs && *result != SCIP_CUTOFF; ++k )
2631 {
2632 SCIP_ROW* row;
2633
2634 /* compute gradient cut */
2635 SCIP_CALL( generateCutSolDisagg(scip, expr, cons, sol, nlhdlrexprdata, k, SCIPgetLPFeastol(scip), rhsval, &row) );
2636
2637 if( row != NULL )
2638 {
2639 SCIP_Real cutefficacy;
2640
2642
2643 SCIPdebugMsg(scip, "generated row for disaggregation %d, efficacy=%g, minefficacy=%g, allowweakcuts=%u\n",
2644 k, cutefficacy, nlhdlrdata->mincutefficacy, allowweakcuts);
2645
2646 /* check whether cut is applicable */
2647 if( SCIPisCutApplicable(scip, row) && (allowweakcuts || cutefficacy >= nlhdlrdata->mincutefficacy) )
2648 {
2649 SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
2650
2651#ifdef SCIP_CONSNONLINEAR_ROWNOTREMOVABLE
2652 /* mark row as not removable from LP for current node, if in enforcement (==addbranchscores) (this can prevent some cycling) */
2653 if( addbranchscores )
2655#endif
2656
2657 if( infeasible )
2659 else
2661 }
2662
2663 /* release row */
2664 SCIP_CALL( SCIPreleaseRow(scip, &row) );
2665 }
2666 }
2667
2668 return SCIP_OKAY;
2669}
2670
2671/*
2672 * nonlinear handler specific interface methods
2673 */
2674
2675/** includes SOC nonlinear handler in nonlinear constraint handler */
2677 SCIP* scip /**< SCIP data structure */
2678 )
2679{
2681 SCIP_NLHDLR* nlhdlr;
2682
2683 assert(scip != NULL);
2684
2685 /* create nonlinear handler data */
2687
2689 assert(nlhdlr != NULL);
2690
2696
2697 /* add soc nlhdlr parameters */
2698 /* TODO should we get rid of this and use separating/mineffiacy(root) instead, which is 1e-4? */
2699 SCIP_CALL( SCIPaddRealParam(scip, "nlhdlr/" NLHDLR_NAME "/mincutefficacy",
2700 "Minimum efficacy which a cut needs in order to be added.",
2701 &nlhdlrdata->mincutefficacy, FALSE, DEFAULT_MINCUTEFFICACY, 0.0, SCIPinfinity(scip), NULL, NULL) );
2702
2703 SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/compeigenvalues",
2704 "Should Eigenvalue computations be done to detect complex cases in quadratic constraints?",
2705 &nlhdlrdata->compeigenvalues, FALSE, DEFAULT_COMPEIGENVALUES, NULL, NULL) );
2706
2707 return SCIP_OKAY;
2708}
2709
2710/** checks whether constraint is SOC representable in original variables and returns the SOC representation
2711 *
2712 * The SOC representation has the form:
2713 * \f$\sqrt{\sum_{i=1}^{n} (v_i^T x + \beta_i)^2} - v_{n+1}^T x - \beta_{n+1} \lessgtr 0\f$,
2714 * where \f$n+1 = \text{nterms}\f$ and the inequality type is given by sidetype (`SCIP_SIDETYPE_RIGHT` if inequality
2715 * is \f$\leq\f$, `SCIP_SIDETYPE_LEFT` if \f$\geq\f$).
2716 *
2717 * 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
2718 * corresponding element `offsets[i-1]` and `termbegins[i-1]` is the starting position of the term in arrays
2719 * `transcoefs` and `transcoefsidx`. The overall number of nonzeros is `termbegins[nterms]`.
2720 *
2721 * Arrays `transcoefs` and `transcoefsidx` have size `termbegins[nterms]` and define the linear expressions \f$v_i^T x\f$
2722 * for each term. For a term \f$i\f$ in the above notation, the nonzeroes are given by elements
2723 * `termbegins[i-1]...termbegins[i]` of `transcoefs` and `transcoefsidx`. There may be no nonzeroes for some term (i.e.,
2724 * constant terms are possible). `transcoefs` contains the coefficients \f$v_i\f$ and `transcoefsidx` contains positions of
2725 * variables in the `vars` array.
2726 *
2727 * The `vars` array has size `nvars` and contains \f$x\f$ variables; each variable is included at most once.
2728 *
2729 * The arrays should be freed by calling SCIPfreeSOCArraysNonlinear().
2730 *
2731 * This function uses the methods that are used in the detection algorithm of the SOC nonlinear handler.
2732 */
2734 SCIP* scip, /**< SCIP data structure */
2735 SCIP_CONS* cons, /**< nonlinear constraint */
2736 SCIP_Bool compeigenvalues, /**< whether eigenvalues should be computed to detect complex cases */
2737 SCIP_Bool* success, /**< pointer to store whether SOC structure has been detected */
2738 SCIP_SIDETYPE* sidetype, /**< pointer to store which side of cons is SOC representable; only
2739 * valid when success is TRUE */
2740 SCIP_VAR*** vars, /**< variables (x) that appear on both sides; no duplicates are allowed */
2741 SCIP_Real** offsets, /**< offsets of both sides (beta_i) */
2742 SCIP_Real** transcoefs, /**< non-zeros of linear transformation vectors (v_i) */
2743 int** transcoefsidx, /**< mapping of transformation coefficients to variable indices in vars */
2744 int** termbegins, /**< starting indices of transcoefs for each term */
2745 int* nvars, /**< total number of variables appearing (i.e. size of vars) */
2746 int* nterms /**< number of summands in the SQRT +1 for RHS (n+1) */
2747 )
2748{
2750 SCIP_NLHDLREXPRDATA *nlhdlrexprdata;
2751 SCIP_Real conslhs;
2752 SCIP_Real consrhs;
2753 SCIP_EXPR* expr;
2754 SCIP_Bool enforcebelow;
2755 int i;
2756
2757 assert(cons != NULL);
2758
2759 expr = SCIPgetExprNonlinear(cons);
2760 assert(expr != NULL);
2761
2762 nlhdlrdata.mincutefficacy = 0.0;
2763 nlhdlrdata.compeigenvalues = compeigenvalues;
2764
2767
2768 SCIP_CALL( detectSOC(scip, &nlhdlrdata, expr, conslhs, consrhs, &nlhdlrexprdata, &enforcebelow, success) );
2769
2770 /* the constraint must be SOC representable in original variables */
2771 if( *success )
2772 {
2773 assert(nlhdlrexprdata != NULL);
2774
2775 for( i = 0; i < nlhdlrexprdata->nvars; ++i )
2776 {
2777 if( !SCIPisExprVar(scip, nlhdlrexprdata->vars[i]) )
2778 {
2779 *success = FALSE;
2780 break;
2781 }
2782 }
2783 }
2784
2785 if( *success )
2786 {
2788 SCIP_CALL( SCIPallocBlockMemoryArray(scip, vars, nlhdlrexprdata->nvars) );
2789
2790 for( i = 0; i < nlhdlrexprdata->nvars; ++i )
2791 {
2792 (*vars)[i] = SCIPgetVarExprVar(nlhdlrexprdata->vars[i]);
2793 assert((*vars)[i] != NULL);
2794 }
2795 SCIPfreeBlockMemoryArray(scip, &nlhdlrexprdata->vars, nlhdlrexprdata->nvars);
2796 *offsets = nlhdlrexprdata->offsets;
2797 *transcoefs = nlhdlrexprdata->transcoefs;
2798 *transcoefsidx = nlhdlrexprdata->transcoefsidx;
2799 *termbegins = nlhdlrexprdata->termbegins;
2800 *nvars = nlhdlrexprdata->nvars;
2801 *nterms = nlhdlrexprdata->nterms;
2802 SCIPfreeBlockMemory(scip, &nlhdlrexprdata);
2803 }
2804 else
2805 {
2806 if( nlhdlrexprdata != NULL )
2807 {
2808 SCIP_CALL( freeNlhdlrExprData(scip, &nlhdlrexprdata) );
2809 }
2810 *vars = NULL;
2811 *offsets = NULL;
2812 *transcoefs = NULL;
2813 *transcoefsidx = NULL;
2814 *termbegins = NULL;
2815 *nvars = 0;
2816 *nterms = 0;
2817 }
2818
2819 return SCIP_OKAY;
2820}
2821
2822/** frees arrays created by SCIPisSOCNonlinear() */
2824 SCIP* scip, /**< SCIP data structure */
2825 SCIP_VAR*** vars, /**< variables that appear on both sides (x) */
2826 SCIP_Real** offsets, /**< offsets of both sides (beta_i) */
2827 SCIP_Real** transcoefs, /**< non-zeros of linear transformation vectors (v_i) */
2828 int** transcoefsidx, /**< mapping of transformation coefficients to variable indices in vars */
2829 int** termbegins, /**< starting indices of transcoefs for each term */
2830 int nvars, /**< total number of variables appearing */
2831 int nterms /**< number of summands in the SQRT +1 for RHS (n+1) */
2832 )
2833{
2834 int ntranscoefs;
2835
2836 if( nvars == 0 )
2837 return;
2838
2839 assert(vars != NULL);
2840 assert(offsets != NULL);
2841 assert(transcoefs != NULL);
2842 assert(transcoefsidx != NULL);
2843 assert(termbegins != NULL);
2844
2845 ntranscoefs = (*termbegins)[nterms];
2846
2847 SCIPfreeBlockMemoryArray(scip, termbegins, nterms + 1);
2848 SCIPfreeBlockMemoryArray(scip, transcoefsidx, ntranscoefs);
2852}
constraint handler for nonlinear constraints specified by algebraic expressions
methods for debugging
#define SCIPdebugGetSolVal(scip, var, val)
Definition debug.h:299
#define SCIPdebugAddSolVal(scip, var, val)
Definition debug.h:298
#define SCIP_MAXSTRLEN
Definition def.h:302
#define SCIP_INVALID
Definition def.h:206
#define TRUE
Definition def.h:95
#define FALSE
Definition def.h:96
#define SCIPABORT()
Definition def.h:360
#define SCIP_CALL(x)
Definition def.h:388
power and signed power expression handlers
sum expression handler
variable expression handler
unsigned int SCIPgetExprNAuxvarUsesNonlinear(SCIP_EXPR *expr)
int SCIPgetExprNLocksPosNonlinear(SCIP_EXPR *expr)
SCIP_VAR * SCIPgetExprAuxVarNonlinear(SCIP_EXPR *expr)
SCIP_EXPR * SCIPgetExprNonlinear(SCIP_CONS *cons)
SCIP_Real SCIPgetRhsNonlinear(SCIP_CONS *cons)
int SCIPgetExprNLocksNegNonlinear(SCIP_EXPR *expr)
SCIP_RETCODE SCIPregisterExprUsageNonlinear(SCIP *scip, SCIP_EXPR *expr, SCIP_Bool useauxvar, SCIP_Bool useactivityforprop, SCIP_Bool useactivityforsepabelow, SCIP_Bool useactivityforsepaabove)
SCIP_Real SCIPgetLhsNonlinear(SCIP_CONS *cons)
SCIP_RETCODE SCIPaddVar(SCIP *scip, SCIP_VAR *var)
Definition scip_prob.c:1668
void SCIPhashmapFree(SCIP_HASHMAP **hashmap)
Definition misc.c:3058
int SCIPhashmapGetImageInt(SCIP_HASHMAP *hashmap, void *origin)
Definition misc.c:3231
int SCIPhashmapGetNElements(SCIP_HASHMAP *hashmap)
Definition misc.c:3483
SCIP_RETCODE SCIPhashmapCreate(SCIP_HASHMAP **hashmap, BMS_BLKMEM *blkmem, int mapsize)
Definition misc.c:3024
SCIP_Bool SCIPhashmapExists(SCIP_HASHMAP *hashmap, void *origin)
Definition misc.c:3373
SCIP_RETCODE SCIPhashmapInsertInt(SCIP_HASHMAP *hashmap, void *origin, int image)
Definition misc.c:3142
void SCIPhashsetFree(SCIP_HASHSET **hashset, BMS_BLKMEM *blkmem)
Definition misc.c:3740
SCIP_Bool SCIPhashsetExists(SCIP_HASHSET *hashset, void *element)
Definition misc.c:3767
int SCIPhashsetGetNElements(SCIP_HASHSET *hashset)
Definition misc.c:3942
SCIP_RETCODE SCIPhashsetInsert(SCIP_HASHSET *hashset, BMS_BLKMEM *blkmem, void *element)
Definition misc.c:3750
SCIP_RETCODE SCIPhashsetCreate(SCIP_HASHSET **hashset, BMS_BLKMEM *blkmem, int size)
Definition misc.c:3709
SCIP_RETCODE SCIPhashsetRemove(SCIP_HASHSET *hashset, void *element)
Definition misc.c:3808
void SCIPinfoMessage(SCIP *scip, FILE *file, const char *formatstr,...)
void SCIPverbMessage(SCIP *scip, SCIP_VERBLEVEL msgverblevel, FILE *file, const char *formatstr,...)
#define SCIPdebugMsg
void SCIPfreeSOCArraysNonlinear(SCIP *scip, SCIP_VAR ***vars, SCIP_Real **offsets, SCIP_Real **transcoefs, int **transcoefsidx, int **termbegins, int nvars, int nterms)
SCIP_RETCODE SCIPisSOCNonlinear(SCIP *scip, SCIP_CONS *cons, SCIP_Bool compeigenvalues, SCIP_Bool *success, SCIP_SIDETYPE *sidetype, SCIP_VAR ***vars, SCIP_Real **offsets, SCIP_Real **transcoefs, int **transcoefsidx, int **termbegins, int *nvars, int *nterms)
SCIP_RETCODE SCIPcallLapackDsyevIpopt(SCIP_Bool computeeigenvectors, int N, SCIP_Real *a, SCIP_Real *w)
SCIP_Bool SCIPisIpoptAvailableIpopt(void)
SCIP_RETCODE SCIPincludeNlhdlrSoc(SCIP *scip)
SCIP_RETCODE SCIPaddRealParam(SCIP *scip, const char *name, const char *desc, SCIP_Real *valueptr, SCIP_Bool isadvanced, SCIP_Real defaultvalue, SCIP_Real minvalue, SCIP_Real maxvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition scip_param.c:139
SCIP_RETCODE SCIPaddBoolParam(SCIP *scip, const char *name, const char *desc, SCIP_Bool *valueptr, SCIP_Bool isadvanced, SCIP_Bool defaultvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition scip_param.c:57
void SCIPswapReals(SCIP_Real *value1, SCIP_Real *value2)
Definition misc.c:10294
SCIP_Real SCIPgetCutEfficacy(SCIP *scip, SCIP_SOL *sol, SCIP_ROW *cut)
Definition scip_cut.c:94
SCIP_Bool SCIPisCutApplicable(SCIP *scip, SCIP_ROW *cut)
Definition scip_cut.c:207
SCIP_RETCODE SCIPaddRow(SCIP *scip, SCIP_ROW *row, SCIP_Bool forcecut, SCIP_Bool *infeasible)
Definition scip_cut.c:250
int SCIPexprGetNChildren(SCIP_EXPR *expr)
Definition expr.c:3801
SCIP_Real SCIPgetExponentExprPow(SCIP_EXPR *expr)
Definition expr_pow.c:3351
SCIP_Bool SCIPisExprProduct(SCIP *scip, SCIP_EXPR *expr)
Definition scip_expr.c:1454
SCIP_Bool SCIPisExprSum(SCIP *scip, SCIP_EXPR *expr)
Definition scip_expr.c:1443
SCIP_Real * SCIPgetCoefsExprSum(SCIP_EXPR *expr)
Definition expr_sum.c:1166
SCIP_Bool SCIPisExprVar(SCIP *scip, SCIP_EXPR *expr)
Definition scip_expr.c:1421
SCIP_RETCODE SCIPprintExpr(SCIP *scip, SCIP_EXPR *expr, FILE *file)
Definition scip_expr.c:1476
SCIP_Bool SCIPisExprPower(SCIP *scip, SCIP_EXPR *expr)
Definition scip_expr.c:1465
SCIP_EXPR ** SCIPexprGetChildren(SCIP_EXPR *expr)
Definition expr.c:3811
SCIP_Real SCIPgetConstantExprSum(SCIP_EXPR *expr)
Definition expr_sum.c:1181
SCIP_VAR * SCIPgetVarExprVar(SCIP_EXPR *expr)
Definition expr_var.c:416
SCIP_INTERVAL SCIPexprGetActivity(SCIP_EXPR *expr)
Definition expr.c:3957
SCIP_RETCODE SCIPdismantleExpr(SCIP *scip, FILE *file, SCIP_EXPR *expr)
Definition scip_expr.c:1598
SCIP_RETCODE SCIPevalExprActivity(SCIP *scip, SCIP_EXPR *expr)
Definition scip_expr.c:1707
SCIP_Real SCIPgetLPFeastol(SCIP *scip)
Definition scip_lp.c:428
#define SCIPfreeBlockMemoryArray(scip, ptr, num)
Definition scip_mem.h:110
#define SCIPallocClearBlockMemory(scip, ptr)
Definition scip_mem.h:91
#define SCIPallocClearBufferArray(scip, ptr, num)
Definition scip_mem.h:126
#define SCIPallocBufferArray(scip, ptr, num)
Definition scip_mem.h:124
#define SCIPfreeBufferArray(scip, ptr)
Definition scip_mem.h:136
#define SCIPduplicateBufferArray(scip, ptr, source, num)
Definition scip_mem.h:132
#define SCIPallocBlockMemoryArray(scip, ptr, num)
Definition scip_mem.h:93
#define SCIPfreeBlockMemory(scip, ptr)
Definition scip_mem.h:108
#define SCIPfreeBlockMemoryArrayNull(scip, ptr, num)
Definition scip_mem.h:111
#define SCIPfreeBufferArrayNull(scip, ptr)
Definition scip_mem.h:137
#define SCIPallocBlockMemory(scip, ptr)
Definition scip_mem.h:89
#define SCIPduplicateBlockMemoryArray(scip, ptr, source, num)
Definition scip_mem.h:105
void SCIPnlhdlrSetInitExit(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRINIT((*init)),)
Definition nlhdlr.c:106
SCIP_NLHDLRDATA * SCIPnlhdlrGetData(SCIP_NLHDLR *nlhdlr)
Definition nlhdlr.c:200
void SCIPnlhdlrSetFreeExprData(SCIP_NLHDLR *nlhdlr,)
Definition nlhdlr.c:94
const char * SCIPnlhdlrGetName(SCIP_NLHDLR *nlhdlr)
Definition nlhdlr.c:150
void SCIPnlhdlrSetSepa(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRINITSEPA((*initsepa)), SCIP_DECL_NLHDLRENFO((*enfo)), SCIP_DECL_NLHDLRESTIMATE((*estimate)),)
Definition nlhdlr.c:132
void SCIPnlhdlrSetFreeHdlrData(SCIP_NLHDLR *nlhdlr,)
Definition nlhdlr.c:83
void SCIPnlhdlrSetCopyHdlr(SCIP_NLHDLR *nlhdlr,)
Definition nlhdlr.c:72
SCIP_RETCODE SCIPincludeNlhdlrNonlinear(SCIP *scip, SCIP_NLHDLR **nlhdlr, const char *name, const char *desc, int detectpriority, int enfopriority, SCIP_DECL_NLHDLRDETECT((*detect)), SCIP_DECL_NLHDLREVALAUX((*evalaux)), SCIP_NLHDLRDATA *nlhdlrdata)
SCIP_RETCODE SCIPcreateEmptyRowConshdlr(SCIP *scip, SCIP_ROW **row, SCIP_CONSHDLR *conshdlr, const char *name, SCIP_Real lhs, SCIP_Real rhs, SCIP_Bool local, SCIP_Bool modifiable, SCIP_Bool removable)
Definition scip_lp.c:1391
SCIP_RETCODE SCIPaddVarToRow(SCIP *scip, SCIP_ROW *row, SCIP_VAR *var, SCIP_Real val)
Definition scip_lp.c:1701
SCIP_Real SCIPgetRowSolFeasibility(SCIP *scip, SCIP_ROW *row, SCIP_SOL *sol)
Definition scip_lp.c:2167
SCIP_RETCODE SCIPreleaseRow(SCIP *scip, SCIP_ROW **row)
Definition scip_lp.c:1562
void SCIPmarkRowNotRemovableLocal(SCIP *scip, SCIP_ROW *row)
Definition scip_lp.c:1868
SCIP_Bool SCIProwIsInLP(SCIP_ROW *row)
Definition lp.c:17523
SCIP_Real SCIPgetSolVal(SCIP *scip, SCIP_SOL *sol, SCIP_VAR *var)
Definition scip_sol.c:1361
SCIP_Longint SCIPgetNLPs(SCIP *scip)
SCIP_Real SCIPinfinity(SCIP *scip)
SCIP_Bool SCIPisPositive(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisInfinity(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisNegative(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPvarIsBinary(SCIP_VAR *var)
Definition var.c:17421
void SCIPvarMarkRelaxationOnly(SCIP_VAR *var)
Definition var.c:17546
SCIP_RETCODE SCIPaddVarLocksType(SCIP *scip, SCIP_VAR *var, SCIP_LOCKTYPE locktype, int nlocksdown, int nlocksup)
Definition scip_var.c:4259
const char * SCIPvarGetName(SCIP_VAR *var)
Definition var.c:17241
SCIP_RETCODE SCIPreleaseVar(SCIP *scip, SCIP_VAR **var)
Definition scip_var.c:1248
SCIP_RETCODE SCIPcreateVarBasic(SCIP *scip, SCIP_VAR **var, const char *name, SCIP_Real lb, SCIP_Real ub, SCIP_Real obj, SCIP_VARTYPE vartype)
Definition scip_var.c:194
SCIP_RETCODE SCIPcleanupRowprep2(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_SOL *sol, SCIP_Real maxcoefbound, SCIP_Bool *success)
SCIP_RETCODE SCIPensureRowprepSize(SCIP *scip, SCIP_ROWPREP *rowprep, int size)
SCIP_Real SCIPgetRowprepViolation(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_SOL *sol, SCIP_Bool *reliable)
char * SCIProwprepGetName(SCIP_ROWPREP *rowprep)
SCIP_RETCODE SCIPaddRowprepTerm(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_VAR *var, SCIP_Real coef)
SCIP_RETCODE SCIPgetRowprepRowCons(SCIP *scip, SCIP_ROW **row, SCIP_ROWPREP *rowprep, SCIP_CONS *cons)
SCIP_RETCODE SCIPcreateRowprep(SCIP *scip, SCIP_ROWPREP **rowprep, SCIP_SIDETYPE sidetype, SCIP_Bool local)
void SCIProwprepAddSide(SCIP_ROWPREP *rowprep, SCIP_Real side)
void SCIPfreeRowprep(SCIP *scip, SCIP_ROWPREP **rowprep)
int SCIPsnprintf(char *t, int len, const char *s,...)
Definition misc.c:10788
return SCIP_OKAY
static SCIP_SOL * sol
assert(minobj< SCIPgetCutoffbound(scip))
int nvars
SCIP_VAR * var
static SCIP_VAR ** vars
static volatile int nterms
Definition interrupt.c:47
#define NULL
Definition lpi_spx1.cpp:161
#define NLHDLR_DETECTPRIORITY
Definition nlhdlr_soc.c:58
static SCIP_RETCODE generateCutSolDisagg(SCIP *scip, SCIP_EXPR *expr, SCIP_CONS *cons, SCIP_SOL *sol, SCIP_NLHDLREXPRDATA *nlhdlrexprdata, int disaggidx, SCIP_Real mincutviolation, SCIP_Real rhsval, SCIP_ROW **cut)
Definition nlhdlr_soc.c:586
#define nlhdlrInitSoc
static SCIP_RETCODE tryFillNlhdlrExprDataQuad(SCIP *scip, SCIP_EXPR **occurringexprs, SCIP_Real *eigvecmatrix, SCIP_Real *eigvals, SCIP_Real *bp, int nvars, int *termbegins, SCIP_Real *transcoefs, int *transcoefsidx, SCIP_Real *offsets, SCIP_Real *lhsconstant, int *nterms, SCIP_Bool *success)
Definition nlhdlr_soc.c:930
#define NLHDLR_ENFOPRIORITY
Definition nlhdlr_soc.c:59
#define nlhdlrExitSoc
static SCIP_RETCODE freeNlhdlrExprData(SCIP *scip, SCIP_NLHDLREXPRDATA **nlhdlrexprdata)
Definition nlhdlr_soc.c:377
static SCIP_Real evalSingleTerm(SCIP *scip, SCIP_NLHDLREXPRDATA *nlhdlrexprdata, SCIP_SOL *sol, int k)
Definition nlhdlr_soc.c:404
static SCIP_RETCODE createDisaggrRow(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_EXPR *expr, SCIP_NLHDLREXPRDATA *nlhdlrexprdata)
Definition nlhdlr_soc.c:278
#define NLHDLR_DESC
Definition nlhdlr_soc.c:57
static SCIP_RETCODE detectSocNorm(SCIP *scip, SCIP_EXPR *expr, SCIP_NLHDLREXPRDATA **nlhdlrexprdata, SCIP_Bool *success)
static SCIP_RETCODE generateCutSolSOC(SCIP *scip, SCIP_EXPR *expr, SCIP_CONS *cons, SCIP_SOL *sol, SCIP_NLHDLREXPRDATA *nlhdlrexprdata, SCIP_Real mincutviolation, SCIP_Real rhsval, SCIP_ROW **cut)
Definition nlhdlr_soc.c:450
static SCIP_RETCODE detectSOC(SCIP *scip, SCIP_NLHDLRDATA *nlhdlrdata, SCIP_EXPR *expr, SCIP_Real conslhs, SCIP_Real consrhs, SCIP_NLHDLREXPRDATA **nlhdlrexprdata, SCIP_Bool *enforcebelow, SCIP_Bool *success)
#define NLHDLR_NAME
Definition nlhdlr_soc.c:56
static SCIP_RETCODE checkAndCollectQuadratic(SCIP *scip, SCIP_EXPR *quadexpr, SCIP_HASHMAP *expr2idx, SCIP_EXPR **occurringexprs, int *nexprs, SCIP_Bool *success)
Definition nlhdlr_soc.c:734
#define DEFAULT_MINCUTEFFICACY
Definition nlhdlr_soc.c:60
#define DEFAULT_COMPEIGENVALUES
Definition nlhdlr_soc.c:61
static SCIP_RETCODE createDisaggrVars(SCIP *scip, SCIP_EXPR *expr, SCIP_NLHDLREXPRDATA *nlhdlrexprdata)
Definition nlhdlr_soc.c:214
static SCIP_RETCODE detectSocQuadraticComplex(SCIP *scip, SCIP_EXPR *expr, SCIP_Real conslhs, SCIP_Real consrhs, SCIP_NLHDLREXPRDATA **nlhdlrexprdata, SCIP_Bool *enforcebelow, SCIP_Bool *success)
static SCIP_RETCODE createNlhdlrExprData(SCIP *scip, SCIP_EXPR **vars, SCIP_Real *offsets, SCIP_Real *transcoefs, int *transcoefsidx, int *termbegins, int nvars, int nterms, SCIP_NLHDLREXPRDATA **nlhdlrexprdata)
Definition nlhdlr_soc.c:331
static void buildQuadExprMatrix(SCIP *scip, SCIP_EXPR *quadexpr, SCIP_HASHMAP *expr2idx, int nexprs, SCIP_Real *quadmatrix, SCIP_Real *linvector)
Definition nlhdlr_soc.c:849
static SCIP_RETCODE freeDisaggrVars(SCIP *scip, SCIP_NLHDLREXPRDATA *nlhdlrexprdata)
Definition nlhdlr_soc.c:248
static SCIP_RETCODE detectSocQuadraticSimple(SCIP *scip, SCIP_EXPR *expr, SCIP_Real conslhs, SCIP_Real consrhs, SCIP_NLHDLREXPRDATA **nlhdlrexprdata, SCIP_Bool *enforcebelow, SCIP_Bool *success)
soc nonlinear handler
BMS_BLKMEM * SCIPblkmem(SCIP *scip)
Definition scip_mem.c:57
Ipopt NLP interface.
#define SCIPerrorMessage
Definition pub_message.h:64
public functions of nonlinear handlers of nonlinear constraints
SCIP_Real sup
SCIP_Real inf
#define MAX(x, y)
Definition tclique_def.h:92
@ SCIP_SIDETYPE_RIGHT
Definition type_lp.h:65
@ SCIP_SIDETYPE_LEFT
Definition type_lp.h:64
enum SCIP_SideType SCIP_SIDETYPE
Definition type_lp.h:67
@ SCIP_VERBLEVEL_FULL
#define SCIP_NLHDLR_METHOD_SEPAABOVE
Definition type_nlhdlr.h:52
#define SCIP_DECL_NLHDLREVALAUX(x)
struct SCIP_NlhdlrData SCIP_NLHDLRDATA
#define SCIP_NLHDLR_METHOD_SEPABOTH
Definition type_nlhdlr.h:53
#define SCIP_DECL_NLHDLRCOPYHDLR(x)
Definition type_nlhdlr.h:70
#define SCIP_DECL_NLHDLRINIT(x)
#define SCIP_DECL_NLHDLREXIT(x)
#define SCIP_DECL_NLHDLRFREEEXPRDATA(x)
Definition type_nlhdlr.h:94
#define SCIP_DECL_NLHDLRDETECT(x)
#define SCIP_DECL_NLHDLREXITSEPA(x)
#define SCIP_DECL_NLHDLRINITSEPA(x)
#define SCIP_DECL_NLHDLRFREEHDLRDATA(x)
Definition type_nlhdlr.h:82
struct SCIP_NlhdlrExprData SCIP_NLHDLREXPRDATA
#define SCIP_DECL_NLHDLRENFO(x)
#define SCIP_NLHDLR_METHOD_SEPABELOW
Definition type_nlhdlr.h:51
@ SCIP_CUTOFF
Definition type_result.h:48
@ SCIP_DIDNOTFIND
Definition type_result.h:44
@ SCIP_SEPARATED
Definition type_result.h:49
enum SCIP_Retcode SCIP_RETCODE
@ SCIP_VARTYPE_CONTINUOUS
Definition type_var.h:71
@ SCIP_LOCKTYPE_MODEL
Definition type_var.h:97