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