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