Scippy

SCIP

Solving Constraint Integer Programs

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