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-2019 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 
3599  /* skip upgrade if fixed (or (multi)aggregated) variables and still in fast presolving
3600  * probably cons_quadratic did not yet had the chance to remove/replace this variable (see also #2352)
3601  */
3602  if( !SCIPvarIsActive(term->var) && (presoltiming & SCIP_PRESOLTIMING_FAST) != 0 )
3603  goto cleanup;
3604 
3605  a[i*nquadvars + i] = term->sqrcoef;
3606  quadvars[i] = term->var;
3607  }
3608 
3609  /* set lower triangular entries of A corresponding to bilinear terms */
3610  for( i = 0; i < nbilinterms; ++i )
3611  {
3612  int idx1;
3613  int idx2;
3614 
3615  bilinterm = &SCIPgetBilinTermsQuadratic(scip, cons)[i];
3616 
3617  SCIP_CALL( SCIPfindQuadVarTermQuadratic(scip, cons, bilinterm->var1, &idx1) );
3618  SCIP_CALL( SCIPfindQuadVarTermQuadratic(scip, cons, bilinterm->var2, &idx2) );
3619 
3620  assert(idx1 >= 0);
3621  assert(idx2 >= 0);
3622  assert(idx1 != idx2);
3623 
3624  a[MIN(idx1,idx2) * nquadvars + MAX(idx1,idx2)] = bilinterm->coef / 2.0;
3625  }
3626 
3627  /* compute eigenvalues and vectors, A = PDP^t
3628  * note: a stores P^t
3629  */
3630  if( LapackDsyev(TRUE, nquadvars, a, eigvals) != SCIP_OKAY )
3631  {
3632  SCIPdebugMsg(scip, "Failed to compute eigenvalues and eigenvectors for constraint <%s>.\n", SCIPconsGetName(cons));
3633  goto cleanup;
3634  }
3635 
3636  /* set small eigenvalues to 0 and compute b*P */
3637  nneg = 0;
3638  npos = 0;
3639  for( i = 0; i < nquadvars; ++i )
3640  {
3641  for( j = 0; j < nquadvars; ++j )
3642  {
3643  term = &SCIPgetQuadVarTermsQuadratic(scip, cons)[j];
3644  bp[i] += term->lincoef * a[i * nquadvars + j];
3645  }
3646 
3647  if( SCIPisZero(scip, eigvals[i]) )
3648  {
3649  /* if there is a purely linear variable, the constraint can't be written as a SOC */
3650  if( !SCIPisZero(scip, bp[i]) )
3651  {
3652  goto cleanup;
3653  }
3654 
3655  bp[i] = 0.0;
3656  eigvals[i] = 0.0;
3657  }
3658  else if( eigvals[i] > 0.0 )
3659  npos++;
3660  else
3661  nneg++;
3662  }
3663 
3664  /* currently, a proper SOC constraint needs at least 3 variables */
3665  if( npos + nneg < 3 )
3666  goto cleanup;
3667 
3668  /* determine whether rhs or lhs of cons is potentially SOC, if any */
3669  rhsissoc = nneg == 1 && !SCIPisInfinity(scip, SCIPgetRhsQuadratic(scip, cons));
3670  lhsissoc = npos == 1 && !SCIPisInfinity(scip, -SCIPgetLhsQuadratic(scip, cons));
3671 
3672  /* if non is potentially SOC, stop */
3673  if( !rhsissoc && !lhsissoc )
3674  goto cleanup;
3675 
3676  assert(rhsissoc != lhsissoc);
3677 
3678  if( lhsissoc )
3679  {
3680  /* lhs is potentially SOC, change signs */
3681  lhsconstant = SCIPgetLhsQuadratic(scip, cons);
3682 
3683  for( i = 0; i < nquadvars; ++i )
3684  {
3685  eigvals[i] = -eigvals[i];
3686  bp[i] = -bp[i];
3687  }
3688  }
3689  else
3690  {
3691  lhsconstant = -SCIPgetRhsQuadratic(scip, cons);
3692  }
3693 
3694  /* we have lhsconstant + x^t A x + b x <= 0 and A has a single negative eigenvalue; try to build soc */
3695  negeigpos = -1;
3696  rhsvarfound = FALSE;
3697  for( i = 0; i < nquadvars; ++i )
3698  {
3699  if( SCIPisZero(scip, eigvals[i]) )
3700  continue;
3701 
3702  if( eigvals[i] > 0.0 )
3703  {
3704  lhscoefs[lhscount] = sqrt(eigvals[i]) * UPGSCALE;
3705  lhsoffsets[lhscount] = bp[i] / (2 * eigvals[i]);
3706  lhsconstant -= bp[i] * bp[i] / (4 * eigvals[i]);
3707 
3708  ++lhscount;
3709  }
3710  else
3711  {
3712  /* should enter here only once */
3713  assert(rhsvarfound == FALSE);
3714 
3715  rhsoffset = bp[i] / (2 * eigvals[i]);
3716 
3717  /* the constraint can only be a soc if the resulting rhs var does not change var; the rhs var is going to be a
3718  * multiaggregated variable, so estimate its bounds
3719  */
3720  rhsvarlb = 0.0;
3721  rhsvarub = 0.0;
3722  for( j = 0; j < nquadvars; ++j )
3723  {
3724  SCIP_Real aux;
3725 
3726  if( a[i * nquadvars + j] > 0.0 )
3727  {
3728  aux = SCIPcomputeVarLbGlobal(scip, quadvars[j]);
3729  assert(!SCIPisInfinity(scip, aux));
3730  }
3731  else
3732  {
3733  aux = SCIPcomputeVarUbGlobal(scip, quadvars[j]);
3734  assert(!SCIPisInfinity(scip, -aux));
3735  }
3736 
3737  if( SCIPisInfinity(scip, aux) || SCIPisInfinity(scip, -aux) )
3738  {
3739  rhsvarlb = -SCIPinfinity(scip);
3740  break;
3741  }
3742  else
3743  rhsvarlb += a[i * nquadvars + j] * aux;
3744  }
3745  rhsvarlb += rhsoffset;
3746 
3747  for( j = 0; j < nquadvars; ++j )
3748  {
3749  SCIP_Real aux;
3750 
3751  if( a[i * nquadvars + j] > 0.0 )
3752  {
3753  aux = SCIPcomputeVarUbGlobal(scip, quadvars[j]);
3754  assert(!SCIPisInfinity(scip, -aux));
3755  }
3756  else
3757  {
3758  aux = SCIPcomputeVarLbGlobal(scip, quadvars[j]);
3759  assert(!SCIPisInfinity(scip, aux));
3760  }
3761 
3762  if( SCIPisInfinity(scip, aux) || SCIPisInfinity(scip, -aux) )
3763  {
3764  rhsvarub = SCIPinfinity(scip);
3765  break;
3766  }
3767  else
3768  rhsvarub += a[i * nquadvars + j] * aux;
3769  }
3770  rhsvarub += rhsoffset;
3771 
3772  /* since we are just interested in obtaining an interval that contains the real bounds and is tight enough so
3773  * that we can identify that the rhsvar does not change sign, we swap the bounds in case of numerical troubles
3774  */
3775  if( rhsvarub < rhsvarlb )
3776  {
3777  assert(SCIPisEQ(scip, rhsvarub, rhsvarlb));
3778  SCIPswapReals(&rhsvarub, &rhsvarlb);
3779  }
3780 
3781  /* check whether rhsvar changes sign */
3782  if( SCIPisGE(scip, rhsvarlb, 0.0) || SCIPisLE(scip, rhsvarub, 0.0) )
3783  {
3784  rhsvarfound = TRUE;
3785  negeigpos = i;
3786  lhsconstant -= rhsoffset * rhsoffset * eigvals[i];
3787  rhscoef = SCIPisLE(scip, rhsvarub, 0.0) ? -sqrt(-eigvals[i]) : sqrt(-eigvals[i]);
3788  }
3789  else
3790  {
3791  SCIPdebugMsg(scip, "Failed because rhsvar [%g, %g] changes sign.\n", rhsvarlb, rhsvarub);
3792  rhsvarfound = FALSE;
3793  break;
3794  }
3795  }
3796  }
3797 
3798  if( rhsvarfound && lhscount >= 2 && !SCIPisNegative(scip, lhsconstant) )
3799  {
3800  /* found SOC constraint, so upgrade to SOC constraint(s) (below) and relax right hand side */
3801  SCIPdebugMsg(scip, "found right hand side of constraint <%s> to be SOC\n", SCIPconsGetName(cons));
3802 
3803  /* check if upgdconss is long enough to store upgrade constraints:
3804  * we need two if we will have a quadratic constraint for the left hand side left */
3805  if( rhsissoc )
3806  *nupgdconss = SCIPisInfinity(scip, -SCIPgetLhsQuadratic(scip, cons)) ? 1 : 2;
3807  else
3808  *nupgdconss = SCIPisInfinity(scip, SCIPgetRhsQuadratic(scip, cons)) ? 1 : 2;
3809  if( *nupgdconss > upgdconsssize )
3810  {
3811  /* signal that we need more memory and return */
3812  *nupgdconss = -*nupgdconss;
3813  goto cleanup;
3814  }
3815 
3816  /* create lhs and rhs vars */
3817  lhsnvars = 0;
3818  for( i = 0; i < nquadvars; ++i )
3819  {
3820  if( eigvals[i] <= 0.0 )
3821  continue;
3822 
3823  SCIP_CALL( SCIPcreateVarBasic(scip, &lhsvars[lhsnvars], NULL,
3824  -SCIPinfinity(scip), SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
3825  SCIP_CALL( SCIPaddVar(scip, lhsvars[lhsnvars]) );
3826 
3827  SCIP_CALL( SCIPmultiaggregateVar(scip, lhsvars[lhsnvars], nquadvars, quadvars, &(a[i * nquadvars]),
3828  lhsoffsets[lhsnvars], &infeas, &success) );
3829 
3830  if( infeas || !success )
3831  {
3832  SCIPdebugMsg(scip, "Problem with aggregation while trying to upgrade <%s>.\n", SCIPconsGetName(cons) );
3833 
3834  /* release all created vars so far */
3835  for( j = 0; j <= lhsnvars; ++j )
3836  SCIP_CALL( SCIPreleaseVar(scip, &lhsvars[j]) );
3837 
3838  *nupgdconss = 0;
3839  goto cleanup;
3840  }
3841  lhsnvars++;
3842  }
3843  assert(lhsnvars == lhscount);
3844  assert(negeigpos >= 0);
3845 
3846  SCIP_CALL( SCIPcreateVarBasic(scip, &rhsvar, NULL,
3847  rhsvarlb, rhsvarub, 0.0, SCIP_VARTYPE_CONTINUOUS) );
3848  SCIP_CALL( SCIPaddVar(scip, rhsvar) );
3849  SCIP_CALL( SCIPmultiaggregateVar(scip, rhsvar, nquadvars, quadvars, &(a[negeigpos * nquadvars]),
3850  rhsoffset, &infeas, &success) );
3851 
3852  if( infeas || !success )
3853  {
3854  SCIPdebugMsg(scip, "Problem with aggregation while trying to upgrade <%s>.\n", SCIPconsGetName(cons) );
3855 
3856  /* release all created vars */
3857  SCIP_CALL( SCIPreleaseVar(scip, &rhsvar) );
3858  for( j = 0; j < lhsnvars; ++j )
3859  {
3860  SCIP_CALL( SCIPreleaseVar(scip, &lhsvars[j]) );
3861  }
3862 
3863  *nupgdconss = 0;
3864  goto cleanup;
3865  }
3866 
3867  SCIP_CALL( SCIPcreateConsSOC(scip, &upgdconss[0], SCIPconsGetName(cons),
3868  lhscount, lhsvars, lhscoefs, NULL, MAX(lhsconstant, 0.0) * UPGSCALE * UPGSCALE,
3869  rhsvar, rhscoef * UPGSCALE, 0.0,
3873  SCIPdebugPrintCons(scip, upgdconss[0], NULL);
3874 
3875  /* release created variables */
3876  SCIP_CALL( SCIPreleaseVar(scip, &rhsvar) );
3877  for( i = 0; i < lhsnvars; ++i )
3878  {
3879  SCIP_CALL( SCIPreleaseVar(scip, &lhsvars[i]) );
3880  }
3881 
3882  /* create constraint that is equal to cons except that rhs/lhs is now +/-infinity */
3883  if( rhsissoc && !SCIPisInfinity(scip, -SCIPgetLhsQuadratic(scip, cons)) )
3884  {
3885  SCIPdebugMsg(scip, "rhs is soc, keep quadratic\n");
3886  SCIP_CALL( SCIPcreateConsQuadratic2(scip, &upgdconss[1], SCIPconsGetName(cons),
3890  SCIPgetLhsQuadratic(scip, cons), SCIPinfinity(scip),
3894  }
3895  else if( lhsissoc && !SCIPisInfinity(scip, SCIPgetRhsQuadratic(scip,cons)) )
3896  {
3897  SCIPdebugMsg(scip, "lhs is soc, keep quadratic\n");
3898  SCIP_CALL( SCIPcreateConsQuadratic2(scip, &upgdconss[1], SCIPconsGetName(cons),
3902  -SCIPinfinity(scip), SCIPgetRhsQuadratic(scip, cons),
3906  }
3907  }
3908 #ifdef SCIP_DEBUG
3909  else
3910  {
3911  if( lhscount < 2 )
3912  {
3913  SCIPdebugMsg(scip, "Failed because there are not enough lhsvars (%d)\n", lhscount);
3914  }
3915  if( SCIPisNegative(scip, lhsconstant) )
3916  {
3917  SCIPdebugMsg(scip, "Failed because lhsconstant is negative (%g)\n", lhsconstant);
3918  }
3919  }
3920 #endif
3921 
3922  cleanup:
3923  SCIPfreeBufferArrayNull(scip, &eigvals);
3924  SCIPfreeBufferArrayNull(scip, &quadvars);
3925  SCIPfreeBufferArrayNull(scip, &bp);
3926  SCIPfreeBufferArrayNull(scip, &a);
3927  SCIPfreeBufferArray(scip, &lhsoffsets);
3928  SCIPfreeBufferArray(scip, &lhscoefs);
3929  SCIPfreeBufferArray(scip, &lhsvars);
3930 
3931  return SCIP_OKAY;
3932 } /*lint !e715*/
3933 #endif
3934 
3935 /*
3936  * Callback methods of constraint handler
3937  */
3938 
3939 /** copy method for constraint handler plugins (called when SCIP copies plugins) */
3940 static
3941 SCIP_DECL_CONSHDLRCOPY(conshdlrCopySOC)
3942 { /*lint --e{715}*/
3943  assert(scip != NULL);
3944  assert(conshdlr != NULL);
3945  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
3946 
3947  /* call inclusion method of constraint handler */
3949 
3950  *valid = TRUE;
3951 
3952  return SCIP_OKAY;
3953 }
3954 
3955 /** destructor of constraint handler to free constraint handler data (called when SCIP is exiting) */
3956 static
3957 SCIP_DECL_CONSFREE(consFreeSOC)
3958 {
3959  SCIP_CONSHDLRDATA* conshdlrdata;
3960 
3961  assert( scip != NULL );
3962  assert( conshdlr != NULL );
3963  assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
3964 
3965  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3966  assert(conshdlrdata != NULL);
3967 
3968  SCIPfreeBlockMemory(scip, &conshdlrdata);
3969 
3970  return SCIP_OKAY;
3971 }
3972 
3973 
3974 /** initialization method of constraint handler (called after problem was transformed) */
3975 static
3976 SCIP_DECL_CONSINIT(consInitSOC)
3977 { /*lint --e{715}*/
3978  SCIP_CONSHDLRDATA* conshdlrdata;
3979  int c;
3980 
3981  assert(scip != NULL);
3982  assert(conshdlr != NULL);
3983  assert(conss != NULL || nconss == 0);
3984 
3985  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3986  assert(conshdlrdata != NULL);
3987 
3988  conshdlrdata->subnlpheur = SCIPfindHeur(scip, "subnlp");
3989  conshdlrdata->trysolheur = SCIPfindHeur(scip, "trysol");
3990  conshdlrdata->haveexprint = (strcmp(SCIPexprintGetName(), "NONE") != 0);
3991 
3992  /* mark constraints for propagation */
3993  for( c = 0; c < nconss; ++c )
3994  {
3995  SCIP_CALL( SCIPmarkConsPropagate(scip, conss[c]) ); /*lint !e613*/
3996  }
3997 
3998  return SCIP_OKAY;
3999 }
4000 
4001 
4002 /** deinitialization method of constraint handler (called before transformed problem is freed) */
4003 static
4004 SCIP_DECL_CONSEXIT(consExitSOC)
4005 { /*lint --e{715}*/
4006  SCIP_CONSHDLRDATA* conshdlrdata;
4007 
4008  assert(scip != NULL);
4009  assert(conshdlr != NULL);
4010 
4011  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4012  assert(conshdlrdata != NULL);
4013 
4014  conshdlrdata->subnlpheur = NULL;
4015  conshdlrdata->trysolheur = NULL;
4016  conshdlrdata->haveexprint = FALSE;
4017 
4018  return SCIP_OKAY;
4019 }
4020 
4021 /** presolving deinitialization method of constraint handler (called after presolving has been finished) */
4022 static
4023 SCIP_DECL_CONSEXITPRE(consExitpreSOC)
4024 { /*lint --e{715}*/
4025  int c;
4026 
4027  assert(scip != NULL);
4028  assert(conshdlr != NULL);
4029  assert(conss != NULL || nconss == 0);
4030 
4031  /* tell SCIP that we have something nonlinear */
4032  for( c = 0; c < nconss; ++c )
4033  {
4034  if( SCIPconsIsAdded(conss[c]) ) /*lint !e613*/
4035  {
4036  SCIPenableNLP(scip);
4037  break;
4038  }
4039  }
4040 
4041  return SCIP_OKAY;
4042 }
4043 
4044 
4045 /** solving process initialization method of constraint handler (called when branch and bound process is about to begin) */
4046 static
4047 SCIP_DECL_CONSINITSOL(consInitsolSOC)
4048 {
4049  SCIP_CONSHDLRDATA* conshdlrdata;
4050  SCIP_CONSDATA* consdata;
4051  int c;
4052 
4053  assert(scip != NULL);
4054  assert(conshdlr != NULL);
4055 
4056  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4057  assert(conshdlrdata != NULL);
4058  assert(conshdlrdata->eventhdlr != NULL);
4059 
4060  /* add nlrow representation to NLP, if NLP has been enabled */
4061  if( SCIPisNLPConstructed(scip) )
4062  {
4063  for( c = 0; c < nconss; ++c )
4064  {
4065  if( SCIPconsIsEnabled(conss[c]) )
4066  {
4067  consdata = SCIPconsGetData(conss[c]);
4068  assert(consdata != NULL);
4069 
4070  if( consdata->nlrow == NULL )
4071  {
4072  SCIP_CALL( createNlRow(scip, conshdlr, conss[c]) );
4073  assert(consdata->nlrow != NULL);
4074  }
4075  SCIP_CALL( SCIPaddNlRow(scip, consdata->nlrow) );
4076  }
4077  }
4078  }
4079 
4080  conshdlrdata->newsoleventfilterpos = -1;
4081  if( nconss != 0 )
4082  {
4083  SCIP_EVENTHDLR* eventhdlr;
4084 
4085  eventhdlr = SCIPfindEventhdlr(scip, CONSHDLR_NAME"_newsolution");
4086  assert(eventhdlr != NULL);
4087 
4088  SCIP_CALL( SCIPcatchEvent(scip, SCIP_EVENTTYPE_SOLFOUND, eventhdlr, (SCIP_EVENTDATA*)conshdlr, &conshdlrdata->newsoleventfilterpos) );
4089  }
4090 
4091  /* reset flags and counters */
4092  conshdlrdata->sepanlp = FALSE;
4093  conshdlrdata->lastenfonode = NULL;
4094  conshdlrdata->nenforounds = 0;
4095 
4096  return SCIP_OKAY;
4097 }
4098 
4099 
4100 /** solving process deinitialization method of constraint handler (called before branch and bound process data is freed) */
4101 static
4102 SCIP_DECL_CONSEXITSOL(consExitsolSOC)
4104  SCIP_CONSHDLRDATA* conshdlrdata;
4105  SCIP_CONSDATA* consdata;
4106  int c;
4107 
4108  assert(scip != NULL);
4109  assert(conshdlr != NULL);
4110 
4111  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4112  assert(conshdlrdata != NULL);
4113  assert(conshdlrdata->eventhdlr != NULL);
4114 
4115  if( conshdlrdata->newsoleventfilterpos >= 0 )
4116  {
4117  SCIP_EVENTHDLR* eventhdlr;
4118 
4119  eventhdlr = SCIPfindEventhdlr(scip, CONSHDLR_NAME"_newsolution");
4120  assert(eventhdlr != NULL);
4121 
4122  SCIP_CALL( SCIPdropEvent(scip, SCIP_EVENTTYPE_SOLFOUND, eventhdlr, (SCIP_EVENTDATA*)conshdlr, conshdlrdata->newsoleventfilterpos) );
4123  conshdlrdata->newsoleventfilterpos = -1;
4124  }
4125 
4126  for( c = 0; c < nconss; ++c )
4127  {
4128  consdata = SCIPconsGetData(conss[c]);
4129  assert(consdata != NULL);
4130 
4131  /* free nonlinear row representation */
4132  if( consdata->nlrow != NULL )
4133  {
4134  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
4135  }
4136  }
4137 
4138  return SCIP_OKAY;
4139 } /*lint !e715*/
4140 
4141 
4142 /** frees specific constraint data */
4143 static
4144 SCIP_DECL_CONSDELETE(consDeleteSOC)
4146  int i;
4147 
4148  assert(scip != NULL);
4149  assert(conshdlr != NULL);
4150  assert(cons != NULL);
4151  assert(consdata != NULL);
4152  assert(*consdata != NULL);
4153  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
4154 
4155  SCIPdebugMsg(scip, "Deleting SOC constraint <%s>.\n", SCIPconsGetName(cons) );
4156 
4157  if( SCIPconsIsTransformed(cons) )
4158  {
4159  SCIP_CONSHDLRDATA* conshdlrdata;
4160 
4161  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4162  assert(conshdlrdata != NULL);
4163 
4164  SCIP_CALL( dropVarEvents(scip, conshdlrdata->eventhdlr, cons) );
4165  }
4166 
4167  for( i = 0; i < (*consdata)->nvars; ++i )
4168  {
4169  SCIP_CALL( SCIPreleaseVar(scip, &(*consdata)->vars[i]) );
4170  }
4171 
4172  SCIPfreeBlockMemoryArray(scip, &(*consdata)->vars, (*consdata)->nvars);
4173  SCIPfreeBlockMemoryArray(scip, &(*consdata)->coefs, (*consdata)->nvars);
4174  SCIPfreeBlockMemoryArray(scip, &(*consdata)->offsets, (*consdata)->nvars);
4175  assert((*consdata)->lhsbndchgeventdata == NULL);
4176 
4177  if( (*consdata)->rhsvar != NULL )
4178  {
4179  SCIP_CALL( SCIPreleaseVar(scip, &(*consdata)->rhsvar) );
4180  }
4181 
4182  /* free nonlinear row representation
4183  * normally released in exitsol, but constraint may be deleted early (e.g., if found redundant)
4184  */
4185  if( (*consdata)->nlrow != NULL )
4186  {
4187  SCIP_CALL( SCIPreleaseNlRow(scip, &(*consdata)->nlrow) );
4188  }
4189 
4190  SCIPfreeBlockMemory(scip, consdata);
4191 
4192  return SCIP_OKAY;
4193 }
4194 
4195 
4196 /** transforms constraint data into data belonging to the transformed problem */
4197 static
4198 SCIP_DECL_CONSTRANS(consTransSOC)
4199 {
4200  SCIP_CONSDATA* consdata;
4201  SCIP_CONSHDLRDATA* conshdlrdata;
4202  SCIP_CONSDATA* sourcedata;
4203  char s[SCIP_MAXSTRLEN];
4204  int i;
4205 
4206  assert(scip != NULL);
4207  assert(conshdlr != NULL);
4208  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
4209  assert(sourcecons != NULL);
4210  assert(targetcons != NULL);
4211 
4212  /* get constraint handler data */
4213  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4214  assert(conshdlrdata != NULL);
4215 
4216  SCIPdebugMsg(scip, "Transforming SOC constraint: <%s>.\n", SCIPconsGetName(sourcecons) );
4217 
4218  /* get data of original constraint */
4219  sourcedata = SCIPconsGetData(sourcecons);
4220  assert(sourcedata != NULL);
4221  assert(sourcedata->vars != NULL);
4222  assert(sourcedata->coefs != NULL);
4223  assert(sourcedata->offsets != NULL);
4224 
4225  /* create constraint data */
4226  SCIP_CALL( SCIPallocBlockMemory(scip, &consdata) );
4227 
4228  consdata->nvars = sourcedata->nvars;
4229  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->vars, consdata->nvars) );
4230  SCIP_CALL( SCIPgetTransformedVars(scip, consdata->nvars, sourcedata->vars, consdata->vars) );
4231  for( i = 0; i < consdata->nvars; ++i )
4232  {
4233  SCIP_CALL( SCIPcaptureVar(scip, consdata->vars[i]) );
4234  }
4235 
4236  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &consdata->coefs, sourcedata->coefs, consdata->nvars) );
4237  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &consdata->offsets, sourcedata->offsets, consdata->nvars) );
4238  consdata->constant = sourcedata->constant;
4239 
4240  SCIP_CALL( SCIPgetTransformedVar(scip, sourcedata->rhsvar, &consdata->rhsvar) );
4241  consdata->rhsoffset = sourcedata->rhsoffset;
4242  consdata->rhscoeff = sourcedata->rhscoeff;
4243 
4244  SCIP_CALL( SCIPcaptureVar(scip, consdata->rhsvar) );
4245 
4246  consdata->nlrow = NULL;
4247  consdata->lhsbndchgeventdata = NULL;
4248  consdata->isapproxadded = FALSE;
4249 
4250  /* create transformed constraint with the same flags */
4251  (void) SCIPsnprintf(s, SCIP_MAXSTRLEN, "t_%s", SCIPconsGetName(sourcecons));
4252  SCIP_CALL( SCIPcreateCons(scip, targetcons, s, conshdlr, consdata,
4253  SCIPconsIsInitial(sourcecons), SCIPconsIsSeparated(sourcecons),
4254  SCIPconsIsEnforced(sourcecons), SCIPconsIsChecked(sourcecons),
4255  SCIPconsIsPropagated(sourcecons), SCIPconsIsLocal(sourcecons),
4256  SCIPconsIsModifiable(sourcecons), SCIPconsIsDynamic(sourcecons),
4257  SCIPconsIsRemovable(sourcecons), SCIPconsIsStickingAtNode(sourcecons)) );
4258 
4259  SCIP_CALL( catchVarEvents(scip, conshdlrdata->eventhdlr, *targetcons) );
4260 
4261  return SCIP_OKAY;
4262 }
4263 
4264 /** separation method of constraint handler for LP solutions */
4265 static
4266 SCIP_DECL_CONSSEPALP(consSepalpSOC)
4267 {
4268  SCIP_CONSHDLRDATA* conshdlrdata;
4269  SCIP_CONS* maxviolcon;
4270  SCIP_Bool sepasuccess;
4271  SCIP_Bool cutoff;
4272 
4273  assert(scip != NULL);
4274  assert(conshdlr != NULL);
4275  assert(conss != NULL || nconss == 0);
4276  assert(result != NULL);
4277 
4278  *result = SCIP_DIDNOTFIND;
4279 
4280  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4281  assert(conshdlrdata != NULL);
4282 
4283  SCIP_CALL( computeViolations(scip, conshdlr, conss, nconss, NULL, &maxviolcon) );
4284  if( maxviolcon == NULL )
4285  return SCIP_OKAY;
4286 
4287  /* at root, check if we want to solve the NLP relaxation and use its solutions as reference point
4288  * if there is something convex, then linearizing in the solution of the NLP relaxation can be very useful
4289  */
4290  if( SCIPgetDepth(scip) == 0 && !conshdlrdata->sepanlp &&
4291  (SCIPgetNContVars(scip) >= conshdlrdata->sepanlpmincont * SCIPgetNVars(scip) || (SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_UNBOUNDEDRAY && conshdlrdata->sepanlpmincont <= 1.0)) &&
4292  SCIPisNLPConstructed(scip) && SCIPgetNNlpis(scip) > 0 )
4293  {
4294  SCIP_NLPSOLSTAT solstat;
4295  SCIP_Bool solvednlp;
4296 
4297  solstat = SCIPgetNLPSolstat(scip);
4298  solvednlp = FALSE;
4299  if( solstat == SCIP_NLPSOLSTAT_UNKNOWN )
4300  {
4301  /* NLP is not solved yet, so we do it now and update solstat */
4302 
4303  /* ensure linear conss are in NLP */
4304  if( conshdlrdata->subnlpheur != NULL )
4305  {
4306  SCIP_CALL( SCIPaddLinearConsToNlpHeurSubNlp(scip, conshdlrdata->subnlpheur, TRUE, TRUE) );
4307  }
4308 
4309  /* set LP solution as starting values, if available */
4311  {
4313  }
4314 
4315  /* SCIP_CALL( SCIPsetNLPIntPar(scip, SCIP_NLPPAR_VERBLEVEL, 1) ); */
4316  SCIP_CALL( SCIPsolveNLP(scip) );
4317 
4318  solstat = SCIPgetNLPSolstat(scip);
4319  SCIPdebugMsg(scip, "solved NLP relax, solution status: %d\n", solstat);
4320 
4321  solvednlp = TRUE;
4322  }
4323 
4324  conshdlrdata->sepanlp = TRUE;
4325 
4326  if( solstat == SCIP_NLPSOLSTAT_GLOBINFEASIBLE )
4327  {
4328  SCIPdebugMsg(scip, "NLP relaxation is globally infeasible, thus can cutoff node\n");
4329  *result = SCIP_CUTOFF;
4330  return SCIP_OKAY;
4331  }
4332 
4333  if( solstat <= SCIP_NLPSOLSTAT_FEASIBLE )
4334  {
4335  /* if we have feasible NLP solution, generate linearization cuts there */
4336  SCIP_Bool lpsolseparated;
4337  SCIP_SOL* nlpsol;
4338 
4339  SCIP_CALL( SCIPcreateNLPSol(scip, &nlpsol, NULL) );
4340  assert(nlpsol != NULL);
4341 
4342  /* if we solved the NLP and solution is integral, then pass it to trysol heuristic */
4343  if( solvednlp && conshdlrdata->trysolheur != NULL )
4344  {
4345  int nfracvars = 0;
4346 
4347  if( SCIPgetNBinVars(scip) > 0 || SCIPgetNIntVars(scip) > 0 )
4348  {
4349  SCIP_CALL( SCIPgetNLPFracVars(scip, NULL, NULL, NULL, &nfracvars, NULL) );
4350  }
4351 
4352  if( nfracvars == 0 )
4353  {
4354  SCIP_CALL( SCIPheurPassSolTrySol(scip, conshdlrdata->trysolheur, nlpsol) );
4355  }
4356  }
4357 
4358  SCIP_CALL( addLinearizationCuts(scip, conshdlr, conss, nconss, nlpsol, &lpsolseparated, SCIPgetSepaMinEfficacy(scip), &cutoff) );
4359 
4360  SCIP_CALL( SCIPfreeSol(scip, &nlpsol) );
4361 
4362  if( cutoff )
4363  {
4364  *result = SCIP_CUTOFF;
4365  return SCIP_OKAY;
4366  }
4367 
4368  /* if a cut that separated the LP solution was added, then return, otherwise continue with usual separation in LP solution */
4369  if( lpsolseparated )
4370  {
4371  SCIPdebugMsg(scip, "linearization cuts separate LP solution\n");
4372 
4373  *result = SCIP_SEPARATED;
4374 
4375  return SCIP_OKAY;
4376  }
4377  }
4378  }
4379  /* if we do not want to try solving the NLP, or have no NLP, or have no NLP solver, or solving the NLP failed,
4380  * or separating with NLP solution as reference point failed, then try (again) with LP solution as reference point
4381  */
4382 
4383  SCIP_CALL( separatePoint(scip, conshdlr, conss, nconss, nusefulconss, NULL, FALSE, &cutoff, &sepasuccess) );
4384  if( cutoff )
4385  *result = SCIP_CUTOFF;
4386  else if ( sepasuccess )
4387  *result = SCIP_SEPARATED;
4388 
4389  return SCIP_OKAY;
4390 }
4391 
4392 
4393 /** separation method of constraint handler for arbitrary primal solutions */
4394 static
4395 SCIP_DECL_CONSSEPASOL(consSepasolSOC)
4396 {
4397  SCIP_CONS* maxviolcon;
4398  SCIP_Bool sepasuccess;
4399  SCIP_Bool cutoff;
4400 
4401  assert(scip != NULL);
4402  assert(conss != NULL || nconss == 0);
4403  assert(result != NULL);
4404  assert(sol != NULL);
4405 
4406  *result = SCIP_DIDNOTFIND;
4407 
4408  SCIP_CALL( computeViolations(scip, conshdlr, conss, nconss, sol, &maxviolcon) );
4409  if( maxviolcon == NULL )
4410  return SCIP_OKAY;
4411 
4412  SCIP_CALL( separatePoint(scip, conshdlr, conss, nconss, nusefulconss, sol, FALSE, &cutoff, &sepasuccess) );
4413  if( cutoff )
4414  *result = SCIP_CUTOFF;
4415  else if ( sepasuccess )
4416  *result = SCIP_SEPARATED;
4417 
4418  return SCIP_OKAY;
4419 }
4420 
4421 
4422 /** constraint enforcing method of constraint handler for LP solutions */
4423 static
4424 SCIP_DECL_CONSENFOLP(consEnfolpSOC)
4425 { /*lint --e{715}*/
4426  SCIP_CALL( enforceConstraint(scip, conshdlr, conss, nconss, nusefulconss, NULL, result) );
4427 
4428  return SCIP_OKAY;
4429 }
4430 
4431 
4432 /** constraint enforcing method of constraint handler for relaxation solutions */
4433 static
4434 SCIP_DECL_CONSENFORELAX(consEnforelaxSOC)
4435 { /*lint --e{715}*/
4436  SCIP_CALL( enforceConstraint(scip, conshdlr, conss, nconss, nusefulconss, sol, result) );
4437 
4438  return SCIP_OKAY;
4439 }
4440 
4441 
4442 /** constraint enforcing method of constraint handler for pseudo solutions */
4443 static
4444 SCIP_DECL_CONSENFOPS(consEnfopsSOC)
4446  SCIP_CONS* maxviolcons;
4447 
4448  assert(scip != NULL);
4449  assert(conss != NULL || nconss == 0);
4450  assert(result != NULL);
4451 
4452  SCIP_CALL( computeViolations(scip, conshdlr, conss, nconss, NULL, &maxviolcons) );
4453 
4454  if( maxviolcons == NULL )
4455  *result = SCIP_FEASIBLE;
4456 
4457  *result = SCIP_INFEASIBLE;
4458 
4459  return SCIP_OKAY;
4460 } /*lint !e715*/
4461 
4462 
4463 /** feasibility check method of constraint handler for integral solutions */
4464 static
4465 SCIP_DECL_CONSCHECK(consCheckSOC)
4467  SCIP_CONSHDLRDATA* conshdlrdata;
4468  SCIP_CONSDATA* consdata;
4469  SCIP_Real maxviol;
4470  SCIP_Bool dolinfeasshift;
4471  SCIP_SOL* polishedsol;
4472  int c;
4473 
4474  assert(scip != NULL);
4475  assert(conshdlr != NULL);
4476  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
4477  assert(conss != NULL || nconss == 0);
4478  assert(result != NULL );
4479 
4480  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4481  assert(conshdlrdata != NULL);
4482 
4483  *result = SCIP_FEASIBLE;
4484  maxviol = 0.0;
4485 
4486  dolinfeasshift = conshdlrdata->linfeasshift && (conshdlrdata->trysolheur != NULL);
4487  polishedsol = NULL;
4488 
4489  for( c = 0; c < nconss; ++c )
4490  {
4491  SCIP_CALL( computeViolation(scip, conshdlr, conss[c], sol) ); /*lint !e613*/
4492 
4493  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
4494  assert(consdata != NULL);
4495 
4496  /* if feasible, just continue */
4497  if( !SCIPisGT(scip, consdata->violation, SCIPfeastol(scip)) )
4498  continue;
4499 
4500  *result = SCIP_INFEASIBLE;
4501 
4502  if( consdata->violation > maxviol )
4503  maxviol = consdata->violation;
4504 
4505  if( printreason )
4506  {
4507  SCIP_CALL( SCIPprintCons(scip, conss[c], NULL) ); /*lint !e613*/
4508  SCIPinfoMessage(scip, NULL, ";\n\tviolation: %g\n", consdata->violation);
4509  }
4510 
4511  /* if we do linear feasibility shifting, then try to adjust solution */
4512  if( dolinfeasshift )
4513  {
4514  if( SCIPvarGetStatus(consdata->rhsvar) != SCIP_VARSTATUS_MULTAGGR &&
4515  !SCIPisInfinity(scip, REALABS(consdata->lhsval)) &&
4516  ( (consdata->rhscoeff > 0.0 && SCIPvarMayRoundUp (consdata->rhsvar)) ||
4517  (consdata->rhscoeff < 0.0 && SCIPvarMayRoundDown(consdata->rhsvar)) ) )
4518  {
4519  SCIP_Bool success;
4520 
4521  if( polishedsol == NULL )
4522  {
4523  if( sol != NULL )
4524  {
4525  SCIP_CALL( SCIPcreateSolCopy(scip, &polishedsol, sol) );
4526  }
4527  else
4528  {
4529  SCIP_CALL( SCIPcreateLPSol(scip, &polishedsol, NULL) );
4530  }
4531  SCIP_CALL( SCIPunlinkSol(scip, polishedsol) );
4532  }
4533  SCIP_CALL( polishSolution(scip, conss[c], polishedsol, &success) ); /*lint !e613*/
4534 
4535  /* disable solution polishing if we failed for this constraint */
4536  dolinfeasshift = success;
4537  }
4538  else /* if locks of the variable are bad or rhs is multi-aggregated, disable solution polishing */
4539  {
4540  dolinfeasshift = FALSE;
4541  }
4542  }
4543 
4544  /* if solution polishing is off and there is no NLP heuristic or we just check the LP solution,
4545  * then there is no need to check remaining constraints (NLP heuristic will pick up LP solution anyway) */
4546  if( !dolinfeasshift && (conshdlrdata->subnlpheur == NULL || sol == NULL) && !completely )
4547  break;
4548  }
4549 
4550  /* if we failed to polish solution, clear the solution */
4551  if( !dolinfeasshift && polishedsol != NULL )
4552  {
4553  SCIP_CALL( SCIPfreeSol(scip, &polishedsol) );
4554  }
4555 
4556  if( polishedsol != NULL )
4557  {
4558  assert(*result == SCIP_INFEASIBLE);
4559  SCIP_CALL( SCIPheurPassSolTrySol(scip, conshdlrdata->trysolheur, polishedsol) );
4560  SCIP_CALL( SCIPfreeSol(scip, &polishedsol) );
4561  }
4562  else if( conshdlrdata->subnlpheur != NULL && sol != NULL && *result == SCIP_INFEASIBLE && !SCIPisInfinity(scip, maxviol) )
4563  {
4564  SCIP_CALL( SCIPupdateStartpointHeurSubNlp(scip, conshdlrdata->subnlpheur, sol, maxviol) );
4565  }
4566 
4567  return SCIP_OKAY;
4568 } /*lint !e715*/
4569 
4570 
4571 /** domain propagation method of constraint handler */
4572 static
4573 SCIP_DECL_CONSPROP(consPropSOC)
4575  SCIP_RESULT propresult;
4576  int c;
4577  int nchgbds;
4578  SCIP_Bool redundant;
4579 
4580  assert(scip != NULL);
4581  assert(conss != NULL || ((nconss == 0) && (nmarkedconss == 0)));
4582  assert(result != NULL);
4583 
4584  *result = SCIP_DIDNOTFIND;
4585  nchgbds = 0;
4586 
4587  for( c = 0; c < nmarkedconss && *result != SCIP_CUTOFF; ++c )
4588  {
4589  SCIP_CALL( propagateBounds(scip, conss[c], &propresult, &nchgbds, &redundant) ); /*lint !e613*/
4590  if( propresult != SCIP_DIDNOTFIND && propresult != SCIP_DIDNOTRUN )
4591  *result = propresult;
4592  }
4593 
4594  return SCIP_OKAY;
4595 } /*lint !e715*/
4596 
4597 
4598 /** presolving method of constraint handler */
4599 static
4600 SCIP_DECL_CONSPRESOL(consPresolSOC)
4602  SCIP_CONSHDLRDATA* conshdlrdata;
4603  SCIP_CONSDATA* consdata;
4604  int c;
4605  SCIP_RESULT propresult;
4606  SCIP_Bool iscutoff;
4607  SCIP_Bool isdeleted;
4608 
4609  assert(scip != NULL);
4610  assert(conss != NULL || nconss == 0);
4611  assert(conshdlr != NULL);
4612  assert(result != NULL);
4613 
4614  *result = SCIP_DIDNOTFIND;
4615 
4616  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4617  assert(conshdlrdata != NULL);
4618 
4619  for( c = 0; c < nconss; ++c )
4620  {
4621  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
4622  assert(consdata != NULL);
4623 
4624  SCIP_CALL( presolveRemoveFixedVariables(scip, conshdlr, conss[c], ndelconss, nupgdconss, nchgbds, nfixedvars, &iscutoff, &isdeleted) ); /*lint !e613*/
4625  if( iscutoff )
4626  {
4627  *result = SCIP_CUTOFF;
4628  return SCIP_OKAY;
4629  }
4630  if( isdeleted )
4631  {
4632  /* conss[c] has been deleted */
4633  *result = SCIP_SUCCESS;
4634  continue;
4635  }
4636 
4637  if( conshdlrdata->nauxvars > 0 && !consdata->isapproxadded && consdata->nvars > 1 )
4638  {
4639  SCIP_CALL( presolveCreateOuterApprox(scip, consdata->nvars, consdata->vars, consdata->coefs, consdata->offsets,
4640  consdata->rhsvar, consdata->rhscoeff, consdata->rhscoeff, consdata->constant, SCIPconsGetName(conss[c]), conss[c],
4641  conshdlrdata->nauxvars, conshdlrdata->glineur, naddconss) ); /*lint !e613*/
4642  consdata->isapproxadded = TRUE;
4643  }
4644 
4645  if( (presoltiming & SCIP_PRESOLTIMING_FAST) != 0 )
4646  {
4647  SCIP_Bool redundant;
4648 
4649  SCIP_CALL( propagateBounds(scip, conss[c], &propresult, nchgbds, &redundant) ); /*lint !e613*/
4650  switch( propresult )
4651  {
4652  case SCIP_DIDNOTRUN:
4653  case SCIP_DIDNOTFIND:
4654  break;
4655  case SCIP_REDUCEDDOM:
4656  *result = SCIP_SUCCESS;
4657  break;
4658  case SCIP_CUTOFF:
4659  *result = SCIP_CUTOFF;
4660  SCIPdebugMsg(scip, "infeasible in presolve due to propagation for constraint %s\n", SCIPconsGetName(conss[c])); /*lint !e613*/
4661  return SCIP_OKAY;
4662  default:
4663  SCIPerrorMessage("unexpected result from propagation: %d\n", propresult);
4664  return SCIP_ERROR;
4665  } /*lint !e788*/
4666  if( redundant )
4667  ++*ndelconss;
4668  }
4669 
4670  /* disaggregate each lhs term to a quadratic constraint by using auxiliary variables */
4671  if( conshdlrdata->disaggregate && (presoltiming & SCIP_PRESOLTIMING_EXHAUSTIVE) != 0 )
4672  {
4673  SCIP_Bool success;
4674 
4675  SCIP_CALL( disaggregate(scip, conss[c], consdata, naddconss, ndelconss, &success) ); /*lint !e613*/
4676 
4677  if( success )
4678  {
4679  SCIPdebugMsg(scip, "disaggregated SOC constraint\n");
4680 
4681  /* conss[c] has been deleted */
4682  *result = SCIP_SUCCESS;
4683  continue;
4684  }
4685  }
4686  }
4687 
4688  return SCIP_OKAY;
4689 } /*lint !e715*/
4690 
4691 
4692 /** variable rounding lock method of constraint handler */
4693 static
4694 SCIP_DECL_CONSLOCK(consLockSOC)
4696  SCIP_CONSDATA* consdata;
4697  int i;
4698 
4699  assert(scip != NULL);
4700  assert(conshdlr != NULL);
4701  assert(cons != NULL);
4702  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
4703  assert(locktype == SCIP_LOCKTYPE_MODEL);
4704 
4705  consdata = SCIPconsGetData(cons);
4706  assert(consdata != NULL);
4707 
4708  SCIPdebugMsg(scip, "Locking constraint <%s>.\n", SCIPconsGetName(cons));
4709 
4710  /* Changing variables x_i, i <= n, in both directions can lead to an infeasible solution. */
4711  for( i = 0; i < consdata->nvars; ++i )
4712  {
4713  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->vars[i], locktype, nlockspos + nlocksneg, nlockspos + nlocksneg) );
4714  }
4715 
4716  /* Rounding x_{n+1} up will not violate a solution. */
4717  if( consdata->rhsvar != NULL )
4718  {
4719  SCIP_CALL( SCIPaddVarLocksType(scip, consdata->rhsvar, locktype, consdata->rhscoeff > 0.0 ? nlockspos : nlocksneg, consdata->rhscoeff > 0.0 ? nlocksneg : nlockspos) );
4720  }
4721 
4722  return SCIP_OKAY;
4723 }
4724 
4725 /** constraint display method of constraint handler */
4726 static
4727 SCIP_DECL_CONSPRINT(consPrintSOC)
4728 {
4729  SCIP_CONSDATA* consdata;
4730  int i;
4731 
4732  assert(scip != NULL);
4733  assert(conshdlr != NULL);
4734  assert(cons != NULL);
4735  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
4736 
4737  consdata = SCIPconsGetData(cons);
4738  assert(consdata != NULL);
4739 
4740  SCIPinfoMessage(scip, file, "sqrt( ");
4741  if( consdata->constant != 0.0 )
4742  {
4743  SCIPinfoMessage(scip, file, "%.15g", consdata->constant);
4744  }
4745 
4746  for( i = 0; i < consdata->nvars; ++i )
4747  {
4748  SCIPinfoMessage(scip, file, "+ (%.15g*(", consdata->coefs[i]);
4749  SCIP_CALL( SCIPwriteVarName(scip, file, consdata->vars[i], TRUE) );
4750  SCIPinfoMessage(scip, file, "%+.15g))^2 ", consdata->offsets[i]);
4751  }
4752 
4753  SCIPinfoMessage(scip, file, ") <= ");
4754  if( consdata->rhsvar != NULL )
4755  {
4756  SCIPinfoMessage(scip, file, "%.15g*(", consdata->rhscoeff);
4757  SCIP_CALL( SCIPwriteVarName(scip, file, consdata->rhsvar, TRUE) );
4758  SCIPinfoMessage(scip, file, "%+.15g)", consdata->rhsoffset);
4759  }
4760  else
4761  {
4762  SCIPinfoMessage(scip, file, "%.15g", consdata->rhscoeff*consdata->rhsoffset);
4763  }
4764 
4765  return SCIP_OKAY;
4766 }
4767 
4768 /** constraint copying method of constraint handler */
4769 static
4770 SCIP_DECL_CONSCOPY(consCopySOC)
4772  SCIP_CONSDATA* consdata;
4773  SCIP_VAR** vars;
4774  SCIP_VAR* rhsvar;
4775  int i;
4776 
4777  assert(scip != NULL);
4778  assert(cons != NULL);
4779  assert(sourcescip != NULL);
4780  assert(sourceconshdlr != NULL);
4781  assert(sourcecons != NULL);
4782  assert(varmap != NULL);
4783  assert(valid != NULL);
4784  assert(stickingatnode == FALSE);
4785 
4786  consdata = SCIPconsGetData(sourcecons);
4787  assert(consdata != NULL);
4788 
4789  *valid = TRUE;
4790 
4791  SCIP_CALL( SCIPallocBufferArray(sourcescip, &vars, consdata->nvars) );
4792 
4793  /* map variables to active variables of the target SCIP */
4794  for( i = 0; i < consdata->nvars && *valid; ++i )
4795  {
4796  SCIP_CALL( SCIPgetVarCopy(sourcescip, scip, consdata->vars[i], &vars[i], varmap, consmap, global, valid) );
4797  assert(!(*valid) || vars[i] != NULL);
4798  }
4799 
4800  /* map rhs variable to active variable of the target SCIP */
4801  if( *valid )
4802  {
4803  SCIP_CALL( SCIPgetVarCopy(sourcescip, scip, consdata->rhsvar, &rhsvar, varmap, consmap, global, valid) );
4804  assert(!(*valid) || rhsvar != NULL);
4805  }
4806 
4807  /* only create the target constraint, if all variables could be copied */
4808  if( *valid )
4809  {
4810  SCIP_CALL( SCIPcreateConsSOC(scip, cons, name ? name : SCIPconsGetName(sourcecons),
4811  consdata->nvars, vars, consdata->coefs, consdata->offsets, consdata->constant,
4812  rhsvar, consdata->rhscoeff, consdata->rhsoffset,
4813  initial, separate, enforce, check, propagate, local, modifiable, dynamic, removable) ); /*lint !e644 */
4814  }
4815 
4816  SCIPfreeBufferArray(sourcescip, &vars);
4817 
4818  return SCIP_OKAY;
4819 }
4820 
4821 
4822 /** constraint parsing method of constraint handler */
4823 static
4824 SCIP_DECL_CONSPARSE(consParseSOC)
4825 { /*lint --e{715}*/
4826  SCIP_VAR* var;
4827  SCIP_VAR** vars;
4828  SCIP_Real* coefs;
4829  SCIP_Real* offsets;
4830  int nvars;
4831  int varssize;
4832  SCIP_VAR* rhsvar;
4833  SCIP_Real rhscoef;
4834  SCIP_Real rhsoffset;
4835  SCIP_Real constant;
4836  SCIP_Real coef;
4837  SCIP_Real offset;
4838  char* endptr;
4839 
4840  assert(scip != NULL);
4841  assert(success != NULL);
4842  assert(str != NULL);
4843  assert(name != NULL);
4844  assert(cons != NULL);
4845 
4846  /* check that string starts with "sqrt( " */
4847  if( strncmp(str, "sqrt( ", 6) != 0 )
4848  {
4849  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected 'sqrt( ' at begin of soc constraint string '%s'\n", str);
4850  *success = FALSE;
4851  return SCIP_OKAY;
4852  }
4853  str += 6;
4854 
4855  *success = TRUE;
4856 
4857  /* check if we have a constant in the beginning */
4858  if( SCIPstrToRealValue(str, &constant, &endptr) )
4859  str = endptr;
4860  else
4861  constant = 0.0;
4862 
4863  nvars = 0;
4864  varssize = 5;
4865  SCIP_CALL( SCIPallocBufferArray(scip, &vars, varssize) );
4866  SCIP_CALL( SCIPallocBufferArray(scip, &coefs, varssize) );
4867  SCIP_CALL( SCIPallocBufferArray(scip, &offsets, varssize) );
4868 
4869  /* read '+ (coef*(var+offset))^2' on lhs, as long as possible */
4870  while( *str != '\0' )
4871  {
4872  /* skip whitespace */
4873  while( isspace((int)*str) )
4874  ++str;
4875 
4876  /* stop if no more coefficients */
4877  if( strncmp(str, "+ (", 3) != 0 )
4878  break;
4879 
4880  str += 3;
4881 
4882  /* parse coef */
4883  if( !SCIPstrToRealValue(str, &coef, &endptr) )
4884  {
4885  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected coefficient at begin of '%s'\n", str);
4886  *success = FALSE;
4887  break;
4888  }
4889  str = endptr;
4890 
4891  if( strncmp(str, "*(", 2) != 0 )
4892  {
4893  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected '*(' at begin of '%s'\n", str);
4894  *success = FALSE;
4895  break;
4896  }
4897  str += 2;
4898 
4899  /* parse variable name */
4900  SCIP_CALL( SCIPparseVarName(scip, str, &var, &endptr) );
4901  if( var == NULL )
4902  {
4903  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "unknown variable name at '%s'\n", str);
4904  *success = FALSE;
4905  break;
4906  }
4907  str = endptr;
4908 
4909  /* parse offset */
4910  if( !SCIPstrToRealValue(str, &offset, &endptr) )
4911  {
4912  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected offset at begin of '%s'\n", str);
4913  *success = FALSE;
4914  break;
4915  }
4916  str = endptr;
4917 
4918  if( strncmp(str, "))^2", 4) != 0 )
4919  {
4920  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected '))^2' at begin of '%s'\n", str);
4921  *success = FALSE;
4922  break;
4923  }
4924  str += 4;
4925 
4926  if( varssize <= nvars )
4927  {
4928  varssize = SCIPcalcMemGrowSize(scip, varssize+1);
4929  SCIP_CALL( SCIPreallocBufferArray(scip, &vars, varssize) );
4930  SCIP_CALL( SCIPreallocBufferArray(scip, &coefs, varssize) );
4931  SCIP_CALL( SCIPreallocBufferArray(scip, &offsets, varssize) );
4932  }
4933  vars[nvars] = var;
4934  coefs[nvars] = coef;
4935  offsets[nvars] = offset;
4936  ++nvars;
4937  }
4938 
4939  if( strncmp(str, ") <=", 4) != 0 )
4940  {
4941  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected ') <=' at begin of '%s'\n", str);
4942  *success = FALSE;
4943  }
4944  str += 4;
4945 
4946  /* read rhs coef*(var+offset) or just a constant */
4947 
4948  /* parse coef */
4949  if( *success )
4950  {
4951  /* skip whitespace */
4952  while( isspace((int)*str) )
4953  ++str;
4954 
4955  if( !SCIPstrToRealValue(str, &rhscoef, &endptr) )
4956  {
4957  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected coefficient at begin of '%s'\n", str);
4958  *success = FALSE;
4959  }
4960  str = endptr;
4961 
4962  /* skip whitespace */
4963  while( isspace((int)*str) )
4964  ++str;
4965  }
4966 
4967  /* parse *(var+offset) */
4968  if( *str != '\0' )
4969  {
4970  if( *success )
4971  {
4972  if( strncmp(str, "*(", 2) != 0 )
4973  {
4974  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected '*(' at begin of '%s'\n", str);
4975  *success = FALSE;
4976  }
4977  else
4978  {
4979  str += 2;
4980  }
4981  }
4982 
4983  /* parse variable name */
4984  if( *success )
4985  {
4986  SCIP_CALL( SCIPparseVarName(scip, str, &rhsvar, &endptr) );
4987  if( rhsvar == NULL )
4988  {
4989  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "unknown variable name at '%s'\n", str);
4990  *success = FALSE;
4991  }
4992  else
4993  {
4994  str = endptr;
4995  }
4996  }
4997 
4998  /* parse offset */
4999  if( *success )
5000  {
5001  if( !SCIPstrToRealValue(str, &rhsoffset, &endptr) )
5002  {
5003  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected offset at begin of '%s'\n", str);
5004  *success = FALSE;
5005  }
5006  else
5007  {
5008  str = endptr;
5009  }
5010  }
5011 
5012  if( *success )
5013  {
5014  if( *str != ')' )
5015  {
5016  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL, "expected ')' at begin of '%s'\n", str);
5017  *success = FALSE;
5018  }
5019  }
5020  }
5021  else if( *success )
5022  {
5023  /* only a constant at right hand side */
5024  rhsoffset = rhscoef; /*lint !e644*/
5025  rhscoef = 1.0;
5026  rhsvar = NULL;
5027  }
5028 
5029  if( *success )
5030  {
5031  assert(!stickingatnode);
5032  SCIP_CALL( SCIPcreateConsSOC(scip, cons, name, nvars, vars, coefs, offsets, constant, rhsvar, rhscoef, rhsoffset,
5033  initial, separate, enforce, check, propagate, local, modifiable, dynamic, removable) ); /*lint !e644 */
5034  }
5035 
5036  SCIPfreeBufferArray(scip, &offsets);
5037  SCIPfreeBufferArray(scip, &coefs);
5038  SCIPfreeBufferArray(scip, &vars);
5039 
5040  return SCIP_OKAY;
5041 }
5042 
5043 /** constraint method of constraint handler which returns the variables (if possible) */
5044 static
5045 SCIP_DECL_CONSGETVARS(consGetVarsSOC)
5046 { /*lint --e{715}*/
5047  SCIP_CONSDATA* consdata;
5048 
5049  consdata = SCIPconsGetData(cons);
5050  assert(consdata != NULL);
5051 
5052  if( varssize < consdata->nvars + 1 )
5053  (*success) = FALSE;
5054  else
5055  {
5056  BMScopyMemoryArray(vars, consdata->vars, consdata->nvars);
5057  vars[consdata->nvars] = consdata->rhsvar;
5058  (*success) = TRUE;
5059  }
5060 
5061  return SCIP_OKAY;
5062 }
5063 
5064 /** constraint method of constraint handler which returns the number of variable (if possible) */
5065 static
5066 SCIP_DECL_CONSGETNVARS(consGetNVarsSOC)
5067 { /*lint --e{715}*/
5068  SCIP_CONSDATA* consdata;
5069 
5070  assert(cons != NULL);
5071 
5072  consdata = SCIPconsGetData(cons);
5073  assert(consdata != NULL);
5074 
5075  (*nvars) = consdata->nvars + 1;
5076  (*success) = TRUE;
5077 
5078  return SCIP_OKAY;
5079 }
5080 
5081 /*
5082  * constraint specific interface methods
5083  */
5084 
5085 /** creates the handler for second order cone constraints and includes it in SCIP */
5087  SCIP* scip /**< SCIP data structure */
5088  )
5089 {
5090  SCIP_CONSHDLRDATA* conshdlrdata;
5091  SCIP_CONSHDLR* conshdlr;
5092  SCIP_EVENTHDLR* eventhdlr;
5093 
5094  /* create constraint handler data */
5095  SCIP_CALL( SCIPallocBlockMemory(scip, &conshdlrdata) );
5096  conshdlrdata->subnlpheur = NULL;
5097  conshdlrdata->trysolheur = NULL;
5098 
5099  SCIP_CALL( SCIPincludeEventhdlrBasic(scip, &eventhdlr, CONSHDLR_NAME"_boundchange",
5100  "signals a bound change to a second order cone constraint",
5101  processVarEvent, NULL) );
5102  conshdlrdata->eventhdlr = eventhdlr;
5103 
5104  SCIP_CALL( SCIPincludeEventhdlrBasic(scip, NULL, CONSHDLR_NAME"_newsolution",
5105  "handles the event that a new primal solution has been found",
5106  processNewSolutionEvent, NULL) );
5107 
5108  /* include constraint handler */
5111  consEnfolpSOC, consEnfopsSOC, consCheckSOC, consLockSOC,
5112  conshdlrdata) );
5113  assert(conshdlr != NULL);
5114 
5115  /* set non-fundamental callbacks via specific setter functions */
5116  SCIP_CALL( SCIPsetConshdlrCopy(scip, conshdlr, conshdlrCopySOC, consCopySOC) );
5117  SCIP_CALL( SCIPsetConshdlrDelete(scip, conshdlr, consDeleteSOC) );
5118  SCIP_CALL( SCIPsetConshdlrExit(scip, conshdlr, consExitSOC) );
5119  SCIP_CALL( SCIPsetConshdlrExitpre(scip, conshdlr, consExitpreSOC) );
5120  SCIP_CALL( SCIPsetConshdlrExitsol(scip, conshdlr, consExitsolSOC) );
5121  SCIP_CALL( SCIPsetConshdlrFree(scip, conshdlr, consFreeSOC) );
5122  SCIP_CALL( SCIPsetConshdlrGetVars(scip, conshdlr, consGetVarsSOC) );
5123  SCIP_CALL( SCIPsetConshdlrGetNVars(scip, conshdlr, consGetNVarsSOC) );
5124  SCIP_CALL( SCIPsetConshdlrInit(scip, conshdlr, consInitSOC) );
5125  SCIP_CALL( SCIPsetConshdlrInitsol(scip, conshdlr, consInitsolSOC) );
5126  SCIP_CALL( SCIPsetConshdlrParse(scip, conshdlr, consParseSOC) );
5127  SCIP_CALL( SCIPsetConshdlrPresol(scip, conshdlr, consPresolSOC, CONSHDLR_MAXPREROUNDS, CONSHDLR_PRESOLTIMING) );
5128  SCIP_CALL( SCIPsetConshdlrPrint(scip, conshdlr, consPrintSOC) );
5130  SCIP_CALL( SCIPsetConshdlrSepa(scip, conshdlr, consSepalpSOC, consSepasolSOC, CONSHDLR_SEPAFREQ,
5132  SCIP_CALL( SCIPsetConshdlrTrans(scip, conshdlr, consTransSOC) );
5133  SCIP_CALL( SCIPsetConshdlrEnforelax(scip, conshdlr, consEnforelaxSOC) );
5134 
5135  if( SCIPfindConshdlr(scip,"quadratic") != NULL )
5136  {
5137  /* notify function that upgrades quadratic constraint to SOC's */
5139  }
5140 
5141  /* add soc constraint handler parameters */
5142  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/projectpoint",
5143  "whether the reference point of a cut should be projected onto the feasible set of the SOC constraint",
5144  &conshdlrdata->projectpoint, TRUE, FALSE, NULL, NULL) );
5145 
5146  SCIP_CALL( SCIPaddIntParam (scip, "constraints/" CONSHDLR_NAME "/nauxvars",
5147  "number of auxiliary variables to use when creating a linear outer approx. of a SOC3 constraint; 0 to turn off",
5148  &conshdlrdata->nauxvars, FALSE, 0, 0, INT_MAX, NULL, NULL) );
5149 
5150  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/glineur",
5151  "whether the Glineur Outer Approximation should be used instead of Ben-Tal Nemirovski",
5152  &conshdlrdata->glineur, FALSE, TRUE, NULL, NULL) );
5153 
5154  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/sparsify",
5155  "whether to sparsify cuts",
5156  &conshdlrdata->sparsify, TRUE, FALSE, NULL, NULL) );
5157 
5158  SCIP_CALL( SCIPaddRealParam(scip, "constraints/" CONSHDLR_NAME "/sparsifymaxloss",
5159  "maximal loss in cut efficacy by sparsification",
5160  &conshdlrdata->sparsifymaxloss, TRUE, 0.2, 0.0, 1.0, NULL, NULL) );
5161 
5162  SCIP_CALL( SCIPaddRealParam(scip, "constraints/" CONSHDLR_NAME "/sparsifynzgrowth",
5163  "growth rate of maximal allowed nonzeros in cuts in sparsification",
5164  &conshdlrdata->sparsifynzgrowth, TRUE, 1.3, 1.000001, SCIPinfinity(scip), NULL, NULL) );
5165 
5166  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/linfeasshift",
5167  "whether to try to make solutions feasible in check by shifting the variable on the right hand side",
5168  &conshdlrdata->linfeasshift, FALSE, TRUE, NULL, NULL) );
5169 
5170  SCIP_CALL( SCIPaddCharParam(scip, "constraints/" CONSHDLR_NAME "/nlpform",
5171  "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)",
5172  &conshdlrdata->nlpform, FALSE, 'a', "aqsed", NULL, NULL) );
5173 
5174  SCIP_CALL( SCIPaddRealParam(scip, "constraints/" CONSHDLR_NAME "/sepanlpmincont",
5175  "minimal required fraction of continuous variables in problem to use solution of NLP relaxation in root for separation",
5176  &conshdlrdata->sepanlpmincont, FALSE, 1.0, 0.0, 2.0, NULL, NULL) );
5177 
5178  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/enfocutsremovable",
5179  "are cuts added during enforcement removable from the LP in the same node?",
5180  &conshdlrdata->enfocutsremovable, TRUE, FALSE, NULL, NULL) );
5181 
5182  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/generalsocupgrade",
5183  "try to upgrade more general quadratics to soc?",
5184  &conshdlrdata->generalsocupg, TRUE, TRUE, NULL, NULL) );
5185 
5186  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/" CONSHDLR_NAME "/disaggregate",
5187  "try to completely disaggregate soc?",
5188  &conshdlrdata->disaggregate, TRUE, TRUE, NULL, NULL) );
5189 
5190  return SCIP_OKAY;
5191 }
5192 
5193 /** creates and captures a second order cone constraint
5194  *
5195  * @note the constraint gets captured, hence at one point you have to release it using the method SCIPreleaseCons()
5196  */
5198  SCIP* scip, /**< SCIP data structure */
5199  SCIP_CONS** cons, /**< pointer to hold the created constraint */
5200  const char* name, /**< name of constraint */
5201  int nvars, /**< number of variables on left hand side of constraint (n) */
5202  SCIP_VAR** vars, /**< array with variables on left hand side (x_i) */
5203  SCIP_Real* coefs, /**< array with coefficients of left hand side variables (alpha_i), or NULL if all 1.0 */
5204  SCIP_Real* offsets, /**< array with offsets of variables (beta_i), or NULL if all 0.0 */
5205  SCIP_Real constant, /**< constant on left hand side (gamma) */
5206  SCIP_VAR* rhsvar, /**< variable on right hand side of constraint (x_{n+1}) */
5207  SCIP_Real rhscoeff, /**< coefficient of variable on right hand side (alpha_{n+1}) */
5208  SCIP_Real rhsoffset, /**< offset of variable on right hand side (beta_{n+1}) */
5209  SCIP_Bool initial, /**< should the LP relaxation of constraint be in the initial LP?
5210  * Usually set to TRUE. Set to FALSE for 'lazy constraints'. */
5211  SCIP_Bool separate, /**< should the constraint be separated during LP processing?
5212  * Usually set to TRUE. */
5213  SCIP_Bool enforce, /**< should the constraint be enforced during node processing?
5214  * TRUE for model constraints, FALSE for additional, redundant constraints. */
5215  SCIP_Bool check, /**< should the constraint be checked for feasibility?
5216  * TRUE for model constraints, FALSE for additional, redundant constraints. */
5217  SCIP_Bool propagate, /**< should the constraint be propagated during node processing?
5218  * Usually set to TRUE. */
5219  SCIP_Bool local, /**< is constraint only valid locally?
5220  * Usually set to FALSE. Has to be set to TRUE, e.g., for branching constraints. */
5221  SCIP_Bool modifiable, /**< is constraint modifiable (subject to column generation)?
5222  * Usually set to FALSE. In column generation applications, set to TRUE if pricing
5223  * adds coefficients to this constraint. */
5224  SCIP_Bool dynamic, /**< is constraint subject to aging?
5225  * Usually set to FALSE. Set to TRUE for own cuts which
5226  * are separated as constraints. */
5227  SCIP_Bool removable /**< should the relaxation be removed from the LP due to aging or cleanup?
5228  * Usually set to FALSE. Set to TRUE for 'lazy constraints' and 'user cuts'. */
5229  )
5230 {
5231  SCIP_CONSHDLR* conshdlr;
5232  SCIP_CONSDATA* consdata;
5233  int i;
5234 
5235  assert(scip != NULL);
5236  assert(cons != NULL);
5237  assert(modifiable == FALSE); /* we do not support column generation */
5238 
5239  /* find the soc constraint handler */
5240  conshdlr = SCIPfindConshdlr(scip, CONSHDLR_NAME);
5241  if( conshdlr == NULL )
5242  {
5243  SCIPerrorMessage("SOC constraint handler not found\n");
5244  return SCIP_PLUGINNOTFOUND;
5245  }
5246 
5247  assert(vars != NULL);
5248  assert(nvars >= 1);
5249  assert(constant >= 0.0);
5250  assert(!SCIPisInfinity(scip, ABS(rhsoffset)));
5251  assert(!SCIPisInfinity(scip, constant));
5252  assert(rhsvar == NULL || rhscoeff <= 0.0 || SCIPisGE(scip, local ? SCIPcomputeVarLbLocal(scip, rhsvar) : SCIPcomputeVarLbGlobal(scip, rhsvar), -rhsoffset));
5253  assert(rhsvar == NULL || rhscoeff >= 0.0 || SCIPisLE(scip, local ? SCIPcomputeVarUbLocal(scip, rhsvar) : SCIPcomputeVarUbGlobal(scip, rhsvar), -rhsoffset));
5254 
5255  /* create constraint data */
5256  SCIP_CALL( SCIPallocBlockMemory(scip, &consdata) );
5257 
5258  consdata->nvars = nvars;
5259  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &consdata->vars, vars, nvars) );
5260  for( i = 0; i < nvars; ++i )
5261  {
5262  SCIP_CALL( SCIPcaptureVar(scip, vars[i]) );
5263  }
5264 
5265  if( coefs != NULL )
5266  {
5267  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &consdata->coefs, coefs, nvars) );
5268  for( i = 0; i < nvars; ++i )
5269  {
5270  if( consdata->coefs[i] < 0.0 )
5271  consdata->coefs[i] = -consdata->coefs[i];
5272  }
5273  }
5274  else
5275  {
5276  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->coefs, nvars) );
5277  for( i = 0; i < nvars; ++i )
5278  consdata->coefs[i] = 1.0;
5279  }
5280 
5281  if( offsets != NULL )
5282  {
5283  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &consdata->offsets, offsets, nvars) );
5284  }
5285  else
5286  {
5287  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->offsets, nvars) );
5288  BMSclearMemoryArray(consdata->offsets, nvars);
5289  }
5290 
5291  consdata->constant = constant;
5292  consdata->rhsvar = rhsvar;
5293  consdata->rhscoeff = rhscoeff;
5294  consdata->rhsoffset = rhsoffset;
5295 
5296  if( rhsvar != NULL )
5297  {
5298  SCIP_CALL( SCIPcaptureVar(scip, rhsvar) );
5299  }
5300 
5301  consdata->nlrow = NULL;
5302 
5303  consdata->lhsbndchgeventdata = NULL;
5304  consdata->isapproxadded = FALSE;
5305 
5306  /* create constraint */
5307  SCIP_CALL( SCIPcreateCons(scip, cons, name, conshdlr, consdata, initial, separate, enforce, check, propagate,
5308  local, modifiable, dynamic, removable, FALSE) );
5309 
5310  if( SCIPisTransformed(scip) )
5311  {
5312  SCIP_CONSHDLRDATA* conshdlrdata = SCIPconshdlrGetData(conshdlr);
5313  assert(conshdlrdata != NULL);
5314 
5315  SCIP_CALL( catchVarEvents(scip, conshdlrdata->eventhdlr, *cons) );
5316  }
5317 
5318  return SCIP_OKAY;
5319 }
5320 
5321 /** creates and captures a second order cone constraint with all its constraint flags
5322  * set to their default values
5323  *
5324  * @note the constraint gets captured, hence at one point you have to release it using the method SCIPreleaseCons()
5325  */
5327  SCIP* scip, /**< SCIP data structure */
5328  SCIP_CONS** cons, /**< pointer to hold the created constraint */
5329  const char* name, /**< name of constraint */
5330  int nvars, /**< number of variables on left hand side of constraint (n) */
5331  SCIP_VAR** vars, /**< array with variables on left hand side (x_i) */
5332  SCIP_Real* coefs, /**< array with coefficients of left hand side variables (alpha_i), or NULL if all 1.0 */
5333  SCIP_Real* offsets, /**< array with offsets of variables (beta_i), or NULL if all 0.0 */
5334  SCIP_Real constant, /**< constant on left hand side (gamma) */
5335  SCIP_VAR* rhsvar, /**< variable on right hand side of constraint (x_{n+1}) */
5336  SCIP_Real rhscoeff, /**< coefficient of variable on right hand side (alpha_{n+1}) */
5337  SCIP_Real rhsoffset /**< offset of variable on right hand side (beta_{n+1}) */
5338  )
5339 {
5340  SCIP_CALL( SCIPcreateConsSOC(scip, cons, name, nvars, vars, coefs, offsets, constant,
5341  rhsvar, rhscoeff, rhsoffset,
5342  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE) );
5343 
5344  return SCIP_OKAY;
5345 }
5346 
5347 /** Gets the SOC constraint as a nonlinear row representation. */
5349  SCIP* scip, /**< SCIP data structure */
5350  SCIP_CONS* cons, /**< constraint */
5351  SCIP_NLROW** nlrow /**< pointer to store nonlinear row */
5352  )
5353 {
5354  SCIP_CONSDATA* consdata;
5355 
5356  assert(cons != NULL);
5357  assert(nlrow != NULL);
5358 
5359  consdata = SCIPconsGetData(cons);
5360  assert(consdata != NULL);
5361 
5362  if( consdata->nlrow == NULL )
5363  {
5364  SCIP_CALL( createNlRow(scip, SCIPconsGetHdlr(cons), cons) );
5365  }
5366  assert(consdata->nlrow != NULL);
5367  *nlrow = consdata->nlrow;
5368 
5369  return SCIP_OKAY;
5370 }
5371 
5372 /** Gets the number of variables on the left hand side of a SOC constraint. */
5373 int SCIPgetNLhsVarsSOC(
5374  SCIP* scip, /**< SCIP data structure */
5375  SCIP_CONS* cons /**< constraint data */
5376  )
5377 {
5378  assert(cons != NULL);
5379  assert(SCIPconsGetData(cons) != NULL);
5380 
5381  return SCIPconsGetData(cons)->nvars;
5382 }
5383 
5384 /** Gets the variables on the left hand side of a SOC constraint. */