Scippy

SCIP

Solving Constraint Integer Programs

cons_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-2021 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file cons_soc.c
17  * @ingroup DEFPLUGINS_CONS
18  * @brief constraint handler for second order cone constraints \f$\sqrt{\gamma + \sum_{i=1}^{n} (\alpha_i\, (x_i + \beta_i))^2} \leq \alpha_{n+1}\, (x_{n+1}+\beta_{n+1})\f$
19  * @author Stefan Vigerske
20  * @author Marc Pfetsch
21  * @author Felipe Serrano
22  * @author Benjamin Mueller
23  *
24  * @todo rhsvar == NULL is supported in some routines, but not everywhere
25  * @todo merge square terms with same variables in presol/exitpre
26  */
27 
28 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
29 
30 #define _USE_MATH_DEFINES /* to get M_PI on Windows */
31 #define SCIP_PRIVATE_ROWPREP
32 
33 #include "blockmemshell/memory.h"
34 #include <ctype.h>
35 #include "nlpi/exprinterpret.h"
36 #include "nlpi/nlpi.h"
37 #include "nlpi/nlpi_ipopt.h"
38 #include "nlpi/pub_expr.h"
39 #include "nlpi/type_expr.h"
40 #include "scip/cons_linear.h"
41 #include "scip/cons_quadratic.h"
42 #include "scip/cons_soc.h"
43 #include "scip/heur_subnlp.h"
44 #include "scip/heur_trysol.h"
45 #include "scip/intervalarith.h"
46 #include "scip/pub_cons.h"
47 #include "scip/pub_event.h"
48 #include "scip/pub_heur.h"
49 #include "scip/pub_message.h"
50 #include "scip/pub_misc.h"
51 #include "scip/pub_misc_sort.h"
52 #include "scip/pub_nlp.h"
53 #include "scip/pub_sol.h"
54 #include "scip/pub_tree.h"
55 #include "scip/pub_var.h"
56 #include "scip/scip_branch.h"
57 #include "scip/scip_cons.h"
58 #include "scip/scip_copy.h"
59 #include "scip/scip_cut.h"
60 #include "scip/scip_event.h"
61 #include "scip/scip_general.h"
62 #include "scip/scip_heur.h"
63 #include "scip/scip_lp.h"
64 #include "scip/scip_mem.h"
65 #include "scip/scip_message.h"
66 #include "scip/scip_nlp.h"
67 #include "scip/scip_numerics.h"
68 #include "scip/scip_param.h"
69 #include "scip/scip_prob.h"
70 #include "scip/scip_sepa.h"
71 #include "scip/scip_sol.h"
72 #include "scip/scip_solvingstats.h"
73 #include "scip/scip_tree.h"
74 #include "scip/scip_var.h"
75 #include <string.h>
76 
77 /* constraint handler properties */
78 #define CONSHDLR_NAME "soc"
79 #define CONSHDLR_DESC "constraint handler for second order cone constraints"
80 #define CONSHDLR_SEPAPRIORITY 10 /**< priority of the constraint handler for separation */
81 #define CONSHDLR_ENFOPRIORITY -40 /**< priority of the constraint handler for constraint enforcing */
82 #define CONSHDLR_CHECKPRIORITY -10 /**< priority of the constraint handler for checking feasibility */
83 #define CONSHDLR_SEPAFREQ 1 /**< frequency for separating cuts; zero means to separate only in the root node */
84 #define CONSHDLR_PROPFREQ 1 /**< frequency for propagating domains; zero means only preprocessing propagation */
85 #define CONSHDLR_EAGERFREQ 100 /**< frequency for using all instead of only the useful constraints in separation,
86  * propagation and enforcement, -1 for no eager evaluations, 0 for first only */
87 #define CONSHDLR_MAXPREROUNDS -1 /**< maximal number of presolving rounds the constraint handler participates in (-1: no limit) */
88 #define CONSHDLR_DELAYSEPA FALSE /**< should separation method be delayed, if other separators found cuts? */
89 #define CONSHDLR_DELAYPROP FALSE /**< should propagation method be delayed, if other propagators found reductions? */
90 #define CONSHDLR_NEEDSCONS TRUE /**< should the constraint handler be skipped, if no constraints are available? */
91 
92 #define CONSHDLR_PROP_TIMING SCIP_PROPTIMING_BEFORELP /**< propagation timing mask of the constraint handler */
93 #define CONSHDLR_PRESOLTIMING SCIP_PRESOLTIMING_ALWAYS /**< presolving timing of the constraint handler (fast, medium, or exhaustive) */
94 
95 #define QUADCONSUPGD_PRIORITY 60000 /**< priority of the constraint handler for upgrading of quadratic constraints */
96 
97 #define UPGSCALE 10 /* scale factor used in general upgrades of quadratic cons to soc */
98 
99 /*
100  * Data structures
101  */
102 
103 /** Eventdata for variable bound change events. */
104 struct VarEventData
105 {
106  SCIP_CONS* cons; /**< the constraint */
107  int varidx; /**< the index of a variable on the left hand side which bound change is caught, or -1 for variable on right hand side */
108  int filterpos; /**< position of corresponding event in event filter */
109 };
110 typedef struct VarEventData VAREVENTDATA;
112 /** constraint data for soc constraints */
113 struct SCIP_ConsData
114 {
115  int nvars; /**< number of variables on left hand side (n) */
116  SCIP_VAR** vars; /**< variables on left hand side (x_i) */
117  SCIP_Real* coefs; /**< coefficients for variables on left hand side (alpha_i) */
118  SCIP_Real* offsets; /**< offsets for variables on left hand side (beta_i) */
119  SCIP_Real constant; /**< constant on left hand side (gamma) */
120 
121  SCIP_VAR* rhsvar; /**< variable on right hand side (x_{n+1}) */
122  SCIP_Real rhscoeff; /**< coefficient of square term on right hand side (alpha_{n+1}) */
123  SCIP_Real rhsoffset; /**< offset for variable on right hand side (beta_{n+1}) */
124 
125  SCIP_NLROW* nlrow; /**< nonlinear row representation of constraint */
126 
127  SCIP_Real lhsval; /**< value of left hand side in current point */
128  SCIP_Real violation; /**< violation of constraint in current point */
129 
130  VAREVENTDATA* lhsbndchgeventdata; /**< eventdata for bound change events on left hand side variables */
131  VAREVENTDATA rhsbndchgeventdata; /**< eventdata for bound change event on right hand side variable */
132  SCIP_Bool isapproxadded; /**< has a linear outer approximation be added? */
133 };
134 
135 /** constraint handler data */
136 struct SCIP_ConshdlrData
137 {
138  SCIP_HEUR* subnlpheur; /**< a pointer to the subNLP heuristic, if available */
139  SCIP_HEUR* trysolheur; /**< a pointer to the trysol heuristic, if available */
140  SCIP_EVENTHDLR* eventhdlr; /**< event handler for bound change events */
141  int newsoleventfilterpos;/**< filter position of new solution event handler, if caught */
142  SCIP_Bool haveexprint; /**< indicates whether an expression interpreter is available */
143  SCIP_Bool sepanlp; /**< where linearization of the NLP relaxation solution added? */
144 
145  SCIP_Bool glineur; /**< is the Glineur outer approx preferred to Ben-Tal Nemirovski? */
146  SCIP_Bool projectpoint; /**< is the point in which a cut is generated projected onto the feasible set? */
147  int nauxvars; /**< number of auxiliary variables to use when creating a linear outer approx. of a SOC3 constraint */
148  SCIP_Bool sparsify; /**< whether to sparsify cuts */
149  SCIP_Real sparsifymaxloss; /**< maximal loss in cut efficacy by sparsification */
150  SCIP_Real sparsifynzgrowth; /**< growth rate of maximal allowed nonzeros in cuts in sparsification */
151  SCIP_Bool linfeasshift; /**< whether to try to make solutions feasible in check by shifting the variable on the right hand side */
152  char nlpform; /**< formulation of SOC constraint in NLP */
153  SCIP_Real sepanlpmincont; /**< minimal required fraction of continuous variables in problem to use solution of NLP relaxation in root for separation */
154  SCIP_Bool enfocutsremovable; /**< are cuts added during enforcement removable from the LP in the same node? */
155  SCIP_Bool generalsocupg; /**< try to upgrade more general quadratics to soc? */
156  SCIP_Bool disaggregate; /**< try to completely disaggregate soc? */
157 
158  SCIP_NODE* lastenfonode; /**< the node for which enforcement was called the last time (and some constraint was violated) */
159  int nenforounds; /**< counter on number of enforcement rounds for the current node */
160 };
161 
162 
163 /*
164  * Local methods
165  */
166 
167 /** catch left hand side variable events */
168 static
170  SCIP* scip, /**< SCIP data structure */
171  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
172  SCIP_CONS* cons, /**< constraint for which to catch bound change events */
173  int varidx /**< index of the variable which events to catch */
174  )
175 {
176  SCIP_CONSDATA* consdata;
177 
178  assert(scip != NULL);
179  assert(cons != NULL);
180  assert(eventhdlr != NULL);
181 
182  consdata = SCIPconsGetData(cons);
183  assert(consdata != NULL);
184  assert(varidx >= 0);
185  assert(varidx < consdata->nvars);
186  assert(consdata->lhsbndchgeventdata != NULL);
187 
188  consdata->lhsbndchgeventdata[varidx].cons = cons;
189  consdata->lhsbndchgeventdata[varidx].varidx = varidx;
190  SCIP_CALL( SCIPcatchVarEvent(scip, consdata->vars[varidx], SCIP_EVENTTYPE_BOUNDTIGHTENED, eventhdlr,
191  (SCIP_EVENTDATA*)&consdata->lhsbndchgeventdata[varidx], &consdata->lhsbndchgeventdata[varidx].filterpos) );
192 
193  /* since bound changes were not catched before, a possibly stored activity may have become outdated */
194  SCIP_CALL( SCIPmarkConsPropagate(scip, cons) );
195 
196  return SCIP_OKAY;
197 }
198 
199 /** catch right hand side variable events */
200 static
202  SCIP* scip, /**< SCIP data structure */
203  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
204  SCIP_CONS* cons /**< constraint for which to catch bound change events */
205  )
206 {
207  SCIP_CONSDATA* consdata;
208 
209  assert(scip != NULL);
210  assert(cons != NULL);
211  assert(eventhdlr != NULL);
212 
213  consdata = SCIPconsGetData(cons);
214  assert(consdata != NULL);
215 
216  consdata->rhsbndchgeventdata.cons = cons;
217  consdata->rhsbndchgeventdata.varidx = -1;
218  SCIP_CALL( SCIPcatchVarEvent(scip, consdata->rhsvar, SCIP_EVENTTYPE_UBTIGHTENED, eventhdlr,
219  (SCIP_EVENTDATA*)&consdata->rhsbndchgeventdata, &consdata->rhsbndchgeventdata.filterpos) );
220 
221  /* since bound changes were not catched before, a possibly stored activity may have become outdated */
222  SCIP_CALL( SCIPmarkConsPropagate(scip, cons) );
223 
224  return SCIP_OKAY;
225 }
226 
227 /** catch variables events */
228 static
230  SCIP* scip, /**< SCIP data structure */
231  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
232  SCIP_CONS* cons /**< constraint for which to catch bound change events */
233  )
234 {
235  SCIP_CONSDATA* consdata;
236  int i;
237 
238  assert(scip != NULL);
239  assert(cons != NULL);
240  assert(eventhdlr != NULL);
241 
242  consdata = SCIPconsGetData(cons);
243  assert(consdata != NULL);
244  assert(consdata->lhsbndchgeventdata == NULL);
245 
246  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->lhsbndchgeventdata, consdata->nvars) );
247 
248  for( i = 0; i < consdata->nvars; ++i )
249  {
250  if( consdata->vars[i] != NULL )
251  {
252  SCIP_CALL( catchLhsVarEvents(scip, eventhdlr, cons, i) );
253  }
254  }
255 
256  if( consdata->rhsvar != NULL )
257  {
258  SCIP_CALL( catchRhsVarEvents(scip, eventhdlr, cons) );
259  }
260 
261  return SCIP_OKAY;
262 }
263 
264 /** drop left hand side variable events */
265 static
267  SCIP* scip, /**< SCIP data structure */
268  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
269  SCIP_CONS* cons, /**< constraint for which to catch bound change events */
270  int varidx /**< index of the variable which events to catch */
271  )
272 {
273  SCIP_CONSDATA* consdata;
274 
275  assert(scip != NULL);
276  assert(cons != NULL);
277  assert(eventhdlr != NULL);
278 
279  consdata = SCIPconsGetData(cons);
280  assert(consdata != NULL);
281  assert(varidx >= 0);
282  assert(varidx < consdata->nvars);
283  assert(consdata->lhsbndchgeventdata != NULL);
284  assert(consdata->lhsbndchgeventdata[varidx].varidx == varidx);
285 
286  SCIP_CALL( SCIPdropVarEvent(scip, consdata->vars[varidx], SCIP_EVENTTYPE_BOUNDTIGHTENED, eventhdlr,
287  (SCIP_EVENTDATA*)&consdata->lhsbndchgeventdata[varidx], consdata->lhsbndchgeventdata[varidx].filterpos) );
288 
289  return SCIP_OKAY;
290 }
291 
292 /** drop right hand side variable events */
293 static
295  SCIP* scip, /**< SCIP data structure */
296  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
297  SCIP_CONS* cons /**< constraint for which to catch bound change events */
298  )
299 {
300  SCIP_CONSDATA* consdata;
301 
302  assert(scip != NULL);
303  assert(cons != NULL);
304  assert(eventhdlr != NULL);
305 
306  consdata = SCIPconsGetData(cons);
307  assert(consdata != NULL);
308  assert(consdata->rhsbndchgeventdata.varidx == -1);
309 
310  SCIP_CALL( SCIPdropVarEvent(scip, consdata->rhsvar, SCIP_EVENTTYPE_UBTIGHTENED, eventhdlr,
311  (SCIP_EVENTDATA*)&consdata->rhsbndchgeventdata, consdata->rhsbndchgeventdata.filterpos) );
312 
313  return SCIP_OKAY;
314 }
315 
316 /** drop variable events */
317 static
319  SCIP* scip, /**< SCIP data structure */
320  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
321  SCIP_CONS* cons /**< constraint for which to catch bound change events */
322  )
323 {
324  SCIP_CONSDATA* consdata;
325  int i;
326 
327  assert(scip != NULL);
328  assert(eventhdlr != NULL);
329  assert(cons != NULL);
330 
331  consdata = SCIPconsGetData(cons);
332  assert(consdata != NULL);
333 
334  for( i = 0; i < consdata->nvars; ++i )
335  {
336  if( consdata->vars[i] != NULL )
337  {
338  SCIP_CALL( dropLhsVarEvents(scip, eventhdlr, cons, i) );
339  }
340  }
341 
342  SCIPfreeBlockMemoryArray(scip, &consdata->lhsbndchgeventdata, consdata->nvars);
343 
344  if( consdata->rhsvar != NULL )
345  {
346  SCIP_CALL( dropRhsVarEvents(scip, eventhdlr, cons) );
347  }
348 
349  return SCIP_OKAY;
350 }
351 
352 /** process variable bound tightening event */
353 static
354 SCIP_DECL_EVENTEXEC(processVarEvent)
355 {
356  SCIP_CONS* cons;
357 
358  assert(scip != NULL);
359  assert(event != NULL);
360  assert(eventdata != NULL);
361  assert(eventhdlr != NULL);
362 
363  cons = ((VAREVENTDATA*)eventdata)->cons;
364  assert(cons != NULL);
365 
366  SCIP_CALL( SCIPmarkConsPropagate(scip, cons) );
367  /* @todo look at bounds on x_i to decide whether propagation makes sense */
368 
369  return SCIP_OKAY;
370 }
371 
372 /** create a nonlinear row representation of the constraint and stores them in consdata */
373 static
375  SCIP* scip, /**< SCIP data structure */
376  SCIP_CONSHDLR* conshdlr, /**< SOC constraint handler */
377  SCIP_CONS* cons /**< SOC constraint */
378  )
379 {
380  SCIP_CONSHDLRDATA* conshdlrdata;
381  SCIP_CONSDATA* consdata;
382  char nlpform;
383  int i;
384 
385  assert(scip != NULL);
386  assert(cons != NULL);
387 
388  conshdlrdata = SCIPconshdlrGetData(conshdlr);
389  assert(conshdlrdata != NULL);
390 
391  consdata = SCIPconsGetData(cons);
392  assert(consdata != NULL);
393 
394  if( consdata->nlrow != NULL )
395  {
396  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
397  }
398 
399  nlpform = conshdlrdata->nlpform;
400  if( nlpform == 'a' )
401  {
402  /* if the user let us choose, then we take 's' for "small" SOC constraints, but 'q' for large ones,
403  * since the 's' form leads to nvars^2 elements in Hessian, while the 'q' form yields only n elements
404  * however, if there is no expression interpreter, then the NLPI may have trouble, so we always use 'q' in this case
405  */
406  if( consdata->nvars < 100 && conshdlrdata->haveexprint )
407  nlpform = 's';
408  else
409  nlpform = 'q';
410  }
411 
412  switch( nlpform )
413  {
414  case 'e':
415  {
416  /* construct expression exp(\sqrt{\gamma + \sum_{i=1}^{n} (\alpha_i\, (x_i + \beta_i))^2} - alpha_{n+1}(x_{n+1} + beta_{n+1})) */
417 
418  if( consdata->nvars > 0 )
419  {
420  SCIP_EXPR* expr;
421  SCIP_EXPR* exprterm;
422  SCIP_EXPR* expr2;
423  SCIP_EXPRTREE* exprtree;
424 
425  if( consdata->constant != 0.0 )
426  {
427  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &exprterm, SCIP_EXPR_CONST, consdata->constant) ); /* gamma */
428  }
429  else
430  {
431  exprterm = NULL;
432  }
433 
434  for( i = 0; i < consdata->nvars; ++i )
435  {
436  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, i) ); /* x_i */
437  if( consdata->offsets[i] != 0.0 )
438  {
439  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_CONST, consdata->offsets[i]) ); /* beta_i */
440  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_PLUS, expr, expr2) ); /* x_i + beta_i */
441  }
442  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_SQUARE, expr) ); /* (x_i + beta_i)^2 */
443  if( consdata->coefs[i] != 1.0 )
444  {
445  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_CONST, SQR(consdata->coefs[i])) ); /* (alpha_i)^2 */
446  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_MUL, expr, expr2) ); /* (alpha_i)^2 * (x_i + beta_i)^2 */
447  }
448  if( exprterm != NULL )
449  {
450  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &exprterm, SCIP_EXPR_PLUS, exprterm, expr) );
451  }
452  else
453  {
454  exprterm = expr;
455  }
456  }
457 
458  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &exprterm, SCIP_EXPR_SQRT, exprterm) ); /* sqrt(gamma + sum_i (...)^2) */
459 
460  if( consdata->rhsvar != NULL )
461  {
462  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, consdata->nvars) ); /* x_{n+1} */
463  if( consdata->rhsoffset != 0.0 )
464  {
465  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_CONST, consdata->rhsoffset) ); /* beta_{n+1} */
466  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_PLUS, expr, expr2) ); /* x_{n+1} + beta_{n+1} */
467  }
468  if( consdata->rhscoeff != 1.0 )
469  {
470  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_CONST, consdata->rhscoeff) ); /* alpha_{n+1} */
471  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_MUL, expr, expr2) ); /* alpha_{n+1} * (x_{n+1} + beta_{n+1}) */
472  }
473  }
474  else
475  {
476  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_CONST, consdata->rhscoeff * consdata->rhsoffset) );
477  }
478  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &exprterm, SCIP_EXPR_MINUS, exprterm, expr) ); /* sqrt(gamma + sum_i (...)^2) - alpha_{n+1} * (x_{n+1} + beta_{n+1}) */
479 
480  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &exprterm, SCIP_EXPR_EXP, exprterm) ); /* exp(sqrt(gamma + sum_i (...)^2) - alpha_{n+1} * (x_{n+1} + beta_{n+1})) */
481 
482  SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &exprtree, exprterm, consdata->nvars+1, 0, NULL) );
483 
484  SCIP_CALL( SCIPexprtreeSetVars(exprtree, consdata->nvars, consdata->vars) );
485  SCIP_CALL( SCIPexprtreeAddVars(exprtree, 1, &consdata->rhsvar) );
486 
487  SCIP_CALL( SCIPcreateNlRow(scip, &consdata->nlrow, SCIPconsGetName(cons),
488  0.0,
489  0, NULL, NULL,
490  0, NULL, 0, NULL,
491  exprtree, -SCIPinfinity(scip), 1.0,
493 
494  SCIP_CALL( SCIPexprtreeFree(&exprtree) );
495 
496  break;
497  }
498  /* if there are no left-hand-side variables, then we let the 's' case handle it */
499  } /*lint -fallthrough */
500 
501  case 's':
502  {
503  /* construct expression \sqrt{\gamma + \sum_{i=1}^{n} (\alpha_i\, (x_i + \beta_i))^2} */
504 
505  SCIP_EXPR* expr;
506  SCIP_EXPR* exprterm;
507  SCIP_EXPR* expr2;
508  SCIP_EXPRTREE* exprtree;
509  SCIP_Real lincoef;
510 
511  if( consdata->constant != 0.0 )
512  {
513  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &exprterm, SCIP_EXPR_CONST, consdata->constant) ); /* gamma */
514  }
515  else
516  {
517  exprterm = NULL;
518  }
519 
520  for( i = 0; i < consdata->nvars; ++i )
521  {
522  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, i) ); /* x_i */
523  if( consdata->offsets[i] != 0.0 )
524  {
525  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_CONST, consdata->offsets[i]) ); /* beta_i */
526  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_PLUS, expr, expr2) ); /* x_i + beta_i */
527  }
528  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_SQUARE, expr) ); /* (x_i + beta_i)^2 */
529  if( consdata->coefs[i] != 1.0 )
530  {
531  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_CONST, SQR(consdata->coefs[i])) ); /* (alpha_i)^2 */
532  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_MUL, expr, expr2) ); /* (alpha_i)^2 * (x_i + beta_i)^2 */
533  }
534  if( exprterm != NULL )
535  {
536  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &exprterm, SCIP_EXPR_PLUS, exprterm, expr) );
537  }
538  else
539  {
540  exprterm = expr;
541  }
542  }
543 
544  if( exprterm != NULL )
545  {
546  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &exprterm, SCIP_EXPR_SQRT, exprterm) ); /* sqrt(gamma + sum_i (...)^2) */
547  SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &exprtree, exprterm, consdata->nvars, 0, NULL) );
548  SCIP_CALL( SCIPexprtreeSetVars(exprtree, consdata->nvars, consdata->vars) );
549  }
550  else
551  {
552  assert(consdata->nvars == 0);
553  assert(consdata->constant == 0.0);
554  exprtree = NULL;
555  }
556 
557  /* linear and constant part is -\alpha_{n+1} (x_{n+1}+\beta_{n+1}) */
558  lincoef = -consdata->rhscoeff;
559  SCIP_CALL( SCIPcreateNlRow(scip, &consdata->nlrow, SCIPconsGetName(cons),
560  -consdata->rhscoeff * consdata->rhsoffset,
561  1, &consdata->rhsvar, &lincoef,
562  0, NULL, 0, NULL,
563  exprtree, -SCIPinfinity(scip), 0.0,
565 
566  SCIP_CALL( SCIPexprtreeFree(&exprtree) );
567 
568  break;
569  }
570 
571  case 'q':
572  {
573  /* construct quadratic form gamma + sum_{i=1}^{n} (alpha_i (x_i + beta_i))^2 <= (alpha_{n+1} (x_{n+1} + beta_{n+1})^2 */
574  SCIP_QUADELEM sqrterm;
575  SCIP_EXPRCURV curvature;
576  SCIP_Real rhs;
577  int rhsvarpos;
578 
579  /* the expression is convex if alpha_{n+1} is zero */
580  curvature = consdata->rhscoeff == 0.0 ? SCIP_EXPRCURV_CONVEX : SCIP_EXPRCURV_UNKNOWN;
581 
582  /* create initial empty row with left hand side variables */
583  SCIP_CALL( SCIPcreateNlRow(scip, &consdata->nlrow, SCIPconsGetName(cons), 0.0,
584  0, NULL, NULL,
585  consdata->nvars, consdata->vars, 0, NULL,
586  NULL, -SCIPinfinity(scip), 0.0,
587  curvature) );
588 
589  /* add gamma + sum_{i=1}^{n} (alpha_i x_i)^2 + 2 alpha_i beta_i x_i + beta_i^2 */
590  rhs = -consdata->constant;
591  for( i = 0; i < consdata->nvars; ++i )
592  {
593  sqrterm.idx1 = i;
594  sqrterm.idx2 = i;
595  sqrterm.coef = consdata->coefs[i] * consdata->coefs[i];
596  SCIP_CALL( SCIPaddQuadElementToNlRow(scip, consdata->nlrow, sqrterm) );
597 
598  if( ! SCIPisZero(scip, consdata->offsets[i]) )
599  {
600  rhs -= consdata->offsets[i] * consdata->offsets[i];
601  SCIP_CALL( SCIPaddLinearCoefToNlRow(scip, consdata->nlrow, consdata->vars[i], 2.0 * consdata->coefs[i] * consdata->offsets[i]) );
602  }
603  }
604 
605  /* add rhsvar to quadvars of nlrow, if not there yet */
606  rhsvarpos = SCIPnlrowSearchQuadVar(consdata->nlrow, consdata->rhsvar);
607  if( rhsvarpos == -1 )
608  {
609  SCIP_CALL( SCIPaddQuadVarToNlRow(scip, consdata->nlrow, consdata->rhsvar) );
610  rhsvarpos = SCIPnlrowSearchQuadVar(consdata->nlrow, consdata->rhsvar);
611  assert(rhsvarpos >= 0);
612  }
613 
614  /* add -(alpha_{n+1} x_{n+1))^2 - 2 alpha_{n+1} beta_{n+1} x_{n+1} - beta_{n+1}^2 */
615  sqrterm.idx1 = rhsvarpos;
616  sqrterm.idx2 = rhsvarpos;
617  sqrterm.coef = -consdata->rhscoeff * consdata->rhscoeff;
618  SCIP_CALL( SCIPaddQuadElementToNlRow(scip, consdata->nlrow, sqrterm) );
619 
620  if( consdata->rhsoffset != 0.0 )
621  {
622  rhs += consdata->rhsoffset * consdata->rhsoffset;
623  SCIP_CALL( SCIPaddLinearCoefToNlRow(scip, consdata->nlrow, consdata->rhsvar, -2.0 * consdata->rhscoeff * consdata->rhsoffset) );
624  }
625 
626  SCIP_CALL( SCIPchgNlRowRhs(scip, consdata->nlrow, rhs) );
627 
628  break;
629  }
630 
631  case 'd':
632  {
633  /* construct division form (gamma + sum_{i=1}^n (alpha_i(x_i+beta_i))^2)/(alpha_{n+1}(x_{n+1}+beta_{n+1})) <= alpha_{n+1}(x_{n+1}+beta_{n+1}) */
634  SCIP_EXPRTREE* exprtree;
635  SCIP_EXPR* expr;
636  SCIP_EXPR* nominator;
637  SCIP_EXPR* denominator;
638  SCIP_EXPR** exprs;
639  SCIP_EXPRDATA_MONOMIAL** monomials;
640  SCIP_Real lincoef;
641  SCIP_Real one;
642  SCIP_Real two;
643 
644  SCIP_CALL( SCIPallocBufferArray(scip, &exprs, consdata->nvars) );
645  SCIP_CALL( SCIPallocBufferArray(scip, &monomials, consdata->nvars) );
646  one = 1.0;
647  two = 2.0;
648 
649  for( i = 0; i < consdata->nvars; ++i )
650  {
651  /* put x_i + beta_i into exprs[i] */
652  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &exprs[i], SCIP_EXPR_VARIDX, i) );
653  if( consdata->offsets[i] != 0.0 )
654  {
655  SCIP_CALL( SCIPexprCreateLinear(SCIPblkmem(scip), &exprs[i], 1, &exprs[i], &one, consdata->offsets[i]) );
656  }
657 
658  /* create monomial alpha_i^2 y_i^2, where y_i will be x_i + beta_i */
659  SCIP_CALL( SCIPexprCreateMonomial(SCIPblkmem(scip), &monomials[i], consdata->coefs[i] * consdata->coefs[i], 1, &i, &two) );
660  }
661 
662  /* setup polynomial expression for gamma + sum_{i=1}^n alpha_i^2 (x_i+beta_i)^2 */
663  SCIP_CALL( SCIPexprCreatePolynomial(SCIPblkmem(scip), &nominator, consdata->nvars, exprs, consdata->nvars, monomials, consdata->constant, FALSE) ); /*lint !e850 */
664 
665  SCIPfreeBufferArray(scip, &monomials);
666  SCIPfreeBufferArray(scip, &exprs);
667 
668  /* setup alpha_{n+1}(x_{n+1}+beta_{n+1})
669  * assert that this term is >= 0.0 (otherwise constraint is infeasible anyway) */
670  assert(consdata->rhsvar != NULL);
671  assert((consdata->rhscoeff >= 0.0 && !SCIPisNegative(scip, SCIPvarGetLbGlobal(consdata->rhsvar) + consdata->rhsoffset)) ||
672  (consdata->rhscoeff <= 0.0 && !SCIPisPositive(scip, SCIPvarGetUbGlobal(consdata->rhsvar) + consdata->rhsoffset)));
673  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &denominator, SCIP_EXPR_VARIDX, consdata->nvars) );
674  if( consdata->rhscoeff != 1.0 || consdata->rhsoffset != 0.0 )
675  {
676  SCIP_CALL( SCIPexprCreateLinear(SCIPblkmem(scip), &denominator, 1, &denominator, &consdata->rhscoeff, consdata->rhscoeff * consdata->rhsoffset) );
677  }
678 
679  /* setup nominator/denominator */
680  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_DIV, nominator, denominator) );
681 
682  SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &exprtree, expr, 0, 0, NULL) );
683  SCIP_CALL( SCIPexprtreeSetVars(exprtree, consdata->nvars, consdata->vars) );
684  SCIP_CALL( SCIPexprtreeAddVars(exprtree, 1, &consdata->rhsvar) );
685 
686  /* linear and constant part is -\alpha_{n+1} (x_{n+1}+\beta_{n+1}) */
687  lincoef = -consdata->rhscoeff;
688  SCIP_CALL( SCIPcreateNlRow(scip, &consdata->nlrow, SCIPconsGetName(cons),
689  -consdata->rhscoeff * consdata->rhsoffset,
690  1, &consdata->rhsvar, &lincoef,
691  0, NULL, 0, NULL,
692  exprtree, -SCIPinfinity(scip), 0.0,
694 
695  SCIP_CALL( SCIPexprtreeFree(&exprtree) );
696 
697  break;
698  }
699 
700  default:
701  SCIPerrorMessage("unknown value for nlp formulation parameter\n");
702  return SCIP_ERROR;
703  }
704 
705  SCIPdebugMsg(scip, "created nonlinear row representation of SOC constraint\n");
706  SCIPdebugPrintCons(scip, cons, NULL);
707  SCIPdebug( SCIP_CALL( SCIPprintNlRow(scip, consdata->nlrow, NULL) ) );
708 
709  return SCIP_OKAY;
710 }
711 
712 /** evaluates the left hand side of a SOC constraint */
713 static
715  SCIP* scip, /**< SCIP data structure */
716  SCIP_CONS* cons, /**< constraint to evaluate */
717  SCIP_SOL* sol /**< solution to evaluate, or NULL if LP solution should be used */
718  )
719 {
720  SCIP_CONSDATA* consdata;
721  SCIP_VAR* var;
722  SCIP_Real val;
723  int i;
724 
725  assert(scip != NULL);
726  assert(cons != NULL);
727 
728  consdata = SCIPconsGetData(cons);
729  assert(consdata != NULL);
730 
731  consdata->lhsval = consdata->constant;
732 
733  for( i = 0; i < consdata->nvars; ++i )
734  {
735  var = consdata->vars[i];
736  val = SCIPgetSolVal(scip, sol, var);
737 
738  if( SCIPisInfinity(scip, val) || SCIPisInfinity(scip, -val) )
739  {
740  consdata->lhsval = SCIPinfinity(scip);
741  return SCIP_OKAY;
742  }
743 
744  val = consdata->coefs[i] * (val + consdata->offsets[i]);
745  consdata->lhsval += val * val;
746  }
747  consdata->lhsval = sqrt(consdata->lhsval);
748 
749  return SCIP_OKAY;
750 }
751 
752 /** computes violation of a SOC constraint */
753 static
755  SCIP* scip, /**< SCIP data structure */
756  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
757  SCIP_CONS* cons, /**< constraint to evaluate */
758  SCIP_SOL* sol /**< solution to evaluate, or NULL if LP solution should be used */
759  )
760 {
761  SCIP_CONSDATA* consdata;
762  SCIP_Real rhsval;
763  SCIP_Real rhs;
764  SCIP_Real absviol;
765  SCIP_Real relviol;
766 
767  assert(scip != NULL);
768  assert(conshdlr != NULL);
769  assert(cons != NULL);
770 
771  consdata = SCIPconsGetData(cons);
772  assert(consdata != NULL);
773 
774  SCIP_CALL( evalLhs(scip, cons, sol) );
775 
776  if( SCIPisInfinity(scip, consdata->lhsval) )
777  {
778  /* infinity <= infinity is feasible
779  * infinity <= finite value is not feasible and has violation infinity
780  */
781  if( (consdata->rhscoeff > 0.0 && SCIPisInfinity(scip, SCIPgetSolVal(scip, sol, consdata->rhsvar))) ||
782  ( consdata->rhscoeff < 0.0 && SCIPisInfinity(scip, -SCIPgetSolVal(scip, sol, consdata->rhsvar))) )
783  consdata->violation = 0.0;
784  else
785  consdata->violation = SCIPinfinity(scip);
786  return SCIP_OKAY;
787  }
788 
789  rhsval = SCIPgetSolVal(scip, sol, consdata->rhsvar);
790  if( SCIPisInfinity(scip, rhsval) )
791  {
792  consdata->violation = consdata->rhscoeff > 0.0 ? 0.0 : SCIPinfinity(scip);
793  return SCIP_OKAY;
794  }
795  if( SCIPisInfinity(scip, -rhsval) )
796  {
797  consdata->violation = consdata->rhscoeff < 0.0 ? 0.0 : SCIPinfinity(scip);
798  return SCIP_OKAY;
799  }
800 
801  rhs = consdata->rhscoeff * (rhsval + consdata->rhsoffset);
802  consdata->violation = consdata->lhsval - rhs;
803  absviol = consdata->violation;
804  relviol = SCIPrelDiff(consdata->lhsval, rhs);
805  if( consdata->violation <= 0.0 )
806  {
807  /* constraint is not violated for sure */
808  consdata->violation = 0.0;
809  return SCIP_OKAY;
810  }
811 
812  if( sol != NULL )
813  SCIPupdateSolConsViolation(scip, sol, absviol, relviol);
814 
815  return SCIP_OKAY;
816 }
817 
818 /** computes violations for a set of constraints */
819 static
821  SCIP* scip, /**< SCIP data structure */
822  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
823  SCIP_CONS** conss, /**< constraints to evaluate */
824  int nconss, /**< number of constraints to evaluate */
825  SCIP_SOL* sol, /**< solution to evaluate, or NULL if LP solution should be used */
826  SCIP_CONS** maxviolcons /**< a buffer to store pointer to maximal violated constraint, or NULL if of no interest */
827  )
828 {
829  SCIP_CONSDATA* consdata;
830  SCIP_Real maxviol = 0.0;
831  int c;
832 
833  assert(scip != NULL);
834  assert(conss != NULL || nconss == 0);
835 
836  if( maxviolcons != NULL )
837  *maxviolcons = NULL;
838 
839  for( c = 0; c < nconss; ++c )
840  {
841  SCIP_CALL( computeViolation(scip, conshdlr, conss[c], sol) ); /*lint !e613*/
842  if( maxviolcons != NULL )
843  {
844  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
845  assert(consdata != NULL);
846  if( consdata->violation > maxviol && SCIPisGT(scip, consdata->violation, SCIPfeastol(scip)) )
847  {
848  maxviol = consdata->violation;
849  *maxviolcons = conss[c]; /*lint !e613*/
850  }
851  }
852  }
853 
854  return SCIP_OKAY;
855 }
856 
857 /** generate supporting hyperplane in a given solution */
858 static
860  SCIP* scip, /**< SCIP pointer */
861  SCIP_CONS* cons, /**< constraint */
862  SCIP_SOL* sol, /**< solution to separate, or NULL for LP solution */
863  SCIP_ROWPREP** rowprep /**< place to store cut */
864  )
865 {
866  SCIP_CONSDATA* consdata;
867  SCIP_Real val;
868  int i;
869 
870  assert(scip != NULL);
871  assert(cons != NULL);
872  assert(rowprep != NULL);
873 
874  consdata = SCIPconsGetData(cons);
875  assert(consdata != NULL);
876 
877  assert(!SCIPisInfinity(scip, consdata->lhsval));
878 
880  SCIP_CALL( SCIPensureRowprepSize(scip, *rowprep, consdata->nvars+1) );
881  (void) SCIPsnprintf((*rowprep)->name, SCIP_MAXSTRLEN, "%s_linearization_%" SCIP_LONGINT_FORMAT, SCIPconsGetName(cons), SCIPgetNLPs(scip));
882 
883  if( SCIPisPositive(scip, consdata->lhsval) )
884  {
885  /* if lhs is 0, then we cannot linearize
886  * but since we are violated, we have rhs < 0, so underestimating lhs by 0 could still give us a useful cut
887  */
888  SCIP_CALL( SCIPensureRowprepSize(scip, *rowprep, consdata->nvars+1) );
889  for( i = 0; i < consdata->nvars; ++i )
890  {
891  val = SCIPgetSolVal(scip, sol, consdata->vars[i]) + consdata->offsets[i];
892  val *= consdata->coefs[i] * consdata->coefs[i];
893 
894  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, consdata->vars[i], val / consdata->lhsval) );
895 
896  val *= SCIPgetSolVal(scip, sol, consdata->vars[i]);
897  SCIPaddRowprepSide(*rowprep, val);
898  }
899  (*rowprep)->side /= consdata->lhsval;
900  (*rowprep)->side -= consdata->lhsval;
901  }
902 
903  /* add linear rhs: rhscoeff * (rhsvar + rhsoffset) */
904  (*rowprep)->side += consdata->rhscoeff * consdata->rhsoffset;
905  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, consdata->rhsvar, -consdata->rhscoeff) );
906 
907  return SCIP_OKAY;
908 }
909 
910 /** generate supporting hyperplane in a given point */
911 static
913  SCIP* scip, /**< SCIP pointer */
914  SCIP_CONS* cons, /**< constraint */
915  SCIP_Real* x, /**< point (lhs-vars) where to generate cut */
916  SCIP_ROWPREP** rowprep /**< place to store cut */
917  )
918 {
919  SCIP_CONSDATA* consdata;
920  SCIP_Real lhsval;
921  SCIP_Real val;
922  int i;
923 
924  assert(scip != NULL);
925  assert(cons != NULL);
926  assert(rowprep != NULL);
927 
928  consdata = SCIPconsGetData(cons);
929  assert(consdata != NULL);
930 
931  lhsval = consdata->constant;
932  for( i = 0; i < consdata->nvars; ++i )
933  {
934  assert(!SCIPisInfinity(scip, ABS(x[i])));
935 
936  val = consdata->coefs[i] * (x[i] + consdata->offsets[i]);
937  lhsval += val * val;
938  }
939  lhsval = sqrt(lhsval);
940 
941  if( SCIPisZero(scip, lhsval) )
942  { /* do not like to linearize in 0 */
943  return SCIP_OKAY;
944  }
945 
947  SCIP_CALL( SCIPensureRowprepSize(scip, *rowprep, consdata->nvars+1) );
948  (void) SCIPsnprintf((*rowprep)->name, SCIP_MAXSTRLEN, "%s_linearization_%" SCIP_LONGINT_FORMAT, SCIPconsGetName(cons), SCIPgetNLPs(scip));
949 
950  for( i = 0; i < consdata->nvars; ++i )
951  {
952  val = x[i] + consdata->offsets[i];
953  if( SCIPisZero(scip, val) )
954  continue;
955 
956  val *= consdata->coefs[i] * consdata->coefs[i];
957 
958  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, consdata->vars[i], val / lhsval) );
959 
960  val *= x[i];
961  SCIPaddRowprepSide(*rowprep, val);
962  }
963  (*rowprep)->side /= lhsval;
964  (*rowprep)->side -= lhsval - consdata->rhscoeff * consdata->rhsoffset;
965 
966  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, consdata->rhsvar, -consdata->rhscoeff) );
967 
968  return SCIP_OKAY;
969 }
970 
971 /** generate supporting hyperplane w.r.t. solution projected on feasible set
972  *
973  * Instead of linearizing the SOC constraint in the given solution point, this function projects the point
974  * first onto the feasible set of the SOC constraint (w.r.t. euclidean norm (scaled by alpha))
975  * and linearizes the SOC constraint there.
976  * The hope is that this produces somewhat tighter cuts.
977  *
978  * The projection has only be computed for the case gamma = 0.
979  * If gamma > 0, generateCut is called.
980  *
981  * Let \f$\hat x\f$ be sol or the LP solution if sol == NULL.
982  * Then the projected point \f$\tilde x\f$ is given by
983  * \f[
984  * \tilde x_i = \frac{\hat x_i + \lambda \beta_i}{1-\lambda}, \quad i=1,\ldots, n; \quad
985  * \tilde x_{n+1} = \frac{\hat x_{n+1} - \lambda \beta_{n+1}}{1+\lambda}
986  * \f]
987  * where
988  * \f[
989  * \lambda = \frac{1-A}{1+A}, \qquad
990  * A = \frac{\alpha_{n+1}(\hat x_{n+1}+\beta_{n+1})}{\sqrt{\sum_{i=1}^n (\alpha_i(\hat x_i+\beta_i))^2}}
991  * \f]
992  *
993  * If lambda is very close to 1, generateCut is called.
994  *
995  * The generated cut is very similar to the unprojected form.
996  * The only difference is in the right hand side, which is (in the case beta = 0) multiplied by 1/(1-lambda).
997  */
998 static
1000  SCIP* scip, /**< SCIP pointer */
1001  SCIP_CONS* cons, /**< constraint */
1002  SCIP_SOL* sol, /**< solution to separate, or NULL for LP solution */
1003  SCIP_ROWPREP** rowprep /**< place to store cut */
1004  )
1005 {
1006  SCIP_CONSDATA* consdata;
1007  SCIP_Real val;
1008  SCIP_Real A, lambda;
1009  int i;
1010 
1011  assert(scip != NULL);
1012  assert(cons != NULL);
1013  assert(rowprep != NULL);
1014 
1015  consdata = SCIPconsGetData(cons);
1016  assert(consdata != NULL);
1017 
1018  assert(!SCIPisInfinity(scip, consdata->lhsval));
1019 
1020  if( !SCIPisZero(scip, consdata->constant) || !SCIPisPositive(scip, consdata->lhsval) )
1021  {
1022  /* have not thought about the constant=0 case yet; if lhsval is 0, also fall back to simple case */
1023  SCIP_CALL( generateCutSol(scip, cons, sol, rowprep) );
1024  return SCIP_OKAY;
1025  }
1026 
1027  A = consdata->rhscoeff * (SCIPgetSolVal(scip, sol, consdata->rhsvar) + consdata->rhsoffset);
1028  A /= consdata->lhsval;
1029 
1030  lambda = (1.0 - A) / (1.0 + A);
1031 
1032  assert(!SCIPisNegative(scip, lambda)); /* otherwise A > 1, so constraint is not violated */
1033 
1034  SCIPdebugMsg(scip, "A = %g \t lambda = %g\n", A, lambda);
1035 
1036  if( SCIPisFeasEQ(scip, lambda, 1.0) )
1037  { /* avoid numerical difficulties when dividing by (1-lambda) below */
1038  SCIP_CALL( generateCutSol(scip, cons, sol, rowprep) );
1039  return SCIP_OKAY;
1040  }
1041 
1043  SCIP_CALL( SCIPensureRowprepSize(scip, *rowprep, consdata->nvars+1) );
1044  (void) SCIPsnprintf((*rowprep)->name, SCIP_MAXSTRLEN, "%s_linearization_%" SCIP_LONGINT_FORMAT, SCIPconsGetName(cons), SCIPgetNLPs(scip));
1045 
1046  for( i = 0; i < consdata->nvars; ++i )
1047  {
1048  val = SCIPgetSolVal(scip, sol, consdata->vars[i]) + consdata->offsets[i];
1049  val *= consdata->coefs[i] * consdata->coefs[i];
1050 
1051  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, consdata->vars[i], val / consdata->lhsval) );
1052 
1053  val *= SCIPgetSolVal(scip, sol, consdata->vars[i]) + lambda * consdata->offsets[i];
1054  SCIPaddRowprepSide(*rowprep, val);
1055  }
1056  (*rowprep)->side /= consdata->lhsval;
1057  (*rowprep)->side-= consdata->lhsval;
1058  (*rowprep)->side /= 1.0 - lambda;
1059  (*rowprep)->side -= consdata->rhscoeff * consdata->rhsoffset;
1060 
1061  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, consdata->rhsvar, -consdata->rhscoeff) );
1062 
1063  return SCIP_OKAY;
1064 }
1065 
1066 /** generates sparsified supporting hyperplane */
1067 static
1069  SCIP* scip, /**< SCIP pointer */
1070  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
1071  SCIP_CONS* cons, /**< constraint */
1072  SCIP_SOL* sol, /**< solution to separate, or NULL for LP solution */
1073  SCIP_ROWPREP** rowprep, /**< place to store cut */
1074  SCIP_Real minefficacy /**< minimal efficacy for a cut to be accepted */
1075  )
1076 {
1077  SCIP_CONSHDLRDATA* conshdlrdata;
1078  SCIP_CONSDATA* consdata;
1079  SCIP_Real* x;
1080  SCIP_Real* dist; /* distance to 0 */
1081  int* ind; /* indices */
1082  int i;
1083  int maxnz, nextmaxnz;
1084  SCIP_Real efficacy;
1085  SCIP_Real goodefficacy;
1086 
1087  assert(scip != NULL);
1088  assert(conshdlr != NULL);
1089  assert(cons != NULL);
1090  assert(rowprep != NULL);
1091 
1092  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1093  assert(conshdlrdata != NULL);
1094 
1095  consdata = SCIPconsGetData(cons);
1096  assert(consdata != NULL);
1097 
1098  assert(!SCIPisInfinity(scip, consdata->lhsval));
1099 
1100  if( consdata->nvars <= 3 || !SCIPisPositive(scip, consdata->lhsval) )
1101  {
1102  SCIP_CALL( generateCutSol(scip, cons, sol, rowprep) );
1103  return SCIP_OKAY;
1104  }
1105 
1106  goodefficacy = MAX((1.0-conshdlrdata->sparsifymaxloss) * consdata->violation, minefficacy);
1107 
1108  SCIP_CALL( SCIPallocBufferArray(scip, &x, consdata->nvars) );
1109  SCIP_CALL( SCIPallocBufferArray(scip, &dist, consdata->nvars) );
1110  SCIP_CALL( SCIPallocBufferArray(scip, &ind, consdata->nvars) );
1111 
1112  SCIP_CALL( SCIPgetSolVals(scip, sol, consdata->nvars, consdata->vars, x) );
1113  /* distance to "-offset" * alpha_i^2 should indicate loss when moving refpoint to x[i] = -offset[i] */
1114  for( i = 0; i < consdata->nvars; ++i )
1115  {
1116  ind[i] = i;
1117  dist[i] = ABS(x[i] + consdata->offsets[i]);
1118  dist[i] *= consdata->coefs[i] * consdata->coefs[i];
1119  }
1120 
1121  /* sort variables according to dist */
1122  SCIPsortRealInt(dist, ind, consdata->nvars);
1123 
1124  maxnz = 2;
1125  /* set all except biggest maxnz entries in x to -offset */
1126  for( i = 0; i < consdata->nvars - maxnz; ++i )
1127  x[ind[i]] = -consdata->offsets[i];
1128 
1129  do
1130  {
1131  /* @todo speed up a bit by computing efficacy of new cut from efficacy of old cut
1132  * generate row only if efficient enough
1133  */
1134  SCIP_CALL( generateCutPoint(scip, cons, x, rowprep) );
1135 
1136  if( *rowprep != NULL )
1137  {
1138  efficacy = SCIPgetRowprepViolation(scip, *rowprep, sol);
1139  if( SCIPisGT(scip, efficacy, goodefficacy) ||
1140  (maxnz >= consdata->nvars && SCIPisGT(scip, efficacy, minefficacy)) )
1141  {
1142  /* cut cuts off solution and is efficient enough */
1143  SCIPdebugMsg(scip, "accepted cut with %d of %d nonzeros, efficacy = %g\n", maxnz, consdata->nvars, efficacy);
1144  break;
1145  }
1146 
1147  SCIPfreeRowprep(scip, rowprep);
1148  }
1149 
1150  /* cut also not efficient enough if generated in original refpoint (that's bad) */
1151  if( maxnz >= consdata->nvars )
1152  break;
1153 
1154  nextmaxnz = (int)(conshdlrdata->sparsifynzgrowth * maxnz);
1155  if( nextmaxnz == consdata->nvars - 1)
1156  nextmaxnz = consdata->nvars;
1157  else if( nextmaxnz == maxnz )
1158  ++nextmaxnz;
1159 
1160  /* restore entries of x that are nonzero in next attempt */
1161  for( i = MAX(0, consdata->nvars - nextmaxnz); i < consdata->nvars - maxnz; ++i )
1162  x[ind[i]] = SCIPgetSolVal(scip, sol, consdata->vars[ind[i]]);
1163 
1164  maxnz = nextmaxnz;
1165  }
1166  while( TRUE ); /*lint !e506*/
1167 
1168  SCIPfreeBufferArray(scip, &x);
1169  SCIPfreeBufferArray(scip, &dist);
1170  SCIPfreeBufferArray(scip, &ind);
1171 
1172  return SCIP_OKAY;
1173 }
1174 
1175 /** separates a point, if possible */
1176 static
1178  SCIP* scip, /**< SCIP data structure */
1179  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
1180  SCIP_CONS** conss, /**< constraints */
1181  int nconss, /**< number of constraints */
1182  int nusefulconss, /**< number of constraints that seem to be useful */
1183  SCIP_SOL* sol, /**< solution to separate, or NULL for LP solution */
1184  SCIP_Bool inenforcement, /**< whether we are in constraint enforcement */
1185  SCIP_Bool* cutoff, /**< pointer to store whether a fixing leads to a cutoff */
1186  SCIP_Bool* success /**< buffer to store whether the point was separated */
1187  )
1188 {
1189  SCIP_CONSHDLRDATA* conshdlrdata;
1190  SCIP_CONSDATA* consdata;
1191  SCIP_Real minefficacy;
1192  int c;
1193  SCIP_ROW* row;
1194  SCIP_ROWPREP* rowprep;
1195 
1196  assert(scip != NULL);
1197  assert(conss != NULL || nconss == 0);
1198  assert(nusefulconss <= nconss);
1199  assert(cutoff != NULL);
1200  assert(success != NULL);
1201 
1202  *cutoff = FALSE;
1203 
1204  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1205  assert(conshdlrdata != NULL);
1206 
1207  *success = FALSE;
1208 
1209  minefficacy = inenforcement ? SCIPgetLPFeastol(scip) : SCIPgetSepaMinEfficacy(scip);
1210 
1211  for( c = 0; c < nconss; ++c )
1212  {
1213  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
1214  assert(consdata != NULL);
1215 
1216  if( SCIPisGT(scip, consdata->violation, SCIPfeastol(scip)) && !SCIPisInfinity(scip, consdata->violation) )
1217  {
1218  SCIP_Real efficacy;
1219 
1220  rowprep = NULL;
1221 
1222  /* generate cut */
1223  if( conshdlrdata->sparsify )
1224  {
1225  SCIP_CALL( generateSparseCut(scip, conshdlr, conss[c], sol, &rowprep, minefficacy) ); /*lint !e613*/
1226  }
1227  else if( conshdlrdata->projectpoint )
1228  {
1229  SCIP_CALL( generateCutProjectedPoint(scip, conss[c], sol, &rowprep) ); /*lint !e613*/
1230  }
1231  else
1232  {
1233  SCIP_CALL( generateCutSol(scip, conss[c], sol, &rowprep) ); /*lint !e613*/
1234  }
1235 
1236  if( rowprep == NULL )
1237  continue;
1238 
1239  /* NOTE: The way that rowprep was constructed, there should be no need to call SCIPmergeRowprep,
1240  * since no variable gets added twice. However, if rowprep were replacing multiaggregated variables
1241  * (as there can exist for soc cons), then SCIPmergeRowprep would be necessary.
1242  */
1243  /* cleanup rowprep (there is no limit on coefrange for cons_soc) TODO add a coefrange limit? */
1244  SCIP_CALL( SCIPcleanupRowprep(scip, rowprep, sol, SCIPinfinity(scip), minefficacy, NULL, &efficacy) );
1245 
1246  if( SCIPisLE(scip, efficacy, minefficacy) )
1247  {
1248  SCIPfreeRowprep(scip, &rowprep);
1249  continue;
1250  }
1251 
1252  /* cut cuts off solution and efficient enough */
1253  SCIP_CALL( SCIPgetRowprepRowConshdlr(scip, &row, rowprep, conshdlr) );
1254  if( SCIPisCutApplicable(scip, row) )
1255  {
1256  SCIP_CALL( SCIPaddRow(scip, row, FALSE, cutoff) );
1257  SCIP_CALL( SCIPresetConsAge(scip, conss[c]) ); /*lint !e613*/
1258 
1259  *success = TRUE;
1260 
1261  SCIPdebugMsg(scip, "added cut with efficacy %g\n", SCIPgetCutEfficacy(scip, sol, row));
1262 
1263  /* mark row as not removable from LP for current node, if in enforcement */
1264  if( inenforcement && !conshdlrdata->enfocutsremovable )
1265  SCIPmarkRowNotRemovableLocal(scip, row);
1266  }
1267 
1268  SCIP_CALL( SCIPreleaseRow (scip, &row) );
1269  SCIPfreeRowprep(scip, &rowprep);
1270  }
1271 
1272  if( *cutoff )
1273  break;
1274 
1275  /* enforce only useful constraints
1276  * others are only checked and enforced if we are still feasible or have not found a separating cut yet
1277  */
1278  if( c >= nusefulconss && *success )
1279  break;
1280  }
1281 
1282  return SCIP_OKAY;
1283 }
1284 
1285 /** adds linearizations cuts for convex constraints w.r.t. a given reference point to cutpool and sepastore
1286  *
1287  * If separatedlpsol is not NULL, then a cut that separates the LP solution is added to the sepastore and is forced to enter the LP.
1288  * If separatedlpsol is not NULL, but cut does not separate the LP solution, then it is added to the cutpool only.
1289  * If separatedlpsol is NULL, then cut is added to cutpool only.
1290  */
1291 static
1293  SCIP* scip, /**< SCIP data structure */
1294  SCIP_CONSHDLR* conshdlr, /**< quadratic constraints handler */
1295  SCIP_CONS** conss, /**< constraints */
1296  int nconss, /**< number of constraints */
1297  SCIP_SOL* ref, /**< reference point where to linearize, or NULL for LP solution */
1298  SCIP_Bool* separatedlpsol, /**< buffer to store whether a cut that separates the current LP solution was found and added to LP, or NULL if adding to cutpool only */
1299  SCIP_Real minefficacy, /**< minimal efficacy of a cut when checking for separation of LP solution */
1300  SCIP_Bool* cutoff /**< pointer to store whether a cutoff was detected */
1301  )
1302 {
1303  SCIP_CONSDATA* consdata;
1304  SCIP_Bool addedtolp;
1305  SCIP_ROWPREP* rowprep;
1306  int c;
1307 
1308  assert(scip != NULL);
1309  assert(conshdlr != NULL);
1310  assert(conss != NULL || nconss == 0);
1311  assert(cutoff != NULL);
1312  *cutoff = FALSE;
1313 
1314  if( separatedlpsol != NULL )
1315  *separatedlpsol = FALSE;
1316 
1317  for( c = 0; c < nconss && !(*cutoff); ++c )
1318  {
1319  assert(conss[c] != NULL); /*lint !e613 */
1320 
1321  if( SCIPconsIsLocal(conss[c]) ) /*lint !e613 */
1322  continue;
1323 
1324  consdata = SCIPconsGetData(conss[c]); /*lint !e613 */
1325  assert(consdata != NULL);
1326 
1327  SCIP_CALL( evalLhs(scip, conss[c], ref) ); /*lint !e613 */
1328  if( !SCIPisPositive(scip, consdata->lhsval) || SCIPisInfinity(scip, consdata->lhsval) )
1329  {
1330  SCIPdebugMsg(scip, "skip adding linearization for <%s> since lhs is %g\n", SCIPconsGetName(conss[c]), consdata->lhsval); /*lint !e613 */
1331  continue;
1332  }
1333 
1334  SCIP_CALL( generateCutSol(scip, conss[c], ref, &rowprep) ); /*lint !e613 */
1335 
1336  if( rowprep == NULL )
1337  continue;
1338 
1339  addedtolp = FALSE;
1340 
1341  /* if caller wants, then check if cut separates LP solution and add to sepastore if so */
1342  if( separatedlpsol != NULL )
1343  {
1344  if( SCIPgetRowprepViolation(scip, rowprep, NULL) >= minefficacy )
1345  {
1346  SCIP_ROW* row;
1347 
1348  SCIP_CALL( SCIPgetRowprepRowConshdlr(scip, &row, rowprep, conshdlr) );
1349  SCIP_CALL( SCIPaddRow(scip, row, TRUE, cutoff) );
1350  SCIP_CALL( SCIPreleaseRow(scip, &row) );
1351 
1352  *separatedlpsol = TRUE;
1353  addedtolp = TRUE;
1354  }
1355  }
1356 
1357  if( !addedtolp && !rowprep->local )
1358  {
1359  SCIP_ROW* row;
1360 
1361  SCIP_CALL( SCIPgetRowprepRowConshdlr(scip, &row, rowprep, conshdlr) );
1362  SCIP_CALL( SCIPaddPoolCut(scip, row) );
1363  SCIP_CALL( SCIPreleaseRow(scip, &row) );
1364  }
1365 
1366  SCIPfreeRowprep(scip, &rowprep);
1367  }
1368 
1369  return SCIP_OKAY;
1370 }
1371 
1372 /** processes the event that a new primal solution has been found */
1373 static
1374 SCIP_DECL_EVENTEXEC(processNewSolutionEvent)
1376  SCIP_CONSHDLRDATA* conshdlrdata;
1377  SCIP_CONSHDLR* conshdlr;
1378  SCIP_CONS** conss;
1379  int nconss;
1380  SCIP_SOL* sol;
1381  SCIP_Bool cutoff;
1382 
1383  assert(scip != NULL);
1384  assert(event != NULL);
1385  assert(eventdata != NULL);
1386  assert(eventhdlr != NULL);
1387 
1388  assert((SCIPeventGetType(event) & SCIP_EVENTTYPE_SOLFOUND) != 0);
1389 
1390  conshdlr = (SCIP_CONSHDLR*)eventdata;
1391 
1392  nconss = SCIPconshdlrGetNConss(conshdlr);
1393 
1394  if( nconss == 0 )
1395  return SCIP_OKAY;
1396 
1397  sol = SCIPeventGetSol(event);
1398  assert(sol != NULL);
1399 
1400  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1401  assert(conshdlrdata != NULL);
1402 
1403  /* we are only interested in solution coming from some heuristic other than trysol, but not from the tree
1404  * the reason for ignoring trysol solutions is that they may come from an NLP solve in sepalp, where we already added linearizations,
1405  * or are from the tree, but postprocessed via proposeFeasibleSolution
1406  */
1407  if( SCIPsolGetHeur(sol) == NULL || SCIPsolGetHeur(sol) == conshdlrdata->trysolheur )
1408  return SCIP_OKAY;
1409 
1410  conss = SCIPconshdlrGetConss(conshdlr);
1411  assert(conss != NULL);
1412 
1413  SCIPdebugMsg(scip, "caught new sol event %" SCIP_EVENTTYPE_FORMAT " from heur <%s>; have %d conss\n", SCIPeventGetType(event),
1414  SCIPheurGetName(SCIPsolGetHeur(sol)), nconss);
1415 
1416  SCIP_CALL( addLinearizationCuts(scip, conshdlr, conss, nconss, sol, NULL, 0.0, &cutoff) );
1417  /* ignore cutoff, cannot return status */
1418 
1419  return SCIP_OKAY;
1420 }
1421 
1422 /** removes fixed variables, replace aggregated and negated variables
1423  *
1424  * repeats replacements until no further change is found;
1425  * takes care of capture/release and locks, but not of variable events (assumes that var events are not caught yet)
1426  */
1427 static
1429  SCIP* scip, /**< SCIP data structure */
1430  SCIP_CONSHDLR* conshdlr, /**< constraint handler for signpower constraints */
1431  SCIP_CONS* cons, /**< constraint */
1432  int* ndelconss, /**< counter for number of deleted constraints */
1433  int* nupgdconss, /**< counter for number of upgraded constraints */
1434  int* nchgbds, /**< counter for number of bound changes */
1435  int* nfixedvars, /**< counter for number of fixed variables */
1436  SCIP_Bool* iscutoff, /**< to store whether constraint cannot be satisfied */
1437  SCIP_Bool* isdeleted /**< to store whether constraint is redundant and can be deleted */
1438  )
1439 {
1440  SCIP_CONSDATA* consdata;
1441  SCIP_CONSHDLRDATA* conshdlrdata;
1442  SCIP_Bool havechange;
1443  SCIP_Bool haveremovedvar;
1444  int i;
1445  SCIP_VAR* x;
1446  SCIP_Real coef;
1447  SCIP_Real offset;
1448 
1449  assert(scip != NULL);
1450  assert(conshdlr != NULL);
1451  assert(cons != NULL);
1452  assert(iscutoff != NULL);
1453  assert(isdeleted != NULL);
1454 
1455  *iscutoff = FALSE;
1456  *isdeleted = FALSE;
1457 
1458  consdata = SCIPconsGetData(cons);
1459  assert(consdata != NULL);
1460 
1461  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1462  assert(conshdlrdata != NULL);
1463 
1464  SCIPdebugMsg(scip, "remove fixed variables from constraint <%s>\n", SCIPconsGetName(cons));
1465  SCIPdebugPrintCons(scip, cons, NULL);
1466 
1467  havechange = FALSE;
1468  haveremovedvar = FALSE;
1469 
1470  /* process variables on left hand side */
1471  for( i = 0; i < consdata->nvars; ++i )
1472  {
1473  x = consdata->vars[i];
1474  assert(x != NULL);
1476 
1478  continue;
1479 
1480  havechange = TRUE;
1481 
1482  /* drop variable event and unlock and release variable */
1483  SCIP_CALL( dropLhsVarEvents(scip, conshdlrdata->eventhdlr, cons, i) );
1484  SCIP_CALL( SCIPunlockVarCons(scip, x, cons, TRUE, TRUE) );
1485  SCIP_CALL( SCIPreleaseVar(scip, &consdata->vars[i]) );
1486 
1487  coef = 1.0;
1488  offset = consdata->offsets[i];
1489  SCIP_CALL( SCIPgetProbvarSum(scip, &x, &coef, &offset) );
1490 
1491  SCIPdebugMsg(scip, " lhs term at position %d is replaced by %g * <%s> + %g\n",
1492  i, coef, SCIPvarGetName(x), offset);
1493 
1494  /* if variable has been fixed, add (alpha*offset)^2 to gamma and continue */
1495  if( coef == 0.0 || x == NULL )
1496  {
1497  consdata->constant += consdata->coefs[i] * consdata->coefs[i] * offset * offset;
1498  consdata->offsets[i] = 0.0;
1499  haveremovedvar = TRUE;
1500  continue;
1501  }
1502 
1504 
1505  /* replace coefs[i] * (vars[i] + offsets[i]) by coefs[i]*coef * (x + offsets[i]/coef) */
1506  consdata->offsets[i] = offset;
1507  if( coef != 1.0 )
1508  {
1509  consdata->coefs[i] = REALABS(coef * consdata->coefs[i]);
1510  consdata->offsets[i] /= coef;
1511  }
1512  consdata->vars[i] = x;
1513 
1514  /* capture and lock new variable, catch variable events */
1515  SCIP_CALL( SCIPcaptureVar(scip, consdata->vars[i]) );
1516  SCIP_CALL( SCIPlockVarCons(scip, consdata->vars[i], cons, TRUE, TRUE) );
1517  SCIP_CALL( catchLhsVarEvents(scip, conshdlrdata->eventhdlr, cons, i) );
1518  }
1519 
1520  /* process variable on right hand side */
1521  x = consdata->rhsvar;
1522  assert(x != NULL);
1524  {
1525  havechange = TRUE;
1526 
1527  /* drop variable event and unlock and release variable */
1528  SCIP_CALL( dropRhsVarEvents(scip, conshdlrdata->eventhdlr, cons) );
1529  SCIP_CALL( SCIPreleaseVar(scip, &consdata->rhsvar) );
1530  SCIP_CALL( SCIPunlockVarCons(scip, x, cons, consdata->rhscoeff > 0.0, consdata->rhscoeff < 0.0) );
1531 
1532  coef = 1.0;
1533  offset = 0.0;
1534  SCIP_CALL( SCIPgetProbvarSum(scip, &x, &coef, &offset) );
1535 
1536  SCIPdebugMsg(scip, " rhs variable is replaced by %g * <%s> + %g\n", coef, SCIPvarGetName(x), offset);
1537 
1538  if( coef == 0.0 || x == NULL )
1539  {
1540  /* if variable has been fixed, add offset to rhsoffset */
1541  consdata->rhsoffset += offset;
1542  }
1543  else
1544  {
1545  /* replace rhscoef * (rhsvar + rhsoffset) by rhscoef*coef * (x + offset/coef + rhsoffset/coef) */
1547 
1548  consdata->rhsoffset = (consdata->rhsoffset + offset) / coef;
1549  consdata->rhscoeff *= coef;
1550  consdata->rhsvar = x;
1551 
1552  /* capture and lock new variable, catch variable events */
1553  SCIP_CALL( SCIPcaptureVar(scip, consdata->rhsvar) );
1554  SCIP_CALL( SCIPlockVarCons(scip, consdata->rhsvar, cons, consdata->rhscoeff > 0.0, consdata->rhscoeff < 0.0) );
1555  SCIP_CALL( catchRhsVarEvents(scip, conshdlrdata->eventhdlr, cons) );
1556  }
1557  }
1558 
1559  if( !havechange )
1560  return SCIP_OKAY;
1561 
1562  /* free nonlinear row representation */
1563  if( consdata->nlrow != NULL )
1564  {
1565  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
1566  }
1567 
1568  /* if a variable has been removed, close gaps in vars array */
1569  if( haveremovedvar )
1570  {
1571  int oldnvars;
1572 
1573  /* due to the realloc of the block memory below and the way we store the eventdata in consdata, we best drop all events here and catch them again below */
1574  SCIP_CALL( dropVarEvents(scip, conshdlrdata->eventhdlr, cons) );
1575 
1576  oldnvars = consdata->nvars;
1577  for( i = 0; i < consdata->nvars; ++i )
1578  {
1579  /* forget about empty places at end of vars array */
1580  while( consdata->nvars && consdata->vars[consdata->nvars-1] == NULL )
1581  --consdata->nvars;
1582 
1583  /* all variables at index >= i have been removed */
1584  if( i == consdata->nvars )
1585  break;
1586 
1587  if( consdata->vars[i] != NULL )
1588  continue;
1589 
1590  /* move variable from position nvars-1 to position i */
1591 
1592  assert(consdata->nvars >= 1);
1593  assert(consdata->vars[consdata->nvars-1] != NULL);
1594 
1595  consdata->vars[i] = consdata->vars[consdata->nvars-1];
1596  consdata->offsets[i] = consdata->offsets[consdata->nvars-1];
1597  consdata->coefs[i] = consdata->coefs[consdata->nvars-1];
1598 
1599  --consdata->nvars;
1600  }
1601 
1602  assert(consdata->nvars < oldnvars);
1603 
1604  /* shrink arrays in consdata */
1605  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->vars, oldnvars, consdata->nvars) );
1606  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->offsets, oldnvars, consdata->nvars) );
1607  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->coefs, oldnvars, consdata->nvars) );
1608 
1609  SCIP_CALL( catchVarEvents(scip, conshdlrdata->eventhdlr, cons) );
1610  }
1611 
1612  SCIPdebugMsg(scip, "\t-> ");
1613  SCIPdebugPrintCons(scip, cons, NULL);
1614 
1615  if( consdata->nvars == 0 )
1616  { /* all variables on left hand size have been removed, remaining constraint is sqrt(gamma) <= ... */
1617  assert(!SCIPisNegative(scip, consdata->constant));
1618  if( consdata->rhsvar == NULL )
1619  { /* also rhsvar has been removed, remaining constraint is sqrt(gamma) <= rhscoeff * rhsoffset */
1620  if( SCIPisFeasLE(scip, sqrt(consdata->constant), consdata->rhscoeff*consdata->rhsoffset) )
1621  {
1622  SCIPdebugMsg(scip, "remove redundant constraint <%s> after fixing all variables\n", SCIPconsGetName(cons));
1623  }
1624  else
1625  {
1626  SCIPdebugMsg(scip, "found problem infeasible after fixing all variables in <%s>\n", SCIPconsGetName(cons));
1627  *iscutoff = TRUE;
1628  }
1629  ++*ndelconss;
1630  }
1631  else if( !SCIPvarIsActive(consdata->rhsvar) )
1632  { /* remaining constraint is sqrt(gamma) - rhscoeff * rhsoffset <= rhscoeff * rhsvar, and rhsvar is probably multi-aggregated */
1633  SCIP_CONS* lincons;
1634 
1635  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, SCIPconsGetName(cons), 1, &consdata->rhsvar, &consdata->rhscoeff,
1636  sqrt(consdata->constant) - consdata->rhscoeff * consdata->rhsoffset, SCIPinfinity(scip),
1640  SCIPconsIsStickingAtNode(cons)) );
1641  SCIP_CALL( SCIPaddCons(scip, lincons) );
1642  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
1643  ++*nupgdconss;
1644  }
1645  else if( consdata->rhscoeff > 0.0 )
1646  { /* remaining constraint is sqrt(gamma) / rhscoeff - rhsoffset <= rhsvar */
1647  SCIP_Bool tightened;
1648  SCIP_CALL( SCIPtightenVarLb(scip, consdata->rhsvar, sqrt(consdata->constant) / consdata->rhscoeff - consdata->rhsoffset, TRUE, iscutoff, &tightened) );
1649  if( *iscutoff )
1650  {
1651  SCIPdebugMsg(scip, "found problem infeasible after fixing all lhs variables in <%s> and tightening lower bound of rhs var\n", SCIPconsGetName(cons));
1652  }
1653  else if( tightened )
1654  {
1655  SCIPdebugMsg(scip, "remove redundant constraint <%s> after fixing all lhs variables and tightening lower bound of rhs var\n", SCIPconsGetName(cons));
1656  ++*nchgbds;
1657  }
1658  else
1659  {
1660  SCIPdebugMsg(scip, "remove redundant constraint <%s> after fixing all lhs variables\n", SCIPconsGetName(cons));
1661  }
1662  ++*ndelconss;
1663  }
1664  else
1665  { /* remaining constraint is sqrt(gamma) / rhscoeff - rhsoffset >= rhsvar */
1666  SCIP_Bool tightened;
1667  SCIP_CALL( SCIPtightenVarUb(scip, consdata->rhsvar, sqrt(consdata->constant) / consdata->rhscoeff - consdata->rhsoffset, TRUE, iscutoff, &tightened) );
1668  if( *iscutoff )
1669  {
1670  SCIPdebugMsg(scip, "found problem infeasible after fixing all lhs variables in <%s> and tightening upper bound of rhs var\n", SCIPconsGetName(cons));
1671  }
1672  else if( tightened )
1673  {
1674  SCIPdebugMsg(scip, "remove redundant constraint <%s> after fixing all lhs variables and tightening upper bound of rhs var\n", SCIPconsGetName(cons));
1675  ++*nchgbds;
1676  }
1677  else
1678  {
1679  SCIPdebugMsg(scip, "remove redundant constraint <%s> after fixing all lhs variables\n", SCIPconsGetName(cons));
1680  }
1681  ++*ndelconss;
1682  }
1683  SCIP_CALL( SCIPdelCons(scip, cons) );
1684  *isdeleted = TRUE;
1685  return SCIP_OKAY;
1686  }
1687 
1688  if( consdata->rhsvar == NULL )
1689  { /* constraint becomes sum_i (alpha_i*(x_i+beta_i))^2 <= (rhscoeff*rhsoffset)^2 - gamma */
1690  if( consdata->nvars > 1 )
1691  { /* upgrade to quadratic constraint */
1692  SCIP_CONS* quadcons;
1693  SCIP_QUADVARTERM* quadvarterms;
1694  SCIP_Real rhs;
1695 
1696  SCIP_CALL( SCIPallocBufferArray(scip, &quadvarterms, consdata->nvars) );
1697  BMSclearMemoryArray(quadvarterms, consdata->nvars);
1698  rhs = consdata->rhscoeff * consdata->rhsoffset;
1699  rhs = rhs * rhs - consdata->constant;
1700 
1701  for( i = 0; i < consdata->nvars; ++i )
1702  {
1703  quadvarterms[i].var = consdata->vars[i];
1704  quadvarterms[i].sqrcoef = consdata->coefs[i] * consdata->coefs[i];
1705  if( consdata->offsets[i] != 0.0 )
1706  {
1707  quadvarterms[i].lincoef = 2 * consdata->offsets[i] * quadvarterms[i].sqrcoef;
1708  rhs -= quadvarterms[i].sqrcoef * consdata->offsets[i]*consdata->offsets[i];
1709  }
1710  }
1711 
1712  assert(!SCIPconsIsStickingAtNode(cons));
1713  SCIP_CALL( SCIPcreateConsQuadratic2(scip, &quadcons, SCIPconsGetName(cons), 0, NULL, NULL,
1714  consdata->nvars, quadvarterms, 0, NULL, -SCIPinfinity(scip), rhs,
1718  SCIP_CALL( SCIPaddCons(scip, quadcons) );
1719  SCIPdebugMsg(scip, "upgraded <%s> to quadratic constraint: ", SCIPconsGetName(cons));
1720  SCIPdebugPrintCons(scip, quadcons, NULL);
1721 
1722  SCIP_CALL( SCIPreleaseCons(scip, &quadcons) );
1723 
1724  SCIPfreeBufferArray(scip, &quadvarterms);
1725 
1726  ++*nupgdconss;
1727  }
1728  else if( !SCIPvarIsActive(consdata->vars[0]) )
1729  { /* constraint is |alpha*(x+beta)| <= sqrt((rhscoeff*rhsoffset)^2 - gamma), but x is probably multaggr. -> turn into ranged linear constraint */
1730  SCIP_CONS* lincons;
1731 
1732  /* create constraint alpha*x <= sqrt((rhscoeff*rhsoffset)^2 - gamma) - alpha*beta
1733  * alpha*x >= -sqrt((rhscoeff*rhsoffset)^2 - gamma) - alpha*beta */
1734  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, SCIPconsGetName(cons), 1, &consdata->vars[0], &consdata->coefs[0],
1735  -sqrt(consdata->rhscoeff * consdata->rhscoeff * consdata->rhsoffset * consdata->rhsoffset - consdata->constant) - consdata->coefs[0] * consdata->offsets[0],
1736  +sqrt(consdata->rhscoeff * consdata->rhscoeff * consdata->rhsoffset * consdata->rhsoffset - consdata->constant) - consdata->coefs[0] * consdata->offsets[0],
1740  SCIPconsIsStickingAtNode(cons)) );
1741  SCIP_CALL( SCIPaddCons(scip, lincons) );
1742  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
1743 
1744  ++*nupgdconss;
1745  }
1746  else
1747  { /* constraint is |alpha*(x+beta)| <= sqrt((rhscoeff*rhsoffset)^2 - gamma) -> propagate bounds */
1748  SCIP_Bool tightened;
1749  SCIP_Real rhs;
1750  assert(consdata->nvars == 1); /* case == 0 handled before */
1751  rhs = consdata->rhscoeff * consdata->rhsoffset;
1752  rhs = rhs * rhs;
1753  if( SCIPisNegative(scip, rhs - consdata->constant) )
1754  { /* take this as infeasible */
1755  SCIPdebugMsg(scip, "found problem infeasible after fixing rhs and all except one lhs variables in <%s>\n", SCIPconsGetName(cons));
1756  *iscutoff = TRUE;
1757  }
1758  else
1759  {
1760  rhs -= consdata->constant;
1761  rhs = rhs < 0.0 ? 0.0 : sqrt(rhs);
1762 
1763  if( SCIPisZero(scip, rhs) )
1764  { /* constraint is x = -beta */
1765  SCIP_CALL( SCIPfixVar(scip, consdata->vars[0], -consdata->offsets[0], iscutoff, &tightened) );
1766  if( *iscutoff )
1767  {
1768  SCIPdebugMsg(scip, "found problem infeasible after fixing rhs and all except one lhs variables and fixing remaining lhs var in <%s>\n", SCIPconsGetName(cons));
1769  }
1770  else if( tightened )
1771  {
1772  SCIPdebugMsg(scip, "remove redundant constraint <%s> after fixing rhs and all except one lhs variables and fixing remaining lhs var\n", SCIPconsGetName(cons));
1773  ++*nfixedvars;
1774  }
1775  else
1776  {
1777  SCIPdebugMsg(scip, "remove redundant constraint <%s> after fixing rhs and all except one lhs variables and fixing remaining lhs var\n", SCIPconsGetName(cons));
1778  }
1779  }
1780  else
1781  { /* constraint is -rhs/|alpha| - beta <= x <= rhs/|alpha| - beta */
1782  rhs /= ABS(consdata->coefs[0]);
1783  SCIP_CALL( SCIPtightenVarLb(scip, consdata->vars[0], -rhs - consdata->offsets[0], TRUE, iscutoff, &tightened) );
1784  if( *iscutoff )
1785  {
1786  SCIPdebugMsg(scip, "found problem infeasible after fixing rhs and all except one lhs variables and tightening lower bound of remaining lhs var in <%s>\n", SCIPconsGetName(cons));
1787  }
1788  else
1789  {
1790  if( tightened )
1791  ++*nchgbds;
1792  SCIP_CALL( SCIPtightenVarUb(scip, consdata->vars[0], rhs - consdata->offsets[0], TRUE, iscutoff, &tightened) );
1793  if( *iscutoff )
1794  {
1795  SCIPdebugMsg(scip, "found problem infeasible after fixing rhs and all except one lhs variables and tightening upper bound of remaining lhs var in <%s>\n", SCIPconsGetName(cons));
1796  }
1797  else if( tightened )
1798  ++*nchgbds;
1799  }
1800  if( !*iscutoff )
1801  {
1802  SCIPdebugMsg(scip, "remove redundant constraint <%s> after fixing rhs and all except one lhs variables and tightening bounds on remaining lhs var\n", SCIPconsGetName(cons));
1803  }
1804  }
1805  }
1806  ++*ndelconss;
1807  }
1808  *isdeleted = TRUE;
1809  SCIP_CALL( SCIPdelCons(scip, cons) );
1810  return SCIP_OKAY;
1811  }
1812 
1813  if( consdata->nvars == 1 && SCIPisZero(scip, consdata->constant) )
1814  { /* one variable on lhs left and no constant, constraint becomes |alpha*(x+beta)| <= rhscoef*(rhsvar+rhsoffset) -> upgrade to two linear constraints */
1815  SCIP_CONS* lincons;
1816  SCIP_VAR* vars[2];
1817  SCIP_Real coefs[2];
1818  SCIP_Real rhs;
1819  assert(consdata->rhsvar != NULL); /* case == NULL has been handled before */
1820 
1821  vars[0] = consdata->vars[0];
1822  vars[1] = consdata->rhsvar;
1823  coefs[0] = consdata->coefs[0];
1824  coefs[1] = -consdata->rhscoeff;
1825  rhs = consdata->rhscoeff * consdata->rhsoffset - coefs[0] * consdata->offsets[0];
1826 
1827  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, SCIPconsGetName(cons), 2, vars, coefs, -SCIPinfinity(scip), rhs,
1831  SCIPconsIsStickingAtNode(cons)) );
1832  SCIP_CALL( SCIPaddCons(scip, lincons) );
1833  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
1834 
1835  coefs[0] = -coefs[0];
1836  rhs = consdata->rhscoeff * consdata->rhsoffset - coefs[0] * consdata->offsets[0];
1837 
1838  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, SCIPconsGetName(cons), 2, vars, coefs, -SCIPinfinity(scip), rhs,
1842  SCIPconsIsStickingAtNode(cons)) );
1843  SCIP_CALL( SCIPaddCons(scip, lincons) );
1844  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
1845 
1846  SCIPdebugMsg(scip, "upgraded <%s> to two linear constraint\n", SCIPconsGetName(cons));
1847 
1848  ++*nupgdconss;
1849  SCIP_CALL( SCIPdelCons(scip, cons) );
1850  *isdeleted = TRUE;
1851  return SCIP_OKAY;
1852  }
1853 
1854  return SCIP_OKAY;
1855 }
1856 
1857 
1858 /** adds the linear outer-approximation of Glineur et.al. for a SOC constraint of dimension 3
1859  *
1860  * Input is the data for a constraint \f$\sqrt{(\alpha_1(x_1+offset1))^2 + (\alpha_2(x_2+offset2))^2)} \leq \alpha_3(x_3+offset3)\f$.
1861  * Here \f$\alpha3 > 0\f$, and the lower bound of \f$x_3 \geq -offset3\f$.
1862  * Also x2 = NULL is allowed, in which case the second term is assumed to be constant, and \f$offset2 \neq 0\f$ is needed.
1863  */
1864 static
1866  SCIP* scip, /**< SCIP data structure */
1867  SCIP_CONS* cons, /**< original constraint */
1868  SCIP_VAR* x1, /**< variable x1 */
1869  SCIP_VAR* x2, /**< variable x2 */
1870  SCIP_VAR* x3, /**< variable x3 */
1871  SCIP_Real alpha1, /**< coefficient of x1 */
1872  SCIP_Real alpha2, /**< coefficient of x2 */
1873  SCIP_Real alpha3, /**< coefficient of x3 */
1874  SCIP_Real offset1, /**< offset of x1 */
1875  SCIP_Real offset2, /**< offset of x2 */
1876  SCIP_Real offset3, /**< offset of x3 */
1877  int N, /**< size of linear approximation, need to be >= 1 */
1878  const char* basename, /**< string to use for building variable and constraint names */
1879  int* naddconss /**< buffer where to add the number of added constraints */
1880  )
1881 {
1882  SCIP_CONS* lincons;
1883  SCIP_VAR* vars[3];
1884  SCIP_Real vals[3];
1885  char varname[255];
1886  char linname[255];
1887  int i;
1888  SCIP_VAR** avars;
1889  SCIP_VAR** bvars;
1890  SCIP_Real val;
1891 
1892  assert(scip != NULL);
1893  assert(cons != NULL);
1894  assert(x1 != NULL);
1895  assert(x2 != NULL || !SCIPisZero(scip, offset2));
1896  assert(x3 != NULL);
1897  assert(SCIPisPositive(scip, alpha3));
1898  assert(SCIPisGE(scip, SCIPconsIsLocal(cons) ? SCIPvarGetLbLocal(x3) : SCIPvarGetLbGlobal(x3), -offset3));
1899  assert(basename != NULL);
1900  assert(N >= 1);
1901  assert(naddconss != NULL);
1902 
1903  SCIPdebugMsg(scip, "Creating linear Glineur outer-approximation for <%s>.\n", basename);
1904  SCIPdebugMsg(scip, "sqr(%g(%s+%g)) + sqr(%g(%s+%g)) <= sqr(%g(%s+%g)).\n",
1905  alpha1, SCIPvarGetName(x1), offset1, alpha2, x2 ? SCIPvarGetName(x2) : "0", offset2, alpha3, SCIPvarGetName(x3), offset3);
1906 
1907  SCIP_CALL( SCIPallocBufferArray(scip, &avars, N+1) );
1908  SCIP_CALL( SCIPallocBufferArray(scip, &bvars, N+1) );
1909 
1910  /* create additional variables; we do not use avars[0] and bvars[0] */
1911  for( i = 1; i <= N; ++i )
1912  {
1913  (void) SCIPsnprintf(varname, 255, "soc#%s_a%d", basename, i);
1914  SCIP_CALL( SCIPcreateVar(scip, &avars[i], varname, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0,
1916  SCIP_CALL( SCIPaddVar(scip, avars[i]) );
1917 
1918  (void) SCIPsnprintf(varname, 255, "soc#%s_b%d", basename, i);
1919  SCIP_CALL( SCIPcreateVar(scip, &bvars[i], varname, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0,
1921  SCIP_CALL( SCIPaddVar(scip, bvars[i]) );
1922  }
1923 
1924  /* create linear constraints for the first case
1925  * cos(pi) = -1, sin(pi) = 0
1926  * -> a_1 = - alpha1 (x1 + offset1) -> -alpha1*x1 - a_1 = alpha1*offset1
1927  * -> b_1 >= | alpha2 (x2 + offset2) | -> alpha2*x2 - b_1 <= -alpha2*offset2
1928  * alpha2*x2 + b_1 >= -alpha2*offset2
1929  */
1930 
1931  vars[0] = x1;
1932  vals[0] = -alpha1;
1933  vars[1] = avars[1];
1934  vals[1] = -1.0;
1935 
1936  (void) SCIPsnprintf(linname, 255, "soc#%s#a%d", basename, 0);
1937  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, alpha1*offset1, alpha1*offset1,
1943  SCIP_CALL( SCIPaddCons(scip, lincons) );
1944  SCIPdebugPrintCons(scip, lincons, NULL);
1945  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
1946  ++(*naddconss);
1947 
1948  if( x2 != NULL )
1949  {
1950  vars[0] = x2;
1951  vals[0] = alpha2;
1952  vars[1] = bvars[1];
1953  vals[1] = -1.0;
1954 
1955  (void) SCIPsnprintf(linname, 255, "soc#%s#b%d", basename, 0);
1956  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, -SCIPinfinity(scip), -alpha2*offset2,
1962  SCIP_CALL( SCIPaddCons(scip, lincons) );
1963  SCIPdebugPrintCons(scip, lincons, NULL);
1964  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
1965  ++(*naddconss);
1966 
1967  vars[0] = x2;
1968  vals[0] = alpha2;
1969  vars[1] = bvars[1];
1970  vals[1] = 1.0;
1971 
1972  (void) SCIPsnprintf(linname, 255, "soc#%s#B%d", basename, 0);
1973  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, -alpha2*offset2, SCIPinfinity(scip),
1979  SCIP_CALL( SCIPaddCons(scip, lincons) );
1980  SCIPdebugPrintCons(scip, lincons, NULL);
1981  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
1982  ++(*naddconss);
1983  }
1984  else
1985  { /* x2 == NULL -> b_1 >= |alpha2*offset2| */
1986  SCIP_Bool infeas;
1987  SCIP_Bool tightened;
1988  SCIP_CALL( SCIPtightenVarLb(scip, bvars[1], ABS(alpha2 * offset2), TRUE, &infeas, &tightened) );
1989  if( infeas == TRUE )
1990  {
1991  SCIPwarningMessage(scip, "creating glineur outer approximation of SOC3 constraint found problem infeasible.\n");
1992  }
1993  }
1994 
1995  /* create intermediate linear constraints */
1996  val = M_PI;
1997  for( i = 1; i < N; ++i )
1998  {
1999  val /= 2.0;
2000 
2001  vars[0] = avars[i];
2002  vals[0] = cos(val);
2003  vars[1] = bvars[i];
2004  vals[1] = sin(val);
2005  vars[2] = avars[i+1]; /*lint !e679*/
2006  vals[2] = -1.0;
2007 
2008  (void) SCIPsnprintf(linname, 255, "soc#%s#a%d", basename, i);
2009  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 3, vars, vals, 0.0, 0.0,
2015  SCIP_CALL( SCIPaddCons(scip, lincons) );
2016  SCIPdebugPrintCons(scip, lincons, NULL);
2017  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2018  ++(*naddconss);
2019 
2020  vars[0] = avars[i];
2021  vals[0] = -sin(val);
2022  vars[1] = bvars[i];
2023  vals[1] = cos(val);
2024  vars[2] = bvars[i+1]; /*lint !e679*/
2025  vals[2] = -1.0;
2026 
2027  (void) SCIPsnprintf(linname, 255, "soc#%s#b%d", basename, i);
2028  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 3, vars, vals, -SCIPinfinity(scip), 0.0,
2034  SCIP_CALL( SCIPaddCons(scip, lincons) );
2035  SCIPdebugPrintCons(scip, lincons, NULL);
2036  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2037  ++(*naddconss);
2038 
2039  vars[0] = avars[i];
2040  vals[0] = -sin(val);
2041  vars[1] = bvars[i];
2042  vals[1] = cos(val);
2043  vars[2] = bvars[i+1]; /*lint !e679*/
2044  vals[2] = 1.0;
2045 
2046  (void) SCIPsnprintf(linname, 255, "soc#%s#B%d", basename, i);
2047  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 3, vars, vals, 0.0, SCIPinfinity(scip),
2053  SCIP_CALL( SCIPaddCons(scip, lincons) );
2054  SCIPdebugPrintCons(scip, lincons, NULL);
2055  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2056  ++(*naddconss);
2057  }
2058 
2059  /* create last linear constraint */
2060  val /= 2.0;
2061  vars[0] = avars[N];
2062  vals[0] = -cos(val);
2063  vars[1] = bvars[N];
2064  vals[1] = -sin(val);
2065  vars[2] = x3;
2066  vals[2] = alpha3;
2067 
2068  (void) SCIPsnprintf(linname, 255, "soc#%s#a%d", basename, N);
2069  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 3, vars, vals, -alpha3*offset3, -alpha3*offset3,
2075  SCIP_CALL( SCIPaddCons(scip, lincons) );
2076  SCIPdebugPrintCons(scip, lincons, NULL);
2077  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2078  ++(*naddconss);
2079 
2080  for( i = 1; i <= N; ++i )
2081  {
2082  SCIP_CALL( SCIPreleaseVar(scip, &avars[i]) );
2083  SCIP_CALL( SCIPreleaseVar(scip, &bvars[i]) );
2084  }
2085  SCIPfreeBufferArray(scip, &avars);
2086  SCIPfreeBufferArray(scip, &bvars);
2087 
2088  return SCIP_OKAY;
2089 }
2090 
2091 /** adds the linear outer-approximation of Ben-Tal and Nemirovski for a SOC constraint of dimension 3
2092  *
2093  * Input is the data for a constraint \f$\sqrt{constant + (\alpha_1(x_1+offset1))^2 + (\alpha_2(x_2+offset2))^2)} \leq \alpha_3(x_3+offset3)\f$.
2094  * Here \f$\alpha3 > 0.0\f$, and the lower bound of \f$x_3 \geq -offset3\f$.
2095  * Also x2 = NULL is allowed, in which case the second term is assumed to be constant, and \f$offset2 \neq 0\f$ is needed.
2096  * */
2097 static
2099  SCIP* scip, /**< SCIP data structure */
2100  SCIP_CONS* cons, /**< original constraint */
2101  SCIP_VAR* x1, /**< variable x1 */
2102  SCIP_VAR* x2, /**< variable x2 */
2103  SCIP_VAR* x3, /**< variable x3 */
2104  SCIP_Real alpha1, /**< coefficient of x1 */
2105  SCIP_Real alpha2, /**< coefficient of x2 */
2106  SCIP_Real alpha3, /**< coefficient of x3 */
2107  SCIP_Real offset1, /**< offset of x1 */
2108  SCIP_Real offset2, /**< offset of x2 */
2109  SCIP_Real offset3, /**< offset of x3 */
2110  int N, /**< size of linear approximation, need to be >= 1 */
2111  const char* basename, /**< string to use for building variable and constraint names */
2112  int* naddconss /**< buffer where to add the number of added constraints */
2113  )
2114 {
2115  SCIP_CONS* lincons;
2116  SCIP_VAR* vars[3];
2117  SCIP_Real vals[3];
2118  char varname[255];
2119  char linname[255];
2120  int i;
2121  SCIP_VAR** avars;
2122  SCIP_VAR** bvars;
2123 
2124  assert(scip != NULL);
2125  assert(cons != NULL);
2126  assert(x1 != NULL);
2127  assert(x2 != NULL || !SCIPisZero(scip, offset2));
2128  assert(x3 != NULL);
2129  assert(SCIPisPositive(scip, alpha3));
2130  assert(SCIPisGE(scip, SCIPconsIsLocal(cons) ? SCIPvarGetLbLocal(x3) : SCIPvarGetLbGlobal(x3), -offset3));
2131  assert(basename != NULL);
2132  assert(N >= 1);
2133  assert(naddconss != NULL);
2134 
2135  SCIPdebugMsg(scip, "Creating linear Ben-Tal Nemirovski outer-approximation for <%s>.\n", basename);
2136 
2137  SCIP_CALL( SCIPallocBufferArray(scip, &avars, N+1) );
2138  SCIP_CALL( SCIPallocBufferArray(scip, &bvars, N+1) );
2139 
2140  /* create additional variables */
2141  for( i = 0; i <= N; ++i )
2142  {
2143  (void) SCIPsnprintf(varname, 255, "soc#%s_a%d", basename, i);
2144  SCIP_CALL( SCIPcreateVar(scip, &avars[i], varname, 0.0, SCIPinfinity(scip), 0.0,
2146  SCIP_CALL( SCIPaddVar(scip, avars[i]) );
2147 
2148  (void) SCIPsnprintf(varname, 255, "soc#%s_b%d", basename, i);
2149  SCIP_CALL( SCIPcreateVar(scip, &bvars[i], varname, 0.0, SCIPinfinity(scip), 0.0,
2151  SCIP_CALL( SCIPaddVar(scip, bvars[i]) );
2152  }
2153 
2154  /* create first linear constraints - split into two because of the absolute value */
2155  vars[0] = avars[0];
2156  vals[0] = 1.0;
2157  vars[1] = x1;
2158  vals[1] = -alpha1;
2159 
2160  (void) SCIPsnprintf(linname, 255, "soc#%s#a%d", basename, 0);
2161  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, alpha1 * offset1, SCIPinfinity(scip),
2166  TRUE /* removable */, SCIPconsIsStickingAtNode(cons)) );
2167  SCIP_CALL( SCIPaddCons(scip, lincons) );
2168  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2169  ++(*naddconss);
2170 
2171  vars[0] = avars[0];
2172  vals[0] = 1.0;
2173  vars[1] = x1;
2174  vals[1] = alpha1;
2175 
2176  (void) SCIPsnprintf(linname, 255, "soc#%s#A%d", basename, 0);
2177  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, -alpha1 * offset1, SCIPinfinity(scip),
2182  TRUE /* removable */, SCIPconsIsStickingAtNode(cons)) );
2183  SCIP_CALL( SCIPaddCons(scip, lincons) );
2184  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2185  ++(*naddconss);
2186 
2187  if( x2 != NULL )
2188  {
2189  vars[0] = bvars[0];
2190  vals[0] = 1.0;
2191  vars[1] = x2;
2192  vals[1] = -alpha2;
2193 
2194  (void) SCIPsnprintf(linname, 255, "soc#%s#b%d", basename, 0);
2195  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, alpha2 * offset2, SCIPinfinity(scip),
2200  TRUE /* removable */, SCIPconsIsStickingAtNode(cons)) );
2201  SCIP_CALL( SCIPaddCons(scip, lincons) );
2202  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2203  ++(*naddconss);
2204 
2205  vars[0] = bvars[0];
2206  vals[0] = 1.0;
2207  vars[1] = x2;
2208  vals[1] = alpha2;
2209 
2210  (void) SCIPsnprintf(linname, 255, "soc#%s#B%d", basename, 0);
2211  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, -alpha2 * offset2, SCIPinfinity(scip),
2216  TRUE /* removable */, SCIPconsIsStickingAtNode(cons)) );
2217  SCIP_CALL( SCIPaddCons(scip, lincons) );
2218  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2219  ++(*naddconss);
2220  }
2221  else
2222  { /* second summand is just a constant */
2223  if( SCIPconsIsLocal(cons) )
2224  {
2225  SCIP_CALL( SCIPchgVarLbNode(scip, NULL, bvars[0], ABS(alpha2 * offset2)) );
2226  }
2227  else
2228  {
2229  SCIP_CALL( SCIPchgVarLbGlobal(scip, bvars[0], ABS(alpha2 * offset2)) );
2230  }
2231  }
2232 
2233  /* create intermediate linear constraints */
2234  for( i = 1; i <= N; ++i )
2235  {
2236  SCIP_Real val;
2237 
2238  val = M_PI / pow(2.0, (double) (i+1));
2239 
2240  vars[0] = avars[i-1];
2241  vals[0] = cos(val);
2242  vars[1] = bvars[i-1];
2243  vals[1] = sin(val);
2244  vars[2] = avars[i];
2245  vals[2] = -1.0;
2246 
2247  (void) SCIPsnprintf(linname, 255, "soc#%s#a%d", basename, i);
2248  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 3, vars, vals, 0.0, 0.0,
2253  TRUE /* removable */, SCIPconsIsStickingAtNode(cons)) );
2254  SCIP_CALL( SCIPaddCons(scip, lincons) );
2255  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2256  ++(*naddconss);
2257 
2258  vars[0] = avars[i-1];
2259  vals[0] = sin(val);
2260  vars[1] = bvars[i-1];
2261  vals[1] = -cos(val);
2262  vars[2] = bvars[i];
2263  vals[2] = 1.0;
2264 
2265  (void) SCIPsnprintf(linname, 255, "soc#%s#b%d", basename, i);
2266  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 3, vars, vals, 0.0, SCIPinfinity(scip),
2271  TRUE /* removable */, SCIPconsIsStickingAtNode(cons)) );
2272  SCIP_CALL( SCIPaddCons(scip, lincons) );
2273  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2274  ++(*naddconss);
2275 
2276  vars[0] = avars[i-1];
2277  vals[0] = -sin(val);
2278  vars[1] = bvars[i-1];
2279  vals[1] = cos(val);
2280  vars[2] = bvars[i];
2281  vals[2] = 1.0;
2282 
2283  (void) SCIPsnprintf(linname, 255, "soc#%s#B%d", basename, i);
2284  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 3, vars, vals, 0.0, SCIPinfinity(scip),
2289  TRUE /* removable */, SCIPconsIsStickingAtNode(cons)) );
2290  SCIP_CALL( SCIPaddCons(scip, lincons) );
2291  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2292  ++(*naddconss);
2293  }
2294 
2295  /* create last linear constraints */
2296  vars[0] = x3;
2297  vals[0] = alpha3;
2298  vars[1] = avars[N];
2299  vals[1] = -1.0;
2300 
2301  (void) SCIPsnprintf(linname, 255, "soc#%s#a%d", basename, N);
2302  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, -alpha3 * offset3, SCIPinfinity(scip),
2308  SCIP_CALL( SCIPaddCons(scip, lincons) );
2309  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2310  ++(*naddconss);
2311 
2312  vars[0] = avars[N];
2313  vals[0] = tan( M_PI / pow(2.0, (double) (N+1)) );
2314  vars[1] = bvars[N];
2315  vals[1] = -1.0;
2316 
2317  (void) SCIPsnprintf(linname, 255, "soc#%s#b%d", basename, i);
2318  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, 0.0, SCIPinfinity(scip),
2323  TRUE /* removable */, SCIPconsIsStickingAtNode(cons)) );
2324  SCIP_CALL( SCIPaddCons(scip, lincons) );
2325  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2326  ++(*naddconss);
2327 
2328  for( i = 0; i <= N; ++i )
2329  {
2330  SCIP_CALL( SCIPreleaseVar(scip, &avars[i]) );
2331  SCIP_CALL( SCIPreleaseVar(scip, &bvars[i]) );
2332  }
2333  SCIPfreeBufferArray(scip, &avars);
2334  SCIPfreeBufferArray(scip, &bvars);
2335 
2336  return SCIP_OKAY;
2337 }
2338 
2339 /** adds a linear outer approx for a three dimensional SOC constraint
2340  *
2341  * chooses between Ben-Tan/Nemirovski and Glineur and calls the corresponding function
2342  */
2343 static
2345  SCIP* scip, /**< SCIP data structure */
2346  SCIP_CONS* cons, /**< original constraint */
2347  SCIP_VAR* x1, /**< variable x1 */
2348  SCIP_VAR* x2, /**< variable x2 */
2349  SCIP_VAR* x3, /**< variable x3 */
2350  SCIP_Real alpha1, /**< coefficient of x1 */
2351  SCIP_Real alpha2, /**< coefficient of x2 */
2352  SCIP_Real alpha3, /**< coefficient of x3 */
2353  SCIP_Real offset1, /**< offset of x1 */
2354  SCIP_Real offset2, /**< offset of x2 */
2355  SCIP_Real offset3, /**< offset of x3 */
2356  int N, /**< size of linear approximation, need to be >= 1 */
2357  SCIP_Bool glineur, /**< whether to prefer Glineur to Ben-Tal Nemirovski */
2358  const char* basename, /**< string to use for building variable and constraint names */
2359  int* naddconss /**< buffer where to add the number of added constraints */
2360  )
2361 {
2362  if( glineur )
2363  {
2364  SCIP_CALL( presolveCreateGlineurApproxDim3(scip, cons, x1, x2, x3, alpha1, alpha2, alpha3, offset1, offset2, offset3, N, basename, naddconss) );
2365  }
2366  else
2367  {
2368  SCIP_CALL( presolveCreateBenTalNemirovskiApproxDim3(scip, cons, x1, x2, x3, alpha1, alpha2, alpha3, offset1, offset2, offset3, N, basename, naddconss) );
2369  }
2370 
2371  return SCIP_OKAY;
2372 }
2373 
2374 /** adds linear outer approximation of Ben-Tal and Nemirovski for a constraint \f$\gamma + \sum_{i=1}^n (\alpha_i (x_i + \beta_i))^2 \leq (\alpha_{n+1} (x_{n+1} + \beta_{n+1}))^2\f$ to the LP
2375  *
2376  * if n > 2, calls same function recursively;
2377  * if n = 2, calls presolveCreateBenTalNemirovskiApproxDim3
2378  */
2379 static
2381  SCIP* scip, /**< SCIP data structure */
2382  int nlhsvars, /**< number of variables on left hand side (n) */
2383  SCIP_VAR** lhsvars, /**< variables on left hand side (x_i) */
2384  SCIP_Real* lhscoefs, /**< coefficients of variables on left hand side (alpha_i) */
2385  SCIP_Real* lhsoffsets, /**< offsets of variable on left hand side (beta_i) */
2386  SCIP_VAR* rhsvar, /**< variable on right hand side (y) */
2387  SCIP_Real rhscoeff, /**< coefficient of variable on right hand side (alpha_{n+1}) */
2388  SCIP_Real rhsoffset, /**< offset of variable on right hand side (beta_{n+1}) */
2389  SCIP_Real constant, /**< constant term (gamma) */
2390  const char* basename, /**< prefix for variable and constraint name */
2391  SCIP_CONS* origcons, /**< original constraint for which this SOC3 set is added */
2392  int soc3_nr_auxvars, /**< number of auxiliary variables to use for a SOC3 constraint, or 0 if automatic */
2393  SCIP_Bool glineur, /**< whether Glineur should be preferred to Ben-Tal Nemirovski */
2394  int* naddconss /**< buffer where to add the number of added constraints */
2395  )
2396 {
2397  char name[255];
2398  SCIP_VAR* auxvar1;
2399  SCIP_VAR* auxvar2;
2400 
2401  assert(scip != NULL);
2402  assert(lhsvars != NULL);
2403  assert(nlhsvars >= 1);
2404  assert(lhscoefs != NULL);
2405  assert(lhsoffsets != NULL);
2406  assert(rhsvar != NULL);
2407  assert(basename != NULL);
2408  assert(!SCIPisNegative(scip, constant));
2409  assert(naddconss != NULL);
2410 
2411  if( nlhsvars == 1 )
2412  { /* end of recursion */
2413  assert(SCIPisPositive(scip, constant));
2414  SCIP_CALL( presolveCreateOuterApproxDim3(scip, origcons,
2415  lhsvars[0], NULL, rhsvar,
2416  lhscoefs[0], 1.0, rhscoeff,
2417  lhsoffsets[0], sqrt(constant), rhsoffset,
2418  soc3_nr_auxvars, glineur, basename, naddconss) );
2419 
2420  return SCIP_OKAY;
2421  }
2422 
2423  if( nlhsvars == 2 && SCIPisZero(scip, constant) )
2424  { /* end of recursion */
2425  assert(lhsvars[0] != NULL);
2426  assert(lhsvars[1] != NULL);
2427  assert(rhsvar != NULL);
2428  SCIP_CALL( presolveCreateOuterApproxDim3(scip, origcons,
2429  lhsvars[0], lhsvars[1], rhsvar,
2430  lhscoefs[0], lhscoefs[1], rhscoeff,
2431  lhsoffsets[0], lhsoffsets[1], rhsoffset,
2432  soc3_nr_auxvars, glineur, basename, naddconss) );
2433 
2434  return SCIP_OKAY;
2435  }
2436 
2437  if( nlhsvars == 3 || (nlhsvars == 2 && !SCIPisZero(scip, constant)) )
2438  {
2439  /* a bit special case too */
2440  /* for first two variables on lhs, create a new aux.var and a new SOC3 */
2441  (void) SCIPsnprintf(name, 255, "%s#z1", basename);
2442  SCIP_CALL( SCIPcreateVar(scip, &auxvar1, name, 0.0, SCIPinfinity(scip), 0.0,
2444  SCIP_CALL( SCIPaddVar(scip, auxvar1) );
2445 
2446  /* constraint alpha_0 (x_0+beta0)^2 + alpha_1 (x_1+beta1)^2 <= auxvar^2 */
2447  SCIP_CALL( presolveCreateOuterApproxDim3(scip, origcons,
2448  lhsvars[0], lhsvars[1], auxvar1,
2449  lhscoefs[0], lhscoefs[1], 1.0,
2450  lhsoffsets[0], lhsoffsets[1], 0.0,
2451  soc3_nr_auxvars, glineur, name, naddconss) );
2452 
2453  (void) SCIPsnprintf(name, 255, "%s_soc3", basename);
2454  if( nlhsvars == 3 )
2455  { /* create new constraint alpha_2 (x_2+beta2)^2 + auxvar^2 <= (rhscoeff * (rhsvar+rhsoffset))^2 */
2456  SCIP_CALL( presolveCreateOuterApproxDim3(scip, origcons,
2457  lhsvars[2], auxvar1, rhsvar,
2458  lhscoefs[2], 1.0, rhscoeff,
2459  lhsoffsets[2], 0.0, rhsoffset,
2460  soc3_nr_auxvars, glineur, name, naddconss) );
2461  }
2462  else
2463  { /* create new constraint auxvar^2 + sqrt(constant)^2 <= (rhscoeff * (rhsvar+rhsoffset))^2 */
2464  SCIP_CALL( presolveCreateOuterApproxDim3(scip, origcons,
2465  auxvar1, NULL, rhsvar,
2466  1.0, 1.0, rhscoeff,
2467  0.0, sqrt(constant), rhsoffset,
2468  soc3_nr_auxvars, glineur, name, naddconss) );
2469  }
2470 
2471  SCIP_CALL( SCIPreleaseVar(scip, &auxvar1) );
2472 
2473  return SCIP_OKAY;
2474  }
2475 
2476  /* nlhsvars >= 4 */
2477 
2478  (void) SCIPsnprintf(name, 255, "%s#z1", basename);
2479  SCIP_CALL( SCIPcreateVar(scip, &auxvar1, name, 0.0, SCIPinfinity(scip), 0.0,
2481  SCIP_CALL( SCIPaddVar(scip, auxvar1) );
2482 
2483  /* approx for left half of lhs */
2485  nlhsvars/2, lhsvars, lhscoefs, lhsoffsets,
2486  auxvar1, 1.0, 0.0,
2487  constant, name, origcons, soc3_nr_auxvars, glineur, naddconss) );
2488 
2489  (void) SCIPsnprintf(name, 255, "%s#z2", basename);
2490  SCIP_CALL( SCIPcreateVar(scip, &auxvar2, name, 0., SCIPinfinity(scip), 0.0,
2492  SCIP_CALL( SCIPaddVar(scip, auxvar2) );
2493 
2494  /* approx for right half of lhs */
2496  nlhsvars-nlhsvars/2, &lhsvars[nlhsvars/2], &lhscoefs[nlhsvars/2], &lhsoffsets[nlhsvars/2],
2497  auxvar2, 1.0, 0.0,
2498  0.0, name, origcons, soc3_nr_auxvars, glineur, naddconss) );
2499 
2500  /* SOC constraint binding both auxvar's */
2501  (void)SCIPsnprintf(name, 255, "%s_soc3", basename);
2502  SCIP_CALL( presolveCreateOuterApproxDim3(scip, origcons,
2503  auxvar1, auxvar2, rhsvar,
2504  1.0, 1.0, rhscoeff,
2505  0.0, 0.0, rhsoffset,
2506  soc3_nr_auxvars, glineur, name, naddconss) );
2507 
2508  SCIP_CALL( SCIPreleaseVar(scip, &auxvar1) );
2509  SCIP_CALL( SCIPreleaseVar(scip, &auxvar2) );
2510 
2511  return SCIP_OKAY;
2512 }
2513 
2514 /** propagates variable bounds */
2515 static
2517  SCIP* scip, /**< SCIP data structure */
2518  SCIP_CONS* cons, /**< constraint */
2519  SCIP_RESULT* result, /**< buffer to store result of propagation */
2520  int* nchgbds, /**< buffer where to add number of tightened bounds */
2521  SCIP_Bool* redundant /**< buffer to indicate whether constraint was marked for deletion because of redundancy */
2522  )
2523 {
2524  SCIP_CONSDATA* consdata;
2525  SCIP_INTERVAL lhsrange;
2526  SCIP_INTERVAL lhsrange_squared;
2527  SCIP_INTERVAL* lhsranges;
2528  SCIP_INTERVAL rhsrange;
2529  SCIP_INTERVAL a, b, c;
2530  SCIP_ROUNDMODE roundmode;
2531  SCIP_Bool infeas, tightened;
2532  int i;
2533  SCIP_Real lb, ub;
2534 
2535  assert(scip != NULL);
2536  assert(cons != NULL);
2537  assert(result != NULL);
2538  assert(nchgbds != NULL);
2539  assert(redundant != NULL);
2540 
2541  consdata = SCIPconsGetData(cons);
2542  assert(consdata != NULL);
2543  assert(consdata->constant >= 0.0);
2544 
2545  *redundant = FALSE;
2546 
2547  if( !SCIPconsIsMarkedPropagate(cons) )
2548  {
2549  SCIPdebugMsg(scip, "skip propagation for constraint %s\n", SCIPconsGetName(cons));
2550  *result = SCIP_DIDNOTRUN;
2551  return SCIP_OKAY;
2552  }
2553  else
2554  {
2555  SCIPdebugMsg(scip, "try propagation for constraint %s\n", SCIPconsGetName(cons));
2556  }
2557 
2558  *result = SCIP_DIDNOTFIND;
2559  SCIP_CALL( SCIPunmarkConsPropagate(scip, cons) );
2560 
2561  /* @todo do something clever to decide whether propagation should be tried */
2562 
2563  /* compute activity on lhs */
2564  SCIPintervalSet(&lhsrange, consdata->constant);
2565  SCIP_CALL( SCIPallocBufferArray(scip, &lhsranges, consdata->nvars) );
2566  for( i = 0; i < consdata->nvars; ++i )
2567  {
2568  lb = SCIPcomputeVarLbLocal(scip, consdata->vars[i]) - SCIPepsilon(scip);
2569  ub = SCIPcomputeVarUbLocal(scip, consdata->vars[i]) + SCIPepsilon(scip);
2570  SCIPintervalSetBounds(&lhsranges[i], MIN(lb, ub), MAX(lb, ub));
2571  if( consdata->offsets[i] != 0.0 )
2572  SCIPintervalAddScalar(SCIPinfinity(scip), &lhsranges[i], lhsranges[i], consdata->offsets[i]);
2573  if( consdata->coefs[i] != 1.0 )
2574  SCIPintervalMulScalar(SCIPinfinity(scip), &lhsranges[i], lhsranges[i], consdata->coefs[i]);
2575  SCIPintervalSquare(SCIPinfinity(scip), &lhsranges[i], lhsranges[i]);
2576 
2577  SCIPintervalAdd(SCIPinfinity(scip), &lhsrange, lhsrange, lhsranges[i]);
2578  }
2579  assert(lhsrange.inf >= 0); /* a sum of squares plus positive constant should be non-negative */
2580  lhsrange_squared = lhsrange; /* we will need the squared version later */
2581  SCIPintervalSquareRoot(SCIPinfinity(scip), &lhsrange, lhsrange);
2582 
2583  /* compute activity on rhs: rhsrange = rhscoeff * (rhsvar + rhsoffset) */
2584  lb = SCIPcomputeVarLbLocal(scip, consdata->rhsvar) - SCIPepsilon(scip);
2585  ub = SCIPcomputeVarUbLocal(scip, consdata->rhsvar) + SCIPepsilon(scip);
2586  SCIPintervalSetBounds(&rhsrange, MIN(lb, ub), MAX(lb, ub));
2587 
2588  if( consdata->rhsoffset != 0.0 )
2589  SCIPintervalAddScalar(SCIPinfinity(scip), &rhsrange, rhsrange, consdata->rhsoffset);
2590  if( consdata->rhscoeff != 1.0 )
2591  SCIPintervalMulScalar(SCIPinfinity(scip), &rhsrange, rhsrange, consdata->rhscoeff);
2592 
2593  /* check for infeasibility */
2594  if( SCIPisGT(scip, lhsrange.inf-SCIPfeastol(scip), rhsrange.sup) )
2595  {
2596  SCIPdebugMsg(scip, "propagation found constraint <%s> infeasible: lhs = [%.15g,%.15g]-feastol-eps > rhs = [%.15g,%.15g]\n",
2597  SCIPconsGetName(cons), lhsrange.inf, lhsrange.sup, rhsrange.inf, rhsrange.sup);
2598  *result = SCIP_CUTOFF;
2599  goto TERMINATE;
2600  }
2601 
2602  /* check for redundancy: max(lhsrange) <= min(rhsrange) */
2603  if( SCIPisLE(scip, lhsrange.sup, rhsrange.inf) )
2604  {
2605  SCIPdebugMsg(scip, "propagation found constraint <%s> redundant: lhs = [%.15g,%.15g] <= rhs = [%.15g,%.15g]\n",
2606  SCIPconsGetName(cons), lhsrange.inf, lhsrange.sup, rhsrange.inf, rhsrange.sup);
2607  SCIP_CALL( SCIPdelConsLocal(scip, cons) );
2608  goto TERMINATE;
2609  }
2610 
2611  /* try to tighten variable on rhs */
2612  if( SCIPvarGetStatus(consdata->rhsvar) != SCIP_VARSTATUS_MULTAGGR )
2613  {
2614  /* we have lhsrange <= rhscoeff * (rhsvar + rhsoffset)
2615  * if rhscoeff > 0, then lhsrange / rhscoeff - rhsoffset <= rhsvar
2616  * if rhscoeff < 0, then lhsrange / rhscoeff - rhsoffset >= rhsvar
2617  */
2618  a = lhsrange;
2619  if( consdata->rhscoeff != 1.0 )
2620  SCIPintervalDivScalar(SCIPinfinity(scip), &a, a, consdata->rhscoeff);
2621  if( consdata->rhsoffset != 0.0 )
2622  SCIPintervalSubScalar(SCIPinfinity(scip), &a, a, consdata->rhsoffset);
2623 
2624  if( consdata->rhscoeff > 0.0 )
2625  {
2626  SCIP_CALL( SCIPtightenVarLb(scip, consdata->rhsvar, SCIPintervalGetInf(a), FALSE, &infeas, &tightened) );
2627  }
2628  else
2629  {
2630  SCIP_CALL( SCIPtightenVarUb(scip, consdata->rhsvar, SCIPintervalGetSup(a), FALSE, &infeas, &tightened) );
2631  }
2632  if( infeas )
2633  {
2634  SCIPdebugMsg(scip, "propagation found constraint <%s> infeasible\n", SCIPconsGetName(cons));
2635  *result = SCIP_CUTOFF;
2636  goto TERMINATE;
2637  }
2638  if( tightened )
2639  {
2640  SCIPdebugMsg(scip, "propagation tightened bounds of rhs variable <%s> in constraint <%s>\n", SCIPvarGetName(consdata->rhsvar), SCIPconsGetName(cons));
2641  *result = SCIP_REDUCEDDOM;
2642  ++*nchgbds;
2643  }
2644  }
2645 
2646  /* try to tighten variables on lhs:
2647  * (coefs[i] * (vars[i] + offset[i]))^2 <= sqr(rhsrange) - (constant + sum_{j != i} (coefs[j] * (vars[j] + offset[j]))^2)
2648  *
2649  * first, set b = sqr(rhsrange) - (constant + sum_i (coefs[i] * (vars[i] + offset[i]))^2)
2650  * then, for each i, we undo the subtraction of (coefs[i] * (vars[i] + offset[i]))^2 in b and take a square root
2651  * thus, we get a = sqrt(sqr(rhsrange) - (constant + sum_{j != i} (coefs[j] * (vars[j] + offset[j]))^2))
2652  * and |coefs[i] * (vars[i] + offset[i])| <= a
2653  * this gives us the two inequalities
2654  * vars[i] <= sqrt(b)/coefs[i] - offset[i]
2655  * vars[i] >= -sqrt(b)/coefs[i] - offset[i]
2656  * (note that coefs[i] >= 0)
2657  */
2658  SCIPintervalSquare(SCIPinfinity(scip), &b, rhsrange);
2659  SCIPintervalSub(SCIPinfinity(scip), &b, b, lhsrange_squared); /*lint !e644 */
2660  for( i = 0; i < consdata->nvars; ++i )
2661  {
2662  if( SCIPvarGetStatus(consdata->vars[i]) == SCIP_VARSTATUS_MULTAGGR )
2663  continue;
2664 
2665  roundmode = SCIPintervalGetRoundingMode();
2666  if( !SCIPisInfinity(scip, b.sup) )
2667  {
2669  a.sup = b.sup + lhsranges[i].inf;
2670  }
2671  else
2672  {
2673  a.sup = SCIPinfinity(scip);
2674  }
2675  if( !SCIPisInfinity(scip, -b.inf) )
2676  {
2678  a.inf = b.inf + lhsranges[i].sup;
2679  }
2680  else
2681  {
2682  a.inf = -SCIPinfinity(scip);
2683  }
2684  SCIPintervalSetRoundingMode(roundmode);
2685  SCIPintervalSquareRoot(SCIPinfinity(scip), &a, a);
2686 
2687  assert(consdata->coefs[i] >= 0.0); /* should be ensured in create and presolveRemoveFixed */
2688 
2689  c = a;
2690  if( consdata->coefs[i] != 1.0 )
2691  SCIPintervalDivScalar(SCIPinfinity(scip), &c, c, consdata->coefs[i]);
2692  if( consdata->offsets[i] != 0.0 )
2693  SCIPintervalSubScalar(SCIPinfinity(scip), &c, c, consdata->offsets[i]);
2694 
2695  SCIP_CALL( SCIPtightenVarUb(scip, consdata->vars[i], SCIPintervalGetSup(c), FALSE, &infeas, &tightened) );
2696  if( infeas )
2697  {
2698  SCIPdebugMsg(scip, "propagation found constraint <%s> infeasible\n", SCIPconsGetName(cons));
2699  *result = SCIP_CUTOFF;
2700  goto TERMINATE;
2701  }
2702  if( tightened )
2703  {
2704  SCIPdebugMsg(scip, "propagation tightened bounds of lhs variable <%s> in constraint <%s>\n", SCIPvarGetName(consdata->vars[i]), SCIPconsGetName(cons));
2705  *result = SCIP_REDUCEDDOM;
2706  ++*nchgbds;
2707  }
2708 
2709  c = a;
2710  SCIPintervalDivScalar(SCIPinfinity(scip), &c, c, -consdata->coefs[i]);
2711  if( consdata->offsets[i] != 0.0 )
2712  SCIPintervalSubScalar(SCIPinfinity(scip), &c, c, consdata->offsets[i]);
2713 
2714  SCIP_CALL( SCIPtightenVarLb(scip, consdata->vars[i], SCIPintervalGetInf(c), FALSE, &infeas, &tightened) );
2715  if( infeas )
2716  {
2717  SCIPdebugMsg(scip, "propagation found constraint <%s> infeasible\n", SCIPconsGetName(cons));
2718  *result = SCIP_CUTOFF;
2719  goto TERMINATE;
2720  }
2721  if( tightened )
2722  {
2723  SCIPdebugMsg(scip, "propagation tightened bounds of lhs variable <%s> in constraint <%s>\n", SCIPvarGetName(consdata->vars[i]), SCIPconsGetName(cons));
2724  *result = SCIP_REDUCEDDOM;
2725  ++*nchgbds;
2726  }
2727  }
2728 
2729 TERMINATE:
2730  SCIPfreeBufferArray(scip, &lhsranges);
2731 
2732  if( *result != SCIP_DIDNOTFIND )
2733  {
2734  SCIP_CALL( SCIPresetConsAge(scip, cons) );
2735  }
2736 
2737  return SCIP_OKAY;
2738 }
2739 
2740 /** tries to adjust a solution such that it satisfies a given constraint by increasing the value for the constraints right hand side variable */
2741 static
2743  SCIP* scip, /**< SCIP data structure */
2744  SCIP_CONS* cons, /**< constraint */
2745  SCIP_SOL* sol, /**< solution to polish */
2746  SCIP_Bool* success /**< buffer to store whether polishing was successful */
2747  )
2748 {
2749  SCIP_CONSDATA* consdata;
2750  SCIP_Real rhsval;
2751 
2752  assert(scip != NULL);
2753  assert(cons != NULL);
2754  assert(sol != NULL);
2755  assert(success != NULL);
2756 
2757  consdata = SCIPconsGetData(cons);
2758  assert(consdata != NULL);
2759  assert(!SCIPisZero(scip, consdata->rhscoeff));
2760 
2761  /* compute minimal rhs variable value so that constraint is satisfied */
2762  if( !SCIPisInfinity(scip, consdata->lhsval) )
2763  rhsval = consdata->lhsval / consdata->rhscoeff - consdata->rhsoffset;
2764  else
2765  rhsval = consdata->rhscoeff > 0.0 ? SCIPinfinity(scip) : -SCIPinfinity(scip);
2766 
2767  if( consdata->rhscoeff > 0.0 )
2768  {
2769  assert(SCIPvarMayRoundUp(consdata->rhsvar));
2770 
2771  /* round rhsval up, if variable is integral */
2772  if( SCIPvarIsIntegral(consdata->rhsvar) && !SCIPisInfinity(scip, rhsval) )
2773  rhsval = SCIPceil(scip, rhsval);
2774 
2775  /* if new value is above upper bound, we are lost */
2776  if( SCIPisGT(scip, rhsval, SCIPvarGetUbGlobal(consdata->rhsvar)) )
2777  {
2778  *success = FALSE;
2779  }
2780  else
2781  {
2782  /* if new value is larger then current one, increase to new value */
2783  if( rhsval > SCIPgetSolVal(scip, sol, consdata->rhsvar) )
2784  {
2785  SCIPdebugMsg(scip, "increase <%s> to %g\n", SCIPvarGetName(consdata->rhsvar), rhsval);
2786  SCIP_CALL( SCIPsetSolVal(scip, sol, consdata->rhsvar, rhsval) );
2787  }
2788 
2789  *success = TRUE;
2790  }
2791  }
2792  else
2793  {
2794  assert(SCIPvarMayRoundDown(consdata->rhsvar));
2795 
2796  /* round rhsval down, if variable is integral */
2797  if( SCIPvarIsIntegral(consdata->rhsvar) )
2798  rhsval = SCIPfloor(scip, rhsval);
2799 
2800  /* if new value is below lower bound, we are lost */
2801  if( SCIPisLT(scip, rhsval, SCIPvarGetLbGlobal(consdata->rhsvar)) )
2802  {
2803  *success = FALSE;
2804  }
2805  else
2806  {
2807  /* if new value is below current one, decrease to new value */
2808  if( rhsval < SCIPgetSolVal(scip, sol, consdata->rhsvar) )
2809  {
2810  SCIPdebugMsg(scip, "decrease <%s> to %g\n", SCIPvarGetName(consdata->rhsvar), rhsval);
2811  SCIP_CALL( SCIPsetSolVal(scip, sol, consdata->rhsvar, rhsval) );
2812  }
2813 
2814  *success = TRUE;
2815  }
2816  }
2817 
2818  SCIPdebugMsg(scip, "polishing solution for constraint <%s> was %ssuccessful\n", SCIPconsGetName(cons), *success ? "" : "not ");
2819 
2820  return SCIP_OKAY;
2821 }
2822 
2823 /** disaggregates a (sufficiently large) SOC constraint into smaller ones; for each term on the lhs we add a quadratic
2824  * constraint \f$(\alpha_i * (x_i + \beta_i))^2 \leq \alpha_{n+1} (x_{n+1} + \beta_{n+1})\, z_i\f$ and a single linear constraint
2825  * \f$\sum_i z_i \leq \alpha_{n+1}\, (x_{n+1} + \beta_{n+1})\f$; each quadratic constraint might be upgraded to a SOC; since the
2826  * violations of all quadratic constraints sum up we scale each constraint by the number of lhs terms + 1
2827  *
2828  * @todo if rhsvar is NULL, then the disaggregation does not produce further cones. Should it then be upgraded
2829  * to a quadratic and let the quadratic desaggregate it?
2830  * The code assumes now that the rhsvar is not NULL in order build the direct SOC -> SOC disaggregation
2831  */
2832 static
2834  SCIP* scip, /**< SCIP data structure */
2835  SCIP_CONS* cons, /**< constraint */
2836  SCIP_CONSDATA* consdata, /**< constraint data */
2837  int* naddconss, /**< pointer to count total number of added constraints */
2838  int* ndelconss, /**< pointer to count total number of deleted constraints */
2839  SCIP_Bool* success /**< pointer to store whether disaggregation was successful */
2840  )
2841 {
2842  SCIP_CONS* discons;
2843  SCIP_VAR** disvars;
2844  SCIP_VAR** sumvars;
2845  SCIP_VAR** difvars;
2846  SCIP_Real* discoefs;
2847  SCIP_VAR* lhsvars[2];
2848  SCIP_VAR* aggvars[2];
2849  SCIP_Real coefs[2];
2850  SCIP_Real offsets[2];
2851  SCIP_Real scalars[2];
2852  char name[SCIP_MAXSTRLEN];
2853  SCIP_Real constant;
2854  SCIP_Real scale;
2855  SCIP_Bool infeas;
2856  int ndisvars;
2857  int i;
2858 
2859  assert(naddconss != NULL);
2860  assert(ndelconss != NULL);
2861  assert(success != NULL);
2862 
2863  *success = FALSE;
2864 
2865  /* disaggregation does not make much sense if there are too few variables */
2866  if( consdata->nvars < 3 )
2867  {
2868  SCIPdebugMsg(scip, "can not disaggregate too small soc constraint %s\n", SCIPconsGetName(cons));
2869  return SCIP_OKAY;
2870  }
2871 
2872  if( consdata->rhsvar == NULL )
2873  {
2874  SCIPdebugMsg(scip, "can not disaggregate directly into a soc without rhs var %s\n", SCIPconsGetName(cons));
2875  return SCIP_OKAY;
2876  }
2877 
2878  /* there are at most n + 2 many linear varibles */
2879  SCIP_CALL( SCIPallocBufferArray(scip, &disvars, consdata->nvars + 2) );
2880  SCIP_CALL( SCIPallocBufferArray(scip, &sumvars, consdata->nvars + 2) );
2881  SCIP_CALL( SCIPallocBufferArray(scip, &difvars, consdata->nvars + 2) );
2882  SCIP_CALL( SCIPallocBufferArray(scip, &discoefs, consdata->nvars + 2) );
2883  ndisvars = 0;
2884 
2885  scale = 1.0 * (consdata->nvars + 1)/4.0;
2886 
2887  /* add (*) (alpha_i * (x_i + beta_i))^2 <= alpha_{n+1} * (x_{n+1} + beta_{n+1}) * z_i:
2888  * create sumvar = alpha_{n+1} * (x_{n+1} + beta_{n+1}) + z_i (multiagg)
2889  * create difvar = alpha_{n+1} * (x_{n+1} + beta_{n+1}) - z_i (multiagg)
2890  * note that (*) is equiv to sqrt( (2 * alpha_i * (x_i + beta_i))^2 + difvar^2) <= sumvar
2891  * scaling give us: sqrt( (2 * scale * alpha_i * (x_i + beta_i))^2 + (scale * difvar)^2) <= scale * sumvar
2892  */
2893  aggvars[0] = consdata->rhsvar;
2894  scalars[0] = consdata->rhscoeff;
2895  for( i = 0; i < consdata->nvars; ++i )
2896  {
2897  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedis_%s_%d", SCIPvarGetName(consdata->vars[i]), i);
2898  SCIP_CALL( SCIPcreateVar(scip, &disvars[i], name, 0.0, SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS, TRUE, FALSE,
2899  NULL, NULL, NULL, NULL, NULL) );
2900  SCIP_CALL( SCIPaddVar(scip, disvars[i]) );
2901 
2902  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedisS_%s_%d", SCIPvarGetName(consdata->vars[i]), i);
2903  SCIP_CALL( SCIPcreateVar(scip, &sumvars[i], name, 0.0, SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS, TRUE, FALSE,
2904  NULL, NULL, NULL, NULL, NULL) );
2905  SCIP_CALL( SCIPaddVar(scip, sumvars[i]) );
2906 
2907  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedisD_%s_%d", SCIPvarGetName(consdata->vars[i]), i);
2908  SCIP_CALL( SCIPcreateVar(scip, &difvars[i], name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0,
2910  SCIP_CALL( SCIPaddVar(scip, difvars[i]) );
2911 
2912  aggvars[1] = disvars[i];
2913  scalars[1] = 1.0;
2914  constant = consdata->rhscoeff * consdata->rhsoffset;
2915  SCIP_CALL( SCIPmultiaggregateVar(scip, sumvars[i], 2, aggvars, scalars, constant, &infeas, success) );
2916  /* @todo what shall we do if multiagg fails? */
2917  assert(!infeas && *success);
2918 
2919  scalars[1] = -1.0;
2920  SCIP_CALL( SCIPmultiaggregateVar(scip, difvars[i], 2, aggvars, scalars, constant, &infeas, success) );
2921  assert(!infeas && *success);
2922 
2923  /* create soc */
2924  lhsvars[0] = difvars[i];
2925  coefs[0] = scale;
2926  offsets[0] = 0.0;
2927  lhsvars[1] = consdata->vars[i];
2928  coefs[1] = scale * 2 * consdata->coefs[i];
2929  offsets[1] = consdata->offsets[i];
2930 
2931  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "consdis_%s_%d", SCIPconsGetName(cons), i);
2932  SCIP_CALL( SCIPcreateConsBasicSOC(scip, &discons, name, 2, lhsvars, coefs, offsets, 0.0, sumvars[i], scale, 0.0) );
2933  SCIP_CALL( SCIPaddCons(scip, discons) );
2934 #ifdef SCIP_DEBUG
2935  SCIP_CALL( SCIPprintCons(scip, discons, NULL) );
2936 #endif
2937  SCIP_CALL( SCIPreleaseCons(scip, &discons) );
2938  ++(*naddconss);
2939 
2940  /* linear coefficient in the linear constraint */
2941  discoefs[ndisvars] = 1.0;
2942  ++ndisvars;
2943  }
2944  assert(ndisvars == consdata->nvars);
2945 
2946  /* add gamma <= alpha_{n+1} * (x_{n+1} + beta_{n+1}) * z_i
2947  * sumvar and difvar are the same as before, but the equivalent soc now is
2948  * sqrt(4 * gamma + difvar^2) <= sumvar
2949  * scaling give us: sqrt( (4 * scale^2 * gamma + (scale * difvar)^2) <= scale * sumvar
2950  */
2951  if( !SCIPisZero(scip, consdata->constant) )
2952  {
2953  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedis_const_%s", SCIPconsGetName(cons));
2954  SCIP_CALL( SCIPcreateVar(scip, &disvars[ndisvars], name, 0.0, SCIPinfinity(scip), 0.0,
2956  SCIP_CALL( SCIPaddVar(scip, disvars[ndisvars]) );
2957 
2958  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedisS_const_%s", SCIPconsGetName(cons));
2959  SCIP_CALL( SCIPcreateVar(scip, &sumvars[ndisvars], name, 0.0, SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS, TRUE, FALSE,
2960  NULL, NULL, NULL, NULL, NULL) );
2961  SCIP_CALL( SCIPaddVar(scip, sumvars[ndisvars]) );
2962 
2963  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedisD_const_%s", SCIPconsGetName(cons));
2964  SCIP_CALL( SCIPcreateVar(scip, &difvars[ndisvars], name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0,
2966  SCIP_CALL( SCIPaddVar(scip, difvars[ndisvars]) );
2967 
2968  aggvars[1] = disvars[i];
2969  scalars[1] = 1.0;
2970  constant = consdata->rhscoeff * consdata->rhsoffset;
2971  SCIP_CALL( SCIPmultiaggregateVar(scip, sumvars[i], 2, aggvars, scalars, constant, &infeas, success) );
2972  assert(!infeas && *success);
2973 
2974  scalars[1] = -1.0;
2975  SCIP_CALL( SCIPmultiaggregateVar(scip, difvars[i], 2, aggvars, scalars, constant, &infeas, success) );
2976  assert(!infeas && *success);
2977 
2978  /* create soc */
2979  lhsvars[0] = difvars[ndisvars];
2980  coefs[0] = scale;
2981  offsets[0] = 0.0;
2982  constant = 4.0 * SQR(scale) * consdata->constant;
2983  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "consdis_%s_constant", SCIPconsGetName(cons));
2984  SCIP_CALL( SCIPcreateConsBasicSOC(scip, &discons, name, 1, lhsvars, coefs, offsets, constant,
2985  sumvars[ndisvars], scale, 0.0) );
2986  SCIP_CALL( SCIPaddCons(scip, discons) );
2987  SCIP_CALL( SCIPreleaseCons(scip, &discons) );
2988  ++(*naddconss);
2989 
2990  /* linear coefficient in the linear constraint */
2991  discoefs[ndisvars] = 1.0;
2992  ++ndisvars;
2993  }
2994 
2995  /* create linear constraint sum z_i <= alpha_{n+1} * (x_{n+1} + beta_{n+1}); first add extra coefficient for the rhs */
2996  discoefs[ndisvars] = -1.0 * consdata->rhscoeff;
2997  disvars[ndisvars] = consdata->rhsvar;
2998 
2999  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "consdis_linear_%s", SCIPconsGetName(cons));
3000  SCIP_CALL( SCIPcreateConsBasicLinear(scip, &discons, name, ndisvars + 1, disvars, discoefs, -SCIPinfinity(scip),
3001  consdata->rhscoeff * consdata->rhsoffset) );
3002 
3003  SCIP_CALL( SCIPaddCons(scip, discons) );
3004  SCIP_CALL( SCIPreleaseCons(scip, &discons) );
3005  ++(*naddconss);
3006 
3007  /* release all variables */
3008  for( i = ndisvars - 1; i >= 0; --i )
3009  {
3010  SCIP_CALL( SCIPreleaseVar(scip, &disvars[i]) );
3011  SCIP_CALL( SCIPreleaseVar(scip, &sumvars[i]) );
3012  SCIP_CALL( SCIPreleaseVar(scip, &difvars[i]) );
3013  }
3014  SCIPfreeBufferArray(scip, &discoefs);
3015  SCIPfreeBufferArray(scip, &difvars);
3016  SCIPfreeBufferArray(scip, &sumvars);
3017  SCIPfreeBufferArray(scip, &disvars);
3018 
3019  /* delete constraint */
3020  SCIP_CALL( SCIPdelCons(scip, cons) );
3021  ++(*ndelconss);
3022 
3023  *success = TRUE;
3024 
3025  return SCIP_OKAY;
3026 }
3027 
3028 
3029 /** helper function to enforce constraints */
3030 static
3032  SCIP* scip, /**< SCIP data structure */
3033  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3034  SCIP_CONS** conss, /**< constraints to process */
3035  int nconss, /**< number of constraints */
3036  int nusefulconss, /**< number of useful (non-obsolete) constraints to process */
3037  SCIP_SOL* sol, /**< solution to enforce (NULL for the LP solution) */
3038  SCIP_RESULT* result /**< pointer to store the result of the enforcing call */
3039  )
3040 {
3041  SCIP_CONSHDLRDATA* conshdlrdata;
3042  SCIP_CONSDATA* consdata;
3043  SCIP_CONS* maxviolcons;
3044  SCIP_Bool success;
3045  SCIP_Bool cutoff;
3046  SCIP_Bool redundant;
3047  int nbndchg;
3048  int c;
3049 
3050  assert(scip != NULL);
3051  assert(conshdlr != NULL);
3052  assert(conss != NULL || nconss == 0);
3053  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
3054  assert(result != NULL);
3055 
3056  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3057  assert(conshdlrdata != NULL);
3058 
3059  SCIP_CALL( computeViolations(scip, conshdlr, conss, nconss, sol, &maxviolcons) );
3060 
3061  if( maxviolcons == NULL )
3062  {
3063  *result = SCIP_FEASIBLE;
3064  return SCIP_OKAY;
3065  }
3066 
3067  /* if we are above the 100'th enforcement round for this node, something is strange
3068  * (maybe the LP does not think that the cuts we add are violated, or we do ECP on a high-dimensional convex function)
3069  * in this case, check if some limit is hit or SCIP should stop for some other reason and terminate enforcement by creating a dummy node
3070  * (in optimized more, returning SCIP_INFEASIBLE in *result would be sufficient, but in debug mode this would give an assert in scip.c)
3071  * the reason to wait for 100 rounds is to avoid calls to SCIPisStopped in normal runs, which may be expensive
3072  * we only increment nenforounds until 101 to avoid an overflow
3073  */
3074  if( conshdlrdata->lastenfonode == SCIPgetCurrentNode(scip) )
3075  {
3076  if( conshdlrdata->nenforounds > 100 )
3077  {
3078  if( SCIPisStopped(scip) )
3079  {
3080  SCIP_NODE* child;
3081 
3082  SCIP_CALL( SCIPcreateChild(scip, &child, 1.0, SCIPnodeGetEstimate(SCIPgetCurrentNode(scip))) );
3083  *result = SCIP_BRANCHED;
3084 
3085  return SCIP_OKAY;
3086  }
3087  }
3088  else
3089  ++conshdlrdata->nenforounds;
3090  }
3091  else
3092  {
3093  conshdlrdata->lastenfonode = SCIPgetCurrentNode(scip);
3094  conshdlrdata->nenforounds = 0;
3095  }
3096 
3097  /* try separation, this should usually work */
3098  SCIP_CALL( separatePoint(scip, conshdlr, conss, nconss, nusefulconss, sol, TRUE, &cutoff, &success) );
3099  if( cutoff )
3100  {
3101  *result = SCIP_CUTOFF;
3102  return SCIP_OKAY;
3103  }
3104  if( success )
3105  {
3106  SCIPdebugMsg(scip, "enforced by separation\n");
3107  *result = SCIP_SEPARATED;
3108  return SCIP_OKAY;
3109  }
3110 
3111  /* try propagation */
3112  for( c = 0; c < nconss; ++c )
3113  {
3114  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
3115  if( !SCIPisGT(scip, consdata->violation, SCIPfeastol(scip)) )
3116  continue;
3117 
3118  nbndchg = 0;
3119  SCIP_CALL( propagateBounds(scip, conss[c], result, &nbndchg, &redundant) ); /*lint !e613*/
3120  assert(!redundant); /* constraint should not be violated and redundant simultaneously (unless solution is far out of bounds) */
3121  if( *result == SCIP_CUTOFF || *result == SCIP_REDUCEDDOM )
3122  {
3123  SCIPdebugMsg(scip, "enforced by %s\n", *result == SCIP_CUTOFF ? "cutting off node" : "reducing domain");
3124  return SCIP_OKAY;
3125  }
3126  }
3127 
3128  SCIPwarningMessage(scip, "could not enforce feasibility by separating or branching; declaring solution with viol %g feasible\n", SCIPconsGetData(maxviolcons)->violation);
3129  *result = SCIP_FEASIBLE;
3130 
3131  return SCIP_OKAY;
3132 }
3133 
3134 /*
3135  * Quadratic constraint upgrading
3136  */
3137 
3138 
3139 #ifdef QUADCONSUPGD_PRIORITY
3140 /** tries to upgrade a quadratic constraint to a SOC constraint
3141  *
3142  * We treat as special cases, quadratic with no bilinear terms and hyperbolic quadratic
3143  * constraints with exactly on bilinear component containing nonnegative variables. For this we use the formula:
3144  * \f[
3145  * x^T x \leq yz,\; y \geq 0,\; z \geq 0
3146  * \qquad\Leftrightarrow\qquad
3147  * \left\| \left(\begin{array}{c} x \\ \frac{1}{2}(y - z)\end{array}\right) \right\| \leq \frac{1}{2}(y + z).
3148  * \f]
3149  *
3150  * @todo implement more general hyperbolic upgrade, e.g., for -x^T x + yz >= 0 or x^T x <= ax + by + cyz
3151  */
3152 static
3153 SCIP_DECL_QUADCONSUPGD(upgradeConsQuadratic)
3155  int nquadvars;
3156  int nbilinterms;
3157  SCIP_QUADVARTERM* quadterms;
3158  SCIP_BILINTERM* bilinterm;
3159  SCIP_VAR** quadvars;
3160  SCIP_VAR** lhsvars;
3161  SCIP_Real* lhscoefs;
3162  SCIP_Real* lhsoffsets;
3163  SCIP_Real lhsconstant;
3164  int lhscount;
3165  int lhsnvars;
3166  SCIP_VAR* rhsvar;
3167  SCIP_Real rhscoef;
3168  SCIP_Real rhsoffset;
3169  SCIP_VAR* bilinvar1;
3170  SCIP_VAR* bilinvar2;
3171  SCIP_Real bilincoef;
3172  int i;
3173  int j;
3174  int negeigpos;
3175  SCIP_Real* a;
3176  SCIP_Real* bp;
3177  SCIP_Real* eigvals;
3178  SCIP_Bool infeas;
3179  SCIP_Bool success;
3180  SCIP_Bool trygeneralupg;
3181  int nneg;
3182  int npos;
3183  SCIP_Bool rhsvarfound;
3184  SCIP_Bool rhsissoc;
3185  SCIP_Bool lhsissoc;
3186  SCIP_Real rhsvarlb;
3187  SCIP_Real rhsvarub;
3188  SCIP_CONSHDLR* conshdlr;
3189  SCIP_CONSHDLRDATA* conshdlrdata;
3190 
3191  assert(scip != NULL);
3192  assert(cons != NULL);
3193  assert(nupgdconss != NULL);
3194  assert(upgdconss != NULL);
3195 
3196  *nupgdconss = 0;
3197 
3198  SCIPdebugMsg(scip, "upgradeConsQuadratic called for constraint <%s>\n", SCIPconsGetName(cons));
3199  SCIPdebugPrintCons(scip, cons, NULL);
3200 
3201  /* currently do not support linear parts in upgrading of SOC constraints; binary vars we can treat as if squared */
3202  if( ncontlin > 0 || nimpllin > 0 || nintlin > 0 )
3203  return SCIP_OKAY;
3204  assert(SCIPgetNLinearVarsQuadratic(scip, cons) == nbinlin);
3205 
3206  nbilinterms = SCIPgetNBilinTermsQuadratic(scip, cons);
3207  nquadvars = SCIPgetNQuadVarTermsQuadratic(scip, cons);
3208 
3209  /* a proper SOC constraint needs at least 2 variables (at least one will be quadratic)
3210  * but performance-wise that doesn't give a clear advantage on product(2), so let's even require 3 vars
3211  */
3212  if( nbinlin + nquadvars < 3 )
3213  return SCIP_OKAY;
3214 
3215  /* reserve space */
3216  SCIP_CALL( SCIPallocBufferArray(scip, &lhsvars, nbinlin + nquadvars - 1) );
3217  SCIP_CALL( SCIPallocBufferArray(scip, &lhscoefs, nbinlin + nquadvars - 1) );
3218  SCIP_CALL( SCIPallocBufferArray(scip, &lhsoffsets, nbinlin + nquadvars - 1) );
3219 
3220  /* initialize data */
3221  a = NULL;
3222  bp = NULL;
3223  eigvals = NULL;
3224  quadvars = NULL;
3225  bilinvar1 = NULL;
3226  bilinvar2 = NULL;
3227  lhsconstant = 0.0;
3228  lhscount = 0;
3229  rhsvar = NULL;
3230  rhscoef = 0.0;
3231  rhsoffset = 0.0;
3232 
3233  trygeneralupg = FALSE;
3234 
3235  /* if more than one bilinear term is present, we are in the general case */
3236  if( nbilinterms > 1 )
3237  {
3238  trygeneralupg = TRUE;
3239  goto GENERALUPG;
3240  }
3241 
3242  /* check hyperbolic part */
3243  if( nbilinterms == 1 && nbinlin == 0 )
3244  {
3245  bilinterm = SCIPgetBilinTermsQuadratic(scip, cons);
3246  bilinvar1 = bilinterm->var1;
3247  bilinvar2 = bilinterm->var2;
3248  bilincoef = bilinterm->coef;
3249 
3250  /* the variables in the bilinear term need to be nonnegative */
3251  if ( SCIPisNegative(scip, SCIPvarGetLbGlobal(bilinvar1)) || SCIPisNegative(scip, SCIPvarGetLbGlobal(bilinvar2)) )
3252  {
3253  trygeneralupg = TRUE;
3254  goto GENERALUPG;
3255  }
3256 
3257  /* we need a rhs */
3258  if ( ! SCIPisZero(scip, SCIPgetRhsQuadratic(scip, cons)) )
3259  {
3260  trygeneralupg = TRUE;
3261  goto GENERALUPG;
3262  }
3263 
3264  /* we only allow for -1.0 bilincoef */
3265  if ( ! SCIPisEQ(scip, bilincoef, -1.0) )
3266  {
3267  trygeneralupg = TRUE;
3268  goto GENERALUPG;
3269  }
3270 
3271  /* check that bilinear terms do not appear in the rest and quadratic terms have positive sqrcoef have no lincoef */
3272  quadterms = SCIPgetQuadVarTermsQuadratic(scip, cons);
3273  for (i = 0; i < nquadvars; ++i)
3274  {
3275  SCIP_QUADVARTERM* term = &quadterms[i];
3276 
3277  if( ! SCIPisZero(scip, term->lincoef) || SCIPisNegative(scip, term->sqrcoef) )
3278  {
3279  trygeneralupg = TRUE;
3280  goto GENERALUPG;
3281  }
3282 
3283  if ( (term->var == bilinvar1 || term->var == bilinvar2) && ! SCIPisZero(scip, term->sqrcoef) )
3284  {
3285  trygeneralupg = TRUE;
3286  goto GENERALUPG;
3287  }
3288  }
3289  }
3290 
3291  quadterms = SCIPgetQuadVarTermsQuadratic(scip, cons);
3292  assert( quadterms != NULL );
3293 
3294  if( !SCIPisInfinity(scip, SCIPgetRhsQuadratic(scip, cons)) )
3295  {
3296  /* try whether constraint on right hand side is SOC */
3297  lhsconstant = -SCIPgetRhsQuadratic(scip, cons);
3298 
3299  for( i = 0; i < nquadvars + nbinlin; ++i )
3300  {
3301  SCIP_Real sqrcoef;
3302  SCIP_Real lincoef;
3303  SCIP_VAR* var;
3304 
3305  if( i < nquadvars )
3306  {
3307  SCIP_QUADVARTERM* term = &quadterms[i];
3308 
3309  sqrcoef = term->sqrcoef;
3310  lincoef = term->lincoef;
3311  var = term->var;
3312  }
3313  else
3314  {
3315  sqrcoef = SCIPgetCoefsLinearVarsQuadratic(scip, cons)[i-nquadvars];
3316  lincoef = 0.0;
3317  var = SCIPgetLinearVarsQuadratic(scip, cons)[i-nquadvars];
3318  }
3319 
3320  /* skip terms with 0 coefficients */
3321  if( sqrcoef == 0.0 )
3322  continue;
3323 
3324  if( sqrcoef > 0.0 )
3325  {
3326  if( lhscount >= nbinlin + nquadvars - 1 )
3327  { /* too many variables on lhs, i.e., all variables seem to have positive coefficient */
3328  rhsvar = NULL;
3329  break;
3330  }
3331 
3332  lhsvars[lhscount] = var;
3333  lhscoefs[lhscount] = sqrt(sqrcoef);
3334 
3335  if( lincoef != 0.0 )
3336  {
3337  lhsoffsets[lhscount] = lincoef / (2.0 * sqrcoef);
3338  lhsconstant -= lincoef * lincoef / (4.0 * sqrcoef);
3339  }
3340  else
3341  {
3342  lhsoffsets[lhscount] = 0.0;
3343  }
3344 
3345  ++lhscount;
3346  }
3347  else if( rhsvar != NULL ||
3348  (SCIPisLT(scip, SCIPcomputeVarLbGlobal(scip, var), -lincoef / (2.0 * sqrcoef)) &&
3349  SCIPisGT(scip, SCIPcomputeVarUbGlobal(scip, var), -lincoef / (2.0 * sqrcoef))) )
3350  { /* second variable with negative coefficient -> cannot be SOC */
3351  /* or rhs side changes sign */
3352  rhsvar = NULL;
3353  break;
3354  }
3355  else if( SCIPvarIsBinary(var) )
3356  {
3357  /* binary variable on rhs? then we originally had a convex quadratic */
3358  break;
3359  }
3360  else
3361  {
3362  rhsvar = var;
3363  rhsoffset = lincoef / (2.0 * sqrcoef);
3364  lhsconstant -= lincoef * lincoef / (4.0 * sqrcoef);
3365 
3366  if( SCIPisGE(scip, SCIPcomputeVarLbGlobal(scip, var), -lincoef / (2.0 * sqrcoef)) )
3367  rhscoef = sqrt(-sqrcoef);
3368  else
3369  rhscoef = -sqrt(-sqrcoef);
3370  }
3371  }
3372  }
3373 
3374  /* treat hyberbolic case */
3375  if( nbilinterms == 1 && nbinlin == 0 )
3376  {
3377  char name[SCIP_MAXSTRLEN];
3378  SCIP_VAR* auxvarsum;
3379  SCIP_VAR* auxvardiff;
3380  SCIP_CONS* couplingcons;
3381  SCIP_VAR* consvars[3];
3382  SCIP_Real consvals[3];
3383 
3384  /* can only upgrade if rhs is 0 */
3385  if ( rhsvar != NULL )
3386  goto cleanup;
3387 
3388  SCIPdebugMsg(scip, "found hyberbolic quadratic constraint <%s> to be SOC\n", SCIPconsGetName(cons));
3389 
3390  /* check if upgdconss is long enough to store upgrade constraints: we need two if we will have a quadratic
3391  * constraint for the left hand side left */
3392  *nupgdconss = SCIPisInfinity(scip, -SCIPgetLhsQuadratic(scip, cons)) ? 1 : 2;
3393  if ( *nupgdconss > upgdconsssize )
3394  {
3395  /* signal that we need more memory and return */
3396  *nupgdconss = -*nupgdconss;
3397  goto cleanup;
3398  }
3399 
3400  /* create auxiliary variable for sum (nonnegative) */
3401  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "soc#%s_s", SCIPconsGetName(cons));
3402  SCIP_CALL( SCIPcreateVar(scip, &auxvarsum, name, 0.0, SCIPinfinity(scip), 0.0,
3404  SCIP_CALL( SCIPaddVar(scip, auxvarsum) );
3405 
3406  /* create auxiliary variable for difference (free) */
3407  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "soc#%s_d", SCIPconsGetName(cons));
3408  SCIP_CALL( SCIPcreateVar(scip, &auxvardiff, name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0,
3410  SCIP_CALL( SCIPaddVar(scip, auxvardiff) );
3411 
3412  /* add coupling constraint for sum */
3413  consvars[0] = bilinvar1;
3414  consvars[1] = bilinvar2;
3415  consvars[2] = auxvarsum;
3416  consvals[0] = 1.0;
3417  consvals[1] = 1.0;
3418  consvals[2] = -1.0;
3419 
3420  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "soc#%s_cs", SCIPconsGetName(cons));
3421  SCIP_CALL( SCIPcreateConsLinear(scip, &couplingcons, name, 3, consvars, consvals, 0.0, 0.0,
3425  SCIPconsIsStickingAtNode(cons)) );
3426  SCIP_CALL( SCIPaddCons(scip, couplingcons) );
3427  SCIP_CALL( SCIPreleaseCons(scip, &couplingcons) );
3428 
3429  /* add coupling constraint for difference */
3430  consvars[0] = bilinvar1;
3431  consvars[1] = bilinvar2;
3432  consvars[2] = auxvardiff;
3433  consvals[0] = 1.0;
3434  consvals[1] = -1.0;
3435  consvals[2] = -1.0;
3436 
3437  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "soc#%s_cd", SCIPconsGetName(cons));
3438  SCIP_CALL( SCIPcreateConsLinear(scip, &couplingcons, name, 3, consvars, consvals, 0.0, 0.0,
3442  SCIPconsIsStickingAtNode(cons)) );
3443  SCIP_CALL( SCIPaddCons(scip, couplingcons) );
3444  SCIP_CALL( SCIPreleaseCons(scip, &couplingcons) );
3445 
3446  /* add difference variable to constraint */
3447  lhsvars[lhscount] = auxvardiff;
3448  lhscoefs[lhscount] = 0.5;
3449  lhsoffsets[lhscount] = 0.0;
3450 
3451  SCIP_CALL( SCIPcreateConsSOC(scip, &upgdconss[0], SCIPconsGetName(cons),
3452  lhscount + 1, lhsvars, lhscoefs, lhsoffsets, MAX(lhsconstant, 0.0),
3453  auxvarsum, 0.5, 0.0,
3457  SCIPdebugPrintCons(scip, upgdconss[0], NULL);
3458 
3459  /* create constraint that is equal to cons except that rhs is now infinity */
3460  if( !SCIPisInfinity(scip, -SCIPgetLhsQuadratic(scip, cons)) )
3461  {
3462  SCIP_CALL( SCIPcreateConsQuadratic2(scip, &upgdconss[1], SCIPconsGetName(cons),
3466  SCIPgetLhsQuadratic(scip, cons), SCIPinfinity(scip),
3470  }
3471  SCIP_CALL( SCIPreleaseVar(scip, &auxvardiff) );
3472  SCIP_CALL( SCIPreleaseVar(scip, &auxvarsum) );
3473 
3474  for (i = 0; i < lhscount; ++i)
3475  {
3476  SCIP_CALL( SCIPmarkDoNotMultaggrVar(scip, lhsvars[i]) );
3477  }
3478 
3479  goto cleanup;
3480  }
3481 
3482  if( rhsvar != NULL && lhscount >= 1 && !SCIPisNegative(scip, lhsconstant) )
3483  { /* found SOC constraint, so upgrade to SOC constraint(s) (below) and relax right hand side */
3484  SCIPdebugMsg(scip, "found right hand side of constraint <%s> to be SOC\n", SCIPconsGetName(cons));
3485 
3486  /* check if upgdconss is long enough to store upgrade constraints:
3487  * we need two if we will have a quadratic constraint for the left hand side left */
3488  *nupgdconss = SCIPisInfinity(scip, -SCIPgetLhsQuadratic(scip, cons)) ? 1 : 2;
3489  if( *nupgdconss > upgdconsssize )
3490  {
3491  /* signal that we need more memory and return */
3492  *nupgdconss = -*nupgdconss;
3493  goto cleanup;
3494  }
3495 
3496  SCIP_CALL( SCIPcreateConsSOC(scip, &upgdconss[0], SCIPconsGetName(cons),
3497  lhscount, lhsvars, lhscoefs, lhsoffsets, MAX(lhsconstant, 0.0),
3498  rhsvar, rhscoef, rhsoffset,
3502  SCIPdebugPrintCons(scip, upgdconss[0], NULL);
3503 
3504  /* create constraint that is equal to cons except that rhs is now infinity */
3505  if( !SCIPisInfinity(scip, -SCIPgetLhsQuadratic(scip, cons)) )
3506  {
3507  SCIP_CALL( SCIPcreateConsQuadratic2(scip, &upgdconss[1], SCIPconsGetName(cons),
3511  SCIPgetLhsQuadratic(scip, cons), SCIPinfinity(scip),
3515  }
3516  }
3517  else if( !SCIPisInfinity(scip, - SCIPgetLhsQuadratic(scip, cons)) )
3518  { /* if the first failed, try if constraint on left hand side is SOC (using negated coefficients) */
3519  lhscount = 0;
3520  rhsvar = NULL;
3521 
3522  lhsconstant = SCIPgetLhsQuadratic(scip, cons);
3523 
3524  for( i = 0; i < nquadvars + nbinlin; ++i )
3525  {
3526  SCIP_Real sqrcoef;
3527  SCIP_Real lincoef;
3528  SCIP_VAR* var;
3529 
3530  if( i < nquadvars )
3531  {
3532  SCIP_QUADVARTERM* term = &quadterms[i];
3533 
3534  sqrcoef = term->sqrcoef;
3535  lincoef = term->lincoef;
3536  var = term->var;
3537  }
3538  else
3539  {
3540  sqrcoef = SCIPgetCoefsLinearVarsQuadratic(scip, cons)[i-nquadvars];
3541  lincoef = 0.0;
3542  var = SCIPgetLinearVarsQuadratic(scip, cons)[i-nquadvars];
3543  }
3544 
3545  /* if there is a linear variable that is still considered as quadratic (constraint probably not presolved yet),
3546  * then give up
3547  */
3548  if( sqrcoef == 0.0 )
3549  goto cleanup;
3550 
3551  if( sqrcoef < 0.0 )
3552  {
3553  if( lhscount >= nquadvars + nbinlin - 1 )
3554  { /* too many variables on lhs, i.e., all variables seem to have negative coefficient */
3555  rhsvar = NULL;
3556  break;
3557  }
3558 
3559  lhsvars[lhscount] = var;
3560  lhscoefs[lhscount] = sqrt(-sqrcoef);
3561 
3562  if( lincoef != 0.0 )
3563  {
3564  lhsoffsets[lhscount] = lincoef / (2.0 * sqrcoef);
3565  lhsconstant += lincoef * lincoef / (4.0 * sqrcoef);
3566  }
3567  else
3568  {
3569  lhsoffsets[lhscount] = 0.0;
3570  }
3571 
3572  ++lhscount;
3573  }
3574  else if( rhsvar != NULL ||
3575  (SCIPisLT(scip, SCIPcomputeVarLbGlobal(scip, var), -lincoef / (2.0 * sqrcoef))
3576  && SCIPisGT(scip, SCIPcomputeVarUbGlobal(scip, var), -lincoef / (2.0 * sqrcoef))) )
3577  { /* second variable with positive coefficient -> cannot be SOC */
3578  /* or rhs side changes sign */
3579  rhsvar = NULL;
3580  break;
3581  }
3582  else if( SCIPvarIsBinary(var) )
3583  {
3584  /* binary variable on rhs? then we originally had a convex quadratic */
3585  break;
3586  }
3587  else
3588  {
3589  rhsvar = var;
3590  rhsoffset = lincoef / (2.0 * sqrcoef);
3591  lhsconstant += lincoef * lincoef / (4.0 * sqrcoef);
3592 
3593  if( SCIPisGE(scip, SCIPcomputeVarLbGlobal(scip, var), -lincoef / (2.0 * sqrcoef)) )
3594  rhscoef = sqrt(sqrcoef);
3595  else
3596  rhscoef = -sqrt(sqrcoef);
3597  }
3598  }
3599 
3600  if( rhsvar != NULL && lhscount >= 1 && !SCIPisNegative(scip, lhsconstant) )
3601  { /* found SOC constraint, so upgrade to SOC constraint(s) (below) and relax left hand side */
3602  SCIPdebugMsg(scip, "found left hand side of constraint <%s> to be SOC\n", SCIPconsGetName(cons));
3603 
3604  /* check if upgdconss is long enough to store upgrade constraints:
3605  * we need two if we will have a quadratic constraint for the right hand side left */
3606  *nupgdconss = SCIPisInfinity(scip, SCIPgetRhsQuadratic(scip, cons)) ? 1 : 2;
3607  if( *nupgdconss > upgdconsssize )
3608  {
3609  /* signal that we need more memory and return */
3610  *nupgdconss = -*nupgdconss;
3611  goto cleanup;
3612  }
3613 
3614  SCIP_CALL( SCIPcreateConsSOC(scip, &upgdconss[0], SCIPconsGetName(cons),
3615  lhscount, lhsvars, lhscoefs, lhsoffsets, MAX(lhsconstant, 0.0),
3616  rhsvar, rhscoef, rhsoffset,
3620  SCIPdebugPrintCons(scip, upgdconss[0], NULL);
3621 
3622  /* create constraint that is equal to cons except that lhs is now -infinity */
3623  if( !SCIPisInfinity(scip, SCIPgetRhsQuadratic(scip, cons)) )
3624  {
3625  SCIP_CALL( SCIPcreateConsQuadratic2(scip, &upgdconss[1], SCIPconsGetName(cons),
3629  -SCIPinfinity(scip), SCIPgetRhsQuadratic(scip, cons),
3633  }
3634  }
3635  }
3636 
3637 GENERALUPG:
3638  if( !trygeneralupg )
3639  goto cleanup;
3640 
3641  /* find the soc constraint handler */
3642  conshdlr = SCIPfindConshdlr(scip, CONSHDLR_NAME);
3643  assert(conshdlr != NULL);
3644  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3645  assert(conshdlrdata != NULL);
3646 
3647  if( !conshdlrdata->generalsocupg )
3648  goto cleanup;
3649 
3650  SCIPdebugMsg(scip, "Trying general method of upgrade to a soc const\n");
3651 
3652  rhsvarlb = 1.0;
3653  rhsvarub = 0.0;
3654 
3655 #ifndef SCIP_STATISTIC
3656  /* skip large matrices (needs to compute eigenvalues/vectors) according to presolve timing */
3657  if( (presoltiming & SCIP_PRESOLTIMING_FAST) != 0 && nquadvars > 10 )
3658  goto cleanup;
3659  if( (presoltiming & SCIP_PRESOLTIMING_MEDIUM) != 0 && nquadvars > 50 )
3660  goto cleanup;
3661 #endif
3662 
3663  /* need routine to compute eigenvalues/eigenvectors */
3664  if( !SCIPisIpoptAvailableIpopt() )
3665  goto cleanup;
3666 
3667  /* build lower triangular part the A matrix */
3668  SCIP_CALL( SCIPallocClearBufferArray(scip, &a, nquadvars*nquadvars) ); /*lint !e647*/
3669  SCIP_CALL( SCIPallocClearBufferArray(scip, &bp, nquadvars) );
3670  SCIP_CALL( SCIPallocBufferArray(scip, &quadvars, nquadvars) );
3671  SCIP_CALL( SCIPallocBufferArray(scip, &eigvals, nquadvars) );
3672 
3673  /* make sure quadratic variables terms are sorted */
3675 
3676  /* set lower triangular entries of A corresponding to square terms */
3677  for( i = 0; i < nquadvars; ++i )
3678  {
3679  SCIP_QUADVARTERM* term = &SCIPgetQuadVarTermsQuadratic(scip, cons)[i];
3680 
3681  /* skip upgrade if fixed (or (multi)aggregated) variables and still in fast presolving
3682  * probably cons_quadratic did not yet had the chance to remove/replace this variable (see also #2352)
3683  */
3684  if( !SCIPvarIsActive(term->var) && (presoltiming & SCIP_PRESOLTIMING_FAST) != 0 )
3685  goto cleanup;
3686 
3687  a[i*nquadvars + i] = term->sqrcoef;
3688  quadvars[i] = term->var;
3689  }
3690 
3691  /* set lower triangular entries of A corresponding to bilinear terms */
3692  for( i = 0; i < nbilinterms; ++i )
3693  {
3694  int idx1;
3695  int idx2;
3696 
3697  bilinterm = &SCIPgetBilinTermsQuadratic(scip, cons)[i];
3698 
3699  SCIP_CALL( SCIPfindQuadVarTermQuadratic(scip, cons, bilinterm->var1, &idx1) );
3700  SCIP_CALL( SCIPfindQuadVarTermQuadratic(scip, cons, bilinterm->var2, &idx2) );
3701 
3702  assert(idx1 >= 0);
3703  assert(idx2 >= 0);
3704  assert(idx1 != idx2);
3705 
3706  a[MIN(idx1,idx2) * nquadvars + MAX(idx1,idx2)] = bilinterm->coef / 2.0;
3707  }
3708 
3709  /* compute eigenvalues and vectors, A = PDP^t
3710  * note: a stores P^t
3711  */
3712  if( LapackDsyev(TRUE, nquadvars, a, eigvals) != SCIP_OKAY )
3713  {
3714  SCIPdebugMsg(scip, "Failed to compute eigenvalues and eigenvectors for constraint <%s>.\n", SCIPconsGetName(cons));
3715  goto cleanup;
3716  }
3717 
3718  /* set small eigenvalues to 0 and compute b*P */
3719  nneg = 0;
3720  npos = 0;
3721  for( i = 0; i < nquadvars; ++i )
3722  {
3723  for( j = 0; j < nquadvars; ++j )
3724  {
3725  SCIP_QUADVARTERM* term = &SCIPgetQuadVarTermsQuadratic(scip, cons)[j];
3726  bp[i] += term->lincoef * a[i * nquadvars + j];
3727  }
3728 
3729  if( SCIPisZero(scip, eigvals[i]) )
3730  {
3731  /* if there is a purely linear variable, the constraint can't be written as a SOC */
3732  if( !SCIPisZero(scip, bp[i]) )
3733  {
3734  goto cleanup;
3735  }
3736 
3737  bp[i] = 0.0;
3738  eigvals[i] = 0.0;
3739  }
3740  else if( eigvals[i] > 0.0 )
3741  npos++;
3742  else
3743  nneg++;
3744  }
3745 
3746  /* a proper SOC constraint needs at least 2 variables
3747  * let's make sure we have at least 3, though, as this upgrade comes with extra (multiaggr.) vars
3748  */
3749  if( npos + nneg < 3 )
3750  goto cleanup;
3751 
3752  /* determine whether rhs or lhs of cons is potentially SOC, if any */
3753  rhsissoc = nneg == 1 && !SCIPisInfinity(scip, SCIPgetRhsQuadratic(scip, cons));
3754  lhsissoc = npos == 1 && !SCIPisInfinity(scip, -SCIPgetLhsQuadratic(scip, cons));
3755 
3756  /* if non is potentially SOC, stop */
3757  if( !rhsissoc && !lhsissoc )
3758  goto cleanup;
3759 
3760  if( rhsissoc && lhsissoc )
3761  {
3762  /* only handle rhs-SOC here for now
3763  * if the upgrade is run again on the remaining lhs-SOC, it would upgrade that side
3764  */
3765  lhsissoc = FALSE;
3766  }
3767 
3768  if( lhsissoc )
3769  {
3770  /* lhs is potentially SOC, change signs */
3771  lhsconstant = SCIPgetLhsQuadratic(scip, cons);
3772 
3773  for( i = 0; i < nquadvars; ++i )
3774  {
3775  eigvals[i] = -eigvals[i];
3776  bp[i] = -bp[i];
3777  }
3778  }
3779  else
3780  {
3781  lhsconstant = -SCIPgetRhsQuadratic(scip, cons);
3782  }
3783 
3784  /* we have lhsconstant + x^t A x + b x <= 0 and A has a single negative eigenvalue; try to build soc */
3785  negeigpos = -1;
3786  rhsvarfound = FALSE;
3787  for( i = 0; i < nquadvars; ++i )
3788  {
3789  if( SCIPisZero(scip, eigvals[i]) )
3790  continue;
3791 
3792  if( eigvals[i] > 0.0 )
3793  {
3794  lhscoefs[lhscount] = sqrt(eigvals[i]) * UPGSCALE;
3795  lhsoffsets[lhscount] = bp[i] / (2 * eigvals[i]);
3796  lhsconstant -= bp[i] * bp[i] / (4 * eigvals[i]);
3797 
3798  ++lhscount;
3799  }
3800  else
3801  {
3802  /* should enter here only once */
3803  assert(rhsvarfound == FALSE);
3804 
3805  rhsoffset = bp[i] / (2 * eigvals[i]);
3806 
3807  /* the constraint can only be a soc if the resulting rhs var does not change var; the rhs var is going to be a
3808  * multiaggregated variable, so estimate its bounds
3809  */
3810  rhsvarlb = 0.0;
3811  rhsvarub = 0.0;
3812  for( j = 0; j < nquadvars; ++j )
3813  {
3814  SCIP_Real aux;
3815 
3816  if( SCIPisZero(scip, a[i * nquadvars + j]) )
3817  continue;
3818 
3819  if( a[i * nquadvars + j] > 0.0 )
3820  {
3821  aux = SCIPcomputeVarLbGlobal(scip, quadvars[j]);
3822  assert(!SCIPisInfinity(scip, aux));
3823  }
3824  else
3825  {
3826  aux = SCIPcomputeVarUbGlobal(scip, quadvars[j]);
3827  assert(!SCIPisInfinity(scip, -aux));
3828  }
3829 
3830  if( SCIPisInfinity(scip, aux) || SCIPisInfinity(scip, -aux) )
3831  {
3832  rhsvarlb = -SCIPinfinity(scip);
3833  break;
3834  }
3835  else
3836  rhsvarlb += a[i * nquadvars + j] * aux;
3837  }
3838  rhsvarlb += rhsoffset;
3839 
3840  for( j = 0; j < nquadvars; ++j )
3841  {
3842  SCIP_Real aux;
3843 
3844  if( SCIPisZero(scip, a[i * nquadvars + j]) )
3845  continue;
3846 
3847  if( a[i * nquadvars + j] > 0.0 )
3848  {
3849  aux = SCIPcomputeVarUbGlobal(scip, quadvars[j]);
3850  assert(!SCIPisInfinity(scip, -aux));
3851  }
3852  else
3853  {
3854  aux = SCIPcomputeVarLbGlobal(scip, quadvars[j]);
3855  assert(!SCIPisInfinity(scip, aux));
3856  }
3857 
3858  if( SCIPisInfinity(scip, aux) || SCIPisInfinity(scip, -aux) )
3859  {
3860  rhsvarub = SCIPinfinity(scip);
3861  break;
3862  }
3863  else
3864  rhsvarub += a[i * nquadvars + j] * aux;
3865  }
3866  rhsvarub += rhsoffset;
3867 
3868  /* since we are just interested in obtaining an interval that contains the real bounds and is tight enough so
3869  * that we can identify that the rhsvar does not change sign, we swap the bounds in case of numerical troubles
3870  */
3871  if( rhsvarub < rhsvarlb )
3872  {
3873  assert(SCIPisEQ(scip, rhsvarub, rhsvarlb));
3874  SCIPswapReals(&rhsvarub, &rhsvarlb);
3875  }
3876 
3877  /* check whether rhsvar changes sign */
3878  if( SCIPisGE(scip, rhsvarlb, 0.0) || SCIPisLE(scip, rhsvarub, 0.0) )
3879  {
3880  rhsvarfound = TRUE;
3881  negeigpos = i;
3882  lhsconstant -= rhsoffset * rhsoffset * eigvals[i];
3883  rhscoef = SCIPisLE(scip, rhsvarub, 0.0) ? -sqrt(-eigvals[i]) : sqrt(-eigvals[i]);
3884  }
3885  else
3886  {
3887  SCIPdebugMsg(scip, "Failed because rhsvar [%g, %g] changes sign.\n", rhsvarlb, rhsvarub);
3888  rhsvarfound = FALSE;
3889  break;
3890  }
3891  }
3892  }
3893 
3894  if( rhsvarfound && lhscount >= 1 && !SCIPisNegative(scip, lhsconstant) )
3895  {
3896  /* found SOC constraint, so upgrade to SOC constraint(s) (below) and relax right hand side */
3897  SCIPdebugMsg(scip, "found right hand side of constraint <%s> to be SOC\n", SCIPconsGetName(cons));
3898 
3899  /* check if upgdconss is long enough to store upgrade constraints:
3900  * we need two if we will have a quadratic constraint for the left hand side left */
3901  if( rhsissoc )
3902  *nupgdconss = SCIPisInfinity(scip, -SCIPgetLhsQuadratic(scip, cons)) ? 1 : 2;
3903  else
3904  *nupgdconss = SCIPisInfinity(scip, SCIPgetRhsQuadratic(scip, cons)) ? 1 : 2;
3905  if( *nupgdconss > upgdconsssize )
3906  {
3907  /* signal that we need more memory and return */
3908  *nupgdconss = -*nupgdconss;
3909  goto cleanup;
3910  }
3911 
3912  /* create lhs and rhs vars */
3913  lhsnvars = 0;
3914  for( i = 0; i < nquadvars; ++i )
3915  {
3916  if( eigvals[i] <= 0.0 )
3917  continue;
3918 
3919  SCIP_CALL( SCIPcreateVarBasic(scip, &lhsvars[lhsnvars], NULL,
3920  -SCIPinfinity(scip), SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
3921  SCIP_CALL( SCIPaddVar(scip, lhsvars[lhsnvars]) );
3922 
3923  SCIP_CALL( SCIPmultiaggregateVar(scip, lhsvars[lhsnvars], nquadvars, quadvars, &(a[i * nquadvars]),
3924  lhsoffsets[lhsnvars], &infeas, &success) );
3925 
3926  if( infeas || !success )
3927  {
3928  SCIPdebugMsg(scip, "Problem with aggregation while trying to upgrade <%s>.\n", SCIPconsGetName(cons) );
3929 
3930  /* release all created vars so far */
3931  for( j = 0; j <= lhsnvars; ++j )
3932  SCIP_CALL( SCIPreleaseVar(scip, &lhsvars[j]) );
3933 
3934  *nupgdconss = 0;
3935  goto cleanup;
3936  }
3937  lhsnvars++;
3938  }
3939  assert(lhsnvars == lhscount);
3940  assert(negeigpos >= 0);
3941 
3942  SCIP_CALL( SCIPcreateVarBasic(scip, &rhsvar, NULL,
3943  rhsvarlb, rhsvarub, 0.0, SCIP_VARTYPE_CONTINUOUS) );
3944  SCIP_CALL( SCIPaddVar(scip, rhsvar) );
3945  SCIP_CALL( SCIPmultiaggregateVar(scip, rhsvar, nquadvars, quadvars, &(a[negeigpos * nquadvars]),
3946  rhsoffset, &infeas, &success) );
3947 
3948  if( infeas || !success )
3949  {
3950  SCIPdebugMsg(scip, "Problem with aggregation while trying to upgrade <%s>.\n", SCIPconsGetName(cons) );
3951 
3952  /* release all created vars */
3953  SCIP_CALL( SCIPreleaseVar(scip, &rhsvar) );
3954  for( j = 0; j < lhsnvars; ++j )
3955  {
3956  SCIP_CALL( SCIPreleaseVar(scip, &lhsvars[j]) );
3957  }
3958 
3959  *nupgdconss = 0;
3960  goto cleanup;
3961  }
3962 
3963  SCIP_CALL( SCIPcreateConsSOC(scip, &upgdconss[0], SCIPconsGetName(cons),
3964  lhscount, lhsvars, lhscoefs, NULL, MAX(lhsconstant, 0.0) * UPGSCALE * UPGSCALE,
3965  rhsvar, rhscoef * UPGSCALE, 0.0,
3969  SCIPdebugPrintCons(scip, upgdconss[0], NULL);
3970 
3971  /* release created variables */
3972  SCIP_CALL( SCIPreleaseVar(scip, &rhsvar) );
3973  for( i = 0; i < lhsnvars; ++i )
3974  {
3975  SCIP_CALL( SCIPreleaseVar(scip, &lhsvars[i]) );
3976  }
3977 
3978  /* create constraint that is equal to cons except that rhs/lhs is now +/-infinity */
3979  if( rhsissoc && !SCIPisInfinity(scip, -SCIPgetLhsQuadratic(scip, cons)) )
3980  {
3981  SCIPdebugMsg(scip, "rhs is soc, keep quadratic\n");
3982  SCIP_CALL( SCIPcreateConsQuadratic2(scip, &upgdconss[1], SCIPconsGetName(cons),
3986  SCIPgetLhsQuadratic(scip, cons), SCIPinfinity(scip),
3990  }
3991  else if( lhsissoc && !SCIPisInfinity(scip, SCIPgetRhsQuadratic(scip,cons)) )
3992  {
3993  SCIPdebugMsg(scip, "lhs is soc, keep quadratic\n");
3994  SCIP_CALL( SCIPcreateConsQuadratic2(scip, &upgdconss[1], SCIPconsGetName(cons),
3998  -SCIPinfinity(scip), SCIPgetRhsQuadratic(scip, cons),
4002  }
4003  }
4004 #ifdef SCIP_DEBUG
4005  else
4006  {
4007  if( lhscount < 1 )
4008  {
4009  SCIPdebugMsg(scip, "Failed because there are not enough lhsvars (%d)\n", lhscount);
4010  }
4011  if( SCIPisNegative(scip, lhsconstant) )
4012  {
4013  SCIPdebugMsg(scip, "Failed because lhsconstant is negative (%g)\n", lhsconstant);
4014  }
4015  }
4016 #endif
4017 
4018  cleanup:
4019  SCIPfreeBufferArrayNull(scip, &eigvals);
4020  SCIPfreeBufferArrayNull(scip, &quadvars);
4021  SCIPfreeBufferArrayNull(scip, &bp);
4022  SCIPfreeBufferArrayNull(scip, &a);
4023  SCIPfreeBufferArray(scip, &lhsoffsets);
4024  SCIPfreeBufferArray(scip, &lhscoefs);
4025  SCIPfreeBufferArray(scip, &lhsvars);
4026 
4027  return SCIP_OKAY;
4028 } /*lint !e715*/
4029 #endif
4030 
4031 /*
4032  * Callback methods of constraint handler
4033  */
4034 
4035 /** copy method for constraint handler plugins (called when SCIP copies plugins) */
4036 static
4037 SCIP_DECL_CONSHDLRCOPY(conshdlrCopySOC)
4038 { /*lint --e{715}*/
4039  assert(scip != NULL);
4040  assert(conshdlr != NULL);
4041  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
4042 
4043  /* call inclusion method of constraint handler */
4045 
4046  *valid = TRUE;
4047 
4048  return SCIP_OKAY;
4049 }
4050 
4051 /** destructor of constraint handler to free constraint handler data (called when SCIP is exiting) */
4052 static
4053 SCIP_DECL_CONSFREE(consFreeSOC)
4054 {
4055  SCIP_CONSHDLRDATA* conshdlrdata;
4056 
4057  assert( scip != NULL );
4058  assert( conshdlr != NULL );
4059  assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
4060 
4061  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4062  assert(conshdlrdata != NULL);
4063 
4064  SCIPfreeBlockMemory(scip, &conshdlrdata);
4065 
4066  return SCIP_OKAY;
4067 }
4068 
4069 
4070 /** initialization method of constraint handler (called after problem was transformed) */
4071 static
4072 SCIP_DECL_CONSINIT(consInitSOC)
4073 { /*lint --e{715}*/
4074  SCIP_CONSHDLRDATA* conshdlrdata;
4075  int c;
4076 
4077  assert(scip != NULL);
4078  assert(conshdlr != NULL);
4079  assert(conss != NULL || nconss == 0);
4080 
4081  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4082  assert(conshdlrdata != NULL);
4083 
4084  conshdlrdata->subnlpheur = SCIPfindHeur(scip, "subnlp");
4085  conshdlrdata->trysolheur = SCIPfindHeur(scip, "trysol");
4086  conshdlrdata->haveexprint = (strcmp(SCIPexprintGetName(), "NONE") != 0);
4087 
4088  /* mark constraints for propagation */
4089  for( c = 0; c < nconss; ++c )
4090  {
4091  SCIP_CALL( SCIPmarkConsPropagate(scip, conss[c]) ); /*lint !e613*/
4092  }
4093 
4094  return SCIP_OKAY;
4095 }
4096 
4097 
4098 /** deinitialization method of constraint handler (called before transformed problem is freed) */
4099 static
4100 SCIP_DECL_CONSEXIT(consExitSOC)
4101 { /*lint --e{715}*/
4102  SCIP_CONSHDLRDATA* conshdlrdata;
4103 
4104  assert(scip != NULL);
4105  assert(conshdlr != NULL);
4106 
4107  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4108  assert(conshdlrdata != NULL);
4109 
4110  conshdlrdata->subnlpheur = NULL;
4111  conshdlrdata->trysolheur = NULL;
4112  conshdlrdata->haveexprint = FALSE;
4113 
4114  return SCIP_OKAY;
4115 }
4116 
4117 /** presolving deinitialization method of constraint handler (called after presolving has been finished) */
4118 static
4119 SCIP_DECL_CONSEXITPRE(consExitpreSOC)
4120 { /*lint --e{715}*/
4121  int c;
4122 
4123  assert(scip != NULL);
4124  assert(conshdlr != NULL);
4125  assert(conss != NULL || nconss == 0);
4126 
4127  /* tell SCIP that we have something nonlinear */
4128  for( c = 0; c < nconss; ++c )
4129  {
4130  if( SCIPconsIsAdded(conss[c]) ) /*lint !e613*/
4131  {
4132  SCIPenableNLP(scip);
4133  break;
4134  }
4135  }
4136 
4137  return SCIP_OKAY;
4138 }
4139 
4140 
4141 /** solving process initialization method of constraint handler (called when branch and bound process is about to begin) */
4142 static
4143 SCIP_DECL_CONSINITSOL(consInitsolSOC)
4144 {
4145  SCIP_CONSHDLRDATA* conshdlrdata;
4146  SCIP_CONSDATA* consdata;
4147  int c;
4148 
4149  assert(scip != NULL);
4150  assert(conshdlr != NULL);
4151 
4152  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4153  assert(conshdlrdata != NULL);
4154  assert(conshdlrdata->eventhdlr != NULL);
4155 
4156  /* add nlrow representation to NLP, if NLP has been enabled */
4157  if( SCIPisNLPConstructed(scip) )
4158  {
4159  for( c = 0; c < nconss; ++c )
4160  {
4161  if( SCIPconsIsEnabled(conss[c]) )
4162  {
4163  consdata = SCIPconsGetData(conss[c]);
4164  assert(consdata != NULL);
4165 
4166  if( consdata->nlrow == NULL )
4167  {
4168  SCIP_CALL( createNlRow(scip, conshdlr, conss[c]) );
4169  assert(consdata->nlrow != NULL);
4170  }
4171  SCIP_CALL( SCIPaddNlRow(scip, consdata->nlrow) );
4172  }
4173  }
4174  }
4175 
4176  conshdlrdata->newsoleventfilterpos = -1;
4177  if( nconss != 0 )
4178  {
4179  SCIP_EVENTHDLR* eventhdlr;
4180 
4181  eventhdlr = SCIPfindEventhdlr(scip, CONSHDLR_NAME"_newsolution");
4182  assert(eventhdlr != NULL);
4183 
4184  SCIP_CALL( SCIPcatchEvent(scip, SCIP_EVENTTYPE_SOLFOUND, eventhdlr, (SCIP_EVENTDATA*)conshdlr, &conshdlrdata->newsoleventfilterpos) );
4185  }
4186 
4187  /* reset flags and counters */
4188  conshdlrdata->sepanlp = FALSE;
4189  conshdlrdata->lastenfonode = NULL;
4190  conshdlrdata->nenforounds = 0;
4191 
4192  return SCIP_OKAY;
4193 }
4194 
4195 
4196 /** solving process deinitialization method of constraint handler (called before branch and bound process data is freed) */
4197 static
4198 SCIP_DECL_CONSEXITSOL(consExitsolSOC)
4200  SCIP_CONSHDLRDATA* conshdlrdata;
4201  SCIP_CONSDATA* consdata;
4202  int c;
4203 
4204  assert(scip != NULL);
4205  assert(conshdlr != NULL);
4206 
4207  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4208  assert(conshdlrdata != NULL);
4209  assert(conshdlrdata->eventhdlr != NULL);
4210 
4211  if( conshdlrdata->newsoleventfilterpos >= 0 )
4212  {
4213  SCIP_EVENTHDLR* eventhdlr;
4214 
4215  eventhdlr = SCIPfindEventhdlr(scip, CONSHDLR_NAME"_newsolution");
4216  assert(eventhdlr != NULL);
4217 
4218  SCIP_CALL( SCIPdropEvent(scip, SCIP_EVENTTYPE_SOLFOUND, eventhdlr, (SCIP_EVENTDATA*)conshdlr, conshdlrdata->newsoleventfilterpos) );
4219  conshdlrdata->newsoleventfilterpos = -1;
4220  }
4221 
4222  for( c = 0; c < nconss; ++c )
4223  {
4224  consdata = SCIPconsGetData(conss[c]);
4225  assert(consdata != NULL);
4226 
4227  /* free nonlinear row representation */
4228  if( consdata->nlrow != NULL )
4229  {
4230  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
4231  }
4232  }
4233 
4234  return SCIP_OKAY;
4235 } /*lint !e715*/
4236 
4237 
4238 /** frees specific constraint data */
4239 static
4240 SCIP_DECL_CONSDELETE(consDeleteSOC)
4242  int i;
4243 
4244  assert(scip != NULL);
4245  assert(conshdlr != NULL);
4246  assert(cons != NULL);
4247  assert(consdata != NULL);
4248  assert(*consdata != NULL);
4249  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
4250 
4251  SCIPdebugMsg(scip, "Deleting SOC constraint <%s>.\n", SCIPconsGetName(cons) );
4252 
4253  if( SCIPconsIsTransformed(cons) )
4254  {
4255  SCIP_CONSHDLRDATA* conshdlrdata;
4256 
4257  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4258  assert(conshdlrdata != NULL);
4259 
4260  SCIP_CALL( dropVarEvents(scip, conshdlrdata->eventhdlr, cons) );
4261  }
4262 
4263  for( i = 0; i < (*consdata)->nvars; ++i )
4264  {
4265  SCIP_CALL( SCIPreleaseVar(scip, &(*consdata)->vars[i]) );
4266  }
4267 
4268  SCIPfreeBlockMemoryArray(scip, &(*consdata)->vars, (*consdata)->nvars);
4269  SCIPfreeBlockMemoryArray(scip, &(*consdata)->coefs, (*consdata)->nvars);
4270  SCIPfreeBlockMemoryArray(scip, &(*consdata)->offsets, (*consdata)->nvars);
4271  assert((*consdata)->lhsbndchgeventdata == NULL);
4272 
4273  if( (*consdata)->rhsvar != NULL )
4274  {
4275  SCIP_CALL( SCIPreleaseVar(scip, &(*consdata)->rhsvar) );
4276  }
4277 
4278  /* free nonlinear row representation
4279  * normally released in exitsol, but constraint may be deleted early (e.g., if found redundant)
4280  */
4281  if( (*consdata)->nlrow != NULL )
4282  {
4283  SCIP_CALL( SCIPreleaseNlRow(scip, &(*consdata)->nlrow) );
4284  }
4285 
4286  SCIPfreeBlockMemory(scip, consdata);
4287 
4288  return SCIP_OKAY;
4289 }
4290 
4291 
4292 /** transforms constraint data into data belonging to the transformed problem */
4293 static
4294 SCIP_DECL_CONSTRANS(consTransSOC)
4295 {
4296  SCIP_CONSDATA* consdata;
4297  SCIP_CONSHDLRDATA* conshdlrdata;
4298  SCIP_CONSDATA* sourcedata;
4299  char s[SCIP_MAXSTRLEN];
4300  int i;
4301 
4302  assert(scip != NULL);
4303  assert(conshdlr != NULL);
4304  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
4305  assert(sourcecons != NULL);
4306  assert(targetcons != NULL);
4307 
4308  /* get constraint handler data */
4309  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4310  assert(conshdlrdata != NULL);
4311 
4312  SCIPdebugMsg(scip, "Transforming SOC constraint: <%s>.\n", SCIPconsGetName(sourcecons) );
4313 
4314  /* get data of original constraint */
4315  sourcedata = SCIPconsGetData(sourcecons);
4316  assert(sourcedata != NULL);
4317  assert(sourcedata->vars != NULL);
4318  assert(sourcedata->coefs != NULL);
4319  assert(sourcedata->offsets != NULL);
4320 
4321  /* create constraint data */
4322  SCIP_CALL( SCIPallocBlockMemory(scip, &consdata) );
4323 
4324  consdata->nvars = sourcedata->nvars;
4325  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->vars, consdata->nvars) );
4326  SCIP_CALL( SCIPgetTransformedVars(scip, consdata->nvars, sourcedata->vars, consdata->vars) );
4327  for( i = 0; i < consdata->nvars; ++i )
4328  {
4329  SCIP_CALL( SCIPcaptureVar(scip, consdata->vars[i]) );
4330  }
4331 
4332  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &consdata->coefs, sourcedata->coefs, consdata->nvars) );
4333  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &consdata->offsets, sourcedata->offsets, consdata->nvars) );
4334  consdata->constant = sourcedata->constant;
4335 
4336  SCIP_CALL( SCIPgetTransformedVar(scip, sourcedata->rhsvar, &consdata->rhsvar) );
4337  consdata->rhsoffset = sourcedata->rhsoffset;
4338  consdata->rhscoeff = sourcedata->rhscoeff;
4339 
4340  SCIP_CALL( SCIPcaptureVar(scip, consdata->rhsvar) );
4341 
4342  consdata->nlrow = NULL;
4343  consdata->lhsbndchgeventdata = NULL;
4344  consdata->isapproxadded = FALSE;
4345 
4346  /* create transformed constraint with the same flags */
4347  (void) SCIPsnprintf(s, SCIP_MAXSTRLEN, "t_%s", SCIPconsGetName(sourcecons));
4348  SCIP_CALL( SCIPcreateCons(scip, targetcons, s, conshdlr, consdata,
4349  SCIPconsIsInitial(sourcecons), SCIPconsIsSeparated(sourcecons),
4350  SCIPconsIsEnforced(sourcecons), SCIPconsIsChecked(sourcecons),
4351  SCIPconsIsPropagated(sourcecons), SCIPconsIsLocal(sourcecons),
4352  SCIPconsIsModifiable(sourcecons), SCIPconsIsDynamic(sourcecons),
4353  SCIPconsIsRemovable(sourcecons), SCIPconsIsStickingAtNode(sourcecons)) );
4354 
4355  SCIP_CALL( catchVarEvents(scip, conshdlrdata->eventhdlr, *targetcons) );
4356 
4357  return SCIP_OKAY;
4358 }
4359 
4360 /** separation method of constraint handler for LP solutions */
4361 static
4362 SCIP_DECL_CONSSEPALP(consSepalpSOC)
4363 {
4364  SCIP_CONSHDLRDATA* conshdlrdata;
4365  SCIP_CONS* maxviolcon;
4366  SCIP_Bool sepasuccess;
4367  SCIP_Bool cutoff;
4368 
4369  assert(scip != NULL);
4370  assert(conshdlr != NULL);
4371  assert(conss != NULL || nconss == 0);
4372  assert(result != NULL);
4373 
4374  *result = SCIP_DIDNOTFIND;
4375 
4376  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4377  assert(conshdlrdata != NULL);
4378 
4379  SCIP_CALL( computeViolations(scip, conshdlr, conss, nconss, NULL, &maxviolcon) );
4380  if( maxviolcon == NULL )
4381  return SCIP_OKAY;
4382 
4383  /* at root, check if we want to solve the NLP relaxation and use its solutions as reference point
4384  * if there is something convex, then linearizing in the solution of the NLP relaxation can be very useful
4385  */
4386  if( SCIPgetDepth(scip) == 0 && !conshdlrdata->sepanlp &&
4387  (SCIPgetNContVars(scip) >= conshdlrdata->sepanlpmincont * SCIPgetNVars(scip) || (SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_UNBOUNDEDRAY && conshdlrdata->sepanlpmincont <= 1.0)) &&
4388  SCIPisNLPConstructed(scip) && SCIPgetNNlpis(scip) > 0 )
4389  {
4390  SCIP_NLPSOLSTAT solstat;
4391  SCIP_Bool solvednlp;
4392 
4393  solstat = SCIPgetNLPSolstat(scip);
4394  solvednlp = FALSE;
4395  if( solstat == SCIP_NLPSOLSTAT_UNKNOWN )
4396  {
4397  /* NLP is not solved yet, so we do it now and update solstat */
4398 
4399  /* ensure linear conss are in NLP */
4400  if( conshdlrdata->subnlpheur != NULL )
4401  {
4402  SCIP_CALL( SCIPaddLinearConsToNlpHeurSubNlp(scip, conshdlrdata->subnlpheur, TRUE, TRUE) );
4403  }
4404 
4405  /* set LP solution as starting values, if available */
4407  {
4409  }
4410 
4411  /* SCIP_CALL( SCIPsetNLPIntPar(scip, SCIP_NLPPAR_VERBLEVEL, 1) ); */
4412  SCIP_CALL( SCIPsolveNLP(scip) );
4413 
4414  solstat = SCIPgetNLPSolstat(scip);
4415  SCIPdebugMsg(scip, "solved NLP relax, solution status: %d\n", solstat);
4416 
4417  solvednlp = TRUE;
4418  }
4419 
4420  conshdlrdata->sepanlp = TRUE;
4421 
4422  if( solstat == SCIP_NLPSOLSTAT_GLOBINFEASIBLE )
4423  {
4424  SCIPdebugMsg(scip, "NLP relaxation is globally infeasible, thus can cutoff node\n");
4425  *result = SCIP_CUTOFF;
4426  return SCIP_OKAY;
4427  }
4428 
4429  if( solstat <= SCIP_NLPSOLSTAT_FEASIBLE )
4430  {
4431  /* if we have feasible NLP solution, generate linearization cuts there */
4432  SCIP_Bool lpsolseparated;
4433  SCIP_SOL* nlpsol;
4434 
4435  SCIP_CALL( SCIPcreateNLPSol(scip, &nlpsol, NULL) );
4436  assert(nlpsol != NULL);
4437 
4438  /* if we solved the NLP and solution is integral, then pass it to trysol heuristic */
4439  if( solvednlp && conshdlrdata->trysolheur != NULL )
4440  {
4441  int nfracvars = 0;
4442 
4443  if( SCIPgetNBinVars(scip) > 0 || SCIPgetNIntVars(scip) > 0 )
4444  {
4445  SCIP_CALL( SCIPgetNLPFracVars(scip, NULL, NULL, NULL, &nfracvars, NULL) );
4446  }
4447 
4448  if( nfracvars == 0 )
4449  {
4450  SCIP_CALL( SCIPheurPassSolTrySol(scip, conshdlrdata->trysolheur, nlpsol) );
4451  }
4452  }
4453 
4454  SCIP_CALL( addLinearizationCuts(scip, conshdlr, conss, nconss, nlpsol, &lpsolseparated, SCIPgetSepaMinEfficacy(scip), &cutoff) );
4455 
4456  SCIP_CALL( SCIPfreeSol(scip, &nlpsol) );
4457 
4458  if( cutoff )
4459  {
4460  *result = SCIP_CUTOFF;
4461  return SCIP_OKAY;
4462  }
4463 
4464  /* if a cut that separated the LP solution was added, then return, otherwise continue with usual separation in LP solution */
4465  if( lpsolseparated )
4466  {
4467  SCIPdebugMsg(scip, "linearization cuts separate LP solution\n");
4468 
4469  *result = SCIP_SEPARATED;
4470 
4471  return SCIP_OKAY;
4472  }
4473  }
4474  }
4475  /* if we do not want to try solving the NLP, or have no NLP, or have no NLP solver, or solving the NLP failed,
4476  * or separating with NLP solution as reference point failed, then try (again) with LP solution as reference point
4477  */
4478 
4479  SCIP_CALL( separatePoint(scip, conshdlr, conss, nconss, nusefulconss, NULL, FALSE, &cutoff, &sepasuccess) );
4480  if( cutoff )
4481  *result = SCIP_CUTOFF;
4482  else if ( sepasuccess )
4483  *result = SCIP_SEPARATED;
4484 
4485  return SCIP_OKAY;
4486 }
4487 
4488 
4489 /** separation method of constraint handler for arbitrary primal solutions */
4490 static
4491 SCIP_DECL_CONSSEPASOL(consSepasolSOC)
4492 {
4493  SCIP_CONS* maxviolcon;
4494  SCIP_Bool sepasuccess;
4495  SCIP_Bool cutoff;
4496 
4497  assert(scip != NULL);
4498  assert(conss != NULL || nconss == 0);
4499  assert(result != NULL);
4500  assert(sol != NULL);
4501 
4502  *result = SCIP_DIDNOTFIND;
4503 
4504  SCIP_CALL( computeViolations(scip, conshdlr, conss, nconss, sol, &maxviolcon) );
4505  if( maxviolcon == NULL )
4506  return SCIP_OKAY;
4507 
4508  SCIP_CALL( separatePoint(scip, conshdlr, conss, nconss, nusefulconss, sol, FALSE, &cutoff, &sepasuccess) );
4509  if( cutoff )
4510  *result = SCIP_CUTOFF;
4511  else if ( sepasuccess )
4512  *result = SCIP_SEPARATED;
4513 
4514  return SCIP_OKAY;
4515 }
4516 
4517 
4518 /** constraint enforcing method of constraint handler for LP solutions */
4519 static
4520 SCIP_DECL_CONSENFOLP(consEnfolpSOC)
4521 { /*lint --e{715}*/
4522  SCIP_CALL( enforceConstraint(scip, conshdlr, conss, nconss, nusefulconss, NULL, result) );
4523 
4524  return SCIP_OKAY;
4525 }
4526 
4527 
4528 /** constraint enforcing method of constraint handler for relaxation solutions */
4529 static
4530 SCIP_DECL_CONSENFORELAX(consEnforelaxSOC)
4531 { /*lint --e{715}*/
4532  SCIP_CALL( enforceConstraint(scip, conshdlr, conss, nconss, nusefulconss, sol, result) );
4533 
4534  return SCIP_OKAY;
4535 }
4536 
4537 
4538 /** constraint enforcing method of constraint handler for pseudo solutions */
4539 static
4540 SCIP_DECL_CONSENFOPS(consEnfopsSOC)
4542  SCIP_CONS* maxviolcons;
4543 
4544  assert(scip != NULL);
4545  assert(conss != NULL || nconss == 0);
4546  assert(result != NULL);
4547 
4548  SCIP_CALL( computeViolations(scip, conshdlr, conss, nconss, NULL, &maxviolcons) );
4549 
4550  if( maxviolcons == NULL )
4551  *result = SCIP_FEASIBLE;
4552 
4553  *result = SCIP_INFEASIBLE;
4554 
4555  return SCIP_OKAY;
4556 } /*lint !e715*/
4557 
4558 
4559 /** feasibility check method of constraint handler for integral solutions */
4560 static
4561 SCIP_DECL_CONSCHECK(consCheckSOC)
4563  SCIP_CONSHDLRDATA* conshdlrdata;
4564  SCIP_CONSDATA* consdata;
4565  SCIP_Real maxviol;
4566  SCIP_Bool dolinfeasshift;
4567  SCIP_SOL* polishedsol;
4568  int c;
4569 
4570  assert(scip != NULL);
4571  assert(conshdlr != NULL);
4572  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
4573  assert(conss != NULL || nconss == 0);
4574  assert(result != NULL );
4575 
4576  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4577  assert(conshdlrdata != NULL);
4578 
4579  *result = SCIP_FEASIBLE;
4580  maxviol = 0.0;
4581 
4582  dolinfeasshift = conshdlrdata->linfeasshift && (conshdlrdata->trysolheur != NULL);
4583  polishedsol = NULL;
4584 
4585  for( c = 0; c < nconss; ++c )
4586  {
4587  SCIP_CALL( computeViolation(scip, conshdlr, conss[c], sol) ); /*lint !e613*/
4588 
4589  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
4590  assert(consdata != NULL);
4591 
4592  /* if feasible, just continue */
4593  if( !SCIPisGT(scip, consdata->violation, SCIPfeastol(scip)) )
4594  continue;
4595 
4596  *result = SCIP_INFEASIBLE;
4597 
4598  if( consdata->violation > maxviol )
4599  maxviol = consdata->violation;
4600 
4601  if( printreason )
4602  {
4603  SCIP_CALL( SCIPprintCons(scip, conss[c], NULL) ); /*lint !e613*/
4604  SCIPinfoMessage(scip, NULL, ";\n\tviolation: %g\n", consdata->violation);
4605  }
4606 
4607  /* if we do linear feasibility shifting, then try to adjust solution */
4608  if( dolinfeasshift )
4609  {
4610  if( SCIPvarGetStatus(consdata->rhsvar) != SCIP_VARSTATUS_MULTAGGR &&
4611  !SCIPisInfinity(scip, REALABS(consdata->lhsval)) &&
4612  ( (consdata->rhscoeff > 0.0 && SCIPvarMayRoundUp (consdata->rhsvar)) ||
4613  (consdata->rhscoeff < 0.0 && SCIPvarMayRoundDown(consdata->rhsvar)) ) )
4614  {
4615  SCIP_Bool success;
4616 
4617  if( polishedsol == NULL )
4618  {
4619  if( sol != NULL )
4620  {
4621  SCIP_CALL( SCIPcreateSolCopy(scip, &polishedsol, sol) );
4622  }
4623  else
4624  {
4625  SCIP_CALL( SCIPcreateLPSol(scip, &polishedsol, NULL) );
4626  }
4627  SCIP_CALL( SCIPunlinkSol(scip, polishedsol) );
4628  }
4629  SCIP_CALL( polishSolution(scip, conss[c], polishedsol, &success) ); /*lint !e613*/
4630 
4631  /* disable solution polishing if we failed for this constraint */
4632  dolinfeasshift = success;
4633  }
4634  else /* if locks of the variable are bad or rhs is multi-aggregated, disable solution polishing */
4635  {
4636  dolinfeasshift = FALSE;
4637  }
4638  }
4639 
4640  /* if solution polishing is off and there is no NLP heuristic or we just check the LP solution,
4641  * then there is no need to check remaining constraints (NLP heuristic will pick up LP solution anyway) */
4642  if( !dolinfeasshift && (conshdlrdata->subnlpheur == NULL || sol == NULL) && !completely )
4643  break;
4644  }
4645 
4646  /* if we failed to polish solution, clear the solution */
4647  if( !dolinfeasshift && polishedsol != NULL )
4648  {
4649  SCIP_CALL( SCIPfreeSol(scip, &polishedsol) );
4650  }
4651 
4652  if( polishedsol != NULL )
4653  {
4654  assert(*result == SCIP_INFEASIBLE);
4655  SCIP_CALL( SCIPheurPassSolTrySol(scip, conshdlrdata->trysolheur, polishedsol) );
4656  SCIP_CALL( SCIPfreeSol(scip, &polishedsol) );
4657  }
4658  else if( conshdlrdata->subnlpheur != NULL && sol != NULL && *result == SCIP_INFEASIBLE && !SCIPisInfinity(scip, maxviol) )
4659  {
4660  SCIP_CALL( SCIPupdateStartpointHeurSubNlp(scip, conshdlrdata->subnlpheur, sol, maxviol) );
4661  }
4662 
4663  return SCIP_OKAY;
4664 } /*lint !e715*/
4665 
4666 
4667 /** domain propagation method of constraint handler */
4668 static
4669 SCIP_DECL_CONSPROP(consPropSOC)
4671  SCIP_RESULT propresult;
4672  int c;
4673  int nchgbds;
4674  SCIP_Bool redundant;
4675 
4676  assert(scip != NULL);
4677  assert(conss != NULL || ((nconss == 0) && (nmarkedconss == 0)));
4678  assert(result != NULL);
4679 
4680  *result = SCIP_DIDNOTFIND;
4681  nchgbds = 0;
4682 
4683  for( c = 0; c < nmarkedconss && *result != SCIP_CUTOFF; ++c )
4684  {
4685  SCIP_CALL( propagateBounds(scip, conss[c], &propresult, &nchgbds, &redundant) ); /*lint !e613*/
4686  if( propresult != SCIP_DIDNOTFIND && propresult != SCIP_DIDNOTRUN )
4687  *result = propresult;
4688  }
4689 
4690  return SCIP_OKAY;
4691 } /*lint !e715*/
4692 
4693 
4694 /** presolving method of constraint handler */
4695 static
4696 SCIP_DECL_CONSPRESOL(consPresolSOC)
4698  SCIP_CONSHDLRDATA* conshdlrdata;
4699  SCIP_CONSDATA* consdata;
4700  int c;
4701  SCIP_RESULT propresult;
4702  SCIP_Bool iscutoff;
4703  SCIP_Bool isdeleted;
4704 
4705  assert(scip != NULL);
4706  assert(conss != NULL || nconss == 0);
4707  assert(conshdlr != NULL);
4708  assert(result != NULL);
4709 
4710  *result = SCIP_DIDNOTFIND;
4711 
4712  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4713  assert(conshdlrdata != NULL);
4714 
4715  for( c = 0; c < nconss; ++c )
4716  {
4717  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
4718  assert(consdata != NULL);
4719 
4720  SCIP_CALL( presolveRemoveFixedVariables(scip, conshdlr, conss[c], ndelconss, nupgdconss, nchgbds, nfixedvars, &iscutoff, &isdeleted) ); /*lint !e613*/
4721  if( iscutoff )
4722  {
4723  *result = SCIP_CUTOFF;
4724  return SCIP_OKAY;
4725  }
4726  if( isdeleted )
4727  {
4728  /* conss[c] has been deleted */
4729  *result = SCIP_SUCCESS;
4730  continue;
4731  }
4732 
4733  if( conshdlrdata->nauxvars > 0 && !consdata->isapproxadded && consdata->nvars > 1 )
4734  {
4735  SCIP_CALL( presolveCreateOuterApprox(scip, consdata->nvars, consdata->vars, consdata->coefs, consdata->offsets,
4736  consdata->rhsvar, consdata->rhscoeff, consdata->rhscoeff, consdata->constant, SCIPconsGetName(conss[c]), conss[c],
4737  conshdlrdata->nauxvars, conshdlrdata->glineur, naddconss) ); /*lint !e613*/
4738  consdata->isapproxadded = TRUE;
4739  }
4740 
4741  if( (presoltiming & SCIP_PRESOLTIMING_FAST) != 0 )
4742  {
4743  SCIP_Bool redundant;
4744 
4745  SCIP_CALL( propagateBounds(scip, conss[c], &propresult, nchgbds, &redundant) ); /*lint !e613*/
4746  switch( propresult )
4747  {
4748  case SCIP_DIDNOTRUN:
4749  case SCIP_DIDNOTFIND:
4750  break;
4751  case SCIP_REDUCEDDOM:
4752  *result = SCIP_SUCCESS;
4753  break;
4754  case SCIP_CUTOFF:
4755  *result = SCIP_CUTOFF;
4756  SCIPdebugMsg(scip, "infeasible in presolve due to propagation for constraint %s\n", SCIPconsGetName(conss[c])); /*lint !e613*/
4757  return SCIP_OKAY;
4758  default:
4759  SCIPerrorMessage("unexpected result from propagation: %d\n", propresult);
4760  return SCIP_ERROR;
4761  } /*lint !e788*/
4762  if( redundant )
4763  ++*ndelconss;
4764  }
4765 
4766  /* disaggregate each lhs term to a quadratic constraint by using auxiliary variables */
4767  if( conshdlrdata->disaggregate && (presoltiming & SCIP_PRESOLTIMING_EXHAUSTIVE) != 0 )
4768  {
4769  SCIP_Bool success;
4770 
4771  SCIP_CALL( disaggregate(scip, conss[c], consdata, naddconss, ndelconss, &success) ); /*lint !e613*/
4772 
4773  if( success )
4774  {
4775  SCIPdebugMsg(scip, "disaggregated SOC constraint\n");
4776 
4777  /* conss[c] has been deleted */
4778  *result = SCIP_SUCCESS;
4779  continue;
4780  }
4781  }
4782  }
4783 
4784  return SCIP_OKAY;
4785 } /*lint !e715*/
4786 
4787 
4788 /** variable rounding lock method of constraint handler */
4789 static
4790 SCIP_DECL_CONSLOCK(consLockSOC)
4792  SCIP_CONSDATA* consdata;
4793  int i;
4794 
4795  assert(scip != NULL);
4796  assert(conshdlr != NULL);
4797  assert(cons != NULL);
4798  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
4799  assert(locktype == SCIP_LOCKTYPE_MODEL);
4800 
4801  consdata = SCIPconsGetData(cons);
4802  assert(consdata != NULL);
4803 
4804  SCIPdebugMsg(scip, "Locking constraint <%s>.\n", SCIPconsGetName(cons));
4805 
4806  /* Changing variables x_i, i <= n, in both directions can lead to an infeasible solution. */
4807  for( i = 0; i < consdata->nvars; ++i )
4808  {
4809  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->vars[i], locktype, nlockspos + nlocksneg, nlockspos + nlocksneg) );
4810  }
4811 
4812  /* Rounding x_{n+1} up will not violate a solution. */
4813  if( consdata->rhsvar != NULL )
4814  {
4815  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->rhsvar, locktype, consdata->rhscoeff > 0.0 ? nlockspos : nlocksneg, consdata->rhscoeff > 0.0 ? nlocksneg : nlockspos) );
4816  }
4817 
4818  return SCIP_OKAY;
4819 }
4820 
4821 /** constraint display method of constraint handler */
4822 static
4823 SCIP_DECL_CONSPRINT(consPrintSOC)
4824 {
4825  SCIP_CONSDATA* consdata;
4826  int i;
4827 
4828  assert(scip != NULL);
4829  assert(conshdlr != NULL);
4830  assert(cons != NULL);
4831  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
4832 
4833  consdata = SCIPconsGetData(cons);
4834  assert(consdata != NULL);
4835 
4836  SCIPinfoMessage(scip, file, "sqrt( ");
4837  if( consdata->constant != 0.0 )
4838  {
4839  SCIPinfoMessage(scip, file, "%.15g", consdata->constant);
4840  }
4841 
4842  for( i = 0; i < consdata->nvars; ++i )
4843  {
4844  SCIPinfoMessage(scip, file, "+ (%.15g*(", consdata->coefs[i]);
4845  SCIP_CALL( SCIPwriteVarName(scip, file, consdata->vars[i], TRUE) );
4846  SCIPinfoMessage(scip, file, "%+.15g))^2 ", consdata->offsets[i]);
4847  }
4848 
4849  SCIPinfoMessage(scip, file, ") <= ");
4850  if( consdata->rhsvar != NULL )
4851  {
4852  SCIPinfoMessage(scip, file, "%.15g*(", consdata->rhscoeff);
4853  SCIP_CALL( SCIPwriteVarName(scip, file, consdata->rhsvar, TRUE) );
4854  SCIPinfoMessage(scip, file, "%+.15g)", consdata->rhsoffset);
4855  }
4856  else
4857  {
4858  SCIPinfoMessage(scip, file, "%.15g", consdata->rhscoeff*consdata->rhsoffset);
4859  }
4860 
4861  return SCIP_OKAY;
4862 }
4863 
4864 /** constraint copying method of constraint handler */
4865 static
4866 SCIP_DECL_CONSCOPY(consCopySOC)
4868  SCIP_CONSDATA* consdata;
4869  SCIP_VAR** vars;
4870  SCIP_VAR* rhsvar;
4871  int i;
4872 
4873  assert(scip != NULL);
4874  assert(cons != NULL);
4875  assert(sourcescip != NULL);
4876  assert(sourceconshdlr != NULL);
4877  assert(sourcecons != NULL);
4878  assert(varmap != NULL);
4879  assert(valid != NULL);
4880  assert(stickingatnode == FALSE);
4881 
4882  consdata = SCIPconsGetData(sourcecons);
4883  assert(consdata != NULL);
4884 
4885  *valid = TRUE;
4886 
4887  SCIP_CALL( SCIPallocBufferArray(sourcescip, &vars, consdata->nvars) );
4888 
4889  /* map variables to active variables of the target SCIP */
4890  for( i = 0; i < consdata->nvars && *valid; ++i )
4891  {
4892  SCIP_CALL( SCIPgetVarCopy(sourcescip, scip, consdata->vars[i], &vars[i], varmap, consmap, global, valid) );
4893  assert(!(*valid) || vars[i] != NULL);
4894  }
4895 
4896  /* map rhs variable to active variable of the target SCIP */
4897  if( *valid )
4898  {
4899  SCIP_CALL( SCIPgetVarCopy(sourcescip, scip, consdata->rhsvar, &rhsvar, varmap, consmap, global, valid) );
4900  assert(!(*valid) || rhsvar != NULL);
4901 
4902  /* only create the target constraint, if all variables could be copied */
4903  if( *valid )
4904  {
4905  SCIP_CALL( SCIPcreateConsSOC(scip, cons, name ? name : SCIPconsGetName(sourcecons),
4906  consdata->nvars, vars, consdata->coefs, consdata->offsets, consdata->constant,
4907  rhsvar, consdata->rhscoeff, consdata->rhsoffset,
4908  initial, separate, enforce, check, propagate, local, modifiable, dynamic, removable) ); /*lint !e644 */
4909  }
4910  }
4911 
4912  SCIPfreeBufferArray(sourcescip, &vars);
4913 
4914  return SCIP_OKAY;
4915 }
4916 
4917 
4918 /** constraint parsing method of constraint handler */
4919 static
4920 SCIP_DECL_CONSPARSE(consParseSOC)
4921 { /*lint --e{715}*/
4922  SCIP_VAR* var;
4923  SCIP_VAR** vars;
4924  SCIP_Real* coefs;
4925  SCIP_Real* offsets;
4926  int nvars;
4927  int varssize;
4928  SCIP_VAR* rhsvar;
4929  SCIP_Real rhscoef;
4930  SCIP_Real rhsoffset;
4931  SCIP_Real constant;
4932  SCIP_Real coef;
4933  SCIP_Real offset;
4934  char* endptr;
4935 
4936  assert(scip != NULL);
4937  assert(success != NULL);
4938  assert(str != NULL);
4939  assert(name != NULL);
4940  assert(cons != NULL);
4941 
4942  /* check that string starts with "sqrt( " */
4943  if( strncmp(str, "sqrt( ", 6) != 0 )
4944  {
4945  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected 'sqrt( ' at begin of soc constraint string '%s'\n", str);
4946  *success = FALSE;
4947  return SCIP_OKAY;
4948  }
4949  str += 6;
4950 
4951  *success = TRUE;
4952 
4953  /* check if we have a constant in the beginning */
4954  if( SCIPstrToRealValue(str, &constant, &endptr) )
4955  str = endptr;
4956  else
4957  constant = 0.0;
4958 
4959  nvars = 0;
4960  varssize = 5;
4961  SCIP_CALL( SCIPallocBufferArray(scip, &vars, varssize) );
4962  SCIP_CALL( SCIPallocBufferArray(scip, &coefs, varssize) );
4963  SCIP_CALL( SCIPallocBufferArray(scip, &offsets, varssize) );
4964 
4965  /* read '+ (coef*(var+offset))^2' on lhs, as long as possible */
4966  while( *str != '\0' )
4967  {
4968  /* skip whitespace */
4969  while( isspace((int)*str) )
4970  ++str;
4971 
4972  /* stop if no more coefficients */
4973  if( strncmp(str, "+ (", 3) != 0 )
4974  break;
4975 
4976  str += 3;
4977 
4978  /* parse coef */
4979  if( !SCIPstrToRealValue(str, &coef, &endptr) )
4980  {
4981  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected coefficient at begin of '%s'\n", str);
4982  *success = FALSE;
4983  break;
4984  }
4985  str = endptr;
4986 
4987  if( strncmp(str, "*(", 2) != 0 )
4988  {
4989  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected '*(' at begin of '%s'\n", str);
4990  *success = FALSE;
4991  break;
4992  }
4993  str += 2;
4994 
4995  /* parse variable name */
4996  SCIP_CALL( SCIPparseVarName(scip, str, &var, &endptr) );
4997  if( var == NULL )
4998  {
4999  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "unknown variable name at '%s'\n", str);
5000  *success = FALSE;
5001  break;
5002  }
5003  str = endptr;
5004 
5005  /* parse offset */
5006  if( !SCIPstrToRealValue(str, &offset, &endptr) )
5007  {
5008  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected offset at begin of '%s'\n", str);
5009  *success = FALSE;
5010  break;
5011  }
5012  str = endptr;
5013 
5014  if( strncmp(str, "))^2", 4) != 0 )
5015  {
5016  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected '))^2' at begin of '%s'\n", str);
5017  *success = FALSE;
5018  break;
5019  }
5020  str += 4;
5021 
5022  if( varssize <= nvars )
5023  {
5024  varssize = SCIPcalcMemGrowSize(scip, varssize+1);
5025  SCIP_CALL( SCIPreallocBufferArray(scip, &vars, varssize) );
5026  SCIP_CALL( SCIPreallocBufferArray(scip, &coefs, varssize) );
5027  SCIP_CALL( SCIPreallocBufferArray(scip, &offsets, varssize) );
5028  }
5029  vars[nvars] = var;
5030  coefs[nvars] = coef;
5031  offsets[nvars] = offset;
5032  ++nvars;
5033  }
5034 
5035  if( strncmp(str, ") <=", 4) != 0 )
5036  {
5037  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected ') <=' at begin of '%s'\n", str);
5038  *success = FALSE;
5039  }
5040  str += 4;
5041 
5042  /* read rhs coef*(var+offset) or just a constant */
5043 
5044  /* parse coef */
5045  if( *success )
5046  {
5047  /* skip whitespace */
5048  while( isspace((int)*str) )
5049  ++str;
5050 
5051  if( !SCIPstrToRealValue(str, &rhscoef, &endptr) )
5052  {
5053  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected coefficient at begin of '%s'\n", str);
5054  *success = FALSE;
5055  }
5056  str = endptr;
5057 
5058  /* skip whitespace */
5059  while( isspace((int)*str) )
5060  ++str;
5061  }
5062 
5063  /* parse *(var+offset) */
5064  if( *str != '\0' )
5065  {
5066  if( *success )
5067  {
5068  if( strncmp(str, "*(", 2) != 0 )
5069  {
5070  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected '*(' at begin of '%s'\n", str);
5071  *success = FALSE;
5072  }
5073  else
5074  {
5075  str += 2;
5076  }
5077  }
5078 
5079  /* parse variable name */
5080  if( *success )
5081  {
5082  SCIP_CALL( SCIPparseVarName(scip, str, &rhsvar, &endptr) );
5083  if( rhsvar == NULL )
5084  {
5085  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "unknown variable name at '%s'\n", str);
5086  *success = FALSE;
5087  }
5088  else
5089  {
5090  str = endptr;
5091  }
5092  }
5093 
5094  /* parse offset */
5095  if( *success )
5096  {
5097  if( !SCIPstrToRealValue(str, &rhsoffset, &endptr) )
5098  {
5099  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected offset at begin of '%s'\n", str);
5100  *success = FALSE;
5101  }
5102  else
5103  {
5104  str = endptr;
5105  }
5106  }
5107 
5108  if( *success )
5109  {
5110  if( *str != ')' )
5111  {
5112  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected ')' at begin of '%s'\n", str);
5113  *success = FALSE;
5114  }
5115  }
5116  }
5117  else if( *success )
5118  {
5119  /* only a constant at right hand side */
5120  rhsoffset = rhscoef; /*lint !e644*/
5121  rhscoef = 1.0;
5122  rhsvar = NULL;
5123  }
5124 
5125  if( *success )
5126  {
5127  assert(!stickingatnode);
5128  SCIP_CALL( SCIPcreateConsSOC(scip, cons, name, nvars, vars, coefs, offsets, constant, rhsvar, rhscoef, rhsoffset,
5129  initial, separate, enforce, check, propagate, local, modifiable, dynamic, removable) ); /*lint !e644 */
5130  }
5131 
5132  SCIPfreeBufferArray(scip, &offsets);
5133  SCIPfreeBufferArray(scip, &coefs);
5134  SCIPfreeBufferArray(scip, &vars);
5135 
5136  return SCIP_OKAY;
5137 }
5138 
5139 /** constraint method of constraint handler which returns the variables (if possible) */
5140 static
5141 SCIP_DECL_CONSGETVARS(consGetVarsSOC)
5142 { /*lint --e{715}*/
5143  SCIP_CONSDATA* consdata;
5144 
5145  consdata = SCIPconsGetData(cons);
5146  assert(consdata != NULL);
5147 
5148  if( varssize < consdata->nvars + 1 )
5149  (*success) = FALSE;
5150  else
5151  {
5152  BMScopyMemoryArray(vars, consdata->vars, consdata->nvars);
5153  vars[consdata->nvars] = consdata->rhsvar;
5154  (*success) = TRUE;
5155  }
5156 
5157  return SCIP_OKAY;
5158 }
5159 
5160 /** constraint method of constraint handler which returns the number of variable (if possible) */
5161 static
5162 SCIP_DECL_CONSGETNVARS(consGetNVarsSOC)
5163 { /*lint --e{715}*/
5164  SCIP_CONSDATA* consdata;
5165 
5166  assert(cons != NULL);
5167 
5168  consdata = SCIPconsGetData(cons);
5169  assert(consdata != NULL);
5170 
5171  (*nvars) = consdata->nvars + 1;
5172  (*success) = TRUE;
5173 
5174  return SCIP_OKAY;
5175 }
5176 
5177 /*
5178  * constraint specific interface methods
5179  */
5180 
5181 /** creates the handler for second order cone constraints and includes it in SCIP */
5183  SCIP* scip /**< SCIP data structure */
5184  )
5185 {
5186  SCIP_CONSHDLRDATA* conshdlrdata;
5187  SCIP_CONSHDLR* conshdlr;
5188  SCIP_EVENTHDLR* eventhdlr;
5189 
5190  /* create constraint handler data */
5191  SCIP_CALL( SCIPallocBlockMemory(scip, &conshdlrdata) );
5192  conshdlrdata->subnlpheur = NULL;
5193  conshdlrdata->trysolheur = NULL;
5194 
5195  SCIP_CALL( SCIPincludeEventhdlrBasic(scip, &eventhdlr, CONSHDLR_NAME"_boundchange",
5196  "signals a bound change to a second order cone constraint",
5197  processVarEvent, NULL) );
5198  conshdlrdata->eventhdlr = eventhdlr;
5199 
5200  SCIP_CALL( SCIPincludeEventhdlrBasic(scip, NULL, CONSHDLR_NAME"_newsolution",
5201  "handles the event that a new primal solution has been found",
5202  processNewSolutionEvent, NULL) );
5203 
5204  /* include constraint handler */
5207  consEnfolpSOC, consEnfopsSOC, consCheckSOC, consLockSOC,
5208  conshdlrdata) );
5209  assert(conshdlr != NULL);
5210 
5211  /* set non-fundamental callbacks via specific setter functions */
5212  SCIP_CALL( SCIPsetConshdlrCopy(scip, conshdlr, conshdlrCopySOC, consCopySOC) );
5213  SCIP_CALL( SCIPsetConshdlrDelete(scip, conshdlr, consDeleteSOC) );
5214  SCIP_CALL( SCIPsetConshdlrExit(scip, conshdlr, consExitSOC) );
5215  SCIP_CALL( SCIPsetConshdlrExitpre(scip, conshdlr, consExitpreSOC) );
5216  SCIP_CALL( SCIPsetConshdlrExitsol(scip, conshdlr, consExitsolSOC) );
5217  SCIP_CALL( SCIPsetConshdlrFree(scip, conshdlr, consFreeSOC) );
5218  SCIP_CALL( SCIPsetConshdlrGetVars(scip, conshdlr, consGetVarsSOC) );
5219  SCIP_CALL( SCIPsetConshdlrGetNVars(scip, conshdlr, consGetNVarsSOC) );
5220  SCIP_CALL( SCIPsetConshdlrInit(scip, conshdlr, consInitSOC) );
5221  SCIP_CALL( SCIPsetConshdlrInitsol(scip, conshdlr, consInitsolSOC) );
5222  SCIP_CALL( SCIPsetConshdlrParse(scip, conshdlr, consParseSOC) );
5223  SCIP_CALL( SCIPsetConshdlrPresol(scip, conshdlr, consPresolSOC, CONSHDLR_MAXPREROUNDS, CONSHDLR_PRESOLTIMING) );
5224  SCIP_CALL( SCIPsetConshdlrPrint(scip, conshdlr, consPrintSOC) );
5226  SCIP_CALL( SCIPsetConshdlrSepa(scip, conshdlr, consSepalpSOC, consSepasolSOC, CONSHDLR_SEPAFREQ,
5228  SCIP_CALL( SCIPsetConshdlrTrans(scip, conshdlr, consTransSOC) );
5229  SCIP_CALL( SCIPsetConshdlrEnforelax(scip, conshdlr, consEnforelaxSOC) );
5230 
5231  if( SCIPfindConshdlr(scip,"quadratic") != NULL )
5232  {
5233  /* notify function that upgrades quadratic constraint to SOC's */
5235  }
5236 
5237  /* add soc constraint handler parameters */
5238  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/projectpoint",
5239  "whether the reference point of a cut should be projected onto the feasible set of the SOC constraint",
5240  &conshdlrdata->projectpoint, TRUE, FALSE, NULL, NULL) );
5241 
5242  SCIP_CALL( SCIPaddIntParam (scip, "constraints/" CONSHDLR_NAME "/nauxvars",
5243  "number of auxiliary variables to use when creating a linear outer approx. of a SOC3 constraint; 0 to turn off",
5244  &conshdlrdata->nauxvars, FALSE, 0, 0, INT_MAX, NULL, NULL) );
5245 
5246  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/glineur",
5247  "whether the Glineur Outer Approximation should be used instead of Ben-Tal Nemirovski",
5248  &conshdlrdata->glineur, FALSE, TRUE, NULL, NULL) );
5249 
5250  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/sparsify",
5251  "whether to sparsify cuts",
5252  &conshdlrdata->sparsify, TRUE, FALSE, NULL, NULL) );
5253 
5254  SCIP_CALL( SCIPaddRealParam(scip, "constraints/" CONSHDLR_NAME "/sparsifymaxloss",
5255  "maximal loss in cut efficacy by sparsification",
5256  &conshdlrdata->sparsifymaxloss, TRUE, 0.2, 0.0, 1.0, NULL, NULL) );
5257 
5258  SCIP_CALL( SCIPaddRealParam(scip, "constraints/" CONSHDLR_NAME "/sparsifynzgrowth",
5259  "growth rate of maximal allowed nonzeros in cuts in sparsification",
5260  &conshdlrdata->sparsifynzgrowth, TRUE, 1.3, 1.000001, SCIPinfinity(scip), NULL, NULL) );
5261 
5262  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/linfeasshift",
5263  "whether to try to make solutions feasible in check by shifting the variable on the right hand side",
5264  &conshdlrdata->linfeasshift, FALSE, TRUE, NULL, NULL) );
5265 
5266  SCIP_CALL( SCIPaddCharParam(scip, "constraints/" CONSHDLR_NAME "/nlpform",
5267  "which formulation to use when adding a SOC constraint to the NLP (a: automatic, q: nonconvex quadratic form, s: convex sqrt form, e: convex exponential-sqrt form, d: convex division form)",
5268  &conshdlrdata->nlpform, FALSE, 'a', "aqsed", NULL, NULL) );
5269 
5270  SCIP_CALL( SCIPaddRealParam(scip, "constraints/" CONSHDLR_NAME "/sepanlpmincont",
5271  "minimal required fraction of continuous variables in problem to use solution of NLP relaxation in root for separation",
5272  &conshdlrdata->sepanlpmincont, FALSE, 1.0, 0.0, 2.0, NULL, NULL) );
5273 
5274  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/enfocutsremovable",
5275  "are cuts added during enforcement removable from the LP in the same node?",
5276  &conshdlrdata->enfocutsremovable, TRUE, FALSE, NULL, NULL) );
5277 
5278  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/generalsocupgrade",
5279  "try to upgrade more general quadratics to soc?",
5280  &conshdlrdata->generalsocupg, TRUE, TRUE, NULL, NULL) );
5281 
5282  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/disaggregate",
5283  "try to completely disaggregate soc?",
5284  &conshdlrdata->disaggregate, TRUE, TRUE, NULL, NULL) );
5285 
5286  return SCIP_OKAY;
5287 }
5288 
5289 /** creates and captures a second order cone constraint
5290  *
5291  * @note the constraint gets captured, hence at one point you have to release it using the method SCIPreleaseCons()
5292  */
5294  SCIP* scip, /**< SCIP data structure */
5295  SCIP_CONS** cons, /**< pointer to hold the created constraint */
5296  const char* name, /**< name of constraint */
5297  int nvars, /**< number of variables on left hand side of constraint (n) */
5298  SCIP_VAR** vars, /**< array with variables on left hand side (x_i) */
5299  SCIP_Real* coefs, /**< array with coefficients of left hand side variables (alpha_i), or NULL if all 1.0 */
5300  SCIP_Real* offsets, /**< array with offsets of variables (beta_i), or NULL if all 0.0 */
5301  SCIP_Real constant, /**< constant on left hand side (gamma) */
5302  SCIP_VAR* rhsvar, /**< variable on right hand side of constraint (x_{n+1}) */
5303  SCIP_Real rhscoeff, /**< coefficient of variable on right hand side (alpha_{n+1}) */
5304  SCIP_Real rhsoffset, /**< offset of variable on right hand side (beta_{n+1}) */
5305  SCIP_Bool initial, /**< should the LP relaxation of constraint be in the initial LP?
5306  * Usually set to TRUE. Set to FALSE for 'lazy constraints'. */
5307  SCIP_Bool separate, /**< should the constraint be separated during LP processing?
5308  * Usually set to TRUE. */
5309  SCIP_Bool enforce, /**< should the constraint be enforced during node processing?
5310  * TRUE for model constraints, FALSE for additional, redundant constraints. */
5311  SCIP_Bool check, /**< should the constraint be checked for feasibility?
5312  * TRUE for model constraints, FALSE for additional, redundant constraints. */
5313  SCIP_Bool propagate, /**< should the constraint be propagated during node processing?
5314  * Usually set to TRUE. */
5315  SCIP_Bool local, /**< is constraint only valid locally?
5316  * Usually set to FALSE. Has to be set to TRUE, e.g., for branching constraints. */
5317  SCIP_Bool modifiable, /**< is constraint modifiable (subject to column generation)?
5318  * Usually set to FALSE. In column generation applications, set to TRUE if pricing
5319  * adds coefficients to this constraint. */
5320  SCIP_Bool dynamic, /**< is constraint subject to aging?
5321  * Usually set to FALSE. Set to TRUE for own cuts which
5322  * are separated as constraints. */
5323  SCIP_Bool removable /**< should the relaxation be removed from the LP due to aging or cleanup?
5324  * Usually set to FALSE. Set to TRUE for 'lazy constraints' and 'user cuts'. */
5325  )
5326 {
5327  SCIP_CONSHDLR* conshdlr;
5328  SCIP_CONSDATA* consdata;
5329  int i;
5330 
5331  assert(scip != NULL);
5332  assert(cons != NULL);
5333  assert(modifiable == FALSE); /* we do not support column generation */
5334 
5335  /* find the soc constraint handler */
5336  conshdlr = SCIPfindConshdlr(scip, CONSHDLR_NAME);
5337  if( conshdlr == NULL )
5338  {
5339  SCIPerrorMessage("SOC constraint handler not found\n");
5340  return SCIP_PLUGINNOTFOUND;
5341  }
5342 
5343  assert(vars != NULL);
5344  assert(nvars >= 1);
5345  assert(constant >= 0.0);
5346  assert(!SCIPisInfinity(scip, ABS(rhsoffset)));
5347  assert(!SCIPisInfinity(scip, constant));
5348  assert(rhsvar == NULL || rhscoeff <= 0.0 || SCIPisGE(scip, local ? SCIPcomputeVarLbLocal(scip, rhsvar) : SCIPcomputeVarLbGlobal(scip, rhsvar), -rhsoffset));
5349  assert(rhsvar == NULL || rhscoeff >= 0.0 || SCIPisLE(scip, local ? SCIPcomputeVarUbLocal(scip, rhsvar) : SCIPcomputeVarUbGlobal(scip, rhsvar), -rhsoffset));
5350 
5351  /* create constraint data */
5352  SCIP_CALL( SCIPallocBlockMemory(scip, &consdata) );
5353 
5354  consdata->nvars = nvars;
5355  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &consdata->vars, vars, nvars) );
5356  for( i = 0; i < nvars; ++i )
5357  {
5358  SCIP_CALL( SCIPcaptureVar(scip, vars[i]) );
5359  }
5360 
5361  if( coefs != NULL )
5362  {
5363  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &consdata->coefs, coefs, nvars) );
5364  for( i = 0; i < nvars; ++i )
5365  {
5366  if( consdata->coefs[i] < 0.0 )
5367  consdata->coefs[i] = -consdata->coefs[i];
5368  }
5369  }
5370  else
5371  {
5372  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->coefs, nvars) );
5373  for( i = 0; i < nvars; ++i )
5374  consdata->coefs[i] = 1.0;
5375  }
5376 
5377  if( offsets != NULL )
5378  {
5379  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &consdata->offsets, offsets, nvars) );
5380  }
5381  else
5382  {
5383  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->offsets, nvars) );
5384  BMSclearMemoryArray(consdata->offsets, nvars);
5385  }
5386 
5387  consdata->constant = constant;
5388  consdata->rhsvar = rhsvar;
5389  consdata->rhscoeff = rhscoeff;
5390  consdata->rhsoffset = rhsoffset;
5391 
5392  if( rhsvar != NULL )
5393  {
5394  SCIP_CALL( SCIPcaptureVar(scip, rhsvar) );
5395  }
5396 
5397  consdata->nlrow = NULL;
5398 
5399  consdata->lhsbndchgeventdata = NULL;
5400  consdata->isapproxadded = FALSE;
5401 
5402  /* create constraint */
5403  SCIP_CALL( SCIPcreateCons(scip, cons, name, conshdlr, consdata, initial, separate, enforce, check, propagate,
5404  local, modifiable, dynamic, removable, FALSE) );
5405 
5406  if( SCIPisTransformed(scip) )
5407  {
5408  SCIP_CONSHDLRDATA* conshdlrdata = SCIPconshdlrGetData(conshdlr);
5409  assert(conshdlrdata != NULL);
5410 
5411  SCIP_CALL( catchVarEvents(scip, conshdlrdata->eventhdlr, *cons) );