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-2020 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(SCIPisPositive(scip, consdata->lhsval)); /* do not like to linearize in 0 */
878  assert(!SCIPisInfinity(scip, consdata->lhsval));
879 
881  SCIP_CALL( SCIPensureRowprepSize(scip, *rowprep, consdata->nvars+1) );
882  (void) SCIPsnprintf((*rowprep)->name, SCIP_MAXSTRLEN, "%s_linearization_%d", SCIPconsGetName(cons), SCIPgetNLPs(scip));
883 
884  for( i = 0; i < consdata->nvars; ++i )
885  {
886  val = SCIPgetSolVal(scip, sol, consdata->vars[i]) + consdata->offsets[i];
887  val *= consdata->coefs[i] * consdata->coefs[i];
888 
889  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, consdata->vars[i], val / consdata->lhsval) );
890 
891  val *= SCIPgetSolVal(scip, sol, consdata->vars[i]);
892  SCIPaddRowprepSide(*rowprep, val);
893  }
894  (*rowprep)->side /= consdata->lhsval;
895  (*rowprep)->side -= consdata->lhsval - consdata->rhscoeff * consdata->rhsoffset;
896 
897  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, consdata->rhsvar, -consdata->rhscoeff) );
898 
899  return SCIP_OKAY;
900 }
901 
902 /** generate supporting hyperplane in a given point */
903 static
905  SCIP* scip, /**< SCIP pointer */
906  SCIP_CONS* cons, /**< constraint */
907  SCIP_Real* x, /**< point (lhs-vars) where to generate cut */
908  SCIP_ROWPREP** rowprep /**< place to store cut */
909  )
910 {
911  SCIP_CONSDATA* consdata;
912  SCIP_Real lhsval;
913  SCIP_Real val;
914  int i;
915 
916  assert(scip != NULL);
917  assert(cons != NULL);
918  assert(rowprep != NULL);
919 
920  consdata = SCIPconsGetData(cons);
921  assert(consdata != NULL);
922 
923  lhsval = consdata->constant;
924  for( i = 0; i < consdata->nvars; ++i )
925  {
926  assert(!SCIPisInfinity(scip, ABS(x[i])));
927 
928  val = consdata->coefs[i] * (x[i] + consdata->offsets[i]);
929  lhsval += val * val;
930  }
931  lhsval = sqrt(lhsval);
932 
933  if( SCIPisZero(scip, lhsval) )
934  { /* do not like to linearize in 0 */
935  return SCIP_OKAY;
936  }
937 
939  SCIP_CALL( SCIPensureRowprepSize(scip, *rowprep, consdata->nvars+1) );
940  (void) SCIPsnprintf((*rowprep)->name, SCIP_MAXSTRLEN, "%s_linearization_%d", SCIPconsGetName(cons), SCIPgetNLPs(scip));
941 
942  for( i = 0; i < consdata->nvars; ++i )
943  {
944  val = x[i] + consdata->offsets[i];
945  if( SCIPisZero(scip, val) )
946  continue;
947 
948  val *= consdata->coefs[i] * consdata->coefs[i];
949 
950  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, consdata->vars[i], val / lhsval) );
951 
952  val *= x[i];
953  SCIPaddRowprepSide(*rowprep, val);
954  }
955  (*rowprep)->side /= lhsval;
956  (*rowprep)->side -= lhsval - consdata->rhscoeff * consdata->rhsoffset;
957 
958  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, consdata->rhsvar, -consdata->rhscoeff) );
959 
960  return SCIP_OKAY;
961 }
962 
963 /** generate supporting hyperplane w.r.t. solution projected on feasible set
964  *
965  * Instead of linearizing the SOC constraint in the given solution point, this function projects the point
966  * first onto the feasible set of the SOC constraint (w.r.t. euclidean norm (scaled by alpha))
967  * and linearizes the SOC constraint there.
968  * The hope is that this produces somewhat tighter cuts.
969  *
970  * The projection has only be computed for the case gamma = 0.
971  * If gamma > 0, generateCut is called.
972  *
973  * Let \f$\hat x\f$ be sol or the LP solution if sol == NULL.
974  * Then the projected point \f$\tilde x\f$ is given by
975  * \f[
976  * \tilde x_i = \frac{\hat x_i + \lambda \beta_i}{1-\lambda}, \quad i=1,\ldots, n; \quad
977  * \tilde x_{n+1} = \frac{\hat x_{n+1} - \lambda \beta_{n+1}}{1+\lambda}
978  * \f]
979  * where
980  * \f[
981  * \lambda = \frac{1-A}{1+A}, \qquad
982  * 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}}
983  * \f]
984  *
985  * If lambda is very close to 1, generateCut is called.
986  *
987  * The generated cut is very similar to the unprojected form.
988  * The only difference is in the right hand side, which is (in the case beta = 0) multiplied by 1/(1-lambda).
989  */
990 static
992  SCIP* scip, /**< SCIP pointer */
993  SCIP_CONS* cons, /**< constraint */
994  SCIP_SOL* sol, /**< solution to separate, or NULL for LP solution */
995  SCIP_ROWPREP** rowprep /**< place to store cut */
996  )
997 {
998  SCIP_CONSDATA* consdata;
999  SCIP_Real val;
1000  SCIP_Real A, lambda;
1001  int i;
1002 
1003  assert(scip != NULL);
1004  assert(cons != NULL);
1005  assert(rowprep != NULL);
1006 
1007  consdata = SCIPconsGetData(cons);
1008  assert(consdata != NULL);
1009 
1010  assert(SCIPisPositive(scip, consdata->lhsval)); /* do not like to linearize in 0 */
1011  assert(!SCIPisInfinity(scip, consdata->lhsval));
1012 
1013  if( !SCIPisZero(scip, consdata->constant) )
1014  { /* have not thought about this case yet */
1015  SCIP_CALL( generateCutSol(scip, cons, sol, rowprep) );
1016  return SCIP_OKAY;
1017  }
1018 
1019  A = consdata->rhscoeff * (SCIPgetSolVal(scip, sol, consdata->rhsvar) + consdata->rhsoffset);
1020  A /= consdata->lhsval;
1021 
1022  lambda = (1.0 - A) / (1.0 + A);
1023 
1024  assert(!SCIPisNegative(scip, lambda)); /* otherwise A > 1, so constraint is not violated */
1025 
1026  SCIPdebugMsg(scip, "A = %g \t lambda = %g\n", A, lambda);
1027 
1028  if( SCIPisFeasEQ(scip, lambda, 1.0) )
1029  { /* avoid numerical difficulties when dividing by (1-lambda) below */
1030  SCIP_CALL( generateCutSol(scip, cons, sol, rowprep) );
1031  return SCIP_OKAY;
1032  }
1033 
1035  SCIP_CALL( SCIPensureRowprepSize(scip, *rowprep, consdata->nvars+1) );
1036  (void) SCIPsnprintf((*rowprep)->name, SCIP_MAXSTRLEN, "%s_linearization_%d", SCIPconsGetName(cons), SCIPgetNLPs(scip));
1037 
1038  for( i = 0; i < consdata->nvars; ++i )
1039  {
1040  val = SCIPgetSolVal(scip, sol, consdata->vars[i]) + consdata->offsets[i];
1041  val *= consdata->coefs[i] * consdata->coefs[i];
1042 
1043  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, consdata->vars[i], val / consdata->lhsval) );
1044 
1045  val *= SCIPgetSolVal(scip, sol, consdata->vars[i]) + lambda * consdata->offsets[i];
1046  SCIPaddRowprepSide(*rowprep, val);
1047  }
1048  (*rowprep)->side /= consdata->lhsval;
1049  (*rowprep)->side-= consdata->lhsval;
1050  (*rowprep)->side /= 1.0 - lambda;
1051  (*rowprep)->side -= consdata->rhscoeff * consdata->rhsoffset;
1052 
1053  SCIP_CALL( SCIPaddRowprepTerm(scip, *rowprep, consdata->rhsvar, -consdata->rhscoeff) );
1054 
1055  return SCIP_OKAY;
1056 }
1057 
1058 /** generates sparsified supporting hyperplane */
1059 static
1061  SCIP* scip, /**< SCIP pointer */
1062  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
1063  SCIP_CONS* cons, /**< constraint */
1064  SCIP_SOL* sol, /**< solution to separate, or NULL for LP solution */
1065  SCIP_ROWPREP** rowprep, /**< place to store cut */
1066  SCIP_Real minefficacy /**< minimal efficacy for a cut to be accepted */
1067  )
1068 {
1069  SCIP_CONSHDLRDATA* conshdlrdata;
1070  SCIP_CONSDATA* consdata;
1071  SCIP_Real* x;
1072  SCIP_Real* dist; /* distance to 0 */
1073  int* ind; /* indices */
1074  int i;
1075  int maxnz, nextmaxnz;
1076  SCIP_Real efficacy;
1077  SCIP_Real goodefficacy;
1078 
1079  assert(scip != NULL);
1080  assert(conshdlr != NULL);
1081  assert(cons != NULL);
1082  assert(rowprep != NULL);
1083 
1084  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1085  assert(conshdlrdata != NULL);
1086 
1087  consdata = SCIPconsGetData(cons);
1088  assert(consdata != NULL);
1089 
1090  assert(SCIPisPositive(scip, consdata->lhsval)); /* do not like to linearize in 0 */
1091  assert(!SCIPisInfinity(scip, consdata->lhsval));
1092 
1093  if( consdata->nvars <= 3 )
1094  {
1095  SCIP_CALL( generateCutSol(scip, cons, sol, rowprep) );
1096  return SCIP_OKAY;
1097  }
1098 
1099  goodefficacy = MAX((1.0-conshdlrdata->sparsifymaxloss) * consdata->violation, minefficacy);
1100 
1101  SCIP_CALL( SCIPallocBufferArray(scip, &x, consdata->nvars) );
1102  SCIP_CALL( SCIPallocBufferArray(scip, &dist, consdata->nvars) );
1103  SCIP_CALL( SCIPallocBufferArray(scip, &ind, consdata->nvars) );
1104 
1105  SCIP_CALL( SCIPgetSolVals(scip, sol, consdata->nvars, consdata->vars, x) );
1106  /* distance to "-offset" * alpha_i^2 should indicate loss when moving refpoint to x[i] = -offset[i] */
1107  for( i = 0; i < consdata->nvars; ++i )
1108  {
1109  ind[i] = i;
1110  dist[i] = ABS(x[i] + consdata->offsets[i]);
1111  dist[i] *= consdata->coefs[i] * consdata->coefs[i];
1112  }
1113 
1114  /* sort variables according to dist */
1115  SCIPsortRealInt(dist, ind, consdata->nvars);
1116 
1117  maxnz = 2;
1118  /* set all except biggest maxnz entries in x to -offset */
1119  for( i = 0; i < consdata->nvars - maxnz; ++i )
1120  x[ind[i]] = -consdata->offsets[i];
1121 
1122  do
1123  {
1124  /* @todo speed up a bit by computing efficacy of new cut from efficacy of old cut
1125  * generate row only if efficient enough
1126  */
1127  SCIP_CALL( generateCutPoint(scip, cons, x, rowprep) );
1128 
1129  if( *rowprep != NULL )
1130  {
1131  efficacy = SCIPgetRowprepViolation(scip, *rowprep, sol);
1132  if( SCIPisGT(scip, efficacy, goodefficacy) ||
1133  (maxnz >= consdata->nvars && SCIPisGT(scip, efficacy, minefficacy)) )
1134  {
1135  /* cut cuts off solution and is efficient enough */
1136  SCIPdebugMsg(scip, "accepted cut with %d of %d nonzeros, efficacy = %g\n", maxnz, consdata->nvars, efficacy);
1137  break;
1138  }
1139 
1140  SCIPfreeRowprep(scip, rowprep);
1141  }
1142 
1143  /* cut also not efficient enough if generated in original refpoint (that's bad) */
1144  if( maxnz >= consdata->nvars )
1145  break;
1146 
1147  nextmaxnz = (int)(conshdlrdata->sparsifynzgrowth * maxnz);
1148  if( nextmaxnz == consdata->nvars - 1)
1149  nextmaxnz = consdata->nvars;
1150  else if( nextmaxnz == maxnz )
1151  ++nextmaxnz;
1152 
1153  /* restore entries of x that are nonzero in next attempt */
1154  for( i = MAX(0, consdata->nvars - nextmaxnz); i < consdata->nvars - maxnz; ++i )
1155  x[ind[i]] = SCIPgetSolVal(scip, sol, consdata->vars[ind[i]]);
1156 
1157  maxnz = nextmaxnz;
1158  }
1159  while( TRUE ); /*lint !e506*/
1160 
1161  SCIPfreeBufferArray(scip, &x);
1162  SCIPfreeBufferArray(scip, &dist);
1163  SCIPfreeBufferArray(scip, &ind);
1164 
1165  return SCIP_OKAY;
1166 }
1167 
1168 /** separates a point, if possible */
1169 static
1171  SCIP* scip, /**< SCIP data structure */
1172  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
1173  SCIP_CONS** conss, /**< constraints */
1174  int nconss, /**< number of constraints */
1175  int nusefulconss, /**< number of constraints that seem to be useful */
1176  SCIP_SOL* sol, /**< solution to separate, or NULL for LP solution */
1177  SCIP_Bool inenforcement, /**< whether we are in constraint enforcement */
1178  SCIP_Bool* cutoff, /**< pointer to store whether a fixing leads to a cutoff */
1179  SCIP_Bool* success /**< buffer to store whether the point was separated */
1180  )
1181 {
1182  SCIP_CONSHDLRDATA* conshdlrdata;
1183  SCIP_CONSDATA* consdata;
1184  SCIP_Real minefficacy;
1185  int c;
1186  SCIP_ROW* row;
1187  SCIP_ROWPREP* rowprep;
1188 
1189  assert(scip != NULL);
1190  assert(conss != NULL || nconss == 0);
1191  assert(nusefulconss <= nconss);
1192  assert(cutoff != NULL);
1193  assert(success != NULL);
1194 
1195  *cutoff = FALSE;
1196 
1197  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1198  assert(conshdlrdata != NULL);
1199 
1200  *success = FALSE;
1201 
1202  minefficacy = inenforcement ? SCIPgetLPFeastol(scip) : SCIPgetSepaMinEfficacy(scip);
1203 
1204  for( c = 0; c < nconss; ++c )
1205  {
1206  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
1207  assert(consdata != NULL);
1208 
1209  if( SCIPisGT(scip, consdata->violation, SCIPfeastol(scip)) && !SCIPisInfinity(scip, consdata->violation) )
1210  {
1211  SCIP_Real efficacy;
1212 
1213  rowprep = NULL;
1214 
1215  /* generate cut */
1216  if( conshdlrdata->sparsify )
1217  {
1218  SCIP_CALL( generateSparseCut(scip, conshdlr, conss[c], sol, &rowprep, minefficacy) ); /*lint !e613*/
1219  }
1220  else if( conshdlrdata->projectpoint )
1221  {
1222  SCIP_CALL( generateCutProjectedPoint(scip, conss[c], sol, &rowprep) ); /*lint !e613*/
1223  }
1224  else
1225  {
1226  SCIP_CALL( generateCutSol(scip, conss[c], sol, &rowprep) ); /*lint !e613*/
1227  }
1228 
1229  if( rowprep == NULL )
1230  continue;
1231 
1232  /* NOTE: The way that rowprep was constructed, there should be no need to call SCIPmergeRowprep,
1233  * since no variable gets added twice. However, if rowprep were replacing multiaggregated variables
1234  * (as there can exist for soc cons), then SCIPmergeRowprep would be necessary.
1235  */
1236  /* cleanup rowprep (there is no limit on coefrange for cons_soc) TODO add a coefrange limit? */
1237  SCIP_CALL( SCIPcleanupRowprep(scip, rowprep, sol, SCIPinfinity(scip), minefficacy, NULL, &efficacy) );
1238 
1239  if( SCIPisLE(scip, efficacy, minefficacy) )
1240  {
1241  SCIPfreeRowprep(scip, &rowprep);
1242  continue;
1243  }
1244 
1245  /* cut cuts off solution and efficient enough */
1246  SCIP_CALL( SCIPgetRowprepRowConshdlr(scip, &row, rowprep, conshdlr) );
1247  if( SCIPisCutApplicable(scip, row) )
1248  {
1249  SCIP_CALL( SCIPaddRow(scip, row, FALSE, cutoff) );
1250  SCIP_CALL( SCIPresetConsAge(scip, conss[c]) ); /*lint !e613*/
1251 
1252  *success = TRUE;
1253 
1254  SCIPdebugMsg(scip, "added cut with efficacy %g\n", SCIPgetCutEfficacy(scip, sol, row));
1255 
1256  /* mark row as not removable from LP for current node, if in enforcement */
1257  if( inenforcement && !conshdlrdata->enfocutsremovable )
1258  SCIPmarkRowNotRemovableLocal(scip, row);
1259  }
1260 
1261  SCIP_CALL( SCIPreleaseRow (scip, &row) );
1262  SCIPfreeRowprep(scip, &rowprep);
1263  }
1264 
1265  if( *cutoff )
1266  break;
1267 
1268  /* enforce only useful constraints
1269  * others are only checked and enforced if we are still feasible or have not found a separating cut yet
1270  */
1271  if( c >= nusefulconss && *success )
1272  break;
1273  }
1274 
1275  return SCIP_OKAY;
1276 }
1277 
1278 /** adds linearizations cuts for convex constraints w.r.t. a given reference point to cutpool and sepastore
1279  *
1280  * 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.
1281  * If separatedlpsol is not NULL, but cut does not separate the LP solution, then it is added to the cutpool only.
1282  * If separatedlpsol is NULL, then cut is added to cutpool only.
1283  */
1284 static
1286  SCIP* scip, /**< SCIP data structure */
1287  SCIP_CONSHDLR* conshdlr, /**< quadratic constraints handler */
1288  SCIP_CONS** conss, /**< constraints */
1289  int nconss, /**< number of constraints */
1290  SCIP_SOL* ref, /**< reference point where to linearize, or NULL for LP solution */
1291  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 */
1292  SCIP_Real minefficacy, /**< minimal efficacy of a cut when checking for separation of LP solution */
1293  SCIP_Bool* cutoff /**< pointer to store whether a cutoff was detected */
1294  )
1295 {
1296  SCIP_CONSDATA* consdata;
1297  SCIP_Bool addedtolp;
1298  SCIP_ROWPREP* rowprep;
1299  int c;
1300 
1301  assert(scip != NULL);
1302  assert(conshdlr != NULL);
1303  assert(conss != NULL || nconss == 0);
1304  assert(cutoff != NULL);
1305  *cutoff = FALSE;
1306 
1307  if( separatedlpsol != NULL )
1308  *separatedlpsol = FALSE;
1309 
1310  for( c = 0; c < nconss && !(*cutoff); ++c )
1311  {
1312  assert(conss[c] != NULL); /*lint !e613 */
1313 
1314  if( SCIPconsIsLocal(conss[c]) ) /*lint !e613 */
1315  continue;
1316 
1317  consdata = SCIPconsGetData(conss[c]); /*lint !e613 */
1318  assert(consdata != NULL);
1319 
1320  SCIP_CALL( evalLhs(scip, conss[c], ref) ); /*lint !e613 */
1321  if( !SCIPisPositive(scip, consdata->lhsval) || SCIPisInfinity(scip, consdata->lhsval) )
1322  {
1323  SCIPdebugMsg(scip, "skip adding linearization for <%s> since lhs is %g\n", SCIPconsGetName(conss[c]), consdata->lhsval); /*lint !e613 */
1324  continue;
1325  }
1326 
1327  SCIP_CALL( generateCutSol(scip, conss[c], ref, &rowprep) ); /*lint !e613 */
1328 
1329  if( rowprep == NULL )
1330  continue;
1331 
1332  addedtolp = FALSE;
1333 
1334  /* if caller wants, then check if cut separates LP solution and add to sepastore if so */
1335  if( separatedlpsol != NULL )
1336  {
1337  if( SCIPgetRowprepViolation(scip, rowprep, NULL) >= minefficacy )
1338  {
1339  SCIP_ROW* row;
1340 
1341  SCIP_CALL( SCIPgetRowprepRowConshdlr(scip, &row, rowprep, conshdlr) );
1342  SCIP_CALL( SCIPaddRow(scip, row, TRUE, cutoff) );
1343  SCIP_CALL( SCIPreleaseRow(scip, &row) );
1344 
1345  *separatedlpsol = TRUE;
1346  addedtolp = TRUE;
1347  }
1348  }
1349 
1350  if( !addedtolp && !rowprep->local )
1351  {
1352  SCIP_ROW* row;
1353 
1354  SCIP_CALL( SCIPgetRowprepRowConshdlr(scip, &row, rowprep, conshdlr) );
1355  SCIP_CALL( SCIPaddPoolCut(scip, row) );
1356  SCIP_CALL( SCIPreleaseRow(scip, &row) );
1357  }
1358 
1359  SCIPfreeRowprep(scip, &rowprep);
1360  }
1361 
1362  return SCIP_OKAY;
1363 }
1364 
1365 /** processes the event that a new primal solution has been found */
1366 static
1367 SCIP_DECL_EVENTEXEC(processNewSolutionEvent)
1369  SCIP_CONSHDLRDATA* conshdlrdata;
1370  SCIP_CONSHDLR* conshdlr;
1371  SCIP_CONS** conss;
1372  int nconss;
1373  SCIP_SOL* sol;
1374  SCIP_Bool cutoff;
1375 
1376  assert(scip != NULL);
1377  assert(event != NULL);
1378  assert(eventdata != NULL);
1379  assert(eventhdlr != NULL);
1380 
1381  assert((SCIPeventGetType(event) & SCIP_EVENTTYPE_SOLFOUND) != 0);
1382 
1383  conshdlr = (SCIP_CONSHDLR*)eventdata;
1384 
1385  nconss = SCIPconshdlrGetNConss(conshdlr);
1386 
1387  if( nconss == 0 )
1388  return SCIP_OKAY;
1389 
1390  sol = SCIPeventGetSol(event);
1391  assert(sol != NULL);
1392 
1393  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1394  assert(conshdlrdata != NULL);
1395 
1396  /* we are only interested in solution coming from some heuristic other than trysol, but not from the tree
1397  * the reason for ignoring trysol solutions is that they may come from an NLP solve in sepalp, where we already added linearizations,
1398  * or are from the tree, but postprocessed via proposeFeasibleSolution
1399  */
1400  if( SCIPsolGetHeur(sol) == NULL || SCIPsolGetHeur(sol) == conshdlrdata->trysolheur )
1401  return SCIP_OKAY;
1402 
1403  conss = SCIPconshdlrGetConss(conshdlr);
1404  assert(conss != NULL);
1405 
1406  SCIPdebugMsg(scip, "caught new sol event %" SCIP_EVENTTYPE_FORMAT " from heur <%s>; have %d conss\n", SCIPeventGetType(event),
1407  SCIPheurGetName(SCIPsolGetHeur(sol)), nconss);
1408 
1409  SCIP_CALL( addLinearizationCuts(scip, conshdlr, conss, nconss, sol, NULL, 0.0, &cutoff) );
1410  /* ignore cutoff, cannot return status */
1411 
1412  return SCIP_OKAY;
1413 }
1414 
1415 /** removes fixed variables, replace aggregated and negated variables
1416  *
1417  * repeats replacements until no further change is found;
1418  * takes care of capture/release and locks, but not of variable events (assumes that var events are not caught yet)
1419  */
1420 static
1422  SCIP* scip, /**< SCIP data structure */
1423  SCIP_CONSHDLR* conshdlr, /**< constraint handler for signpower constraints */
1424  SCIP_CONS* cons, /**< constraint */
1425  int* ndelconss, /**< counter for number of deleted constraints */
1426  int* nupgdconss, /**< counter for number of upgraded constraints */
1427  int* nchgbds, /**< counter for number of bound changes */
1428  int* nfixedvars, /**< counter for number of fixed variables */
1429  SCIP_Bool* iscutoff, /**< to store whether constraint cannot be satisfied */
1430  SCIP_Bool* isdeleted /**< to store whether constraint is redundant and can be deleted */
1431  )
1432 {
1433  SCIP_CONSDATA* consdata;
1434  SCIP_CONSHDLRDATA* conshdlrdata;
1435  SCIP_Bool havechange;
1436  SCIP_Bool haveremovedvar;
1437  int i;
1438  SCIP_VAR* x;
1439  SCIP_Real coef;
1440  SCIP_Real offset;
1441 
1442  assert(scip != NULL);
1443  assert(conshdlr != NULL);
1444  assert(cons != NULL);
1445  assert(iscutoff != NULL);
1446  assert(isdeleted != NULL);
1447 
1448  *iscutoff = FALSE;
1449  *isdeleted = FALSE;
1450 
1451  consdata = SCIPconsGetData(cons);
1452  assert(consdata != NULL);
1453 
1454  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1455  assert(conshdlrdata != NULL);
1456 
1457  SCIPdebugMsg(scip, "remove fixed variables from constraint <%s>\n", SCIPconsGetName(cons));
1458  SCIPdebugPrintCons(scip, cons, NULL);
1459 
1460  havechange = FALSE;
1461  haveremovedvar = FALSE;
1462 
1463  /* process variables on left hand side */
1464  for( i = 0; i < consdata->nvars; ++i )
1465  {
1466  x = consdata->vars[i];
1467  assert(x != NULL);
1469 
1471  continue;
1472 
1473  havechange = TRUE;
1474 
1475  /* drop variable event and unlock and release variable */
1476  SCIP_CALL( dropLhsVarEvents(scip, conshdlrdata->eventhdlr, cons, i) );
1477  SCIP_CALL( SCIPunlockVarCons(scip, x, cons, TRUE, TRUE) );
1478  SCIP_CALL( SCIPreleaseVar(scip, &consdata->vars[i]) );
1479 
1480  coef = 1.0;
1481  offset = consdata->offsets[i];
1482  SCIP_CALL( SCIPgetProbvarSum(scip, &x, &coef, &offset) );
1483 
1484  SCIPdebugMsg(scip, " lhs term at position %d is replaced by %g * <%s> + %g\n",
1485  i, coef, SCIPvarGetName(x), offset);
1486 
1487  /* if variable has been fixed, add (alpha*offset)^2 to gamma and continue */
1488  if( coef == 0.0 || x == NULL )
1489  {
1490  consdata->constant += consdata->coefs[i] * consdata->coefs[i] * offset * offset;
1491  consdata->offsets[i] = 0.0;
1492  haveremovedvar = TRUE;
1493  continue;
1494  }
1495 
1497 
1498  /* replace coefs[i] * (vars[i] + offsets[i]) by coefs[i]*coef * (x + offsets[i]/coef) */
1499  consdata->offsets[i] = offset;
1500  if( coef != 1.0 )
1501  {
1502  consdata->coefs[i] = REALABS(coef * consdata->coefs[i]);
1503  consdata->offsets[i] /= coef;
1504  }
1505  consdata->vars[i] = x;
1506 
1507  /* capture and lock new variable, catch variable events */
1508  SCIP_CALL( SCIPcaptureVar(scip, consdata->vars[i]) );
1509  SCIP_CALL( SCIPlockVarCons(scip, consdata->vars[i], cons, TRUE, TRUE) );
1510  SCIP_CALL( catchLhsVarEvents(scip, conshdlrdata->eventhdlr, cons, i) );
1511  }
1512 
1513  /* process variable on right hand side */
1514  x = consdata->rhsvar;
1515  assert(x != NULL);
1517  {
1518  havechange = TRUE;
1519 
1520  /* drop variable event and unlock and release variable */
1521  SCIP_CALL( dropRhsVarEvents(scip, conshdlrdata->eventhdlr, cons) );
1522  SCIP_CALL( SCIPreleaseVar(scip, &consdata->rhsvar) );
1523  SCIP_CALL( SCIPunlockVarCons(scip, x, cons, consdata->rhscoeff > 0.0, consdata->rhscoeff < 0.0) );
1524 
1525  coef = 1.0;
1526  offset = 0.0;
1527  SCIP_CALL( SCIPgetProbvarSum(scip, &x, &coef, &offset) );
1528 
1529  SCIPdebugMsg(scip, " rhs variable is replaced by %g * <%s> + %g\n", coef, SCIPvarGetName(x), offset);
1530 
1531  if( coef == 0.0 || x == NULL )
1532  {
1533  /* if variable has been fixed, add offset to rhsoffset */
1534  consdata->rhsoffset += offset;
1535  }
1536  else
1537  {
1538  /* replace rhscoef * (rhsvar + rhsoffset) by rhscoef*coef * (x + offset/coef + rhsoffset/coef) */
1540 
1541  consdata->rhsoffset = (consdata->rhsoffset + offset) / coef;
1542  consdata->rhscoeff *= coef;
1543  consdata->rhsvar = x;
1544 
1545  /* capture and lock new variable, catch variable events */
1546  SCIP_CALL( SCIPcaptureVar(scip, consdata->rhsvar) );
1547  SCIP_CALL( SCIPlockVarCons(scip, consdata->rhsvar, cons, consdata->rhscoeff > 0.0, consdata->rhscoeff < 0.0) );
1548  SCIP_CALL( catchRhsVarEvents(scip, conshdlrdata->eventhdlr, cons) );
1549  }
1550  }
1551 
1552  if( !havechange )
1553  return SCIP_OKAY;
1554 
1555  /* free nonlinear row representation */
1556  if( consdata->nlrow != NULL )
1557  {
1558  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
1559  }
1560 
1561  /* if a variable has been removed, close gaps in vars array */
1562  if( haveremovedvar )
1563  {
1564  int oldnvars;
1565 
1566  /* 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 */
1567  SCIP_CALL( dropVarEvents(scip, conshdlrdata->eventhdlr, cons) );
1568 
1569  oldnvars = consdata->nvars;
1570  for( i = 0; i < consdata->nvars; ++i )
1571  {
1572  /* forget about empty places at end of vars array */
1573  while( consdata->nvars && consdata->vars[consdata->nvars-1] == NULL )
1574  --consdata->nvars;
1575 
1576  /* all variables at index >= i have been removed */
1577  if( i == consdata->nvars )
1578  break;
1579 
1580  if( consdata->vars[i] != NULL )
1581  continue;
1582 
1583  /* move variable from position nvars-1 to position i */
1584 
1585  assert(consdata->nvars >= 1);
1586  assert(consdata->vars[consdata->nvars-1] != NULL);
1587 
1588  consdata->vars[i] = consdata->vars[consdata->nvars-1];
1589  consdata->offsets[i] = consdata->offsets[consdata->nvars-1];
1590  consdata->coefs[i] = consdata->coefs[consdata->nvars-1];
1591 
1592  --consdata->nvars;
1593  }
1594 
1595  assert(consdata->nvars < oldnvars);
1596 
1597  /* shrink arrays in consdata */
1598  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->vars, oldnvars, consdata->nvars) );
1599  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->offsets, oldnvars, consdata->nvars) );
1600  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->coefs, oldnvars, consdata->nvars) );
1601 
1602  SCIP_CALL( catchVarEvents(scip, conshdlrdata->eventhdlr, cons) );
1603  }
1604 
1605  SCIPdebugMsg(scip, "\t-> ");
1606  SCIPdebugPrintCons(scip, cons, NULL);
1607 
1608  if( consdata->nvars == 0 )
1609  { /* all variables on left hand size have been removed, remaining constraint is sqrt(gamma) <= ... */
1610  assert(!SCIPisNegative(scip, consdata->constant));
1611  if( consdata->rhsvar == NULL )
1612  { /* also rhsvar has been removed, remaining constraint is sqrt(gamma) <= rhscoeff * rhsoffset */
1613  if( SCIPisFeasLE(scip, sqrt(consdata->constant), consdata->rhscoeff*consdata->rhsoffset) )
1614  {
1615  SCIPdebugMsg(scip, "remove redundant constraint <%s> after fixing all variables\n", SCIPconsGetName(cons));
1616  }
1617  else
1618  {
1619  SCIPdebugMsg(scip, "found problem infeasible after fixing all variables in <%s>\n", SCIPconsGetName(cons));
1620  *iscutoff = TRUE;
1621  }
1622  ++*ndelconss;
1623  }
1624  else if( !SCIPvarIsActive(consdata->rhsvar) )
1625  { /* remaining constraint is sqrt(gamma) - rhscoeff * rhsoffset <= rhscoeff * rhsvar, and rhsvar is probably multi-aggregated */
1626  SCIP_CONS* lincons;
1627 
1628  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, SCIPconsGetName(cons), 1, &consdata->rhsvar, &consdata->rhscoeff,
1629  sqrt(consdata->constant) - consdata->rhscoeff * consdata->rhsoffset, SCIPinfinity(scip),
1633  SCIPconsIsStickingAtNode(cons)) );
1634  SCIP_CALL( SCIPaddCons(scip, lincons) );
1635  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
1636  ++*nupgdconss;
1637  }
1638  else if( consdata->rhscoeff > 0.0 )
1639  { /* remaining constraint is sqrt(gamma) / rhscoeff - rhsoffset <= rhsvar */
1640  SCIP_Bool tightened;
1641  SCIP_CALL( SCIPtightenVarLb(scip, consdata->rhsvar, sqrt(consdata->constant) / consdata->rhscoeff - consdata->rhsoffset, TRUE, iscutoff, &tightened) );
1642  if( *iscutoff )
1643  {
1644  SCIPdebugMsg(scip, "found problem infeasible after fixing all lhs variables in <%s> and tightening lower bound of rhs var\n", SCIPconsGetName(cons));
1645  }
1646  else if( tightened )
1647  {
1648  SCIPdebugMsg(scip, "remove redundant constraint <%s> after fixing all lhs variables and tightening lower bound of rhs var\n", SCIPconsGetName(cons));
1649  ++*nchgbds;
1650  }
1651  else
1652  {
1653  SCIPdebugMsg(scip, "remove redundant constraint <%s> after fixing all lhs variables\n", SCIPconsGetName(cons));
1654  }
1655  ++*ndelconss;
1656  }
1657  else
1658  { /* remaining constraint is sqrt(gamma) / rhscoeff - rhsoffset >= rhsvar */
1659  SCIP_Bool tightened;
1660  SCIP_CALL( SCIPtightenVarUb(scip, consdata->rhsvar, sqrt(consdata->constant) / consdata->rhscoeff - consdata->rhsoffset, TRUE, iscutoff, &tightened) );
1661  if( *iscutoff )
1662  {
1663  SCIPdebugMsg(scip, "found problem infeasible after fixing all lhs variables in <%s> and tightening upper bound of rhs var\n", SCIPconsGetName(cons));
1664  }
1665  else if( tightened )
1666  {
1667  SCIPdebugMsg(scip, "remove redundant constraint <%s> after fixing all lhs variables and tightening upper bound of rhs var\n", SCIPconsGetName(cons));
1668  ++*nchgbds;
1669  }
1670  else
1671  {
1672  SCIPdebugMsg(scip, "remove redundant constraint <%s> after fixing all lhs variables\n", SCIPconsGetName(cons));
1673  }
1674  ++*ndelconss;
1675  }
1676  SCIP_CALL( SCIPdelCons(scip, cons) );
1677  *isdeleted = TRUE;
1678  return SCIP_OKAY;
1679  }
1680 
1681  if( consdata->rhsvar == NULL )
1682  { /* constraint becomes sum_i (alpha_i*(x_i+beta_i))^2 <= (rhscoeff*rhsoffset)^2 - gamma */
1683  if( consdata->nvars > 1 )
1684  { /* upgrade to quadratic constraint */
1685  SCIP_CONS* quadcons;
1686  SCIP_QUADVARTERM* quadvarterms;
1687  SCIP_Real rhs;
1688 
1689  SCIP_CALL( SCIPallocBufferArray(scip, &quadvarterms, consdata->nvars) );
1690  BMSclearMemoryArray(quadvarterms, consdata->nvars);
1691  rhs = consdata->rhscoeff * consdata->rhsoffset;
1692  rhs = rhs * rhs - consdata->constant;
1693 
1694  for( i = 0; i < consdata->nvars; ++i )
1695  {
1696  quadvarterms[i].var = consdata->vars[i];
1697  quadvarterms[i].sqrcoef = consdata->coefs[i] * consdata->coefs[i];
1698  if( consdata->offsets[i] != 0.0 )
1699  {
1700  quadvarterms[i].lincoef = 2 * consdata->offsets[i] * quadvarterms[i].sqrcoef;
1701  rhs -= quadvarterms[i].sqrcoef * consdata->offsets[i]*consdata->offsets[i];
1702  }
1703  }
1704 
1705  assert(!SCIPconsIsStickingAtNode(cons));
1706  SCIP_CALL( SCIPcreateConsQuadratic2(scip, &quadcons, SCIPconsGetName(cons), 0, NULL, NULL,
1707  consdata->nvars, quadvarterms, 0, NULL, -SCIPinfinity(scip), rhs,
1711  SCIP_CALL( SCIPaddCons(scip, quadcons) );
1712  SCIPdebugMsg(scip, "upgraded <%s> to quadratic constraint: ", SCIPconsGetName(cons));
1713  SCIPdebugPrintCons(scip, quadcons, NULL);
1714 
1715  SCIP_CALL( SCIPreleaseCons(scip, &quadcons) );
1716 
1717  SCIPfreeBufferArray(scip, &quadvarterms);
1718 
1719  ++*nupgdconss;
1720  }
1721  else if( !SCIPvarIsActive(consdata->vars[0]) )
1722  { /* constraint is |alpha*(x+beta)| <= sqrt((rhscoeff*rhsoffset)^2 - gamma), but x is probably multaggr. -> turn into ranged linear constraint */
1723  SCIP_CONS* lincons;
1724 
1725  /* create constraint alpha*x <= sqrt((rhscoeff*rhsoffset)^2 - gamma) - alpha*beta
1726  * alpha*x >= -sqrt((rhscoeff*rhsoffset)^2 - gamma) - alpha*beta */
1727  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, SCIPconsGetName(cons), 1, &consdata->vars[0], &consdata->coefs[0],
1728  -sqrt(consdata->rhscoeff * consdata->rhscoeff * consdata->rhsoffset * consdata->rhsoffset - consdata->constant) - consdata->coefs[0] * consdata->offsets[0],
1729  +sqrt(consdata->rhscoeff * consdata->rhscoeff * consdata->rhsoffset * consdata->rhsoffset - consdata->constant) - consdata->coefs[0] * consdata->offsets[0],
1733  SCIPconsIsStickingAtNode(cons)) );
1734  SCIP_CALL( SCIPaddCons(scip, lincons) );
1735  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
1736 
1737  ++*nupgdconss;
1738  }
1739  else
1740  { /* constraint is |alpha*(x+beta)| <= sqrt((rhscoeff*rhsoffset)^2 - gamma) -> propagate bounds */
1741  SCIP_Bool tightened;
1742  SCIP_Real rhs;
1743  assert(consdata->nvars == 1); /* case == 0 handled before */
1744  rhs = consdata->rhscoeff * consdata->rhsoffset;
1745  rhs = rhs * rhs;
1746  if( SCIPisNegative(scip, rhs - consdata->constant) )
1747  { /* take this as infeasible */
1748  SCIPdebugMsg(scip, "found problem infeasible after fixing rhs and all except one lhs variables in <%s>\n", SCIPconsGetName(cons));
1749  *iscutoff = TRUE;
1750  }
1751  else
1752  {
1753  rhs -= consdata->constant;
1754  rhs = rhs < 0.0 ? 0.0 : sqrt(rhs);
1755 
1756  if( SCIPisZero(scip, rhs) )
1757  { /* constraint is x = -beta */
1758  SCIP_CALL( SCIPfixVar(scip, consdata->vars[0], -consdata->offsets[0], iscutoff, &tightened) );
1759  if( *iscutoff )
1760  {
1761  SCIPdebugMsg(scip, "found problem infeasible after fixing rhs and all except one lhs variables and fixing remaining lhs var in <%s>\n", SCIPconsGetName(cons));
1762  }
1763  else if( tightened )
1764  {
1765  SCIPdebugMsg(scip, "remove redundant constraint <%s> after fixing rhs and all except one lhs variables and fixing remaining lhs var\n", SCIPconsGetName(cons));
1766  ++*nfixedvars;
1767  }
1768  else
1769  {
1770  SCIPdebugMsg(scip, "remove redundant constraint <%s> after fixing rhs and all except one lhs variables and fixing remaining lhs var\n", SCIPconsGetName(cons));
1771  }
1772  }
1773  else
1774  { /* constraint is -rhs/|alpha| - beta <= x <= rhs/|alpha| - beta */
1775  rhs /= ABS(consdata->coefs[0]);
1776  SCIP_CALL( SCIPtightenVarLb(scip, consdata->vars[0], -rhs - consdata->offsets[0], TRUE, iscutoff, &tightened) );
1777  if( *iscutoff )
1778  {
1779  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));
1780  }
1781  else
1782  {
1783  if( tightened )
1784  ++*nchgbds;
1785  SCIP_CALL( SCIPtightenVarUb(scip, consdata->vars[0], rhs - consdata->offsets[0], TRUE, iscutoff, &tightened) );
1786  if( *iscutoff )
1787  {
1788  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));
1789  }
1790  else if( tightened )
1791  ++*nchgbds;
1792  }
1793  if( !*iscutoff )
1794  {
1795  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));
1796  }
1797  }
1798  }
1799  ++*ndelconss;
1800  }
1801  *isdeleted = TRUE;
1802  SCIP_CALL( SCIPdelCons(scip, cons) );
1803  return SCIP_OKAY;
1804  }
1805 
1806  if( consdata->nvars == 1 && SCIPisZero(scip, consdata->constant) )
1807  { /* one variable on lhs left and no constant, constraint becomes |alpha*(x+beta)| <= rhscoef*(rhsvar+rhsoffset) -> upgrade to two linear constraints */
1808  SCIP_CONS* lincons;
1809  SCIP_VAR* vars[2];
1810  SCIP_Real coefs[2];
1811  SCIP_Real rhs;
1812  assert(consdata->rhsvar != NULL); /* case == NULL has been handled before */
1813 
1814  vars[0] = consdata->vars[0];
1815  vars[1] = consdata->rhsvar;
1816  coefs[0] = consdata->coefs[0];
1817  coefs[1] = -consdata->rhscoeff;
1818  rhs = consdata->rhscoeff * consdata->rhsoffset - coefs[0] * consdata->offsets[0];
1819 
1820  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, SCIPconsGetName(cons), 2, vars, coefs, -SCIPinfinity(scip), rhs,
1824  SCIPconsIsStickingAtNode(cons)) );
1825  SCIP_CALL( SCIPaddCons(scip, lincons) );
1826  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
1827 
1828  coefs[0] = -coefs[0];
1829  rhs = consdata->rhscoeff * consdata->rhsoffset - coefs[0] * consdata->offsets[0];
1830 
1831  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, SCIPconsGetName(cons), 2, vars, coefs, -SCIPinfinity(scip), rhs,
1835  SCIPconsIsStickingAtNode(cons)) );
1836  SCIP_CALL( SCIPaddCons(scip, lincons) );
1837  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
1838 
1839  SCIPdebugMsg(scip, "upgraded <%s> to two linear constraint\n", SCIPconsGetName(cons));
1840 
1841  ++*nupgdconss;
1842  SCIP_CALL( SCIPdelCons(scip, cons) );
1843  *isdeleted = TRUE;
1844  return SCIP_OKAY;
1845  }
1846 
1847  return SCIP_OKAY;
1848 }
1849 
1850 
1851 /** adds the linear outer-approximation of Glineur et.al. for a SOC constraint of dimension 3
1852  *
1853  * 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$.
1854  * Here \f$\alpha3 > 0\f$, and the lower bound of \f$x_3 \geq -offset3\f$.
1855  * Also x2 = NULL is allowed, in which case the second term is assumed to be constant, and \f$offset2 \neq 0\f$ is needed.
1856  */
1857 static
1859  SCIP* scip, /**< SCIP data structure */
1860  SCIP_CONS* cons, /**< original constraint */
1861  SCIP_VAR* x1, /**< variable x1 */
1862  SCIP_VAR* x2, /**< variable x2 */
1863  SCIP_VAR* x3, /**< variable x3 */
1864  SCIP_Real alpha1, /**< coefficient of x1 */
1865  SCIP_Real alpha2, /**< coefficient of x2 */
1866  SCIP_Real alpha3, /**< coefficient of x3 */
1867  SCIP_Real offset1, /**< offset of x1 */
1868  SCIP_Real offset2, /**< offset of x2 */
1869  SCIP_Real offset3, /**< offset of x3 */
1870  int N, /**< size of linear approximation, need to be >= 1 */
1871  const char* basename, /**< string to use for building variable and constraint names */
1872  int* naddconss /**< buffer where to add the number of added constraints */
1873  )
1874 {
1875  SCIP_CONS* lincons;
1876  SCIP_VAR* vars[3];
1877  SCIP_Real vals[3];
1878  char varname[255];
1879  char linname[255];
1880  int i;
1881  SCIP_VAR** avars;
1882  SCIP_VAR** bvars;
1883  SCIP_Real val;
1884 
1885  assert(scip != NULL);
1886  assert(cons != NULL);
1887  assert(x1 != NULL);
1888  assert(x2 != NULL || !SCIPisZero(scip, offset2));
1889  assert(x3 != NULL);
1890  assert(SCIPisPositive(scip, alpha3));
1891  assert(SCIPisGE(scip, SCIPconsIsLocal(cons) ? SCIPvarGetLbLocal(x3) : SCIPvarGetLbGlobal(x3), -offset3));
1892  assert(basename != NULL);
1893  assert(N >= 1);
1894  assert(naddconss != NULL);
1895 
1896  SCIPdebugMsg(scip, "Creating linear Glineur outer-approximation for <%s>.\n", basename);
1897  SCIPdebugMsg(scip, "sqr(%g(%s+%g)) + sqr(%g(%s+%g)) <= sqr(%g(%s+%g)).\n",
1898  alpha1, SCIPvarGetName(x1), offset1, alpha2, x2 ? SCIPvarGetName(x2) : "0", offset2, alpha3, SCIPvarGetName(x3), offset3);
1899 
1900  SCIP_CALL( SCIPallocBufferArray(scip, &avars, N+1) );
1901  SCIP_CALL( SCIPallocBufferArray(scip, &bvars, N+1) );
1902 
1903  /* create additional variables; we do not use avars[0] and bvars[0] */
1904  for( i = 1; i <= N; ++i )
1905  {
1906  (void) SCIPsnprintf(varname, 255, "soc#%s_a%d", basename, i);
1907  SCIP_CALL( SCIPcreateVar(scip, &avars[i], varname, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0,
1909  SCIP_CALL( SCIPaddVar(scip, avars[i]) );
1910 
1911  (void) SCIPsnprintf(varname, 255, "soc#%s_b%d", basename, i);
1912  SCIP_CALL( SCIPcreateVar(scip, &bvars[i], varname, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0,
1914  SCIP_CALL( SCIPaddVar(scip, bvars[i]) );
1915  }
1916 
1917  /* create linear constraints for the first case
1918  * cos(pi) = -1, sin(pi) = 0
1919  * -> a_1 = - alpha1 (x1 + offset1) -> -alpha1*x1 - a_1 = alpha1*offset1
1920  * -> b_1 >= | alpha2 (x2 + offset2) | -> alpha2*x2 - b_1 <= -alpha2*offset2
1921  * alpha2*x2 + b_1 >= -alpha2*offset2
1922  */
1923 
1924  vars[0] = x1;
1925  vals[0] = -alpha1;
1926  vars[1] = avars[1];
1927  vals[1] = -1.0;
1928 
1929  (void) SCIPsnprintf(linname, 255, "soc#%s#a%d", basename, 0);
1930  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, alpha1*offset1, alpha1*offset1,
1936  SCIP_CALL( SCIPaddCons(scip, lincons) );
1937  SCIPdebugPrintCons(scip, lincons, NULL);
1938  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
1939  ++(*naddconss);
1940 
1941  if( x2 != NULL )
1942  {
1943  vars[0] = x2;
1944  vals[0] = alpha2;
1945  vars[1] = bvars[1];
1946  vals[1] = -1.0;
1947 
1948  (void) SCIPsnprintf(linname, 255, "soc#%s#b%d", basename, 0);
1949  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, -SCIPinfinity(scip), -alpha2*offset2,
1955  SCIP_CALL( SCIPaddCons(scip, lincons) );
1956  SCIPdebugPrintCons(scip, lincons, NULL);
1957  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
1958  ++(*naddconss);
1959 
1960  vars[0] = x2;
1961  vals[0] = alpha2;
1962  vars[1] = bvars[1];
1963  vals[1] = 1.0;
1964 
1965  (void) SCIPsnprintf(linname, 255, "soc#%s#B%d", basename, 0);
1966  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, -alpha2*offset2, SCIPinfinity(scip),
1972  SCIP_CALL( SCIPaddCons(scip, lincons) );
1973  SCIPdebugPrintCons(scip, lincons, NULL);
1974  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
1975  ++(*naddconss);
1976  }
1977  else
1978  { /* x2 == NULL -> b_1 >= |alpha2*offset2| */
1979  SCIP_Bool infeas;
1980  SCIP_Bool tightened;
1981  SCIP_CALL( SCIPtightenVarLb(scip, bvars[1], ABS(alpha2 * offset2), TRUE, &infeas, &tightened) );
1982  if( infeas == TRUE )
1983  {
1984  SCIPwarningMessage(scip, "creating glineur outer approximation of SOC3 constraint found problem infeasible.\n");
1985  }
1986  }
1987 
1988  /* create intermediate linear constraints */
1989  val = M_PI;
1990  for( i = 1; i < N; ++i )
1991  {
1992  val /= 2.0;
1993 
1994  vars[0] = avars[i];
1995  vals[0] = cos(val);
1996  vars[1] = bvars[i];
1997  vals[1] = sin(val);
1998  vars[2] = avars[i+1]; /*lint !e679*/
1999  vals[2] = -1.0;
2000 
2001  (void) SCIPsnprintf(linname, 255, "soc#%s#a%d", basename, i);
2002  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 3, vars, vals, 0.0, 0.0,
2008  SCIP_CALL( SCIPaddCons(scip, lincons) );
2009  SCIPdebugPrintCons(scip, lincons, NULL);
2010  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2011  ++(*naddconss);
2012 
2013  vars[0] = avars[i];
2014  vals[0] = -sin(val);
2015  vars[1] = bvars[i];
2016  vals[1] = cos(val);
2017  vars[2] = bvars[i+1]; /*lint !e679*/
2018  vals[2] = -1.0;
2019 
2020  (void) SCIPsnprintf(linname, 255, "soc#%s#b%d", basename, i);
2021  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 3, vars, vals, -SCIPinfinity(scip), 0.0,
2027  SCIP_CALL( SCIPaddCons(scip, lincons) );
2028  SCIPdebugPrintCons(scip, lincons, NULL);
2029  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2030  ++(*naddconss);
2031 
2032  vars[0] = avars[i];
2033  vals[0] = -sin(val);
2034  vars[1] = bvars[i];
2035  vals[1] = cos(val);
2036  vars[2] = bvars[i+1]; /*lint !e679*/
2037  vals[2] = 1.0;
2038 
2039  (void) SCIPsnprintf(linname, 255, "soc#%s#B%d", basename, i);
2040  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 3, vars, vals, 0.0, SCIPinfinity(scip),
2046  SCIP_CALL( SCIPaddCons(scip, lincons) );
2047  SCIPdebugPrintCons(scip, lincons, NULL);
2048  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2049  ++(*naddconss);
2050  }
2051 
2052  /* create last linear constraint */
2053  val /= 2.0;
2054  vars[0] = avars[N];
2055  vals[0] = -cos(val);
2056  vars[1] = bvars[N];
2057  vals[1] = -sin(val);
2058  vars[2] = x3;
2059  vals[2] = alpha3;
2060 
2061  (void) SCIPsnprintf(linname, 255, "soc#%s#a%d", basename, N);
2062  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 3, vars, vals, -alpha3*offset3, -alpha3*offset3,
2068  SCIP_CALL( SCIPaddCons(scip, lincons) );
2069  SCIPdebugPrintCons(scip, lincons, NULL);
2070  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2071  ++(*naddconss);
2072 
2073  for( i = 1; i <= N; ++i )
2074  {
2075  SCIP_CALL( SCIPreleaseVar(scip, &avars[i]) );
2076  SCIP_CALL( SCIPreleaseVar(scip, &bvars[i]) );
2077  }
2078  SCIPfreeBufferArray(scip, &avars);
2079  SCIPfreeBufferArray(scip, &bvars);
2080 
2081  return SCIP_OKAY;
2082 }
2083 
2084 /** adds the linear outer-approximation of Ben-Tal and Nemirovski for a SOC constraint of dimension 3
2085  *
2086  * 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$.
2087  * Here \f$\alpha3 > 0.0\f$, and the lower bound of \f$x_3 \geq -offset3\f$.
2088  * Also x2 = NULL is allowed, in which case the second term is assumed to be constant, and \f$offset2 \neq 0\f$ is needed.
2089  * */
2090 static
2092  SCIP* scip, /**< SCIP data structure */
2093  SCIP_CONS* cons, /**< original constraint */
2094  SCIP_VAR* x1, /**< variable x1 */
2095  SCIP_VAR* x2, /**< variable x2 */
2096  SCIP_VAR* x3, /**< variable x3 */
2097  SCIP_Real alpha1, /**< coefficient of x1 */
2098  SCIP_Real alpha2, /**< coefficient of x2 */
2099  SCIP_Real alpha3, /**< coefficient of x3 */
2100  SCIP_Real offset1, /**< offset of x1 */
2101  SCIP_Real offset2, /**< offset of x2 */
2102  SCIP_Real offset3, /**< offset of x3 */
2103  int N, /**< size of linear approximation, need to be >= 1 */
2104  const char* basename, /**< string to use for building variable and constraint names */
2105  int* naddconss /**< buffer where to add the number of added constraints */
2106  )
2107 {
2108  SCIP_CONS* lincons;
2109  SCIP_VAR* vars[3];
2110  SCIP_Real vals[3];
2111  char varname[255];
2112  char linname[255];
2113  int i;
2114  SCIP_VAR** avars;
2115  SCIP_VAR** bvars;
2116 
2117  assert(scip != NULL);
2118  assert(cons != NULL);
2119  assert(x1 != NULL);
2120  assert(x2 != NULL || !SCIPisZero(scip, offset2));
2121  assert(x3 != NULL);
2122  assert(SCIPisPositive(scip, alpha3));
2123  assert(SCIPisGE(scip, SCIPconsIsLocal(cons) ? SCIPvarGetLbLocal(x3) : SCIPvarGetLbGlobal(x3), -offset3));
2124  assert(basename != NULL);
2125  assert(N >= 1);
2126  assert(naddconss != NULL);
2127 
2128  SCIPdebugMsg(scip, "Creating linear Ben-Tal Nemirovski outer-approximation for <%s>.\n", basename);
2129 
2130  SCIP_CALL( SCIPallocBufferArray(scip, &avars, N+1) );
2131  SCIP_CALL( SCIPallocBufferArray(scip, &bvars, N+1) );
2132 
2133  /* create additional variables */
2134  for( i = 0; i <= N; ++i )
2135  {
2136  (void) SCIPsnprintf(varname, 255, "soc#%s_a%d", basename, i);
2137  SCIP_CALL( SCIPcreateVar(scip, &avars[i], varname, 0.0, SCIPinfinity(scip), 0.0,
2139  SCIP_CALL( SCIPaddVar(scip, avars[i]) );
2140 
2141  (void) SCIPsnprintf(varname, 255, "soc#%s_b%d", basename, i);
2142  SCIP_CALL( SCIPcreateVar(scip, &bvars[i], varname, 0.0, SCIPinfinity(scip), 0.0,
2144  SCIP_CALL( SCIPaddVar(scip, bvars[i]) );
2145  }
2146 
2147  /* create first linear constraints - split into two because of the absolute value */
2148  vars[0] = avars[0];
2149  vals[0] = 1.0;
2150  vars[1] = x1;
2151  vals[1] = -alpha1;
2152 
2153  (void) SCIPsnprintf(linname, 255, "soc#%s#a%d", basename, 0);
2154  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, alpha1 * offset1, SCIPinfinity(scip),
2159  TRUE /* removable */, SCIPconsIsStickingAtNode(cons)) );
2160  SCIP_CALL( SCIPaddCons(scip, lincons) );
2161  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2162  ++(*naddconss);
2163 
2164  vars[0] = avars[0];
2165  vals[0] = 1.0;
2166  vars[1] = x1;
2167  vals[1] = alpha1;
2168 
2169  (void) SCIPsnprintf(linname, 255, "soc#%s#A%d", basename, 0);
2170  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, -alpha1 * offset1, SCIPinfinity(scip),
2175  TRUE /* removable */, SCIPconsIsStickingAtNode(cons)) );
2176  SCIP_CALL( SCIPaddCons(scip, lincons) );
2177  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2178  ++(*naddconss);
2179 
2180  if( x2 != NULL )
2181  {
2182  vars[0] = bvars[0];
2183  vals[0] = 1.0;
2184  vars[1] = x2;
2185  vals[1] = -alpha2;
2186 
2187  (void) SCIPsnprintf(linname, 255, "soc#%s#b%d", basename, 0);
2188  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, alpha2 * offset2, SCIPinfinity(scip),
2193  TRUE /* removable */, SCIPconsIsStickingAtNode(cons)) );
2194  SCIP_CALL( SCIPaddCons(scip, lincons) );
2195  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2196  ++(*naddconss);
2197 
2198  vars[0] = bvars[0];
2199  vals[0] = 1.0;
2200  vars[1] = x2;
2201  vals[1] = alpha2;
2202 
2203  (void) SCIPsnprintf(linname, 255, "soc#%s#B%d", basename, 0);
2204  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, -alpha2 * offset2, SCIPinfinity(scip),
2209  TRUE /* removable */, SCIPconsIsStickingAtNode(cons)) );
2210  SCIP_CALL( SCIPaddCons(scip, lincons) );
2211  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2212  ++(*naddconss);
2213  }
2214  else
2215  { /* second summand is just a constant */
2216  if( SCIPconsIsLocal(cons) )
2217  {
2218  SCIP_CALL( SCIPchgVarLbNode(scip, NULL, bvars[0], ABS(alpha2 * offset2)) );
2219  }
2220  else
2221  {
2222  SCIP_CALL( SCIPchgVarLbGlobal(scip, bvars[0], ABS(alpha2 * offset2)) );
2223  }
2224  }
2225 
2226  /* create intermediate linear constraints */
2227  for( i = 1; i <= N; ++i )
2228  {
2229  SCIP_Real val;
2230 
2231  val = M_PI / pow(2.0, (double) (i+1));
2232 
2233  vars[0] = avars[i-1];
2234  vals[0] = cos(val);
2235  vars[1] = bvars[i-1];
2236  vals[1] = sin(val);
2237  vars[2] = avars[i];
2238  vals[2] = -1.0;
2239 
2240  (void) SCIPsnprintf(linname, 255, "soc#%s#a%d", basename, i);
2241  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 3, vars, vals, 0.0, 0.0,
2246  TRUE /* removable */, SCIPconsIsStickingAtNode(cons)) );
2247  SCIP_CALL( SCIPaddCons(scip, lincons) );
2248  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2249  ++(*naddconss);
2250 
2251  vars[0] = avars[i-1];
2252  vals[0] = sin(val);
2253  vars[1] = bvars[i-1];
2254  vals[1] = -cos(val);
2255  vars[2] = bvars[i];
2256  vals[2] = 1.0;
2257 
2258  (void) SCIPsnprintf(linname, 255, "soc#%s#b%d", basename, i);
2259  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 3, vars, vals, 0.0, SCIPinfinity(scip),
2264  TRUE /* removable */, SCIPconsIsStickingAtNode(cons)) );
2265  SCIP_CALL( SCIPaddCons(scip, lincons) );
2266  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2267  ++(*naddconss);
2268 
2269  vars[0] = avars[i-1];
2270  vals[0] = -sin(val);
2271  vars[1] = bvars[i-1];
2272  vals[1] = cos(val);
2273  vars[2] = bvars[i];
2274  vals[2] = 1.0;
2275 
2276  (void) SCIPsnprintf(linname, 255, "soc#%s#B%d", basename, i);
2277  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 3, vars, vals, 0.0, SCIPinfinity(scip),
2282  TRUE /* removable */, SCIPconsIsStickingAtNode(cons)) );
2283  SCIP_CALL( SCIPaddCons(scip, lincons) );
2284  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2285  ++(*naddconss);
2286  }
2287 
2288  /* create last linear constraints */
2289  vars[0] = x3;
2290  vals[0] = alpha3;
2291  vars[1] = avars[N];
2292  vals[1] = -1.0;
2293 
2294  (void) SCIPsnprintf(linname, 255, "soc#%s#a%d", basename, N);
2295  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, -alpha3 * offset3, SCIPinfinity(scip),
2301  SCIP_CALL( SCIPaddCons(scip, lincons) );
2302  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2303  ++(*naddconss);
2304 
2305  vars[0] = avars[N];
2306  vals[0] = tan( M_PI / pow(2.0, (double) (N+1)) );
2307  vars[1] = bvars[N];
2308  vals[1] = -1.0;
2309 
2310  (void) SCIPsnprintf(linname, 255, "soc#%s#b%d", basename, i);
2311  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, linname, 2, vars, vals, 0.0, SCIPinfinity(scip),
2316  TRUE /* removable */, SCIPconsIsStickingAtNode(cons)) );
2317  SCIP_CALL( SCIPaddCons(scip, lincons) );
2318  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
2319  ++(*naddconss);
2320 
2321  for( i = 0; i <= N; ++i )
2322  {
2323  SCIP_CALL( SCIPreleaseVar(scip, &avars[i]) );
2324  SCIP_CALL( SCIPreleaseVar(scip, &bvars[i]) );
2325  }
2326  SCIPfreeBufferArray(scip, &avars);
2327  SCIPfreeBufferArray(scip, &bvars);
2328 
2329  return SCIP_OKAY;
2330 }
2331 
2332 /** adds a linear outer approx for a three dimensional SOC constraint
2333  *
2334  * chooses between Ben-Tan/Nemirovski and Glineur and calls the corresponding function
2335  */
2336 static
2338  SCIP* scip, /**< SCIP data structure */
2339  SCIP_CONS* cons, /**< original constraint */
2340  SCIP_VAR* x1, /**< variable x1 */
2341  SCIP_VAR* x2, /**< variable x2 */
2342  SCIP_VAR* x3, /**< variable x3 */
2343  SCIP_Real alpha1, /**< coefficient of x1 */
2344  SCIP_Real alpha2, /**< coefficient of x2 */
2345  SCIP_Real alpha3, /**< coefficient of x3 */
2346  SCIP_Real offset1, /**< offset of x1 */
2347  SCIP_Real offset2, /**< offset of x2 */
2348  SCIP_Real offset3, /**< offset of x3 */
2349  int N, /**< size of linear approximation, need to be >= 1 */
2350  SCIP_Bool glineur, /**< whether to prefer Glineur to Ben-Tal Nemirovski */
2351  const char* basename, /**< string to use for building variable and constraint names */
2352  int* naddconss /**< buffer where to add the number of added constraints */
2353  )
2354 {
2355  if( glineur )
2356  {
2357  SCIP_CALL( presolveCreateGlineurApproxDim3(scip, cons, x1, x2, x3, alpha1, alpha2, alpha3, offset1, offset2, offset3, N, basename, naddconss) );
2358  }
2359  else
2360  {
2361  SCIP_CALL( presolveCreateBenTalNemirovskiApproxDim3(scip, cons, x1, x2, x3, alpha1, alpha2, alpha3, offset1, offset2, offset3, N, basename, naddconss) );
2362  }
2363 
2364  return SCIP_OKAY;
2365 }
2366 
2367 /** 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
2368  *
2369  * if n > 2, calls same function recursively;
2370  * if n = 2, calls presolveCreateBenTalNemirovskiApproxDim3
2371  */
2372 static
2374  SCIP* scip, /**< SCIP data structure */
2375  int nlhsvars, /**< number of variables on left hand side (n) */
2376  SCIP_VAR** lhsvars, /**< variables on left hand side (x_i) */
2377  SCIP_Real* lhscoefs, /**< coefficients of variables on left hand side (alpha_i) */
2378  SCIP_Real* lhsoffsets, /**< offsets of variable on left hand side (beta_i) */
2379  SCIP_VAR* rhsvar, /**< variable on right hand side (y) */
2380  SCIP_Real rhscoeff, /**< coefficient of variable on right hand side (alpha_{n+1}) */
2381  SCIP_Real rhsoffset, /**< offset of variable on right hand side (beta_{n+1}) */
2382  SCIP_Real constant, /**< constant term (gamma) */
2383  const char* basename, /**< prefix for variable and constraint name */
2384  SCIP_CONS* origcons, /**< original constraint for which this SOC3 set is added */
2385  int soc3_nr_auxvars, /**< number of auxiliary variables to use for a SOC3 constraint, or 0 if automatic */
2386  SCIP_Bool glineur, /**< whether Glineur should be preferred to Ben-Tal Nemirovski */
2387  int* naddconss /**< buffer where to add the number of added constraints */
2388  )
2389 {
2390  char name[255];
2391  SCIP_VAR* auxvar1;
2392  SCIP_VAR* auxvar2;
2393 
2394  assert(scip != NULL);
2395  assert(lhsvars != NULL);
2396  assert(nlhsvars >= 1);
2397  assert(lhscoefs != NULL);
2398  assert(lhsoffsets != NULL);
2399  assert(rhsvar != NULL);
2400  assert(basename != NULL);
2401  assert(!SCIPisNegative(scip, constant));
2402  assert(naddconss != NULL);
2403 
2404  if( nlhsvars == 1 )
2405  { /* end of recursion */
2406  assert(SCIPisPositive(scip, constant));
2407  SCIP_CALL( presolveCreateOuterApproxDim3(scip, origcons,
2408  lhsvars[0], NULL, rhsvar,
2409  lhscoefs[0], 1.0, rhscoeff,
2410  lhsoffsets[0], sqrt(constant), rhsoffset,
2411  soc3_nr_auxvars, glineur, basename, naddconss) );
2412 
2413  return SCIP_OKAY;
2414  }
2415 
2416  if( nlhsvars == 2 && SCIPisZero(scip, constant) )
2417  { /* end of recursion */
2418  assert(lhsvars[0] != NULL);
2419  assert(lhsvars[1] != NULL);
2420  assert(rhsvar != NULL);
2421  SCIP_CALL( presolveCreateOuterApproxDim3(scip, origcons,
2422  lhsvars[0], lhsvars[1], rhsvar,
2423  lhscoefs[0], lhscoefs[1], rhscoeff,
2424  lhsoffsets[0], lhsoffsets[1], rhsoffset,
2425  soc3_nr_auxvars, glineur, basename, naddconss) );
2426 
2427  return SCIP_OKAY;
2428  }
2429 
2430  if( nlhsvars == 3 || (nlhsvars == 2 && !SCIPisZero(scip, constant)) )
2431  {
2432  /* a bit special case too */
2433  /* for first two variables on lhs, create a new aux.var and a new SOC3 */
2434  (void) SCIPsnprintf(name, 255, "%s#z1", basename);
2435  SCIP_CALL( SCIPcreateVar(scip, &auxvar1, name, 0.0, SCIPinfinity(scip), 0.0,
2437  SCIP_CALL( SCIPaddVar(scip, auxvar1) );
2438 
2439  /* constraint alpha_0 (x_0+beta0)^2 + alpha_1 (x_1+beta1)^2 <= auxvar^2 */
2440  SCIP_CALL( presolveCreateOuterApproxDim3(scip, origcons,
2441  lhsvars[0], lhsvars[1], auxvar1,
2442  lhscoefs[0], lhscoefs[1], 1.0,
2443  lhsoffsets[0], lhsoffsets[1], 0.0,
2444  soc3_nr_auxvars, glineur, name, naddconss) );
2445 
2446  (void) SCIPsnprintf(name, 255, "%s_soc3", basename);
2447  if( nlhsvars == 3 )
2448  { /* create new constraint alpha_2 (x_2+beta2)^2 + auxvar^2 <= (rhscoeff * (rhsvar+rhsoffset))^2 */
2449  SCIP_CALL( presolveCreateOuterApproxDim3(scip, origcons,
2450  lhsvars[2], auxvar1, rhsvar,
2451  lhscoefs[2], 1.0, rhscoeff,
2452  lhsoffsets[2], 0.0, rhsoffset,
2453  soc3_nr_auxvars, glineur, name, naddconss) );
2454  }
2455  else
2456  { /* create new constraint auxvar^2 + sqrt(constant)^2 <= (rhscoeff * (rhsvar+rhsoffset))^2 */
2457  SCIP_CALL( presolveCreateOuterApproxDim3(scip, origcons,
2458  auxvar1, NULL, rhsvar,
2459  1.0, 1.0, rhscoeff,
2460  0.0, sqrt(constant), rhsoffset,
2461  soc3_nr_auxvars, glineur, name, naddconss) );
2462  }
2463 
2464  SCIP_CALL( SCIPreleaseVar(scip, &auxvar1) );
2465 
2466  return SCIP_OKAY;
2467  }
2468 
2469  /* nlhsvars >= 4 */
2470 
2471  (void) SCIPsnprintf(name, 255, "%s#z1", basename);
2472  SCIP_CALL( SCIPcreateVar(scip, &auxvar1, name, 0.0, SCIPinfinity(scip), 0.0,
2474  SCIP_CALL( SCIPaddVar(scip, auxvar1) );
2475 
2476  /* approx for left half of lhs */
2478  nlhsvars/2, lhsvars, lhscoefs, lhsoffsets,
2479  auxvar1, 1.0, 0.0,
2480  constant, name, origcons, soc3_nr_auxvars, glineur, naddconss) );
2481 
2482  (void) SCIPsnprintf(name, 255, "%s#z2", basename);
2483  SCIP_CALL( SCIPcreateVar(scip, &auxvar2, name, 0., SCIPinfinity(scip), 0.0,
2485  SCIP_CALL( SCIPaddVar(scip, auxvar2) );
2486 
2487  /* approx for right half of lhs */
2489  nlhsvars-nlhsvars/2, &lhsvars[nlhsvars/2], &lhscoefs[nlhsvars/2], &lhsoffsets[nlhsvars/2],
2490  auxvar2, 1.0, 0.0,
2491  0.0, name, origcons, soc3_nr_auxvars, glineur, naddconss) );
2492 
2493  /* SOC constraint binding both auxvar's */
2494  (void)SCIPsnprintf(name, 255, "%s_soc3", basename);
2495  SCIP_CALL( presolveCreateOuterApproxDim3(scip, origcons,
2496  auxvar1, auxvar2, rhsvar,
2497  1.0, 1.0, rhscoeff,
2498  0.0, 0.0, rhsoffset,
2499  soc3_nr_auxvars, glineur, name, naddconss) );
2500 
2501  SCIP_CALL( SCIPreleaseVar(scip, &auxvar1) );
2502  SCIP_CALL( SCIPreleaseVar(scip, &auxvar2) );
2503 
2504  return SCIP_OKAY;
2505 }
2506 
2507 /** propagates variable bounds */
2508 static
2510  SCIP* scip, /**< SCIP data structure */
2511  SCIP_CONS* cons, /**< constraint */
2512  SCIP_RESULT* result, /**< buffer to store result of propagation */
2513  int* nchgbds, /**< buffer where to add number of tightened bounds */
2514  SCIP_Bool* redundant /**< buffer to indicate whether constraint was marked for deletion because of redundancy */
2515  )
2516 {
2517  SCIP_CONSDATA* consdata;
2518  SCIP_INTERVAL lhsrange;
2519  SCIP_INTERVAL* lhsranges;
2520  SCIP_INTERVAL rhsrange;
2521  SCIP_INTERVAL a, b, c;
2522  SCIP_ROUNDMODE roundmode;
2523  SCIP_Bool infeas, tightened;
2524  int i;
2525  SCIP_Real lb, ub;
2526 
2527  assert(scip != NULL);
2528  assert(cons != NULL);
2529  assert(result != NULL);
2530  assert(nchgbds != NULL);
2531  assert(redundant != NULL);
2532 
2533  consdata = SCIPconsGetData(cons);
2534  assert(consdata != NULL);
2535 
2536  *redundant = FALSE;
2537 
2538  if( !SCIPconsIsMarkedPropagate(cons) )
2539  {
2540  SCIPdebugMsg(scip, "skip propagation for constraint %s\n", SCIPconsGetName(cons));
2541  *result = SCIP_DIDNOTRUN;
2542  return SCIP_OKAY;
2543  }
2544  else
2545  {
2546  SCIPdebugMsg(scip, "try propagation for constraint %s\n", SCIPconsGetName(cons));
2547  }
2548 
2549  *result = SCIP_DIDNOTFIND;
2550  SCIP_CALL( SCIPunmarkConsPropagate(scip, cons) );
2551 
2552  /* @todo do something clever to decide whether propagation should be tried */
2553 
2554  /* compute activity on lhs */
2555  SCIPintervalSet(&lhsrange, consdata->constant);
2556  SCIP_CALL( SCIPallocBufferArray(scip, &lhsranges, consdata->nvars) );
2557  for( i = 0; i < consdata->nvars; ++i )
2558  {
2559  lb = SCIPcomputeVarLbLocal(scip, consdata->vars[i]) - SCIPepsilon(scip);
2560  ub = SCIPcomputeVarUbLocal(scip, consdata->vars[i]) + SCIPepsilon(scip);
2561  SCIPintervalSetBounds(&lhsranges[i], MIN(lb, ub), MAX(lb, ub));
2562  if( consdata->offsets[i] != 0.0 )
2563  SCIPintervalAddScalar(SCIPinfinity(scip), &lhsranges[i], lhsranges[i], consdata->offsets[i]);
2564  if( consdata->coefs[i] != 1.0 )
2565  SCIPintervalMulScalar(SCIPinfinity(scip), &lhsranges[i], lhsranges[i], consdata->coefs[i]);
2566  SCIPintervalSquare(SCIPinfinity(scip), &lhsranges[i], lhsranges[i]);
2567 
2568  SCIPintervalAdd(SCIPinfinity(scip), &lhsrange, lhsrange, lhsranges[i]);
2569  }
2570 
2571  /* compute activity on rhs: rhsrange = sqr(rhscoeff * (rhsvar + rhsoffset) ) */
2572  lb = SCIPcomputeVarLbLocal(scip, consdata->rhsvar) - SCIPepsilon(scip);
2573  ub = SCIPcomputeVarUbLocal(scip, consdata->rhsvar) + SCIPepsilon(scip);
2574  SCIPintervalSetBounds(&rhsrange, MIN(lb, ub), MAX(lb, ub));
2575 
2576  if( consdata->rhsoffset != 0.0 )
2577  SCIPintervalAddScalar(SCIPinfinity(scip), &rhsrange, rhsrange, consdata->rhsoffset);
2578  if( consdata->rhscoeff != 1.0 )
2579  SCIPintervalMulScalar(SCIPinfinity(scip), &rhsrange, rhsrange, consdata->rhscoeff);
2580  SCIPintervalSquare(SCIPinfinity(scip), &rhsrange, rhsrange);
2581 
2582  /* check for infeasibility */
2583  if( SCIPisGT(scip, lhsrange.inf-SCIPfeastol(scip), rhsrange.sup) )
2584  {
2585  SCIPdebugMsg(scip, "propagation found constraint <%s> infeasible: lhs = [%.15g,%.15g]-feastol-eps > rhs = [%.15g,%.15g]\n",
2586  SCIPconsGetName(cons), lhsrange.inf, lhsrange.sup, rhsrange.inf, rhsrange.sup);
2587  *result = SCIP_CUTOFF;
2588  goto TERMINATE;
2589  }
2590 
2591  /* check for redundancy: max(lhsrange) <= min(rhsrange) */
2592  if( SCIPisLE(scip, lhsrange.sup, rhsrange.inf) )
2593  {
2594  SCIPdebugMsg(scip, "propagation found constraint <%s> redundant: lhs = [%.15g,%.15g] <= rhs = [%.15g,%.15g]\n",
2595  SCIPconsGetName(cons), lhsrange.inf, lhsrange.sup, rhsrange.inf, rhsrange.sup);
2596  SCIP_CALL( SCIPdelConsLocal(scip, cons) );
2597  goto TERMINATE;
2598  }
2599 
2600  /* try to tighten variable on rhs */
2601  if( SCIPvarGetStatus(consdata->rhsvar) != SCIP_VARSTATUS_MULTAGGR )
2602  {
2603  SCIPintervalSquareRoot(SCIPinfinity(scip), &a, lhsrange);
2604  if( consdata->rhscoeff != 1.0 )
2605  SCIPintervalDivScalar(SCIPinfinity(scip), &a, a, consdata->rhscoeff);
2606  if( consdata->rhsoffset != 0.0 )
2607  SCIPintervalSubScalar(SCIPinfinity(scip), &a, a, consdata->rhsoffset);
2608  SCIP_CALL( SCIPtightenVarLb(scip, consdata->rhsvar, SCIPintervalGetInf(a), FALSE, &infeas, &tightened) );
2609  if( infeas )
2610  {
2611  SCIPdebugMsg(scip, "propagation found constraint <%s> infeasible\n", SCIPconsGetName(cons));
2612  *result = SCIP_CUTOFF;
2613  goto TERMINATE;
2614  }
2615  if( tightened )
2616  {
2617  SCIPdebugMsg(scip, "propagation tightened bounds of rhs variable <%s> in constraint <%s>\n", SCIPvarGetName(consdata->rhsvar), SCIPconsGetName(cons));
2618  *result = SCIP_REDUCEDDOM;
2619  ++*nchgbds;
2620  }
2621  }
2622 
2623  /* try to tighten variables on lhs */
2624  SCIPintervalSub(SCIPinfinity(scip), &b, rhsrange, lhsrange); /*lint !e644 */
2625  for( i = 0; i < consdata->nvars; ++i )
2626  {
2627  if( SCIPvarGetStatus(consdata->vars[i]) == SCIP_VARSTATUS_MULTAGGR )
2628  continue;
2629 
2630  roundmode = SCIPintervalGetRoundingMode();
2631  if( !SCIPisInfinity(scip, b.sup) )
2632  {
2634  a.sup = b.sup + lhsranges[i].inf;
2635  }
2636  else
2637  {
2638  a.sup = SCIPinfinity(scip);
2639  }
2640  if( !SCIPisInfinity(scip, -b.inf) )
2641  {
2643  a.inf = b.inf + lhsranges[i].sup;
2644  }
2645  else
2646  {
2647  a.inf = -SCIPinfinity(scip);
2648  }
2649  SCIPintervalSetRoundingMode(roundmode);
2650  SCIPintervalSquareRoot(SCIPinfinity(scip), &a, a);
2651 
2652  assert(consdata->coefs[i] >= 0.0); /* should be ensured in create and presolveRemoveFixed */
2653 
2654  c = a;
2655  if( consdata->coefs[i] != 1.0 )
2656  SCIPintervalDivScalar(SCIPinfinity(scip), &c, c, consdata->coefs[i]);
2657  if( consdata->offsets[i] != 0.0 )
2658  SCIPintervalSubScalar(SCIPinfinity(scip), &c, c, consdata->offsets[i]);
2659 
2660  SCIP_CALL( SCIPtightenVarUb(scip, consdata->vars[i], SCIPintervalGetSup(c), FALSE, &infeas, &tightened) );
2661  if( infeas )
2662  {
2663  SCIPdebugMsg(scip, "propagation found constraint <%s> infeasible\n", SCIPconsGetName(cons));
2664  *result = SCIP_CUTOFF;
2665  goto TERMINATE;
2666  }
2667  if( tightened )
2668  {
2669  SCIPdebugMsg(scip, "propagation tightened bounds of lhs variable <%s> in constraint <%s>\n", SCIPvarGetName(consdata->vars[i]), SCIPconsGetName(cons));
2670  *result = SCIP_REDUCEDDOM;
2671  ++*nchgbds;
2672  }
2673 
2674  c = a;
2675  SCIPintervalDivScalar(SCIPinfinity(scip), &c, c, -consdata->coefs[i]);
2676  if( consdata->offsets[i] != 0.0 )
2677  SCIPintervalSubScalar(SCIPinfinity(scip), &c, c, consdata->offsets[i]);
2678 
2679  SCIP_CALL( SCIPtightenVarLb(scip, consdata->vars[i], SCIPintervalGetInf(c), FALSE, &infeas, &tightened) );
2680  if( infeas )
2681  {
2682  SCIPdebugMsg(scip, "propagation found constraint <%s> infeasible\n", SCIPconsGetName(cons));
2683  *result = SCIP_CUTOFF;
2684  goto TERMINATE;
2685  }
2686  if( tightened )
2687  {
2688  SCIPdebugMsg(scip, "propagation tightened bounds of lhs variable <%s> in constraint <%s>\n", SCIPvarGetName(consdata->vars[i]), SCIPconsGetName(cons));
2689  *result = SCIP_REDUCEDDOM;
2690  ++*nchgbds;
2691  }
2692  }
2693 
2694 TERMINATE:
2695  SCIPfreeBufferArray(scip, &lhsranges);
2696 
2697  if( *result != SCIP_DIDNOTFIND )
2698  {
2699  SCIP_CALL( SCIPresetConsAge(scip, cons) );
2700  }
2701 
2702  return SCIP_OKAY;
2703 }
2704 
2705 /** tries to adjust a solution such that it satisfies a given constraint by increasing the value for the constraints right hand side variable */
2706 static
2708  SCIP* scip, /**< SCIP data structure */
2709  SCIP_CONS* cons, /**< constraint */
2710  SCIP_SOL* sol, /**< solution to polish */
2711  SCIP_Bool* success /**< buffer to store whether polishing was successful */
2712  )
2713 {
2714  SCIP_CONSDATA* consdata;
2715  SCIP_Real rhsval;
2716 
2717  assert(scip != NULL);
2718  assert(cons != NULL);
2719  assert(sol != NULL);
2720  assert(success != NULL);
2721 
2722  consdata = SCIPconsGetData(cons);
2723  assert(consdata != NULL);
2724  assert(!SCIPisZero(scip, consdata->rhscoeff));
2725 
2726  /* compute minimal rhs variable value so that constraint is satisfied */
2727  if( !SCIPisInfinity(scip, consdata->lhsval) )
2728  rhsval = consdata->lhsval / consdata->rhscoeff - consdata->rhsoffset;
2729  else
2730  rhsval = consdata->rhscoeff > 0.0 ? SCIPinfinity(scip) : -SCIPinfinity(scip);
2731 
2732  if( consdata->rhscoeff > 0.0 )
2733  {
2734  assert(SCIPvarMayRoundUp(consdata->rhsvar));
2735 
2736  /* round rhsval up, if variable is integral */
2737  if( SCIPvarIsIntegral(consdata->rhsvar) && !SCIPisInfinity(scip, rhsval) )
2738  rhsval = SCIPceil(scip, rhsval);
2739 
2740  /* if new value is above upper bound, we are lost */
2741  if( SCIPisGT(scip, rhsval, SCIPvarGetUbGlobal(consdata->rhsvar)) )
2742  {
2743  *success = FALSE;
2744  }
2745  else
2746  {
2747  /* if new value is larger then current one, increase to new value */
2748  if( rhsval > SCIPgetSolVal(scip, sol, consdata->rhsvar) )
2749  {
2750  SCIPdebugMsg(scip, "increase <%s> to %g\n", SCIPvarGetName(consdata->rhsvar), rhsval);
2751  SCIP_CALL( SCIPsetSolVal(scip, sol, consdata->rhsvar, rhsval) );
2752  }
2753 
2754  *success = TRUE;
2755  }
2756  }
2757  else
2758  {
2759  assert(SCIPvarMayRoundDown(consdata->rhsvar));
2760 
2761  /* round rhsval down, if variable is integral */
2762  if( SCIPvarIsIntegral(consdata->rhsvar) )
2763  rhsval = SCIPfloor(scip, rhsval);
2764 
2765  /* if new value is below lower bound, we are lost */
2766  if( SCIPisLT(scip, rhsval, SCIPvarGetLbGlobal(consdata->rhsvar)) )
2767  {
2768  *success = FALSE;
2769  }
2770  else
2771  {
2772  /* if new value is below current one, decrease to new value */
2773  if( rhsval < SCIPgetSolVal(scip, sol, consdata->rhsvar) )
2774  {
2775  SCIPdebugMsg(scip, "decrease <%s> to %g\n", SCIPvarGetName(consdata->rhsvar), rhsval);
2776  SCIP_CALL( SCIPsetSolVal(scip, sol, consdata->rhsvar, rhsval) );
2777  }
2778 
2779  *success = TRUE;
2780  }
2781  }
2782 
2783  SCIPdebugMsg(scip, "polishing solution for constraint <%s> was %ssuccessful\n", SCIPconsGetName(cons), *success ? "" : "not ");
2784 
2785  return SCIP_OKAY;
2786 }
2787 
2788 /** disaggregates a (sufficiently large) SOC constraint into smaller ones; for each term on the lhs we add a quadratic
2789  * 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
2790  * \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
2791  * violations of all quadratic constraints sum up we scale each constraint by the number of lhs terms + 1
2792  *
2793  * @todo if rhsvar is NULL, then the disaggregation does not produce further cones. Should it then be upgraded
2794  * to a quadratic and let the quadratic desaggregate it?
2795  * The code assumes now that the rhsvar is not NULL in order build the direct SOC -> SOC disaggregation
2796  */
2797 static
2799  SCIP* scip, /**< SCIP data structure */
2800  SCIP_CONS* cons, /**< constraint */
2801  SCIP_CONSDATA* consdata, /**< constraint data */
2802  int* naddconss, /**< pointer to count total number of added constraints */
2803  int* ndelconss, /**< pointer to count total number of deleted constraints */
2804  SCIP_Bool* success /**< pointer to store whether disaggregation was successful */
2805  )
2806 {
2807  SCIP_CONS* discons;
2808  SCIP_VAR** disvars;
2809  SCIP_VAR** sumvars;
2810  SCIP_VAR** difvars;
2811  SCIP_Real* discoefs;
2812  SCIP_VAR* lhsvars[2];
2813  SCIP_VAR* aggvars[2];
2814  SCIP_Real coefs[2];
2815  SCIP_Real offsets[2];
2816  SCIP_Real scalars[2];
2817  char name[SCIP_MAXSTRLEN];
2818  SCIP_Real constant;
2819  SCIP_Real scale;
2820  SCIP_Bool infeas;
2821  int ndisvars;
2822  int i;
2823 
2824  assert(naddconss != NULL);
2825  assert(ndelconss != NULL);
2826  assert(success != NULL);
2827 
2828  *success = FALSE;
2829 
2830  /* disaggregation does not make much sense if there are too few variables */
2831  if( consdata->nvars < 3 )
2832  {
2833  SCIPdebugMsg(scip, "can not disaggregate too small soc constraint %s\n", SCIPconsGetName(cons));
2834  return SCIP_OKAY;
2835  }
2836 
2837  if( consdata->rhsvar == NULL )
2838  {
2839  SCIPdebugMsg(scip, "can not disaggregate directly into a soc without rhs var %s\n", SCIPconsGetName(cons));
2840  return SCIP_OKAY;
2841  }
2842 
2843  /* there are at most n + 2 many linear varibles */
2844  SCIP_CALL( SCIPallocBufferArray(scip, &disvars, consdata->nvars + 2) );
2845  SCIP_CALL( SCIPallocBufferArray(scip, &sumvars, consdata->nvars + 2) );
2846  SCIP_CALL( SCIPallocBufferArray(scip, &difvars, consdata->nvars + 2) );
2847  SCIP_CALL( SCIPallocBufferArray(scip, &discoefs, consdata->nvars + 2) );
2848  ndisvars = 0;
2849 
2850  scale = 1.0 * (consdata->nvars + 1)/4.0;
2851 
2852  /* add (*) (alpha_i * (x_i + beta_i))^2 <= alpha_{n+1} * (x_{n+1} + beta_{n+1}) * z_i:
2853  * create sumvar = alpha_{n+1} * (x_{n+1} + beta_{n+1}) + z_i (multiagg)
2854  * create difvar = alpha_{n+1} * (x_{n+1} + beta_{n+1}) - z_i (multiagg)
2855  * note that (*) is equiv to sqrt( (2 * alpha_i * (x_i + beta_i))^2 + difvar^2) <= sumvar
2856  * scaling give us: sqrt( (2 * scale * alpha_i * (x_i + beta_i))^2 + (scale * difvar)^2) <= scale * sumvar
2857  */
2858  aggvars[0] = consdata->rhsvar;
2859  scalars[0] = consdata->rhscoeff;
2860  for( i = 0; i < consdata->nvars; ++i )
2861  {
2862  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedis_%s_%d", SCIPvarGetName(consdata->vars[i]), i);
2863  SCIP_CALL( SCIPcreateVar(scip, &disvars[i], name, 0.0, SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS, TRUE, FALSE,
2864  NULL, NULL, NULL, NULL, NULL) );
2865  SCIP_CALL( SCIPaddVar(scip, disvars[i]) );
2866 
2867  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedisS_%s_%d", SCIPvarGetName(consdata->vars[i]), i);
2868  SCIP_CALL( SCIPcreateVar(scip, &sumvars[i], name, 0.0, SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS, TRUE, FALSE,
2869  NULL, NULL, NULL, NULL, NULL) );
2870  SCIP_CALL( SCIPaddVar(scip, sumvars[i]) );
2871 
2872  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedisD_%s_%d", SCIPvarGetName(consdata->vars[i]), i);
2873  SCIP_CALL( SCIPcreateVar(scip, &difvars[i], name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0,
2875  SCIP_CALL( SCIPaddVar(scip, difvars[i]) );
2876 
2877  aggvars[1] = disvars[i];
2878  scalars[1] = 1.0;
2879  constant = consdata->rhscoeff * consdata->rhsoffset;
2880  SCIP_CALL( SCIPmultiaggregateVar(scip, sumvars[i], 2, aggvars, scalars, constant, &infeas, success) );
2881  /* @todo what shall we do if multiagg fails? */
2882  assert(!infeas && *success);
2883 
2884  scalars[1] = -1.0;
2885  SCIP_CALL( SCIPmultiaggregateVar(scip, difvars[i], 2, aggvars, scalars, constant, &infeas, success) );
2886  assert(!infeas && *success);
2887 
2888  /* create soc */
2889  lhsvars[0] = difvars[i];
2890  coefs[0] = scale;
2891  offsets[0] = 0.0;
2892  lhsvars[1] = consdata->vars[i];
2893  coefs[1] = scale * 2 * consdata->coefs[i];
2894  offsets[1] = consdata->offsets[i];
2895 
2896  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "consdis_%s_%d", SCIPconsGetName(cons), i);
2897  SCIP_CALL( SCIPcreateConsBasicSOC(scip, &discons, name, 2, lhsvars, coefs, offsets, 0.0, sumvars[i], scale, 0.0) );
2898  SCIP_CALL( SCIPaddCons(scip, discons) );
2899 #ifdef SCIP_DEBUG
2900  SCIP_CALL( SCIPprintCons(scip, discons, NULL) );
2901 #endif
2902  SCIP_CALL( SCIPreleaseCons(scip, &discons) );
2903  ++(*naddconss);
2904 
2905  /* linear coefficient in the linear constraint */
2906  discoefs[ndisvars] = 1.0;
2907  ++ndisvars;
2908  }
2909  assert(ndisvars == consdata->nvars);
2910 
2911  /* add gamma <= alpha_{n+1} * (x_{n+1} + beta_{n+1}) * z_i
2912  * sumvar and difvar are the same as before, but the equivalent soc now is
2913  * sqrt(4 * gamma + difvar^2) <= sumvar
2914  * scaling give us: sqrt( (4 * scale^2 * gamma + (scale * difvar)^2) <= scale * sumvar
2915  */
2916  if( !SCIPisZero(scip, consdata->constant) )
2917  {
2918  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedis_const_%s", SCIPconsGetName(cons));
2919  SCIP_CALL( SCIPcreateVar(scip, &disvars[ndisvars], name, 0.0, SCIPinfinity(scip), 0.0,
2921  SCIP_CALL( SCIPaddVar(scip, disvars[ndisvars]) );
2922 
2923  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedisS_const_%s", SCIPconsGetName(cons));
2924  SCIP_CALL( SCIPcreateVar(scip, &sumvars[ndisvars], name, 0.0, SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS, TRUE, FALSE,
2925  NULL, NULL, NULL, NULL, NULL) );
2926  SCIP_CALL( SCIPaddVar(scip, sumvars[ndisvars]) );
2927 
2928  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "conedisD_const_%s", SCIPconsGetName(cons));
2929  SCIP_CALL( SCIPcreateVar(scip, &difvars[ndisvars], name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0,
2931  SCIP_CALL( SCIPaddVar(scip, difvars[ndisvars]) );
2932 
2933  aggvars[1] = disvars[i];
2934  scalars[1] = 1.0;
2935  constant = consdata->rhscoeff * consdata->rhsoffset;
2936  SCIP_CALL( SCIPmultiaggregateVar(scip, sumvars[i], 2, aggvars, scalars, constant, &infeas, success) );
2937  assert(!infeas && *success);
2938 
2939  scalars[1] = -1.0;
2940  SCIP_CALL( SCIPmultiaggregateVar(scip, difvars[i], 2, aggvars, scalars, constant, &infeas, success) );
2941  assert(!infeas && *success);
2942 
2943  /* create soc */
2944  lhsvars[0] = difvars[ndisvars];
2945  coefs[0] = scale;
2946  offsets[0] = 0.0;
2947  constant = 4.0 * SQR(scale) * consdata->constant;
2948  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "consdis_%s_constant", SCIPconsGetName(cons));
2949  SCIP_CALL( SCIPcreateConsBasicSOC(scip, &discons, name, 1, lhsvars, coefs, offsets, constant,
2950  sumvars[ndisvars], scale, 0.0) );
2951  SCIP_CALL( SCIPaddCons(scip, discons) );
2952  SCIP_CALL( SCIPreleaseCons(scip, &discons) );
2953  ++(*naddconss);
2954 
2955  /* linear coefficient in the linear constraint */
2956  discoefs[ndisvars] = 1.0;
2957  ++ndisvars;
2958  }
2959 
2960  /* create linear constraint sum z_i <= alpha_{n+1} * (x_{n+1} + beta_{n+1}); first add extra coefficient for the rhs */
2961  discoefs[ndisvars] = -1.0 * consdata->rhscoeff;
2962  disvars[ndisvars] = consdata->rhsvar;
2963 
2964  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "consdis_linear_%s", SCIPconsGetName(cons));
2965  SCIP_CALL( SCIPcreateConsBasicLinear(scip, &discons, name, ndisvars + 1, disvars, discoefs, -SCIPinfinity(scip),
2966  consdata->rhscoeff * consdata->rhsoffset) );
2967 
2968  SCIP_CALL( SCIPaddCons(scip, discons) );
2969  SCIP_CALL( SCIPreleaseCons(scip, &discons) );
2970  ++(*naddconss);
2971 
2972  /* release all variables */
2973  for( i = ndisvars - 1; i >= 0; --i )
2974  {
2975  SCIP_CALL( SCIPreleaseVar(scip, &disvars[i]) );
2976  SCIP_CALL( SCIPreleaseVar(scip, &sumvars[i]) );
2977  SCIP_CALL( SCIPreleaseVar(scip, &difvars[i]) );
2978  }
2979  SCIPfreeBufferArray(scip, &discoefs);
2980  SCIPfreeBufferArray(scip, &difvars);
2981  SCIPfreeBufferArray(scip, &sumvars);
2982  SCIPfreeBufferArray(scip, &disvars);
2983 
2984  /* delete constraint */
2985  SCIP_CALL( SCIPdelCons(scip, cons) );
2986  ++(*ndelconss);
2987 
2988  *success = TRUE;
2989 
2990  return SCIP_OKAY;
2991 }
2992 
2993 
2994 /** helper function to enforce constraints */
2995 static
2997  SCIP* scip, /**< SCIP data structure */
2998  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
2999  SCIP_CONS** conss, /**< constraints to process */
3000  int nconss, /**< number of constraints */
3001  int nusefulconss, /**< number of useful (non-obsolete) constraints to process */
3002  SCIP_SOL* sol, /**< solution to enforce (NULL for the LP solution) */
3003  SCIP_RESULT* result /**< pointer to store the result of the enforcing call */
3004  )
3005 {
3006  SCIP_CONSHDLRDATA* conshdlrdata;
3007  SCIP_CONSDATA* consdata;
3008  SCIP_CONS* maxviolcons;
3009  SCIP_Bool success;
3010  SCIP_Bool cutoff;
3011  SCIP_Bool redundant;
3012  int nbndchg;
3013  int c;
3014 
3015  assert(scip != NULL);
3016  assert(conshdlr != NULL);
3017  assert(conss != NULL || nconss == 0);
3018  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
3019  assert(result != NULL);
3020 
3021  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3022  assert(conshdlrdata != NULL);
3023 
3024  SCIP_CALL( computeViolations(scip, conshdlr, conss, nconss, sol, &maxviolcons) );
3025 
3026  if( maxviolcons == NULL )
3027  {
3028  *result = SCIP_FEASIBLE;
3029  return SCIP_OKAY;
3030  }
3031 
3032  /* if we are above the 100'th enforcement round for this node, something is strange
3033  * (maybe the LP does not think that the cuts we add are violated, or we do ECP on a high-dimensional convex function)
3034  * 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
3035  * (in optimized more, returning SCIP_INFEASIBLE in *result would be sufficient, but in debug mode this would give an assert in scip.c)
3036  * the reason to wait for 100 rounds is to avoid calls to SCIPisStopped in normal runs, which may be expensive
3037  * we only increment nenforounds until 101 to avoid an overflow
3038  */
3039  if( conshdlrdata->lastenfonode == SCIPgetCurrentNode(scip) )
3040  {
3041  if( conshdlrdata->nenforounds > 100 )
3042  {
3043  if( SCIPisStopped(scip) )
3044  {
3045  SCIP_NODE* child;
3046 
3047  SCIP_CALL( SCIPcreateChild(scip, &child, 1.0, SCIPnodeGetEstimate(SCIPgetCurrentNode(scip))) );
3048  *result = SCIP_BRANCHED;
3049 
3050  return SCIP_OKAY;
3051  }
3052  }
3053  else
3054  ++conshdlrdata->nenforounds;
3055  }
3056  else
3057  {
3058  conshdlrdata->lastenfonode = SCIPgetCurrentNode(scip);
3059  conshdlrdata->nenforounds = 0;
3060  }
3061 
3062  /* try separation, this should usually work */
3063  SCIP_CALL( separatePoint(scip, conshdlr, conss, nconss, nusefulconss, sol, TRUE, &cutoff, &success) );
3064  if( cutoff )
3065  {
3066  *result = SCIP_CUTOFF;
3067  return SCIP_OKAY;
3068  }
3069  if( success )
3070  {
3071  SCIPdebugMsg(scip, "enforced by separation\n");
3072  *result = SCIP_SEPARATED;
3073  return SCIP_OKAY;
3074  }
3075 
3076  /* try propagation */
3077  for( c = 0; c < nconss; ++c )
3078  {
3079  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
3080  if( !SCIPisGT(scip, consdata->violation, SCIPfeastol(scip)) )
3081  continue;
3082 
3083  nbndchg = 0;
3084  SCIP_CALL( propagateBounds(scip, conss[c], result, &nbndchg, &redundant) ); /*lint !e613*/
3085  assert(!redundant); /* constraint should not be violated and redundant simultaneously (unless solution is far out of bounds) */
3086  if( *result == SCIP_CUTOFF || *result == SCIP_REDUCEDDOM )
3087  {
3088  SCIPdebugMsg(scip, "enforced by %s\n", *result == SCIP_CUTOFF ? "cutting off node" : "reducing domain");
3089  return SCIP_OKAY;
3090  }
3091  }
3092 
3093  SCIPwarningMessage(scip, "could not enforce feasibility by separating or branching; declaring solution with viol %g feasible\n", SCIPconsGetData(maxviolcons)->violation);
3094  *result = SCIP_FEASIBLE;
3095 
3096  return SCIP_OKAY;
3097 }
3098 
3099 /*
3100  * Quadratic constraint upgrading
3101  */
3102 
3103 
3104 #ifdef QUADCONSUPGD_PRIORITY
3105 /** tries to upgrade a quadratic constraint to a SOC constraint
3106  *
3107  * We treat as special cases, quadratic with no bilinear terms and hyperbolic quadratic
3108  * constraints with exactly on bilinear component containing nonnegative variables. For this we use the formula:
3109  * \f[
3110  * x^T x \leq yz,\; y \geq 0,\; z \geq 0
3111  * \qquad\Leftrightarrow\qquad
3112  * \left\| \left(\begin{array}{c} x \\ \frac{1}{2}(y - z)\end{array}\right) \right\| \leq \frac{1}{2}(y + z).
3113  * \f]
3114  *
3115  * @todo implement more general hyperbolic upgrade, e.g., for -x^T x + yz >= 0 or x^T x <= ax + by + cyz
3116  */
3117 static
3118 SCIP_DECL_QUADCONSUPGD(upgradeConsQuadratic)
3120  int nquadvars;
3121  int nbilinterms;
3122  SCIP_QUADVARTERM* quadterms;
3123  SCIP_BILINTERM* bilinterm;
3124  SCIP_VAR** quadvars;
3125  SCIP_VAR** lhsvars;
3126  SCIP_Real* lhscoefs;
3127  SCIP_Real* lhsoffsets;
3128  SCIP_Real lhsconstant;
3129  int lhscount;
3130  int lhsnvars;
3131  SCIP_VAR* rhsvar;
3132  SCIP_Real rhscoef;
3133  SCIP_Real rhsoffset;
3134  SCIP_VAR* bilinvar1;
3135  SCIP_VAR* bilinvar2;
3136  SCIP_Real bilincoef;
3137  int i;
3138  int j;
3139  int negeigpos;
3140  SCIP_Real* a;
3141  SCIP_Real* bp;
3142  SCIP_Real* eigvals;
3143  SCIP_Bool infeas;
3144  SCIP_Bool success;
3145  SCIP_Bool trygeneralupg;
3146  int nneg;
3147  int npos;
3148  SCIP_Bool rhsvarfound;
3149  SCIP_Bool rhsissoc;
3150  SCIP_Bool lhsissoc;
3151  SCIP_Real rhsvarlb;
3152  SCIP_Real rhsvarub;
3153  SCIP_CONSHDLR* conshdlr;
3154  SCIP_CONSHDLRDATA* conshdlrdata;
3155 
3156  assert(scip != NULL);
3157  assert(cons != NULL);
3158  assert(nupgdconss != NULL);
3159  assert(upgdconss != NULL);
3160 
3161  *nupgdconss = 0;
3162 
3163  SCIPdebugMsg(scip, "upgradeConsQuadratic called for constraint <%s>\n", SCIPconsGetName(cons));
3164  SCIPdebugPrintCons(scip, cons, NULL);
3165 
3166  /* currently do not support linear parts in upgrading of SOC constraints; binary vars we can treat as if squared */
3167  if( ncontlin > 0 || nimpllin > 0 || nintlin > 0 )
3168  return SCIP_OKAY;
3169  assert(SCIPgetNLinearVarsQuadratic(scip, cons) == nbinlin);
3170 
3171  nbilinterms = SCIPgetNBilinTermsQuadratic(scip, cons);
3172  nquadvars = SCIPgetNQuadVarTermsQuadratic(scip, cons);
3173 
3174  /* a proper SOC constraint needs at least 2 variables (at least one will be quadratic)
3175  * but performance-wise that doesn't give a clear advantage on product(2), so let's even require 3 vars
3176  */
3177  if( nbinlin + nquadvars < 3 )
3178  return SCIP_OKAY;
3179 
3180  /* reserve space */
3181  SCIP_CALL( SCIPallocBufferArray(scip, &lhsvars, nbinlin + nquadvars - 1) );
3182  SCIP_CALL( SCIPallocBufferArray(scip, &lhscoefs, nbinlin + nquadvars - 1) );
3183  SCIP_CALL( SCIPallocBufferArray(scip, &lhsoffsets, nbinlin + nquadvars - 1) );
3184 
3185  /* initialize data */
3186  a = NULL;
3187  bp = NULL;
3188  eigvals = NULL;
3189  quadvars = NULL;
3190  bilinvar1 = NULL;
3191  bilinvar2 = NULL;
3192  lhsconstant = 0.0;
3193  lhscount = 0;
3194  rhsvar = NULL;
3195  rhscoef = 0.0;
3196  rhsoffset = 0.0;
3197 
3198  trygeneralupg = FALSE;
3199 
3200  /* if more than one bilinear term is present, we are in the general case */
3201  if( nbilinterms > 1 )
3202  {
3203  trygeneralupg = TRUE;
3204  goto GENERALUPG;
3205  }
3206 
3207  /* check hyperbolic part */
3208  if( nbilinterms == 1 && nbinlin == 0 )
3209  {
3210  bilinterm = SCIPgetBilinTermsQuadratic(scip, cons);
3211  bilinvar1 = bilinterm->var1;
3212  bilinvar2 = bilinterm->var2;
3213  bilincoef = bilinterm->coef;
3214 
3215  /* the variables in the bilinear term need to be nonnegative */
3216  if ( SCIPisNegative(scip, SCIPvarGetLbGlobal(bilinvar1)) || SCIPisNegative(scip, SCIPvarGetLbGlobal(bilinvar2)) )
3217  {
3218  trygeneralupg = TRUE;
3219  goto GENERALUPG;
3220  }
3221 
3222  /* we need a rhs */
3223  if ( ! SCIPisZero(scip, SCIPgetRhsQuadratic(scip, cons)) )
3224  {
3225  trygeneralupg = TRUE;
3226  goto GENERALUPG;
3227  }
3228 
3229  /* we only allow for -1.0 bilincoef */
3230  if ( ! SCIPisEQ(scip, bilincoef, -1.0) )
3231  {
3232  trygeneralupg = TRUE;
3233  goto GENERALUPG;
3234  }
3235 
3236  /* check that bilinear terms do not appear in the rest and quadratic terms have positive sqrcoef have no lincoef */
3237  quadterms = SCIPgetQuadVarTermsQuadratic(scip, cons);
3238  for (i = 0; i < nquadvars; ++i)
3239  {
3240  SCIP_QUADVARTERM* term = &quadterms[i];
3241 
3242  if( ! SCIPisZero(scip, term->lincoef) || SCIPisNegative(scip, term->sqrcoef) )
3243  {
3244  trygeneralupg = TRUE;
3245  goto GENERALUPG;
3246  }
3247 
3248  if ( (term->var == bilinvar1 || term->var == bilinvar2) && ! SCIPisZero(scip, term->sqrcoef) )
3249  {
3250  trygeneralupg = TRUE;
3251  goto GENERALUPG;
3252  }
3253  }
3254  }
3255 
3256  quadterms = SCIPgetQuadVarTermsQuadratic(scip, cons);
3257  assert( quadterms != NULL );
3258 
3259  if( !SCIPisInfinity(scip, SCIPgetRhsQuadratic(scip, cons)) )
3260  {
3261  /* try whether constraint on right hand side is SOC */
3262  lhsconstant = -SCIPgetRhsQuadratic(scip, cons);
3263 
3264  for( i = 0; i < nquadvars + nbinlin; ++i )
3265  {
3266  SCIP_Real sqrcoef;
3267  SCIP_Real lincoef;
3268  SCIP_VAR* var;
3269 
3270  if( i < nquadvars )
3271  {
3272  SCIP_QUADVARTERM* term = &quadterms[i];
3273 
3274  sqrcoef = term->sqrcoef;
3275  lincoef = term->lincoef;
3276  var = term->var;
3277  }
3278  else
3279  {
3280  sqrcoef = SCIPgetCoefsLinearVarsQuadratic(scip, cons)[i-nquadvars];
3281  lincoef = 0.0;
3282  var = SCIPgetLinearVarsQuadratic(scip, cons)[i-nquadvars];
3283  }
3284 
3285  /* skip terms with 0 coefficients */
3286  if( sqrcoef == 0.0 )
3287  continue;
3288 
3289  if( sqrcoef > 0.0 )
3290  {
3291  if( lhscount >= nbinlin + nquadvars - 1 )
3292  { /* too many variables on lhs, i.e., all variables seem to have positive coefficient */
3293  rhsvar = NULL;
3294  break;
3295  }
3296 
3297  lhsvars[lhscount] = var;
3298  lhscoefs[lhscount] = sqrt(sqrcoef);
3299 
3300  if( lincoef != 0.0 )
3301  {
3302  lhsoffsets[lhscount] = lincoef / (2.0 * sqrcoef);
3303  lhsconstant -= lincoef * lincoef / (4.0 * sqrcoef);
3304  }
3305  else
3306  {
3307  lhsoffsets[lhscount] = 0.0;
3308  }
3309 
3310  ++lhscount;
3311  }
3312  else if( rhsvar != NULL ||
3313  (SCIPisLT(scip, SCIPcomputeVarLbGlobal(scip, var), -lincoef / (2.0 * sqrcoef)) &&
3314  SCIPisGT(scip, SCIPcomputeVarUbGlobal(scip, var), -lincoef / (2.0 * sqrcoef))) )
3315  { /* second variable with negative coefficient -> cannot be SOC */
3316  /* or rhs side changes sign */
3317  rhsvar = NULL;
3318  break;
3319  }
3320  else if( SCIPvarIsBinary(var) )
3321  {
3322  /* binary variable on rhs? then we originally had a convex quadratic */
3323  break;
3324  }
3325  else
3326  {
3327  rhsvar = var;
3328  rhsoffset = lincoef / (2.0 * sqrcoef);
3329  lhsconstant -= lincoef * lincoef / (4.0 * sqrcoef);
3330 
3331  if( SCIPisGE(scip, SCIPcomputeVarLbGlobal(scip, var), -lincoef / (2.0 * sqrcoef)) )
3332  rhscoef = sqrt(-sqrcoef);
3333  else
3334  rhscoef = -sqrt(-sqrcoef);
3335  }
3336  }
3337  }
3338 
3339  /* treat hyberbolic case */
3340  if( nbilinterms == 1 && nbinlin == 0 )
3341  {
3342  char name[SCIP_MAXSTRLEN];
3343  SCIP_VAR* auxvarsum;
3344  SCIP_VAR* auxvardiff;
3345  SCIP_CONS* couplingcons;
3346  SCIP_VAR* consvars[3];
3347  SCIP_Real consvals[3];
3348 
3349  /* can only upgrade if rhs is 0 */
3350  if ( rhsvar != NULL )
3351  goto cleanup;
3352 
3353  SCIPdebugMsg(scip, "found hyberbolic quadratic constraint <%s> to be SOC\n", SCIPconsGetName(cons));
3354 
3355  /* check if upgdconss is long enough to store upgrade constraints: we need two if we will have a quadratic
3356  * constraint for the left hand side left */
3357  *nupgdconss = SCIPisInfinity(scip, -SCIPgetLhsQuadratic(scip, cons)) ? 1 : 2;
3358  if ( *nupgdconss > upgdconsssize )
3359  {
3360  /* signal that we need more memory and return */
3361  *nupgdconss = -*nupgdconss;
3362  goto cleanup;
3363  }
3364 
3365  /* create auxiliary variable for sum (nonnegative) */
3366  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "soc#%s_s", SCIPconsGetName(cons));
3367  SCIP_CALL( SCIPcreateVar(scip, &auxvarsum, name, 0.0, SCIPinfinity(scip), 0.0,
3369  SCIP_CALL( SCIPaddVar(scip, auxvarsum) );
3370 
3371  /* create auxiliary variable for difference (free) */
3372  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "soc#%s_d", SCIPconsGetName(cons));
3373  SCIP_CALL( SCIPcreateVar(scip, &auxvardiff, name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0,
3375  SCIP_CALL( SCIPaddVar(scip, auxvardiff) );
3376 
3377  /* add coupling constraint for sum */
3378  consvars[0] = bilinvar1;
3379  consvars[1] = bilinvar2;
3380  consvars[2] = auxvarsum;
3381  consvals[0] = 1.0;
3382  consvals[1] = 1.0;
3383  consvals[2] = -1.0;
3384 
3385  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "soc#%s_cs", SCIPconsGetName(cons));
3386  SCIP_CALL( SCIPcreateConsLinear(scip, &couplingcons, name, 3, consvars, consvals, 0.0, 0.0,
3390  SCIPconsIsStickingAtNode(cons)) );
3391  SCIP_CALL( SCIPaddCons(scip, couplingcons) );
3392  SCIP_CALL( SCIPreleaseCons(scip, &couplingcons) );
3393 
3394  /* add coupling constraint for difference */
3395  consvars[0] = bilinvar1;
3396  consvars[1] = bilinvar2;
3397  consvars[2] = auxvardiff;
3398  consvals[0] = 1.0;
3399  consvals[1] = -1.0;
3400  consvals[2] = -1.0;
3401 
3402  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "soc#%s_cd", SCIPconsGetName(cons));
3403  SCIP_CALL( SCIPcreateConsLinear(scip, &couplingcons, name, 3, consvars, consvals, 0.0, 0.0,
3407  SCIPconsIsStickingAtNode(cons)) );
3408  SCIP_CALL( SCIPaddCons(scip, couplingcons) );
3409  SCIP_CALL( SCIPreleaseCons(scip, &couplingcons) );
3410 
3411  /* add difference variable to constraint */
3412  lhsvars[lhscount] = auxvardiff;
3413  lhscoefs[lhscount] = 0.5;
3414  lhsoffsets[lhscount] = 0.0;
3415 
3416  SCIP_CALL( SCIPcreateConsSOC(scip, &upgdconss[0], SCIPconsGetName(cons),
3417  lhscount + 1, lhsvars, lhscoefs, lhsoffsets, MAX(lhsconstant, 0.0),
3418  auxvarsum, 0.5, 0.0,
3422  SCIPdebugPrintCons(scip, upgdconss[0], NULL);
3423 
3424  /* create constraint that is equal to cons except that rhs is now infinity */
3425  if( !SCIPisInfinity(scip, -SCIPgetLhsQuadratic(scip, cons)) )
3426  {
3427  SCIP_CALL( SCIPcreateConsQuadratic2(scip, &upgdconss[1], SCIPconsGetName(cons),
3431  SCIPgetLhsQuadratic(scip, cons), SCIPinfinity(scip),
3435  }
3436  SCIP_CALL( SCIPreleaseVar(scip, &auxvardiff) );
3437  SCIP_CALL( SCIPreleaseVar(scip, &auxvarsum) );
3438 
3439  for (i = 0; i < lhscount; ++i)
3440  {
3441  SCIP_CALL( SCIPmarkDoNotMultaggrVar(scip, lhsvars[i]) );
3442  }
3443 
3444  goto cleanup;
3445  }
3446 
3447  if( rhsvar != NULL && lhscount >= 1 && !SCIPisNegative(scip, lhsconstant) )
3448  { /* found SOC constraint, so upgrade to SOC constraint(s) (below) and relax right hand side */
3449  SCIPdebugMsg(scip, "found right hand side of constraint <%s> to be SOC\n", SCIPconsGetName(cons));
3450 
3451  /* check if upgdconss is long enough to store upgrade constraints:
3452  * we need two if we will have a quadratic constraint for the left hand side left */
3453  *nupgdconss = SCIPisInfinity(scip, -SCIPgetLhsQuadratic(scip, cons)) ? 1 : 2;
3454  if( *nupgdconss > upgdconsssize )
3455  {
3456  /* signal that we need more memory and return */
3457  *nupgdconss = -*nupgdconss;
3458  goto cleanup;
3459  }
3460 
3461  SCIP_CALL( SCIPcreateConsSOC(scip, &upgdconss[0], SCIPconsGetName(cons),
3462  lhscount, lhsvars, lhscoefs, lhsoffsets, MAX(lhsconstant, 0.0),
3463  rhsvar, rhscoef, rhsoffset,
3467  SCIPdebugPrintCons(scip, upgdconss[0], NULL);
3468 
3469  /* create constraint that is equal to cons except that rhs is now infinity */
3470  if( !SCIPisInfinity(scip, -SCIPgetLhsQuadratic(scip, cons)) )
3471  {
3472  SCIP_CALL( SCIPcreateConsQuadratic2(scip, &upgdconss[1], SCIPconsGetName(cons),
3476  SCIPgetLhsQuadratic(scip, cons), SCIPinfinity(scip),
3480  }
3481  }
3482  else if( !SCIPisInfinity(scip, - SCIPgetLhsQuadratic(scip, cons)) )
3483  { /* if the first failed, try if constraint on left hand side is SOC (using negated coefficients) */
3484  lhscount = 0;
3485  rhsvar = NULL;
3486 
3487  lhsconstant = SCIPgetLhsQuadratic(scip, cons);
3488 
3489  for( i = 0; i < nquadvars + nbinlin; ++i )
3490  {
3491  SCIP_Real sqrcoef;
3492  SCIP_Real lincoef;
3493  SCIP_VAR* var;
3494 
3495  if( i < nquadvars )
3496  {
3497  SCIP_QUADVARTERM* term = &quadterms[i];
3498 
3499  sqrcoef = term->sqrcoef;
3500  lincoef = term->lincoef;
3501  var = term->var;
3502  }
3503  else
3504  {
3505  sqrcoef = SCIPgetCoefsLinearVarsQuadratic(scip, cons)[i-nquadvars];
3506  lincoef = 0.0;
3507  var = SCIPgetLinearVarsQuadratic(scip, cons)[i-nquadvars];
3508  }
3509 
3510  /* if there is a linear variable that is still considered as quadratic (constraint probably not presolved yet),
3511  * then give up
3512  */
3513  if( sqrcoef == 0.0 )
3514  goto cleanup;
3515 
3516  if( sqrcoef < 0.0 )
3517  {
3518  if( lhscount >= nquadvars + nbinlin - 1 )
3519  { /* too many variables on lhs, i.e., all variables seem to have negative coefficient */
3520  rhsvar = NULL;
3521  break;
3522  }
3523 
3524  lhsvars[lhscount] = var;
3525  lhscoefs[lhscount] = sqrt(-sqrcoef);
3526 
3527  if( lincoef != 0.0 )
3528  {
3529  lhsoffsets[lhscount] = lincoef / (2.0 * sqrcoef);
3530  lhsconstant += lincoef * lincoef / (4.0 * sqrcoef);
3531  }
3532  else
3533  {
3534  lhsoffsets[lhscount] = 0.0;
3535  }
3536 
3537  ++lhscount;
3538  }
3539  else if( rhsvar != NULL ||
3540  (SCIPisLT(scip, SCIPcomputeVarLbGlobal(scip, var), -lincoef / (2.0 * sqrcoef))
3541  && SCIPisGT(scip, SCIPcomputeVarUbGlobal(scip, var), -lincoef / (2.0 * sqrcoef))) )
3542  { /* second variable with positive coefficient -> cannot be SOC */
3543  /* or rhs side changes sign */
3544  rhsvar = NULL;
3545  break;
3546  }
3547  else if( SCIPvarIsBinary(var) )
3548  {
3549  /* binary variable on rhs? then we originally had a convex quadratic */
3550  break;
3551  }
3552  else
3553  {
3554  rhsvar = var;
3555  rhsoffset = lincoef / (2.0 * sqrcoef);
3556  lhsconstant += lincoef * lincoef / (4.0 * sqrcoef);
3557 
3558  if( SCIPisGE(scip, SCIPcomputeVarLbGlobal(scip, var), -lincoef / (2.0 * sqrcoef)) )
3559  rhscoef = sqrt(sqrcoef);
3560  else
3561  rhscoef = -sqrt(sqrcoef);
3562  }
3563  }
3564 
3565  if( rhsvar != NULL && lhscount >= 1 && !SCIPisNegative(scip, lhsconstant) )
3566  { /* found SOC constraint, so upgrade to SOC constraint(s) (below) and relax left hand side */
3567  SCIPdebugMsg(scip, "found left hand side of constraint <%s> to be SOC\n", SCIPconsGetName(cons));
3568 
3569  /* check if upgdconss is long enough to store upgrade constraints:
3570  * we need two if we will have a quadratic constraint for the right hand side left */
3571  *nupgdconss = SCIPisInfinity(scip, SCIPgetRhsQuadratic(scip, cons)) ? 1 : 2;
3572  if( *nupgdconss > upgdconsssize )
3573  {
3574  /* signal that we need more memory and return */
3575  *nupgdconss = -*nupgdconss;
3576  goto cleanup;
3577  }
3578 
3579  SCIP_CALL( SCIPcreateConsSOC(scip, &upgdconss[0], SCIPconsGetName(cons),
3580  lhscount, lhsvars, lhscoefs, lhsoffsets, MAX(lhsconstant, 0.0),
3581  rhsvar, rhscoef, rhsoffset,
3585  SCIPdebugPrintCons(scip, upgdconss[0], NULL);
3586 
3587  /* create constraint that is equal to cons except that lhs is now -infinity */
3588  if( !SCIPisInfinity(scip, SCIPgetRhsQuadratic(scip, cons)) )
3589  {
3590  SCIP_CALL( SCIPcreateConsQuadratic2(scip, &upgdconss[1], SCIPconsGetName(cons),
3594  -SCIPinfinity(scip), SCIPgetRhsQuadratic(scip, cons),
3598  }
3599  }
3600  }
3601 
3602 GENERALUPG:
3603  if( !trygeneralupg )
3604  goto cleanup;
3605 
3606  /* find the soc constraint handler */
3607  conshdlr = SCIPfindConshdlr(scip, CONSHDLR_NAME);
3608  assert(conshdlr != NULL);
3609  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3610  assert(conshdlrdata != NULL);
3611 
3612  if( !conshdlrdata->generalsocupg )
3613  goto cleanup;
3614 
3615  SCIPdebugMsg(scip, "Trying general method of upgrade to a soc const\n");
3616 
3617  rhsvarlb = 1.0;
3618  rhsvarub = 0.0;
3619 
3620 #ifndef SCIP_STATISTIC
3621  /* skip large matrices (needs to compute eigenvalues/vectors) according to presolve timing */
3622  if( (presoltiming & SCIP_PRESOLTIMING_FAST) != 0 && nquadvars > 10 )
3623  goto cleanup;
3624  if( (presoltiming & SCIP_PRESOLTIMING_MEDIUM) != 0 && nquadvars > 50 )
3625  goto cleanup;
3626 #endif
3627 
3628  /* need routine to compute eigenvalues/eigenvectors */
3629  if( !SCIPisIpoptAvailableIpopt() )
3630  goto cleanup;
3631 
3632  /* build lower triangular part the A matrix */
3633  SCIP_CALL( SCIPallocClearBufferArray(scip, &a, nquadvars*nquadvars) ); /*lint !e647*/
3634  SCIP_CALL( SCIPallocClearBufferArray(scip, &bp, nquadvars) );
3635  SCIP_CALL( SCIPallocBufferArray(scip, &quadvars, nquadvars) );
3636  SCIP_CALL( SCIPallocBufferArray(scip, &eigvals, nquadvars) );
3637 
3638  /* make sure quadratic variables terms are sorted */
3640 
3641  /* set lower triangular entries of A corresponding to square terms */
3642  for( i = 0; i < nquadvars; ++i )
3643  {
3644  SCIP_QUADVARTERM* term = &SCIPgetQuadVarTermsQuadratic(scip, cons)[i];
3645 
3646  /* skip upgrade if fixed (or (multi)aggregated) variables and still in fast presolving
3647  * probably cons_quadratic did not yet had the chance to remove/replace this variable (see also #2352)
3648  */
3649  if( !SCIPvarIsActive(term->var) && (presoltiming & SCIP_PRESOLTIMING_FAST) != 0 )
3650  goto cleanup;
3651 
3652  a[i*nquadvars + i] = term->sqrcoef;
3653  quadvars[i] = term->var;
3654  }
3655 
3656  /* set lower triangular entries of A corresponding to bilinear terms */
3657  for( i = 0; i < nbilinterms; ++i )
3658  {
3659  int idx1;
3660  int idx2;
3661 
3662  bilinterm = &SCIPgetBilinTermsQuadratic(scip, cons)[i];
3663 
3664  SCIP_CALL( SCIPfindQuadVarTermQuadratic(scip, cons, bilinterm->var1, &idx1) );
3665  SCIP_CALL( SCIPfindQuadVarTermQuadratic(scip, cons, bilinterm->var2, &idx2) );
3666 
3667  assert(idx1 >= 0);
3668  assert(idx2 >= 0);
3669  assert(idx1 != idx2);
3670 
3671  a[MIN(idx1,idx2) * nquadvars + MAX(idx1,idx2)] = bilinterm->coef / 2.0;
3672  }
3673 
3674  /* compute eigenvalues and vectors, A = PDP^t
3675  * note: a stores P^t
3676  */
3677  if( LapackDsyev(TRUE, nquadvars, a, eigvals) != SCIP_OKAY )
3678  {
3679  SCIPdebugMsg(scip, "Failed to compute eigenvalues and eigenvectors for constraint <%s>.\n", SCIPconsGetName(cons));
3680  goto cleanup;
3681  }
3682 
3683  /* set small eigenvalues to 0 and compute b*P */
3684  nneg = 0;
3685  npos = 0;
3686  for( i = 0; i < nquadvars; ++i )
3687  {
3688  for( j = 0; j < nquadvars; ++j )
3689  {
3690  SCIP_QUADVARTERM* term = &SCIPgetQuadVarTermsQuadratic(scip, cons)[j];
3691  bp[i] += term->lincoef * a[i * nquadvars + j];
3692  }
3693 
3694  if( SCIPisZero(scip, eigvals[i]) )
3695  {
3696  /* if there is a purely linear variable, the constraint can't be written as a SOC */
3697  if( !SCIPisZero(scip, bp[i]) )
3698  {
3699  goto cleanup;
3700  }
3701 
3702  bp[i] = 0.0;
3703  eigvals[i] = 0.0;
3704  }
3705  else if( eigvals[i] > 0.0 )
3706  npos++;
3707  else
3708  nneg++;
3709  }
3710 
3711  /* a proper SOC constraint needs at least 2 variables
3712  * let's make sure we have at least 3, though, as this upgrade comes with extra (multiaggr.) vars
3713  */
3714  if( npos + nneg < 3 )
3715  goto cleanup;
3716 
3717  /* determine whether rhs or lhs of cons is potentially SOC, if any */
3718  rhsissoc = nneg == 1 && !SCIPisInfinity(scip, SCIPgetRhsQuadratic(scip, cons));
3719  lhsissoc = npos == 1 && !SCIPisInfinity(scip, -SCIPgetLhsQuadratic(scip, cons));
3720 
3721  /* if non is potentially SOC, stop */
3722  if( !rhsissoc && !lhsissoc )
3723  goto cleanup;
3724 
3725  if( rhsissoc && lhsissoc )
3726  {
3727  /* only handle rhs-SOC here for now
3728  * if the upgrade is run again on the remaining lhs-SOC, it would upgrade that side
3729  */
3730  lhsissoc = FALSE;
3731  }
3732 
3733  if( lhsissoc )
3734  {
3735  /* lhs is potentially SOC, change signs */
3736  lhsconstant = SCIPgetLhsQuadratic(scip, cons);
3737 
3738  for( i = 0; i < nquadvars; ++i )
3739  {
3740  eigvals[i] = -eigvals[i];
3741  bp[i] = -bp[i];
3742  }
3743  }
3744  else
3745  {
3746  lhsconstant = -SCIPgetRhsQuadratic(scip, cons);
3747  }
3748 
3749  /* we have lhsconstant + x^t A x + b x <= 0 and A has a single negative eigenvalue; try to build soc */
3750  negeigpos = -1;
3751  rhsvarfound = FALSE;
3752  for( i = 0; i < nquadvars; ++i )
3753  {
3754  if( SCIPisZero(scip, eigvals[i]) )
3755  continue;
3756 
3757  if( eigvals[i] > 0.0 )
3758  {
3759  lhscoefs[lhscount] = sqrt(eigvals[i]) * UPGSCALE;
3760  lhsoffsets[lhscount] = bp[i] / (2 * eigvals[i]);
3761  lhsconstant -= bp[i] * bp[i] / (4 * eigvals[i]);
3762 
3763  ++lhscount;
3764  }
3765  else
3766  {
3767  /* should enter here only once */
3768  assert(rhsvarfound == FALSE);
3769 
3770  rhsoffset = bp[i] / (2 * eigvals[i]);
3771 
3772  /* the constraint can only be a soc if the resulting rhs var does not change var; the rhs var is going to be a
3773  * multiaggregated variable, so estimate its bounds
3774  */
3775  rhsvarlb = 0.0;
3776  rhsvarub = 0.0;
3777  for( j = 0; j < nquadvars; ++j )
3778  {
3779  SCIP_Real aux;
3780 
3781  if( SCIPisZero(scip, a[i * nquadvars + j]) )
3782  continue;
3783 
3784  if( a[i * nquadvars + j] > 0.0 )
3785  {
3786  aux = SCIPcomputeVarLbGlobal(scip, quadvars[j]);
3787  assert(!SCIPisInfinity(scip, aux));
3788  }
3789  else
3790  {
3791  aux = SCIPcomputeVarUbGlobal(scip, quadvars[j]);
3792  assert(!SCIPisInfinity(scip, -aux));
3793  }
3794 
3795  if( SCIPisInfinity(scip, aux) || SCIPisInfinity(scip, -aux) )
3796  {
3797  rhsvarlb = -SCIPinfinity(scip);
3798  break;
3799  }
3800  else
3801  rhsvarlb += a[i * nquadvars + j] * aux;
3802  }
3803  rhsvarlb += rhsoffset;
3804 
3805  for( j = 0; j < nquadvars; ++j )
3806  {
3807  SCIP_Real aux;
3808 
3809  if( SCIPisZero(scip, a[i * nquadvars + j]) )
3810  continue;
3811 
3812  if( a[i * nquadvars + j] > 0.0 )
3813  {
3814  aux = SCIPcomputeVarUbGlobal(scip, quadvars[j]);
3815  assert(!SCIPisInfinity(scip, -aux));
3816  }
3817  else
3818  {
3819  aux = SCIPcomputeVarLbGlobal(scip, quadvars[j]);
3820  assert(!SCIPisInfinity(scip, aux));
3821  }
3822 
3823  if( SCIPisInfinity(scip, aux) || SCIPisInfinity(scip, -aux) )
3824  {
3825  rhsvarub = SCIPinfinity(scip);
3826  break;
3827  }
3828  else
3829  rhsvarub += a[i * nquadvars + j] * aux;
3830  }
3831  rhsvarub += rhsoffset;
3832 
3833  /* since we are just interested in obtaining an interval that contains the real bounds and is tight enough so
3834  * that we can identify that the rhsvar does not change sign, we swap the bounds in case of numerical troubles
3835  */
3836  if( rhsvarub < rhsvarlb )
3837  {
3838  assert(SCIPisEQ(scip, rhsvarub, rhsvarlb));
3839  SCIPswapReals(&rhsvarub, &rhsvarlb);
3840  }
3841 
3842  /* check whether rhsvar changes sign */
3843  if( SCIPisGE(scip, rhsvarlb, 0.0) || SCIPisLE(scip, rhsvarub, 0.0) )
3844  {
3845  rhsvarfound = TRUE;
3846  negeigpos = i;
3847  lhsconstant -= rhsoffset * rhsoffset * eigvals[i];
3848  rhscoef = SCIPisLE(scip, rhsvarub, 0.0) ? -sqrt(-eigvals[i]) : sqrt(-eigvals[i]);
3849  }
3850  else
3851  {
3852  SCIPdebugMsg(scip, "Failed because rhsvar [%g, %g] changes sign.\n", rhsvarlb, rhsvarub);
3853  rhsvarfound = FALSE;
3854  break;
3855  }
3856  }
3857  }
3858 
3859  if( rhsvarfound && lhscount >= 1 && !SCIPisNegative(scip, lhsconstant) )
3860  {
3861  /* found SOC constraint, so upgrade to SOC constraint(s) (below) and relax right hand side */
3862  SCIPdebugMsg(scip, "found right hand side of constraint <%s> to be SOC\n", SCIPconsGetName(cons));
3863 
3864  /* check if upgdconss is long enough to store upgrade constraints:
3865  * we need two if we will have a quadratic constraint for the left hand side left */
3866  if( rhsissoc )
3867  *nupgdconss = SCIPisInfinity(scip, -SCIPgetLhsQuadratic(scip, cons)) ? 1 : 2;
3868  else
3869  *nupgdconss = SCIPisInfinity(scip, SCIPgetRhsQuadratic(scip, cons)) ? 1 : 2;
3870  if( *nupgdconss > upgdconsssize )
3871  {
3872  /* signal that we need more memory and return */
3873  *nupgdconss = -*nupgdconss;
3874  goto cleanup;
3875  }
3876 
3877  /* create lhs and rhs vars */
3878  lhsnvars = 0;
3879  for( i = 0; i < nquadvars; ++i )
3880  {
3881  if( eigvals[i] <= 0.0 )
3882  continue;
3883 
3884  SCIP_CALL( SCIPcreateVarBasic(scip, &lhsvars[lhsnvars], NULL,
3885  -SCIPinfinity(scip), SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
3886  SCIP_CALL( SCIPaddVar(scip, lhsvars[lhsnvars]) );
3887 
3888  SCIP_CALL( SCIPmultiaggregateVar(scip, lhsvars[lhsnvars], nquadvars, quadvars, &(a[i * nquadvars]),
3889  lhsoffsets[lhsnvars], &infeas, &success) );
3890 
3891  if( infeas || !success )
3892  {
3893  SCIPdebugMsg(scip, "Problem with aggregation while trying to upgrade <%s>.\n", SCIPconsGetName(cons) );
3894 
3895  /* release all created vars so far */
3896  for( j = 0; j <= lhsnvars; ++j )
3897  SCIP_CALL( SCIPreleaseVar(scip, &lhsvars[j]) );
3898 
3899  *nupgdconss = 0;
3900  goto cleanup;
3901  }
3902  lhsnvars++;
3903  }
3904  assert(lhsnvars == lhscount);
3905  assert(negeigpos >= 0);
3906 
3907  SCIP_CALL( SCIPcreateVarBasic(scip, &rhsvar, NULL,
3908  rhsvarlb, rhsvarub, 0.0, SCIP_VARTYPE_CONTINUOUS) );
3909  SCIP_CALL( SCIPaddVar(scip, rhsvar) );
3910  SCIP_CALL( SCIPmultiaggregateVar(scip, rhsvar, nquadvars, quadvars, &(a[negeigpos * nquadvars]),
3911  rhsoffset, &infeas, &success) );
3912 
3913  if( infeas || !success )
3914  {
3915  SCIPdebugMsg(scip, "Problem with aggregation while trying to upgrade <%s>.\n", SCIPconsGetName(cons) );
3916 
3917  /* release all created vars */
3918  SCIP_CALL( SCIPreleaseVar(scip, &rhsvar) );
3919  for( j = 0; j < lhsnvars; ++j )
3920  {
3921  SCIP_CALL( SCIPreleaseVar(scip, &lhsvars[j]) );
3922  }
3923 
3924  *nupgdconss = 0;
3925  goto cleanup;
3926  }
3927 
3928  SCIP_CALL( SCIPcreateConsSOC(scip, &upgdconss[0], SCIPconsGetName(cons),
3929  lhscount, lhsvars, lhscoefs, NULL, MAX(lhsconstant, 0.0) * UPGSCALE * UPGSCALE,
3930  rhsvar, rhscoef * UPGSCALE, 0.0,
3934  SCIPdebugPrintCons(scip, upgdconss[0], NULL);
3935 
3936  /* release created variables */
3937  SCIP_CALL( SCIPreleaseVar(scip, &rhsvar) );
3938  for( i = 0; i < lhsnvars; ++i )
3939  {
3940  SCIP_CALL( SCIPreleaseVar(scip, &lhsvars[i]) );
3941  }
3942 
3943  /* create constraint that is equal to cons except that rhs/lhs is now +/-infinity */
3944  if( rhsissoc && !SCIPisInfinity(scip, -SCIPgetLhsQuadratic(scip, cons)) )
3945  {
3946  SCIPdebugMsg(scip, "rhs is soc, keep quadratic\n");
3947  SCIP_CALL( SCIPcreateConsQuadratic2(scip, &upgdconss[1], SCIPconsGetName(cons),
3951  SCIPgetLhsQuadratic(scip, cons), SCIPinfinity(scip),
3955  }
3956  else if( lhsissoc && !SCIPisInfinity(scip, SCIPgetRhsQuadratic(scip,cons)) )
3957  {
3958  SCIPdebugMsg(scip, "lhs is soc, keep quadratic\n");
3959  SCIP_CALL( SCIPcreateConsQuadratic2(scip, &upgdconss[1], SCIPconsGetName(cons),
3963  -SCIPinfinity(scip), SCIPgetRhsQuadratic(scip, cons),
3967  }
3968  }
3969 #ifdef SCIP_DEBUG
3970  else
3971  {
3972  if( lhscount < 1 )
3973  {
3974  SCIPdebugMsg(scip, "Failed because there are not enough lhsvars (%d)\n", lhscount);
3975  }
3976  if( SCIPisNegative(scip, lhsconstant) )
3977  {
3978  SCIPdebugMsg(scip, "Failed because lhsconstant is negative (%g)\n", lhsconstant);
3979  }
3980  }
3981 #endif
3982 
3983  cleanup:
3984  SCIPfreeBufferArrayNull(scip, &eigvals);
3985  SCIPfreeBufferArrayNull(scip, &quadvars);
3986  SCIPfreeBufferArrayNull(scip, &bp);
3987  SCIPfreeBufferArrayNull(scip, &a);
3988  SCIPfreeBufferArray(scip, &lhsoffsets);
3989  SCIPfreeBufferArray(scip, &lhscoefs);
3990  SCIPfreeBufferArray(scip, &lhsvars);
3991 
3992  return SCIP_OKAY;
3993 } /*lint !e715*/
3994 #endif
3995 
3996 /*
3997  * Callback methods of constraint handler
3998  */
3999 
4000 /** copy method for constraint handler plugins (called when SCIP copies plugins) */
4001 static
4002 SCIP_DECL_CONSHDLRCOPY(conshdlrCopySOC)
4003 { /*lint --e{715}*/
4004  assert(scip != NULL);
4005  assert(conshdlr != NULL);
4006  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
4007 
4008  /* call inclusion method of constraint handler */
4010 
4011  *valid = TRUE;
4012 
4013  return SCIP_OKAY;
4014 }
4015 
4016 /** destructor of constraint handler to free constraint handler data (called when SCIP is exiting) */
4017 static
4018 SCIP_DECL_CONSFREE(consFreeSOC)
4019 {
4020  SCIP_CONSHDLRDATA* conshdlrdata;
4021 
4022  assert( scip != NULL );
4023  assert( conshdlr != NULL );
4024  assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
4025 
4026  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4027  assert(conshdlrdata != NULL);
4028 
4029  SCIPfreeBlockMemory(scip, &conshdlrdata);
4030 
4031  return SCIP_OKAY;
4032 }
4033 
4034 
4035 /** initialization method of constraint handler (called after problem was transformed) */
4036 static
4037 SCIP_DECL_CONSINIT(consInitSOC)
4038 { /*lint --e{715}*/
4039  SCIP_CONSHDLRDATA* conshdlrdata;
4040  int c;
4041 
4042  assert(scip != NULL);
4043  assert(conshdlr != NULL);
4044  assert(conss != NULL || nconss == 0);
4045 
4046  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4047  assert(conshdlrdata != NULL);
4048 
4049  conshdlrdata->subnlpheur = SCIPfindHeur(scip, "subnlp");
4050  conshdlrdata->trysolheur = SCIPfindHeur(scip, "trysol");
4051  conshdlrdata->haveexprint = (strcmp(SCIPexprintGetName(), "NONE") != 0);
4052 
4053  /* mark constraints for propagation */
4054  for( c = 0; c < nconss; ++c )
4055  {
4056  SCIP_CALL( SCIPmarkConsPropagate(scip, conss[c]) ); /*lint !e613*/
4057  }
4058 
4059  return SCIP_OKAY;
4060 }
4061 
4062 
4063 /** deinitialization method of constraint handler (called before transformed problem is freed) */
4064 static
4065 SCIP_DECL_CONSEXIT(consExitSOC)
4066 { /*lint --e{715}*/
4067  SCIP_CONSHDLRDATA* conshdlrdata;
4068 
4069  assert(scip != NULL);
4070  assert(conshdlr != NULL);
4071 
4072  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4073  assert(conshdlrdata != NULL);
4074 
4075  conshdlrdata->subnlpheur = NULL;
4076  conshdlrdata->trysolheur = NULL;
4077  conshdlrdata->haveexprint = FALSE;
4078 
4079  return SCIP_OKAY;
4080 }
4081 
4082 /** presolving deinitialization method of constraint handler (called after presolving has been finished) */
4083 static
4084 SCIP_DECL_CONSEXITPRE(consExitpreSOC)
4085 { /*lint --e{715}*/
4086  int c;
4087 
4088  assert(scip != NULL);
4089  assert(conshdlr != NULL);
4090  assert(conss != NULL || nconss == 0);
4091 
4092  /* tell SCIP that we have something nonlinear */
4093  for( c = 0; c < nconss; ++c )
4094  {
4095  if( SCIPconsIsAdded(conss[c]) ) /*lint !e613*/
4096  {
4097  SCIPenableNLP(scip);
4098  break;
4099  }
4100  }
4101 
4102  return SCIP_OKAY;
4103 }
4104 
4105 
4106 /** solving process initialization method of constraint handler (called when branch and bound process is about to begin) */
4107 static
4108 SCIP_DECL_CONSINITSOL(consInitsolSOC)
4109 {
4110  SCIP_CONSHDLRDATA* conshdlrdata;
4111  SCIP_CONSDATA* consdata;
4112  int c;
4113 
4114  assert(scip != NULL);
4115  assert(conshdlr != NULL);
4116 
4117  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4118  assert(conshdlrdata != NULL);
4119  assert(conshdlrdata->eventhdlr != NULL);
4120 
4121  /* add nlrow representation to NLP, if NLP has been enabled */
4122  if( SCIPisNLPConstructed(scip) )
4123  {
4124  for( c = 0; c < nconss; ++c )
4125  {
4126  if( SCIPconsIsEnabled(conss[c]) )
4127  {
4128  consdata = SCIPconsGetData(conss[c]);
4129  assert(consdata != NULL);
4130 
4131  if( consdata->nlrow == NULL )
4132  {
4133  SCIP_CALL( createNlRow(scip, conshdlr, conss[c]) );
4134  assert(consdata->nlrow != NULL);
4135  }
4136  SCIP_CALL( SCIPaddNlRow(scip, consdata->nlrow) );
4137  }
4138  }
4139  }
4140 
4141  conshdlrdata->newsoleventfilterpos = -1;
4142  if( nconss != 0 )
4143  {
4144  SCIP_EVENTHDLR* eventhdlr;
4145 
4146  eventhdlr = SCIPfindEventhdlr(scip, CONSHDLR_NAME"_newsolution");
4147  assert(eventhdlr != NULL);
4148 
4149  SCIP_CALL( SCIPcatchEvent(scip, SCIP_EVENTTYPE_SOLFOUND, eventhdlr, (SCIP_EVENTDATA*)conshdlr, &conshdlrdata->newsoleventfilterpos) );
4150  }
4151 
4152  /* reset flags and counters */
4153  conshdlrdata->sepanlp = FALSE;
4154  conshdlrdata->lastenfonode = NULL;
4155  conshdlrdata->nenforounds = 0;
4156 
4157  return SCIP_OKAY;
4158 }
4159 
4160 
4161 /** solving process deinitialization method of constraint handler (called before branch and bound process data is freed) */
4162 static
4163 SCIP_DECL_CONSEXITSOL(consExitsolSOC)
4165  SCIP_CONSHDLRDATA* conshdlrdata;
4166  SCIP_CONSDATA* consdata;
4167  int c;
4168 
4169  assert(scip != NULL);
4170  assert(conshdlr != NULL);
4171 
4172  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4173  assert(conshdlrdata != NULL);
4174  assert(conshdlrdata->eventhdlr != NULL);
4175 
4176  if( conshdlrdata->newsoleventfilterpos >= 0 )
4177  {
4178  SCIP_EVENTHDLR* eventhdlr;
4179 
4180  eventhdlr = SCIPfindEventhdlr(scip, CONSHDLR_NAME"_newsolution");
4181  assert(eventhdlr != NULL);
4182 
4183  SCIP_CALL( SCIPdropEvent(scip, SCIP_EVENTTYPE_SOLFOUND, eventhdlr, (SCIP_EVENTDATA*)conshdlr, conshdlrdata->newsoleventfilterpos) );
4184  conshdlrdata->newsoleventfilterpos = -1;
4185  }
4186 
4187  for( c = 0; c < nconss; ++c )
4188  {
4189  consdata = SCIPconsGetData(conss[c]);
4190  assert(consdata != NULL);
4191 
4192  /* free nonlinear row representation */
4193  if( consdata->nlrow != NULL )
4194  {
4195  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
4196  }
4197  }
4198 
4199  return SCIP_OKAY;
4200 } /*lint !e715*/
4201 
4202 
4203 /** frees specific constraint data */
4204 static
4205 SCIP_DECL_CONSDELETE(consDeleteSOC)
4207  int i;
4208 
4209  assert(scip != NULL);
4210  assert(conshdlr != NULL);
4211  assert(cons != NULL);
4212  assert(consdata != NULL);
4213  assert(*consdata != NULL);
4214  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
4215 
4216  SCIPdebugMsg(scip, "Deleting SOC constraint <%s>.\n", SCIPconsGetName(cons) );
4217 
4218  if( SCIPconsIsTransformed(cons) )
4219  {
4220  SCIP_CONSHDLRDATA* conshdlrdata;
4221 
4222  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4223  assert(conshdlrdata != NULL);
4224 
4225  SCIP_CALL( dropVarEvents(scip, conshdlrdata->eventhdlr, cons) );
4226  }
4227 
4228  for( i = 0; i < (*consdata)->nvars; ++i )
4229  {
4230  SCIP_CALL( SCIPreleaseVar(scip, &(*consdata)->vars[i]) );
4231  }
4232 
4233  SCIPfreeBlockMemoryArray(scip, &(*consdata)->vars, (*consdata)->nvars);
4234  SCIPfreeBlockMemoryArray(scip, &(*consdata)->coefs, (*consdata)->nvars);
4235  SCIPfreeBlockMemoryArray(scip, &(*consdata)->offsets, (*consdata)->nvars);
4236  assert((*consdata)->lhsbndchgeventdata == NULL);
4237 
4238  if( (*consdata)->rhsvar != NULL )
4239  {
4240  SCIP_CALL( SCIPreleaseVar(scip, &(*consdata)->rhsvar) );
4241  }
4242 
4243  /* free nonlinear row representation
4244  * normally released in exitsol, but constraint may be deleted early (e.g., if found redundant)
4245  */
4246  if( (*consdata)->nlrow != NULL )
4247  {
4248  SCIP_CALL( SCIPreleaseNlRow(scip, &(*consdata)->nlrow) );
4249  }
4250 
4251  SCIPfreeBlockMemory(scip, consdata);
4252 
4253  return SCIP_OKAY;
4254 }
4255 
4256 
4257 /** transforms constraint data into data belonging to the transformed problem */
4258 static
4259 SCIP_DECL_CONSTRANS(consTransSOC)
4260 {
4261  SCIP_CONSDATA* consdata;
4262  SCIP_CONSHDLRDATA* conshdlrdata;
4263  SCIP_CONSDATA* sourcedata;
4264  char s[SCIP_MAXSTRLEN];
4265  int i;
4266 
4267  assert(scip != NULL);
4268  assert(conshdlr != NULL);
4269  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
4270  assert(sourcecons != NULL);
4271  assert(targetcons != NULL);
4272 
4273  /* get constraint handler data */
4274  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4275  assert(conshdlrdata != NULL);
4276 
4277  SCIPdebugMsg(scip, "Transforming SOC constraint: <%s>.\n", SCIPconsGetName(sourcecons) );
4278 
4279  /* get data of original constraint */
4280  sourcedata = SCIPconsGetData(sourcecons);
4281  assert(sourcedata != NULL);
4282  assert(sourcedata->vars != NULL);
4283  assert(sourcedata->coefs != NULL);
4284  assert(sourcedata->offsets != NULL);
4285 
4286  /* create constraint data */
4287  SCIP_CALL( SCIPallocBlockMemory(scip, &consdata) );
4288 
4289  consdata->nvars = sourcedata->nvars;
4290  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->vars, consdata->nvars) );
4291  SCIP_CALL( SCIPgetTransformedVars(scip, consdata->nvars, sourcedata->vars, consdata->vars) );
4292  for( i = 0; i < consdata->nvars; ++i )
4293  {
4294  SCIP_CALL( SCIPcaptureVar(scip, consdata->vars[i]) );
4295  }
4296 
4297  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &consdata->coefs, sourcedata->coefs, consdata->nvars) );
4298  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &consdata->offsets, sourcedata->offsets, consdata->nvars) );
4299  consdata->constant = sourcedata->constant;
4300 
4301  SCIP_CALL( SCIPgetTransformedVar(scip, sourcedata->rhsvar, &consdata->rhsvar) );
4302  consdata->rhsoffset = sourcedata->rhsoffset;
4303  consdata->rhscoeff = sourcedata->rhscoeff;
4304 
4305  SCIP_CALL( SCIPcaptureVar(scip, consdata->rhsvar) );
4306 
4307  consdata->nlrow = NULL;
4308  consdata->lhsbndchgeventdata = NULL;
4309  consdata->isapproxadded = FALSE;
4310 
4311  /* create transformed constraint with the same flags */
4312  (void) SCIPsnprintf(s, SCIP_MAXSTRLEN, "t_%s", SCIPconsGetName(sourcecons));
4313  SCIP_CALL( SCIPcreateCons(scip, targetcons, s, conshdlr, consdata,
4314  SCIPconsIsInitial(sourcecons), SCIPconsIsSeparated(sourcecons),
4315  SCIPconsIsEnforced(sourcecons), SCIPconsIsChecked(sourcecons),
4316  SCIPconsIsPropagated(sourcecons), SCIPconsIsLocal(sourcecons),
4317  SCIPconsIsModifiable(sourcecons), SCIPconsIsDynamic(sourcecons),
4318  SCIPconsIsRemovable(sourcecons), SCIPconsIsStickingAtNode(sourcecons)) );
4319 
4320  SCIP_CALL( catchVarEvents(scip, conshdlrdata->eventhdlr, *targetcons) );
4321 
4322  return SCIP_OKAY;
4323 }
4324 
4325 /** separation method of constraint handler for LP solutions */
4326 static
4327 SCIP_DECL_CONSSEPALP(consSepalpSOC)
4328 {
4329  SCIP_CONSHDLRDATA* conshdlrdata;
4330  SCIP_CONS* maxviolcon;
4331  SCIP_Bool sepasuccess;
4332  SCIP_Bool cutoff;
4333 
4334  assert(scip != NULL);
4335  assert(conshdlr != NULL);
4336  assert(conss != NULL || nconss == 0);
4337  assert(result != NULL);
4338 
4339  *result = SCIP_DIDNOTFIND;
4340 
4341  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4342  assert(conshdlrdata != NULL);
4343 
4344  SCIP_CALL( computeViolations(scip, conshdlr, conss, nconss, NULL, &maxviolcon) );
4345  if( maxviolcon == NULL )
4346  return SCIP_OKAY;
4347 
4348  /* at root, check if we want to solve the NLP relaxation and use its solutions as reference point
4349  * if there is something convex, then linearizing in the solution of the NLP relaxation can be very useful
4350  */
4351  if( SCIPgetDepth(scip) == 0 && !conshdlrdata->sepanlp &&
4352  (SCIPgetNContVars(scip) >= conshdlrdata->sepanlpmincont * SCIPgetNVars(scip) || (SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_UNBOUNDEDRAY && conshdlrdata->sepanlpmincont <= 1.0)) &&
4353  SCIPisNLPConstructed(scip) && SCIPgetNNlpis(scip) > 0 )
4354  {
4355  SCIP_NLPSOLSTAT solstat;
4356  SCIP_Bool solvednlp;
4357 
4358  solstat = SCIPgetNLPSolstat(scip);
4359  solvednlp = FALSE;
4360  if( solstat == SCIP_NLPSOLSTAT_UNKNOWN )
4361  {
4362  /* NLP is not solved yet, so we do it now and update solstat */
4363 
4364  /* ensure linear conss are in NLP */
4365  if( conshdlrdata->subnlpheur != NULL )
4366  {
4367  SCIP_CALL( SCIPaddLinearConsToNlpHeurSubNlp(scip, conshdlrdata->subnlpheur, TRUE, TRUE) );
4368  }
4369 
4370  /* set LP solution as starting values, if available */
4372  {
4374  }
4375 
4376  /* SCIP_CALL( SCIPsetNLPIntPar(scip, SCIP_NLPPAR_VERBLEVEL, 1) ); */
4377  SCIP_CALL( SCIPsolveNLP(scip) );
4378 
4379  solstat = SCIPgetNLPSolstat(scip);
4380  SCIPdebugMsg(scip, "solved NLP relax, solution status: %d\n", solstat);
4381 
4382  solvednlp = TRUE;
4383  }
4384 
4385  conshdlrdata->sepanlp = TRUE;
4386 
4387  if( solstat == SCIP_NLPSOLSTAT_GLOBINFEASIBLE )
4388  {
4389  SCIPdebugMsg(scip, "NLP relaxation is globally infeasible, thus can cutoff node\n");
4390  *result = SCIP_CUTOFF;
4391  return SCIP_OKAY;
4392  }
4393 
4394  if( solstat <= SCIP_NLPSOLSTAT_FEASIBLE )
4395  {
4396  /* if we have feasible NLP solution, generate linearization cuts there */
4397  SCIP_Bool lpsolseparated;
4398  SCIP_SOL* nlpsol;
4399 
4400  SCIP_CALL( SCIPcreateNLPSol(scip, &nlpsol, NULL) );
4401  assert(nlpsol != NULL);
4402 
4403  /* if we solved the NLP and solution is integral, then pass it to trysol heuristic */
4404  if( solvednlp && conshdlrdata->trysolheur != NULL )
4405  {
4406  int nfracvars = 0;
4407 
4408  if( SCIPgetNBinVars(scip) > 0 || SCIPgetNIntVars(scip) > 0 )
4409  {
4410  SCIP_CALL( SCIPgetNLPFracVars(scip, NULL, NULL, NULL, &nfracvars, NULL) );
4411  }
4412 
4413  if( nfracvars == 0 )
4414  {
4415  SCIP_CALL( SCIPheurPassSolTrySol(scip, conshdlrdata->trysolheur, nlpsol) );
4416  }
4417  }
4418 
4419  SCIP_CALL( addLinearizationCuts(scip, conshdlr, conss, nconss, nlpsol, &lpsolseparated, SCIPgetSepaMinEfficacy(scip), &cutoff) );
4420 
4421  SCIP_CALL( SCIPfreeSol(scip, &nlpsol) );
4422 
4423  if( cutoff )
4424  {
4425  *result = SCIP_CUTOFF;
4426  return SCIP_OKAY;
4427  }
4428 
4429  /* if a cut that separated the LP solution was added, then return, otherwise continue with usual separation in LP solution */
4430  if( lpsolseparated )
4431  {
4432  SCIPdebugMsg(scip, "linearization cuts separate LP solution\n");
4433 
4434  *result = SCIP_SEPARATED;
4435 
4436  return SCIP_OKAY;
4437  }
4438  }
4439  }
4440  /* if we do not want to try solving the NLP, or have no NLP, or have no NLP solver, or solving the NLP failed,
4441  * or separating with NLP solution as reference point failed, then try (again) with LP solution as reference point
4442  */
4443 
4444  SCIP_CALL( separatePoint(scip, conshdlr, conss, nconss, nusefulconss, NULL, FALSE, &cutoff, &sepasuccess) );
4445  if( cutoff )
4446  *result = SCIP_CUTOFF;
4447  else if ( sepasuccess )
4448  *result = SCIP_SEPARATED;
4449 
4450  return SCIP_OKAY;
4451 }
4452 
4453 
4454 /** separation method of constraint handler for arbitrary primal solutions */
4455 static
4456 SCIP_DECL_CONSSEPASOL(consSepasolSOC)
4457 {
4458  SCIP_CONS* maxviolcon;
4459  SCIP_Bool sepasuccess;
4460  SCIP_Bool cutoff;
4461 
4462  assert(scip != NULL);
4463  assert(conss != NULL || nconss == 0);
4464  assert(result != NULL);
4465  assert(sol != NULL);
4466 
4467  *result = SCIP_DIDNOTFIND;
4468 
4469  SCIP_CALL( computeViolations(scip, conshdlr, conss, nconss, sol, &maxviolcon) );
4470  if( maxviolcon == NULL )
4471  return SCIP_OKAY;
4472 
4473  SCIP_CALL( separatePoint(scip, conshdlr, conss, nconss, nusefulconss, sol, FALSE, &cutoff, &sepasuccess) );
4474  if( cutoff )
4475  *result = SCIP_CUTOFF;
4476  else if ( sepasuccess )
4477  *result = SCIP_SEPARATED;
4478 
4479  return SCIP_OKAY;
4480 }
4481 
4482 
4483 /** constraint enforcing method of constraint handler for LP solutions */
4484 static
4485 SCIP_DECL_CONSENFOLP(consEnfolpSOC)
4486 { /*lint --e{715}*/
4487  SCIP_CALL( enforceConstraint(scip, conshdlr, conss, nconss, nusefulconss, NULL, result) );
4488 
4489  return SCIP_OKAY;
4490 }
4491 
4492 
4493 /** constraint enforcing method of constraint handler for relaxation solutions */
4494 static
4495 SCIP_DECL_CONSENFORELAX(consEnforelaxSOC)
4496 { /*lint --e{715}*/
4497  SCIP_CALL( enforceConstraint(scip, conshdlr, conss, nconss, nusefulconss, sol, result) );
4498 
4499  return SCIP_OKAY;
4500 }
4501 
4502 
4503 /** constraint enforcing method of constraint handler for pseudo solutions */
4504 static
4505 SCIP_DECL_CONSENFOPS(consEnfopsSOC)
4507  SCIP_CONS* maxviolcons;
4508 
4509  assert(scip != NULL);
4510  assert(conss != NULL || nconss == 0);
4511  assert(result != NULL);
4512 
4513  SCIP_CALL( computeViolations(scip, conshdlr, conss, nconss, NULL, &maxviolcons) );
4514 
4515  if( maxviolcons == NULL )
4516  *result = SCIP_FEASIBLE;
4517 
4518  *result = SCIP_INFEASIBLE;
4519 
4520  return SCIP_OKAY;
4521 } /*lint !e715*/
4522 
4523 
4524 /** feasibility check method of constraint handler for integral solutions */
4525 static
4526 SCIP_DECL_CONSCHECK(consCheckSOC)
4528  SCIP_CONSHDLRDATA* conshdlrdata;
4529  SCIP_CONSDATA* consdata;
4530  SCIP_Real maxviol;
4531  SCIP_Bool dolinfeasshift;
4532  SCIP_SOL* polishedsol;
4533  int c;
4534 
4535  assert(scip != NULL);
4536  assert(conshdlr != NULL);
4537  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
4538  assert(conss != NULL || nconss == 0);
4539  assert(result != NULL );
4540 
4541  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4542  assert(conshdlrdata != NULL);
4543 
4544  *result = SCIP_FEASIBLE;
4545  maxviol = 0.0;
4546 
4547  dolinfeasshift = conshdlrdata->linfeasshift && (conshdlrdata->trysolheur != NULL);
4548  polishedsol = NULL;
4549 
4550  for( c = 0; c < nconss; ++c )
4551  {
4552  SCIP_CALL( computeViolation(scip, conshdlr, conss[c], sol) ); /*lint !e613*/
4553 
4554  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
4555  assert(consdata != NULL);
4556 
4557  /* if feasible, just continue */
4558  if( !SCIPisGT(scip, consdata->violation, SCIPfeastol(scip)) )
4559  continue;
4560 
4561  *result = SCIP_INFEASIBLE;
4562 
4563  if( consdata->violation > maxviol )
4564  maxviol = consdata->violation;
4565 
4566  if( printreason )
4567  {
4568  SCIP_CALL( SCIPprintCons(scip, conss[c], NULL) ); /*lint !e613*/
4569  SCIPinfoMessage(scip, NULL, ";\n\tviolation: %g\n", consdata->violation);
4570  }
4571 
4572  /* if we do linear feasibility shifting, then try to adjust solution */
4573  if( dolinfeasshift )
4574  {
4575  if( SCIPvarGetStatus(consdata->rhsvar) != SCIP_VARSTATUS_MULTAGGR &&
4576  !SCIPisInfinity(scip, REALABS(consdata->lhsval)) &&
4577  ( (consdata->rhscoeff > 0.0 && SCIPvarMayRoundUp (consdata->rhsvar)) ||
4578  (consdata->rhscoeff < 0.0 && SCIPvarMayRoundDown(consdata->rhsvar)) ) )
4579  {
4580  SCIP_Bool success;
4581 
4582  if( polishedsol == NULL )
4583  {
4584  if( sol != NULL )
4585  {
4586  SCIP_CALL( SCIPcreateSolCopy(scip, &polishedsol, sol) );
4587  }
4588  else
4589  {
4590  SCIP_CALL( SCIPcreateLPSol(scip, &polishedsol, NULL) );
4591  }
4592  SCIP_CALL( SCIPunlinkSol(scip, polishedsol) );
4593  }
4594  SCIP_CALL( polishSolution(scip, conss[c], polishedsol, &success) ); /*lint !e613*/
4595 
4596  /* disable solution polishing if we failed for this constraint */
4597  dolinfeasshift = success;
4598  }
4599  else /* if locks of the variable are bad or rhs is multi-aggregated, disable solution polishing */
4600  {
4601  dolinfeasshift = FALSE;
4602  }
4603  }
4604 
4605  /* if solution polishing is off and there is no NLP heuristic or we just check the LP solution,
4606  * then there is no need to check remaining constraints (NLP heuristic will pick up LP solution anyway) */
4607  if( !dolinfeasshift && (conshdlrdata->subnlpheur == NULL || sol == NULL) && !completely )
4608  break;
4609  }
4610 
4611  /* if we failed to polish solution, clear the solution */
4612  if( !dolinfeasshift && polishedsol != NULL )
4613  {
4614  SCIP_CALL( SCIPfreeSol(scip, &polishedsol) );
4615  }
4616 
4617  if( polishedsol != NULL )
4618  {
4619  assert(*result == SCIP_INFEASIBLE);
4620  SCIP_CALL( SCIPheurPassSolTrySol(scip, conshdlrdata->trysolheur, polishedsol) );
4621  SCIP_CALL( SCIPfreeSol(scip, &polishedsol) );
4622  }
4623  else if( conshdlrdata->subnlpheur != NULL && sol != NULL && *result == SCIP_INFEASIBLE && !SCIPisInfinity(scip, maxviol) )
4624  {
4625  SCIP_CALL( SCIPupdateStartpointHeurSubNlp(scip, conshdlrdata->subnlpheur, sol, maxviol) );
4626  }
4627 
4628  return SCIP_OKAY;
4629 } /*lint !e715*/
4630 
4631 
4632 /** domain propagation method of constraint handler */
4633 static
4634 SCIP_DECL_CONSPROP(consPropSOC)
4636  SCIP_RESULT propresult;
4637  int c;
4638  int nchgbds;
4639  SCIP_Bool redundant;
4640 
4641  assert(scip != NULL);
4642  assert(conss != NULL || ((nconss == 0) && (nmarkedconss == 0)));
4643  assert(result != NULL);
4644 
4645  *result = SCIP_DIDNOTFIND;
4646  nchgbds = 0;
4647 
4648  for( c = 0; c < nmarkedconss && *result != SCIP_CUTOFF; ++c )
4649  {
4650  SCIP_CALL( propagateBounds(scip, conss[c], &propresult, &nchgbds, &redundant) ); /*lint !e613*/
4651  if( propresult != SCIP_DIDNOTFIND && propresult != SCIP_DIDNOTRUN )
4652  *result = propresult;
4653  }
4654 
4655  return SCIP_OKAY;
4656 } /*lint !e715*/
4657 
4658 
4659 /** presolving method of constraint handler */
4660 static
4661 SCIP_DECL_CONSPRESOL(consPresolSOC)
4663  SCIP_CONSHDLRDATA* conshdlrdata;
4664  SCIP_CONSDATA* consdata;
4665  int c;
4666  SCIP_RESULT propresult;
4667  SCIP_Bool iscutoff;
4668  SCIP_Bool isdeleted;
4669 
4670  assert(scip != NULL);
4671  assert(conss != NULL || nconss == 0);
4672  assert(conshdlr != NULL);
4673  assert(result != NULL);
4674 
4675  *result = SCIP_DIDNOTFIND;
4676 
4677  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4678  assert(conshdlrdata != NULL);
4679 
4680  for( c = 0; c < nconss; ++c )
4681  {
4682  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
4683  assert(consdata != NULL);
4684 
4685  SCIP_CALL( presolveRemoveFixedVariables(scip, conshdlr, conss[c], ndelconss, nupgdconss, nchgbds, nfixedvars, &iscutoff, &isdeleted) ); /*lint !e613*/
4686  if( iscutoff )
4687  {
4688  *result = SCIP_CUTOFF;
4689  return SCIP_OKAY;
4690  }
4691  if( isdeleted )
4692  {
4693  /* conss[c] has been deleted */
4694  *result = SCIP_SUCCESS;
4695  continue;
4696  }
4697 
4698  if( conshdlrdata->nauxvars > 0 && !consdata->isapproxadded && consdata->nvars > 1 )
4699  {
4700  SCIP_CALL( presolveCreateOuterApprox(scip, consdata->nvars, consdata->vars, consdata->coefs, consdata->offsets,
4701  consdata->rhsvar, consdata->rhscoeff, consdata->rhscoeff, consdata->constant, SCIPconsGetName(conss[c]), conss[c],
4702  conshdlrdata->nauxvars, conshdlrdata->glineur, naddconss) ); /*lint !e613*/
4703  consdata->isapproxadded = TRUE;
4704  }
4705 
4706  if( (presoltiming & SCIP_PRESOLTIMING_FAST) != 0 )
4707  {
4708  SCIP_Bool redundant;
4709 
4710  SCIP_CALL( propagateBounds(scip, conss[c], &propresult, nchgbds, &redundant) ); /*lint !e613*/
4711  switch( propresult )
4712  {
4713  case SCIP_DIDNOTRUN:
4714  case SCIP_DIDNOTFIND:
4715  break;
4716  case SCIP_REDUCEDDOM:
4717  *result = SCIP_SUCCESS;
4718  break;
4719  case SCIP_CUTOFF:
4720  *result = SCIP_CUTOFF;
4721  SCIPdebugMsg(scip, "infeasible in presolve due to propagation for constraint %s\n", SCIPconsGetName(conss[c])); /*lint !e613*/
4722  return SCIP_OKAY;
4723  default:
4724  SCIPerrorMessage("unexpected result from propagation: %d\n", propresult);
4725  return SCIP_ERROR;
4726  } /*lint !e788*/
4727  if( redundant )
4728  ++*ndelconss;
4729  }
4730 
4731  /* disaggregate each lhs term to a quadratic constraint by using auxiliary variables */
4732  if( conshdlrdata->disaggregate && (presoltiming & SCIP_PRESOLTIMING_EXHAUSTIVE) != 0 )
4733  {
4734  SCIP_Bool success;
4735 
4736  SCIP_CALL( disaggregate(scip, conss[c], consdata, naddconss, ndelconss, &success) ); /*lint !e613*/
4737 
4738  if( success )
4739  {
4740  SCIPdebugMsg(scip, "disaggregated SOC constraint\n");
4741 
4742  /* conss[c] has been deleted */
4743  *result = SCIP_SUCCESS;
4744  continue;
4745  }
4746  }
4747  }
4748 
4749  return SCIP_OKAY;
4750 } /*lint !e715*/
4751 
4752 
4753 /** variable rounding lock method of constraint handler */
4754 static
4755 SCIP_DECL_CONSLOCK(consLockSOC)
4757  SCIP_CONSDATA* consdata;
4758  int i;
4759 
4760  assert(scip != NULL);
4761  assert(conshdlr != NULL);
4762  assert(cons != NULL);
4763  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
4764  assert(locktype == SCIP_LOCKTYPE_MODEL);
4765 
4766  consdata = SCIPconsGetData(cons);
4767  assert(consdata != NULL);
4768 
4769  SCIPdebugMsg(scip, "Locking constraint <%s>.\n", SCIPconsGetName(cons));
4770 
4771  /* Changing variables x_i, i <= n, in both directions can lead to an infeasible solution. */
4772  for( i = 0; i < consdata->nvars; ++i )
4773  {
4774  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->vars[i], locktype, nlockspos + nlocksneg, nlockspos + nlocksneg) );
4775  }
4776 
4777  /* Rounding x_{n+1} up will not violate a solution. */
4778  if( consdata->rhsvar != NULL )
4779  {
4780  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->rhsvar, locktype, consdata->rhscoeff > 0.0 ? nlockspos : nlocksneg, consdata->rhscoeff > 0.0 ? nlocksneg : nlockspos) );
4781  }
4782 
4783  return SCIP_OKAY;
4784 }
4785 
4786 /** constraint display method of constraint handler */
4787 static
4788 SCIP_DECL_CONSPRINT(consPrintSOC)
4789 {
4790  SCIP_CONSDATA* consdata;
4791  int i;
4792 
4793  assert(scip != NULL);
4794  assert(conshdlr != NULL);
4795  assert(cons != NULL);
4796  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
4797 
4798  consdata = SCIPconsGetData(cons);
4799  assert(consdata != NULL);
4800 
4801  SCIPinfoMessage(scip, file, "sqrt( ");
4802  if( consdata->constant != 0.0 )
4803  {
4804  SCIPinfoMessage(scip, file, "%.15g", consdata->constant);
4805  }
4806 
4807  for( i = 0; i < consdata->nvars; ++i )
4808  {
4809  SCIPinfoMessage(scip, file, "+ (%.15g*(", consdata->coefs[i]);
4810  SCIP_CALL( SCIPwriteVarName(scip, file, consdata->vars[i], TRUE) );
4811  SCIPinfoMessage(scip, file, "%+.15g))^2 ", consdata->offsets[i]);
4812  }
4813 
4814  SCIPinfoMessage(scip, file, ") <= ");
4815  if( consdata->rhsvar != NULL )
4816  {
4817  SCIPinfoMessage(scip, file, "%.15g*(", consdata->rhscoeff);
4818  SCIP_CALL( SCIPwriteVarName(scip, file, consdata->rhsvar, TRUE) );
4819  SCIPinfoMessage(scip, file, "%+.15g)", consdata->rhsoffset);
4820  }
4821  else
4822  {
4823  SCIPinfoMessage(scip, file, "%.15g", consdata->rhscoeff*consdata->rhsoffset);
4824  }
4825 
4826  return SCIP_OKAY;
4827 }
4828 
4829 /** constraint copying method of constraint handler */
4830 static
4831 SCIP_DECL_CONSCOPY(consCopySOC)
4833  SCIP_CONSDATA* consdata;
4834  SCIP_VAR** vars;
4835  SCIP_VAR* rhsvar;
4836  int i;
4837 
4838  assert(scip != NULL);
4839  assert(cons != NULL);
4840  assert(sourcescip != NULL);
4841  assert(sourceconshdlr != NULL);
4842  assert(sourcecons != NULL);
4843  assert(varmap != NULL);
4844  assert(valid != NULL);
4845  assert(stickingatnode == FALSE);
4846 
4847  consdata = SCIPconsGetData(sourcecons);
4848  assert(consdata != NULL);
4849 
4850  *valid = TRUE;
4851 
4852  SCIP_CALL( SCIPallocBufferArray(sourcescip, &vars, consdata->nvars) );
4853 
4854  /* map variables to active variables of the target SCIP */
4855  for( i = 0; i < consdata->nvars && *valid; ++i )
4856  {
4857  SCIP_CALL( SCIPgetVarCopy(sourcescip, scip, consdata->vars[i], &vars[i], varmap, consmap, global, valid) );
4858  assert(!(*valid) || vars[i] != NULL);
4859  }
4860 
4861  /* map rhs variable to active variable of the target SCIP */
4862  if( *valid )
4863  {
4864  SCIP_CALL( SCIPgetVarCopy(sourcescip, scip, consdata->rhsvar, &rhsvar, varmap, consmap, global, valid) );
4865  assert(!(*valid) || rhsvar != NULL);
4866 
4867  /* only create the target constraint, if all variables could be copied */
4868  if( *valid )
4869  {
4870  SCIP_CALL( SCIPcreateConsSOC(scip, cons, name ? name : SCIPconsGetName(sourcecons),
4871  consdata->nvars, vars, consdata->coefs, consdata->offsets, consdata->constant,
4872  rhsvar, consdata->rhscoeff, consdata->rhsoffset,
4873  initial, separate, enforce, check, propagate, local, modifiable, dynamic, removable) ); /*lint !e644 */
4874  }
4875  }
4876 
4877  SCIPfreeBufferArray(sourcescip, &vars);
4878 
4879  return SCIP_OKAY;
4880 }
4881 
4882 
4883 /** constraint parsing method of constraint handler */
4884 static
4885 SCIP_DECL_CONSPARSE(consParseSOC)
4886 { /*lint --e{715}*/
4887  SCIP_VAR* var;
4888  SCIP_VAR** vars;
4889  SCIP_Real* coefs;
4890  SCIP_Real* offsets;
4891  int nvars;
4892  int varssize;
4893  SCIP_VAR* rhsvar;
4894  SCIP_Real rhscoef;
4895  SCIP_Real rhsoffset;
4896  SCIP_Real constant;
4897  SCIP_Real coef;
4898  SCIP_Real offset;
4899  char* endptr;
4900 
4901  assert(scip != NULL);
4902  assert(success != NULL);
4903  assert(str != NULL);
4904  assert(name != NULL);
4905  assert(cons != NULL);
4906 
4907  /* check that string starts with "sqrt( " */
4908  if( strncmp(str, "sqrt( ", 6) != 0 )
4909  {
4910  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected 'sqrt( ' at begin of soc constraint string '%s'\n", str);
4911  *success = FALSE;
4912  return SCIP_OKAY;
4913  }
4914  str += 6;
4915 
4916  *success = TRUE;
4917 
4918  /* check if we have a constant in the beginning */
4919  if( SCIPstrToRealValue(str, &constant, &endptr) )
4920  str = endptr;
4921  else
4922  constant = 0.0;
4923 
4924  nvars = 0;
4925  varssize = 5;
4926  SCIP_CALL( SCIPallocBufferArray(scip, &vars, varssize) );
4927  SCIP_CALL( SCIPallocBufferArray(scip, &coefs, varssize) );
4928  SCIP_CALL( SCIPallocBufferArray(scip, &offsets, varssize) );
4929 
4930  /* read '+ (coef*(var+offset))^2' on lhs, as long as possible */
4931  while( *str != '\0' )
4932  {
4933  /* skip whitespace */
4934  while( isspace((int)*str) )
4935  ++str;
4936 
4937  /* stop if no more coefficients */
4938  if( strncmp(str, "+ (", 3) != 0 )
4939  break;
4940 
4941  str += 3;
4942 
4943  /* parse coef */
4944  if( !SCIPstrToRealValue(str, &coef, &endptr) )
4945  {
4946  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected coefficient at begin of '%s'\n", str);
4947  *success = FALSE;
4948  break;
4949  }
4950  str = endptr;
4951 
4952  if( strncmp(str, "*(", 2) != 0 )
4953  {
4954  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected '*(' at begin of '%s'\n", str);
4955  *success = FALSE;
4956  break;
4957  }
4958  str += 2;
4959 
4960  /* parse variable name */
4961  SCIP_CALL( SCIPparseVarName(scip, str, &var, &endptr) );
4962  if( var == NULL )
4963  {
4964  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "unknown variable name at '%s'\n", str);
4965  *success = FALSE;
4966  break;
4967  }
4968  str = endptr;
4969 
4970  /* parse offset */
4971  if( !SCIPstrToRealValue(str, &offset, &endptr) )
4972  {
4973  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected offset at begin of '%s'\n", str);
4974  *success = FALSE;
4975  break;
4976  }
4977  str = endptr;
4978 
4979  if( strncmp(str, "))^2", 4) != 0 )
4980  {
4981  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected '))^2' at begin of '%s'\n", str);
4982  *success = FALSE;
4983  break;
4984  }
4985  str += 4;
4986 
4987  if( varssize <= nvars )
4988  {
4989  varssize = SCIPcalcMemGrowSize(scip, varssize+1);
4990  SCIP_CALL( SCIPreallocBufferArray(scip, &vars, varssize) );
4991  SCIP_CALL( SCIPreallocBufferArray(scip, &coefs, varssize) );
4992  SCIP_CALL( SCIPreallocBufferArray(scip, &offsets, varssize) );
4993  }
4994  vars[nvars] = var;
4995  coefs[nvars] = coef;
4996  offsets[nvars] = offset;
4997  ++nvars;
4998  }
4999 
5000  if( strncmp(str, ") <=", 4) != 0 )
5001  {
5002  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected ') <=' at begin of '%s'\n", str);
5003  *success = FALSE;
5004  }
5005  str += 4;
5006 
5007  /* read rhs coef*(var+offset) or just a constant */
5008 
5009  /* parse coef */
5010  if( *success )
5011  {
5012  /* skip whitespace */
5013  while( isspace((int)*str) )
5014  ++str;
5015 
5016  if( !SCIPstrToRealValue(str, &rhscoef, &endptr) )
5017  {
5018  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected coefficient at begin of '%s'\n", str);
5019  *success = FALSE;
5020  }
5021  str = endptr;
5022 
5023  /* skip whitespace */
5024  while( isspace((int)*str) )
5025  ++str;
5026  }
5027 
5028  /* parse *(var+offset) */
5029  if( *str != '\0' )
5030  {
5031  if( *success )
5032  {
5033  if( strncmp(str, "*(", 2) != 0 )
5034  {
5035  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected '*(' at begin of '%s'\n", str);
5036  *success = FALSE;
5037  }
5038  else
5039  {
5040  str += 2;
5041  }
5042  }
5043 
5044  /* parse variable name */
5045  if( *success )
5046  {
5047  SCIP_CALL( SCIPparseVarName(scip, str, &rhsvar, &endptr) );
5048  if( rhsvar == NULL )
5049  {
5050  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "unknown variable name at '%s'\n", str);
5051  *success = FALSE;
5052  }
5053  else
5054  {
5055  str = endptr;
5056  }
5057  }
5058 
5059  /* parse offset */
5060  if( *success )
5061  {
5062  if( !SCIPstrToRealValue(str, &rhsoffset, &endptr) )
5063  {
5064  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected offset at begin of '%s'\n", str);
5065  *success = FALSE;
5066  }
5067  else
5068  {
5069  str = endptr;
5070  }
5071  }
5072 
5073  if( *success )
5074  {
5075  if( *str != ')' )
5076  {
5077  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected ')' at begin of '%s'\n", str);
5078  *success = FALSE;
5079  }
5080  }
5081  }
5082  else if( *success )
5083  {
5084  /* only a constant at right hand side */
5085  rhsoffset = rhscoef; /*lint !e644*/
5086  rhscoef = 1.0;
5087  rhsvar = NULL;
5088  }
5089 
5090  if( *success )
5091  {
5092  assert(!stickingatnode);
5093  SCIP_CALL( SCIPcreateConsSOC(scip, cons, name, nvars, vars, coefs, offsets, constant, rhsvar, rhscoef, rhsoffset,
5094  initial, separate, enforce, check, propagate, local, modifiable, dynamic, removable) ); /*lint !e644 */
5095  }
5096 
5097  SCIPfreeBufferArray(scip, &offsets);
5098  SCIPfreeBufferArray(scip, &coefs);
5099  SCIPfreeBufferArray(scip, &vars);
5100 
5101  return SCIP_OKAY;
5102 }
5103 
5104 /** constraint method of constraint handler which returns the variables (if possible) */
5105 static
5106 SCIP_DECL_CONSGETVARS(consGetVarsSOC)
5107 { /*lint --e{715}*/
5108  SCIP_CONSDATA* consdata;
5109 
5110  consdata = SCIPconsGetData(cons);
5111  assert(consdata != NULL);
5112 
5113  if( varssize < consdata->nvars + 1 )
5114  (*success) = FALSE;
5115  else
5116  {
5117  BMScopyMemoryArray(vars, consdata->vars, consdata->nvars);
5118  vars[consdata->nvars] = consdata->rhsvar;
5119  (*success) = TRUE;
5120  }
5121 
5122  return SCIP_OKAY;
5123 }
5124 
5125 /** constraint method of constraint handler which returns the number of variable (if possible) */
5126 static
5127 SCIP_DECL_CONSGETNVARS(consGetNVarsSOC)
5128 { /*lint --e{715}*/
5129  SCIP_CONSDATA* consdata;
5130 
5131  assert(cons != NULL);
5132 
5133  consdata = SCIPconsGetData(cons);
5134  assert(consdata != NULL);
5135 
5136  (*nvars) = consdata->nvars + 1;
5137  (*success) = TRUE;
5138 
5139  return SCIP_OKAY;
5140 }
5141 
5142 /*
5143  * constraint specific interface methods
5144  */
5145 
5146 /** creates the handler for second order cone constraints and includes it in SCIP */
5148  SCIP* scip /**< SCIP data structure */
5149  )
5150 {
5151  SCIP_CONSHDLRDATA* conshdlrdata;
5152  SCIP_CONSHDLR* conshdlr;
5153  SCIP_EVENTHDLR* eventhdlr;
5154 
5155  /* create constraint handler data */
5156  SCIP_CALL( SCIPallocBlockMemory(scip, &conshdlrdata) );
5157  conshdlrdata->subnlpheur = NULL;
5158  conshdlrdata->trysolheur = NULL;
5159 
5160  SCIP_CALL( SCIPincludeEventhdlrBasic(scip, &eventhdlr, CONSHDLR_NAME"_boundchange",
5161  "signals a bound change to a second order cone constraint",
5162  processVarEvent, NULL) );
5163  conshdlrdata->eventhdlr = eventhdlr;
5164 
5165  SCIP_CALL( SCIPincludeEventhdlrBasic(scip, NULL, CONSHDLR_NAME"_newsolution",
5166  "handles the event that a new primal solution has been found",
5167  processNewSolutionEvent, NULL) );
5168 
5169  /* include constraint handler */
5172  consEnfolpSOC, consEnfopsSOC, consCheckSOC, consLockSOC,
5173  conshdlrdata) );
5174  assert(conshdlr != NULL);
5175 
5176  /* set non-fundamental callbacks via specific setter functions */
5177  SCIP_CALL( SCIPsetConshdlrCopy(scip, conshdlr, conshdlrCopySOC, consCopySOC) );
5178  SCIP_CALL( SCIPsetConshdlrDelete(scip, conshdlr, consDeleteSOC) );
5179  SCIP_CALL( SCIPsetConshdlrExit(scip, conshdlr, consExitSOC) );
5180  SCIP_CALL( SCIPsetConshdlrExitpre(scip, conshdlr, consExitpreSOC) );
5181  SCIP_CALL( SCIPsetConshdlrExitsol(scip, conshdlr, consExitsolSOC) );
5182  SCIP_CALL( SCIPsetConshdlrFree(scip, conshdlr, consFreeSOC) );
5183  SCIP_CALL( SCIPsetConshdlrGetVars(scip, conshdlr, consGetVarsSOC) );
5184  SCIP_CALL( SCIPsetConshdlrGetNVars(scip, conshdlr, consGetNVarsSOC) );
5185  SCIP_CALL( SCIPsetConshdlrInit(scip, conshdlr, consInitSOC) );
5186  SCIP_CALL( SCIPsetConshdlrInitsol(scip, conshdlr, consInitsolSOC) );
5187  SCIP_CALL( SCIPsetConshdlrParse(scip, conshdlr, consParseSOC) );
5188  SCIP_CALL( SCIPsetConshdlrPresol(scip, conshdlr, consPresolSOC, CONSHDLR_MAXPREROUNDS, CONSHDLR_PRESOLTIMING) );
5189  SCIP_CALL( SCIPsetConshdlrPrint(scip, conshdlr, consPrintSOC) );
5191  SCIP_CALL( SCIPsetConshdlrSepa(scip, conshdlr, consSepalpSOC, consSepasolSOC, CONSHDLR_SEPAFREQ,
5193  SCIP_CALL( SCIPsetConshdlrTrans(scip, conshdlr, consTransSOC) );
5194  SCIP_CALL( SCIPsetConshdlrEnforelax(scip, conshdlr, consEnforelaxSOC) );
5195 
5196  if( SCIPfindConshdlr(scip,"quadratic") != NULL )
5197  {
5198  /* notify function that upgrades quadratic constraint to SOC's */
5200  }
5201 
5202  /* add soc constraint handler parameters */
5203  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/projectpoint",
5204  "whether the reference point of a cut should be projected onto the feasible set of the SOC constraint",
5205  &conshdlrdata->projectpoint, TRUE, FALSE, NULL, NULL) );
5206 
5207  SCIP_CALL( SCIPaddIntParam (scip, "constraints/" CONSHDLR_NAME "/nauxvars",
5208  "number of auxiliary variables to use when creating a linear outer approx. of a SOC3 constraint; 0 to turn off",
5209  &conshdlrdata->nauxvars, FALSE, 0, 0, INT_MAX, NULL, NULL) );
5210 
5211  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/glineur",
5212  "whether the Glineur Outer Approximation should be used instead of Ben-Tal Nemirovski",
5213  &conshdlrdata->glineur, FALSE, TRUE, NULL, NULL) );
5214 
5215  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/sparsify",
5216  "whether to sparsify cuts",
5217  &conshdlrdata->sparsify, TRUE, FALSE, NULL, NULL) );
5218 
5219  SCIP_CALL( SCIPaddRealParam(scip, "constraints/" CONSHDLR_NAME "/sparsifymaxloss",
5220  "maximal loss in cut efficacy by sparsification",
5221  &conshdlrdata->sparsifymaxloss, TRUE, 0.2, 0.0, 1.0, NULL, NULL) );
5222 
5223  SCIP_CALL( SCIPaddRealParam(scip, "constraints/" CONSHDLR_NAME "/sparsifynzgrowth",
5224  "growth rate of maximal allowed nonzeros in cuts in sparsification",
5225  &conshdlrdata->sparsifynzgrowth, TRUE, 1.3, 1.000001, SCIPinfinity(scip), NULL, NULL) );
5226 
5227  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/linfeasshift",
5228  "whether to try to make solutions feasible in check by shifting the variable on the right hand side",
5229  &conshdlrdata->linfeasshift, FALSE, TRUE, NULL, NULL) );
5230 
5231  SCIP_CALL( SCIPaddCharParam(scip, "constraints/" CONSHDLR_NAME "/nlpform",
5232  "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)",
5233  &conshdlrdata->nlpform, FALSE, 'a', "aqsed", NULL, NULL) );
5234 
5235  SCIP_CALL( SCIPaddRealParam(scip, "constraints/" CONSHDLR_NAME "/sepanlpmincont",
5236  "minimal required fraction of continuous variables in problem to use solution of NLP relaxation in root for separation",
5237  &conshdlrdata->sepanlpmincont, FALSE, 1.0, 0.0, 2.0, NULL, NULL) );
5238 
5239  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/enfocutsremovable",
5240  "are cuts added during enforcement removable from the LP in the same node?",
5241  &conshdlrdata->enfocutsremovable, TRUE, FALSE, NULL, NULL) );
5242 
5243  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/generalsocupgrade",
5244  "try to upgrade more general quadratics to soc?",
5245  &conshdlrdata->generalsocupg, TRUE, TRUE, NULL, NULL) );
5246 
5247  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/disaggregate",
5248  "try to completely disaggregate soc?",
5249  &conshdlrdata->disaggregate, TRUE, TRUE, NULL, NULL) );
5250 
5251  return SCIP_OKAY;
5252 }
5253 
5254 /** creates and captures a second order cone constraint
5255  *
5256  * @note the constraint gets captured, hence at one point you have to release it using the method SCIPreleaseCons()
5257  */
5259  SCIP* scip, /**< SCIP data structure */
5260  SCIP_CONS** cons, /**< pointer to hold the created constraint */
5261  const char* name, /**< name of constraint */
5262  int nvars, /**< number of variables on left hand side of constraint (n) */
5263  SCIP_VAR** vars, /**< array with variables on left hand side (x_i) */
5264  SCIP_Real* coefs, /**< array with coefficients of left hand side variables (alpha_i), or NULL if all 1.0 */
5265  SCIP_Real* offsets, /**< array with offsets of variables (beta_i), or NULL if all 0.0 */
5266  SCIP_Real constant, /**< constant on left hand side (gamma) */
5267  SCIP_VAR* rhsvar, /**< variable on right hand side of constraint (x_{n+1}) */
5268  SCIP_Real rhscoeff, /**< coefficient of variable on right hand side (alpha_{n+1}) */
5269  SCIP_Real rhsoffset, /**< offset of variable on right hand side (beta_{n+1}) */
5270  SCIP_Bool initial, /**< should the LP relaxation of constraint be in the initial LP?
5271  * Usually set to TRUE. Set to FALSE for 'lazy constraints'. */
5272  SCIP_Bool separate, /**< should the constraint be separated during LP processing?
5273  * Usually set to TRUE. */
5274  SCIP_Bool enforce, /**< should the constraint be enforced during node processing?
5275  * TRUE for model constraints, FALSE for additional, redundant constraints. */
5276  SCIP_Bool check, /**< should the constraint be checked for feasibility?
5277  * TRUE for model constraints, FALSE for additional, redundant constraints. */
5278  SCIP_Bool propagate, /**< should the constraint be propagated during node processing?
5279  * Usually set to TRUE. */
5280  SCIP_Bool local, /**< is constraint only valid locally?
5281  * Usually set to FALSE. Has to be set to TRUE, e.g., for branching constraints. */
5282  SCIP_Bool modifiable, /**< is constraint modifiable (subject to column generation)?
5283  * Usually set to FALSE. In column generation applications, set to TRUE if pricing
5284  * adds coefficients to this constraint. */
5285  SCIP_Bool dynamic, /**< is constraint subject to aging?
5286  * Usually set to FALSE. Set to TRUE for own cuts which
5287  * are separated as constraints. */
5288  SCIP_Bool removable /**< should the relaxation be removed from the LP due to aging or cleanup?
5289  * Usually set to FALSE. Set to TRUE for 'lazy constraints' and 'user cuts'. */
5290  )
5291 {
5292  SCIP_CONSHDLR* conshdlr;
5293  SCIP_CONSDATA* consdata;
5294  int i;
5295 
5296  assert(scip != NULL);
5297  assert(cons != NULL);
5298  assert(modifiable == FALSE); /* we do not support column generation */
5299 
5300  /* find the soc constraint handler */
5301  conshdlr = SCIPfindConshdlr(scip, CONSHDLR_NAME);
5302  if( conshdlr == NULL )
5303  {
5304  SCIPerrorMessage("SOC constraint handler not found\n");
5305  return SCIP_PLUGINNOTFOUND;
5306  }
5307 
5308  assert(vars != NULL);
5309  assert(nvars >= 1);
5310  assert(constant >= 0.0);
5311  assert(!SCIPisInfinity(scip, ABS(rhsoffset)));
5312  assert(!SCIPisInfinity(scip, constant));
5313  assert(rhsvar == NULL || rhscoeff <= 0.0 || SCIPisGE(scip, local ? SCIPcomputeVarLbLocal(scip, rhsvar) : SCIPcomputeVarLbGlobal(scip, rhsvar), -rhsoffset));
5314  assert(rhsvar == NULL || rhscoeff >= 0.0 || SCIPisLE(scip, local ? SCIPcomputeVarUbLocal(scip, rhsvar) : SCIPcomputeVarUbGlobal(scip, rhsvar), -rhsoffset));
5315 
5316  /* create constraint data */
5317  SCIP_CALL( SCIPallocBlockMemory(scip, &consdata) );
5318 
5319  consdata->nvars = nvars;
5320  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &consdata->vars, vars, nvars) );
5321  for( i = 0; i < nvars; ++i )
5322  {
5323  SCIP_CALL( SCIPcaptureVar(scip, vars[i]) );
5324  }
5325 
5326  if( coefs != NULL )
5327  {
5328  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &consdata->coefs, coefs, nvars) );
5329  for( i = 0; i < nvars; ++i )
5330  {
5331  if( consdata->coefs[i] < 0.0 )
5332  consdata->coefs[i] = -consdata->coefs[i];
5333  }
5334  }
5335  else
5336  {
5337  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->coefs, nvars) );
5338  for( i = 0; i < nvars; ++i )
5339  consdata->coefs[i] = 1.0;
5340  }
5341 
5342  if( offsets != NULL )
5343  {
5344  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &consdata->offsets, offsets, nvars) );
5345  }
5346  else
5347  {
5348  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->offsets, nvars) );
5349  BMSclearMemoryArray(consdata->offsets, nvars);
5350  }
5351 
5352  consdata->constant = constant;
5353  consdata->rhsvar = rhsvar;
5354  consdata->rhscoeff = rhscoeff;
5355  consdata->rhsoffset = rhsoffset;
5356 
5357  if( rhsvar != NULL )
5358  {
5359  SCIP_CALL( SCIPcaptureVar(scip, rhsvar) );
5360  }
5361 
5362  consdata->nlrow = NULL;
5363 
5364  consdata->lhsbndchgeventdata = NULL;
5365  consdata->isapproxadded = FALSE;
5366 
5367  /* create constraint */
5368  SCIP_CALL( SCIPcreateCons(scip, cons, name, conshdlr, consdata, initial, separate, enforce, check, propagate,
5369  local, modifiable, dynamic, removable, FALSE) );
5370 
5371  if( SCIPisTransformed(scip) )
5372  {
5373  SCIP_CONSHDLRDATA* conshdlrdata = SCIPconshdlrGetData(conshdlr);
5374  assert(conshdlrdata != NULL);
5375 
5376  SCIP_CALL( catchVarEvents(scip, conshdlrdata->eventhdlr, *cons) );
5377  }
5378 
5379  return SCIP_OKAY;
5380 }
5381 
5382 /** creates and captures a second order cone constraint with all its constraint flags
5383  * set to their default values
5384  *
5385  * @note the constraint gets captured, hence at one point you have to release it using the method SCIPreleaseCons()
5386  */
5388  SCIP* scip, /**< SCIP data structure */
5389  SCIP_CONS** cons, /**< pointer to hold the created constraint */
5390  const char* name, /**< name of constraint */
5391  int nvars, /**< number of variables on left hand side of constraint (n) */
5392  SCIP_VAR** vars, /**< array with variables on left hand side (x_i) */
5393  SCIP_Real* coefs, /**< array with coefficients of left hand side variables (alpha_i), or NULL if all 1.0 */
5394  SCIP_Real* offsets, /**< array with offsets of variables (beta_i), or NULL if all 0.0 */
5395  SCIP_Real constant, /**< constant on left hand side (gamma) */
5396  SCIP_VAR* rhsvar, /**< variable on right hand side of constraint (x_{n+1}) */
5397  SCIP_Real rhscoeff, /**< coefficient of variable on right hand side (alpha_{n+1}) */
5398  SCIP_Real rhsoffset /**< offset of variable on right hand side (beta_{n+1}) */
5399  )
5400 {
5401  SCIP_CALL( SCIPcreateConsSOC(scip, cons, name, nvars, vars, coefs, offsets, constant,
5402  rhsvar, rhscoeff, rhsoffset,
5403  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE,