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