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 sepa_convexproj.c
17 * @ingroup DEFPLUGINS_SEPA
18 * @brief convexproj separator
19 * @author Felipe Serrano
20 *
21 * @todo should separator only be run when SCIPallColsInLP is true?
22 * @todo check if it makes sense to implement the copy callback
23 * @todo add SCIPisStopped(scip) to the condition of time consuming loops
24 */
25 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
26
27 #include <assert.h>
28 #include <string.h>
29
30 #include "blockmemshell/memory.h"
31 #include "nlpi/exprinterpret.h"
32 #include "nlpi/nlpi.h"
33 #include "nlpi/pub_expr.h"
34 #include "scip/pub_message.h"
35 #include "scip/pub_misc.h"
36 #include "scip/pub_nlp.h"
37 #include "scip/pub_sepa.h"
38 #include "scip/pub_var.h"
39 #include "scip/scip_cut.h"
40 #include "scip/scip_general.h"
41 #include "scip/scip_lp.h"
42 #include "scip/scip_mem.h"
43 #include "scip/scip_message.h"
44 #include "scip/scip_nlp.h"
45 #include "scip/scip_nonlinear.h"
46 #include "scip/scip_numerics.h"
47 #include "scip/scip_param.h"
48 #include "scip/scip_prob.h"
49 #include "scip/scip_sepa.h"
50 #include "scip/scip_sol.h"
51 #include "scip/scip_solvingstats.h"
52 #include "scip/scip_timing.h"
53 #include "scip/scip_tree.h"
54 #include "scip/sepa_convexproj.h"
55 #include <string.h>
56
57
58 #define SEPA_NAME "convexproj"
59 #define SEPA_DESC "separate at projection of point onto convex region"
60 #define SEPA_PRIORITY 0
61 #define SEPA_FREQ -1
62 #define SEPA_MAXBOUNDDIST 1.0
63 #define SEPA_USESSUBSCIP FALSE /**< does the separator use a secondary SCIP instance? */
64 #define SEPA_DELAY TRUE /**< should separation method be delayed, if other separators found cuts? */
65
66 #define DEFAULT_MAXDEPTH -1 /* maximum depth at which the separator is applied; -1 means no limit */
67 #define DEFAULT_NLPTIMELIMIT 0.0 /**< default time limit of NLP solver; 0.0 for no limit */
68 #define DEFAULT_NLPITERLIM 250 /**< default NLP iteration limit */
69
70 #define VIOLATIONFAC 100 /* points regarded violated if max violation > VIOLATIONFAC*SCIPfeastol */
71
72 #define NLPVERBOSITY 0 /**< NLP solver verbosity */
73
74 /*
75 * Data structures
76 */
77
78 /** side that makes an nlrow convex */
79 enum ConvexSide
80 {
81 LHS = 0, /**< left hand side */
82 RHS = 1 /**< right hand side */
83 };
84 typedef enum ConvexSide CONVEXSIDE;
85
86 /** separator data
87 * it keeps the nlpi which represents the projection problem (see sepa_convexproj.h); it also keeps the convex nlrows
88 * and the side which actually makes them convex; when separating, we use the nlpi to compute the projection and then
89 * the convex nlrows to compute the actual gradient cuts */
90 struct SCIP_SepaData
91 {
92 SCIP_NLPI* nlpi; /**< nlpi used to create the nlpi problem */
93 SCIP_NLPIPROBLEM* nlpiprob; /**< nlpi problem representing the convex NLP relaxation */
94 SCIP_VAR** nlpivars; /**< array containing all variables of the nlpi */
95 SCIP_HASHMAP* var2nlpiidx; /**< mapping between variables and nlpi indices */
96 int nlpinvars; /**< total number of nlpi variables */
97
98 SCIP_Bool skipsepa; /**< should separator be skipped? */
99
100 SCIP_NLROW** nlrows; /**< convex nlrows */
101 CONVEXSIDE* convexsides; /**< which sides make the nlrows convex */
102 SCIP_Real* constraintviolation;/**< array storing the violation of constraint by current solution; 0.0 if it is not violated */
103 int nnlrows; /**< total number of nlrows */
104 int nlrowssize; /**< memory allocated for nlrows, convexsides and nlrowsidx */
105
106 SCIP_EXPRINT* exprinterpreter; /**< expression interpreter to compute gradients */
107
108 /* parameter */
109 SCIP_Real nlptimelimit; /**< time limit of NLP solver; 0.0 for no limit */
110 int nlpiterlimit; /**< iteration limit of NLP solver; 0 for no limit */
111 int maxdepth; /**< maximal depth at which the separator is applied */
112
113 int ncuts; /**< number of cuts generated */
114 };
115
116
117 /*
118 * Local methods
119 */
120
121 /** clears the sepadata data */
122 static
sepadataClear(SCIP * scip,SCIP_SEPADATA * sepadata)123 SCIP_RETCODE sepadataClear(
124 SCIP* scip, /**< SCIP data structure */
125 SCIP_SEPADATA* sepadata /**< separator data */
126 )
127 {
128 assert(sepadata != NULL);
129
130 /* nlrowssize gets allocated first and then its decided whether to create the nlpiprob */
131 if( sepadata->nlrowssize > 0 )
132 {
133 SCIPfreeBlockMemoryArray(scip, &sepadata->constraintviolation, sepadata->nlrowssize);
134 SCIPfreeBlockMemoryArray(scip, &sepadata->convexsides, sepadata->nlrowssize);
135 SCIPfreeBlockMemoryArray(scip, &sepadata->nlrows, sepadata->nlrowssize);
136 sepadata->nlrowssize = 0;
137 }
138
139 if( sepadata->nlpiprob != NULL )
140 {
141 assert(sepadata->nlpi != NULL);
142
143 SCIPfreeBlockMemoryArray(scip, &sepadata->nlpivars, sepadata->nlpinvars);
144
145 SCIPhashmapFree(&sepadata->var2nlpiidx);
146 SCIP_CALL( SCIPnlpiFreeProblem(sepadata->nlpi, &sepadata->nlpiprob) );
147 SCIP_CALL( SCIPexprintFree(&sepadata->exprinterpreter) );
148
149 sepadata->nlpinvars = 0;
150 sepadata->nnlrows = 0;
151 }
152 assert(sepadata->nlpinvars == 0);
153 assert(sepadata->nnlrows == 0);
154 assert(sepadata->nlrowssize == 0);
155
156 sepadata->skipsepa = FALSE;
157
158 return SCIP_OKAY;
159 }
160
161 /** computes gradient of exprtree at projection */
162 static
computeGradient(SCIP * scip,SCIP_EXPRINT * exprint,SCIP_SOL * projection,SCIP_EXPRTREE * exprtree,SCIP_Real * grad)163 SCIP_RETCODE computeGradient(
164 SCIP* scip, /**< SCIP data structure */
165 SCIP_EXPRINT* exprint, /**< expressions interpreter */
166 SCIP_SOL* projection, /**< point where we compute gradient */
167 SCIP_EXPRTREE* exprtree, /**< exprtree for which we compute the gradient */
168 SCIP_Real* grad /**< buffer to store the gradient */
169 )
170 {
171 SCIP_Real* x;
172 SCIP_Real val;
173 int nvars;
174 int i;
175
176 assert(scip != NULL);
177 assert(exprint != NULL);
178 assert(projection != NULL);
179 assert(exprtree != NULL);
180 assert(grad != NULL);
181
182 nvars = SCIPexprtreeGetNVars(exprtree);
183 assert(nvars > 0);
184
185 SCIP_CALL( SCIPallocBufferArray(scip, &x, nvars) );
186
187 /* compile expression exprtree, if not done before */
188 if( SCIPexprtreeGetInterpreterData(exprtree) == NULL )
189 {
190 SCIP_CALL( SCIPexprintCompile(exprint, exprtree) );
191 }
192
193 for( i = 0; i < nvars; ++i )
194 {
195 x[i] = SCIPgetSolVal(scip, projection, SCIPexprtreeGetVars(exprtree)[i]);
196 }
197
198 SCIP_CALL( SCIPexprintGrad(exprint, exprtree, x, TRUE, &val, grad) );
199
200 /*SCIPdebug( for( i = 0; i < nvars; ++i ) printf("%e [%s]\n", grad[i], SCIPvarGetName(SCIPexprtreeGetVars(exprtree)[i])) );*/
201
202 SCIPfreeBufferArray(scip, &x);
203
204 return SCIP_OKAY;
205 }
206
207 /** computes gradient cut (linearization) of nlrow at projection */
208 static
generateCut(SCIP * scip,SCIP_SEPA * sepa,SCIP_EXPRINT * exprint,SCIP_SOL * projection,SCIP_NLROW * nlrow,CONVEXSIDE convexside,SCIP_Real activity,SCIP_ROW ** row)209 SCIP_RETCODE generateCut(
210 SCIP* scip, /**< SCIP data structure */
211 SCIP_SEPA* sepa, /**< the cut separator itself */
212 SCIP_EXPRINT* exprint, /**< expression interpreter */
213 SCIP_SOL* projection, /**< point where we compute gradient cut */
214 SCIP_NLROW* nlrow, /**< constraint for which we generate gradient cut */
215 CONVEXSIDE convexside, /**< which side makes the nlrow convex */
216 SCIP_Real activity, /**< activity of constraint at projection */
217 SCIP_ROW** row /**< storage for cut */
218 )
219 {
220 char rowname[SCIP_MAXSTRLEN];
221 SCIP_SEPADATA* sepadata;
222 SCIP_Real gradx0; /* <grad f(x_0), x_0> */
223 int i;
224
225 assert(scip != NULL);
226 assert(sepa != NULL);
227 assert(exprint != NULL);
228 assert(nlrow != NULL);
229 assert(row != NULL);
230
231 sepadata = SCIPsepaGetData(sepa);
232
233 assert(sepadata != NULL);
234
235 gradx0 = 0.0;
236
237 /* an nlrow has a linear part, quadratic part and expression tree; ideally one would just build the gradient but we
238 * do not know if the different parts share variables or not, so we can't just build the gradient; for this reason
239 * we create the row right away and compute the gradients of each part independently and add them to the row; the
240 * row takes care to add coeffs corresponding to the same variable when they appear in different parts of the nlrow
241 * NOTE: a gradient cut is globally valid whenever the constraint from which it is deduced is globally valid; since
242 * we build the convex relaxation using only globally valid constraints, the cuts are globally valid
243 */
244 (void) SCIPsnprintf(rowname, SCIP_MAXSTRLEN, "proj_cut_%s_%u", SCIPnlrowGetName(nlrow), ++(sepadata->ncuts));
245 SCIP_CALL( SCIPcreateEmptyRowSepa(scip, row, sepa, rowname, -SCIPinfinity(scip), SCIPinfinity(scip), TRUE, FALSE ,
246 TRUE) );
247
248 SCIP_CALL( SCIPcacheRowExtensions(scip, *row) );
249
250 /* linear part */
251 for( i = 0; i < SCIPnlrowGetNLinearVars(nlrow); i++ )
252 {
253 gradx0 += SCIPgetSolVal(scip, projection, SCIPnlrowGetLinearVars(nlrow)[i]) * SCIPnlrowGetLinearCoefs(nlrow)[i];
254 SCIP_CALL( SCIPaddVarToRow(scip, *row, SCIPnlrowGetLinearVars(nlrow)[i], SCIPnlrowGetLinearCoefs(nlrow)[i]) );
255 }
256
257 /* quadratic part */
258 for( i = 0; i < SCIPnlrowGetNQuadElems(nlrow); i++ )
259 {
260 SCIP_VAR* var1;
261 SCIP_VAR* var2;
262 SCIP_Real grad1;
263 SCIP_Real grad2;
264
265 var1 = SCIPnlrowGetQuadVars(nlrow)[SCIPnlrowGetQuadElems(nlrow)[i].idx1];
266 var2 = SCIPnlrowGetQuadVars(nlrow)[SCIPnlrowGetQuadElems(nlrow)[i].idx2];
267 grad1 = SCIPnlrowGetQuadElems(nlrow)[i].coef * SCIPgetSolVal(scip, projection, var2);
268 grad2 = SCIPnlrowGetQuadElems(nlrow)[i].coef * SCIPgetSolVal(scip, projection, var1);
269
270 SCIP_CALL( SCIPaddVarToRow(scip, *row, var1, grad1) );
271 SCIP_CALL( SCIPaddVarToRow(scip, *row, var2, grad2) );
272
273 gradx0 += grad1 * SCIPgetSolVal(scip, projection, var1) + grad2 * SCIPgetSolVal(scip, projection, var2);
274 }
275
276 /* expression tree part */
277 {
278 SCIP_Real* grad;
279 SCIP_EXPRTREE* tree;
280
281 tree = SCIPnlrowGetExprtree(nlrow);
282
283 if( tree != NULL && SCIPexprtreeGetNVars(tree) > 0 )
284 {
285 SCIP_CALL( SCIPallocBufferArray(scip, &grad, SCIPexprtreeGetNVars(tree)) );
286
287 SCIP_CALL( computeGradient(scip, sepadata->exprinterpreter, projection, tree, grad) );
288
289 for( i = 0; i < SCIPexprtreeGetNVars(tree); i++ )
290 {
291 gradx0 += grad[i] * SCIPgetSolVal(scip, projection, SCIPexprtreeGetVars(tree)[i]);
292 SCIP_CALL( SCIPaddVarToRow(scip, *row, SCIPexprtreeGetVars(tree)[i], grad[i]) );
293 }
294
295 SCIPfreeBufferArray(scip, &grad);
296 }
297 }
298
299 SCIP_CALL( SCIPflushRowExtensions(scip, *row) );
300
301 SCIPdebugPrintf("gradient: ");
302 SCIPdebug( SCIP_CALL( SCIPprintRow(scip, *row, NULL) ) );
303 SCIPdebugPrintf("gradient dot x_0: %g\n", gradx0);
304
305 /* gradient cut is f(x_0) - <grad f(x_0), x_0> + <grad f(x_0), x> <= rhs or >= lhs */
306 if( convexside == RHS )
307 {
308 assert(!SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrow)));
309 SCIP_CALL( SCIPchgRowRhs(scip, *row, SCIPnlrowGetRhs(nlrow) - activity + gradx0) );
310 }
311 else
312 {
313 assert(convexside == LHS);
314 assert(!SCIPisInfinity(scip, -SCIPnlrowGetLhs(nlrow)));
315 SCIP_CALL( SCIPchgRowLhs(scip, *row, SCIPnlrowGetLhs(nlrow) - activity + gradx0) );
316 }
317
318 SCIPdebugPrintf("gradient cut: ");
319 SCIPdebug( SCIP_CALL( SCIPprintRow(scip, *row, NULL) ) );
320
321 return SCIP_OKAY;
322 }
323
324 /** set quadratic part of objective function: \f$ \sum_i x_i^2 \f$; the objective function is \f$ ||x - x_0||^2 \f$,
325 * where \f$ x_0 \f$ is the point to separate; the only part that changes is the term \f$ -2 \langle x_0, x \rangle \f$
326 * which is linear and is set every time we want to separate a point, see separateCuts
327 */
328 static
setQuadraticObj(SCIP * scip,SCIP_SEPADATA * sepadata)329 SCIP_RETCODE setQuadraticObj(
330 SCIP* scip, /**< SCIP data structure */
331 SCIP_SEPADATA* sepadata /**< the cut separator data */
332 )
333 {
334 SCIP_QUADELEM* quadelems;
335 int i;
336
337 assert(scip != NULL);
338 assert(sepadata != NULL);
339 assert(sepadata->nlpi != NULL);
340 assert(sepadata->nlpiprob != NULL);
341 assert(sepadata->var2nlpiidx != NULL);
342 assert(sepadata->nlpinvars > 0);
343
344 SCIP_CALL( SCIPallocBufferArray(scip, &quadelems, sepadata->nlpinvars) );
345 for( i = 0; i < sepadata->nlpinvars; i++ )
346 {
347 SCIP_VAR* var;
348
349 var = sepadata->nlpivars[i];
350 assert(SCIPhashmapExists(sepadata->var2nlpiidx, (void*)var) );
351
352 quadelems[i].idx1 = SCIPhashmapGetImageInt(sepadata->var2nlpiidx, (void*)var);
353 quadelems[i].idx2 = quadelems[i].idx1;
354 quadelems[i].coef = 1.0;
355 }
356
357 /* set quadratic part of objective function */
358 SCIP_CALL( SCIPnlpiSetObjective(sepadata->nlpi, sepadata->nlpiprob,
359 0, NULL, NULL, sepadata->nlpinvars, quadelems, NULL, NULL, 0.0) );
360
361 /* free memory */
362 SCIPfreeBufferArray(scip, &quadelems);
363
364 return SCIP_OKAY;
365 }
366
367 /** projects sol onto convex relaxation (stored in sepadata) and tries to generate gradient cuts at the projection
368 * it generates cuts only for the constraints that were violated by the LP solution and are now active or still
369 * violated (in case we don't solve to optimality).
370 * @todo: store a feasible solution if one is found to use as warmstart
371 */
372 static
separateCuts(SCIP * scip,SCIP_SEPA * sepa,SCIP_SOL * sol,SCIP_RESULT * result)373 SCIP_RETCODE separateCuts(
374 SCIP* scip, /**< SCIP data structure */
375 SCIP_SEPA* sepa, /**< the cut separator itself */
376 SCIP_SOL* sol, /**< solution that should be separated */
377 SCIP_RESULT* result /**< pointer to store the result of the separation call */
378 )
379 {
380 SCIP_SEPADATA* sepadata;
381 SCIP_SOL* projection;
382 SCIP_Real* linvals;
383 SCIP_Real* nlpisol;
384 SCIP_Real timelimit;
385 int nlpinvars;
386 int i;
387 int iterlimit;
388 int* lininds;
389 SCIP_Bool nlpunstable;
390
391 nlpunstable = FALSE;
392
393 assert(sepa != NULL);
394
395 sepadata = SCIPsepaGetData(sepa);
396
397 assert(result != NULL);
398 assert(sepadata != NULL);
399 assert(sepadata->nnlrows > 0);
400 assert(sepadata->nlpi != NULL);
401 assert(sepadata->nlpinvars > 0);
402 assert(sepadata->nlrows != NULL);
403 assert(sepadata->nlpiprob != NULL);
404 assert(sepadata->var2nlpiidx != NULL);
405 assert(sepadata->convexsides != NULL);
406 assert(sepadata->constraintviolation != NULL);
407
408 nlpinvars = sepadata->nlpinvars;
409 /* set linear part of objective function: \norm(x - x^0)^2 = \norm(x)^2 - \sum 2 * x_i * x^0_i + const
410 * we ignore the constant; x0 is `sol`
411 */
412 SCIP_CALL( SCIPallocBufferArray(scip, &linvals, nlpinvars) );
413 SCIP_CALL( SCIPallocBufferArray(scip, &lininds, nlpinvars) );
414 for( i = 0; i < nlpinvars; i++ )
415 {
416 SCIP_VAR* var;
417
418 var = sepadata->nlpivars[i];
419 assert(SCIPhashmapExists(sepadata->var2nlpiidx, (void*)var) );
420
421 lininds[i] = SCIPhashmapGetImageInt(sepadata->var2nlpiidx, (void*)var);
422 linvals[i] = - 2.0 * SCIPgetSolVal(scip, sol, var);
423
424 /* if coefficient is too large, don't separate */
425 if( SCIPisHugeValue(scip, REALABS(linvals[i])) )
426 {
427 SCIPdebugMsg(scip, "Don't separate points too close to infinity\n");
428 goto CLEANUP;
429 }
430 }
431
432 /* set linear part of objective function */
433 SCIP_CALL( SCIPnlpiChgLinearCoefs(sepadata->nlpi, sepadata->nlpiprob, -1, nlpinvars, lininds, linvals) );
434
435 /* set parameters in nlpi; time and iterations limit, tolerance, verbosity; for time limit, get time limit of scip;
436 * if scip doesn't have much time left, don't run separator. otherwise, timelimit is the minimum between whats left
437 * for scip and the timelimit setting
438 */
439 SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) );
440 if( !SCIPisInfinity(scip, timelimit) )
441 {
442 timelimit -= SCIPgetSolvingTime(scip);
443 if( timelimit <= 1.0 )
444 {
445 SCIPdebugMsg(scip, "skip NLP solve; no time left\n");
446 goto CLEANUP;
447 }
448 }
449 if( sepadata->nlptimelimit > 0.0 )
450 timelimit = MIN(sepadata->nlptimelimit, timelimit);
451 SCIP_CALL( SCIPnlpiSetRealPar(sepadata->nlpi, sepadata->nlpiprob, SCIP_NLPPAR_TILIM, timelimit) );
452
453 iterlimit = sepadata->nlpiterlimit > 0 ? sepadata->nlpiterlimit : INT_MAX;
454 SCIP_CALL( SCIPnlpiSetIntPar(sepadata->nlpi, sepadata->nlpiprob, SCIP_NLPPAR_ITLIM, iterlimit) );
455 SCIP_CALL( SCIPnlpiSetRealPar(sepadata->nlpi, sepadata->nlpiprob, SCIP_NLPPAR_FEASTOL, SCIPfeastol(scip) / 10.0) ); /* use tighter tolerances for the NLP solver */
456 SCIP_CALL( SCIPnlpiSetRealPar(sepadata->nlpi, sepadata->nlpiprob, SCIP_NLPPAR_RELOBJTOL, MAX(SCIPfeastol(scip), SCIPdualfeastol(scip))) ); /*lint !e666*/
457 SCIP_CALL( SCIPnlpiSetIntPar(sepadata->nlpi, sepadata->nlpiprob, SCIP_NLPPAR_VERBLEVEL, NLPVERBOSITY) );
458
459 /* compute the projection onto the convex NLP relaxation */
460 SCIP_CALL( SCIPnlpiSolve(sepadata->nlpi, sepadata->nlpiprob) );
461 SCIPdebugMsg(scip, "NLP solstat = %d\n", SCIPnlpiGetSolstat(sepadata->nlpi, sepadata->nlpiprob));
462
463 /* if solution is feasible, add cuts */
464 switch( SCIPnlpiGetSolstat(sepadata->nlpi, sepadata->nlpiprob) )
465 {
466 case SCIP_NLPSOLSTAT_GLOBOPT:
467 case SCIP_NLPSOLSTAT_LOCOPT:
468 /* @todo: if solution is optimal, we might as well add the cut <x - P(x_0), x_0 - P(x_0)> <= 0
469 * even though this cut is implied by all the gradient cuts of the rows active at the projection,
470 * we do not add them all (only the gradient cuts of constraints that violated the LP solution */
471 case SCIP_NLPSOLSTAT_FEASIBLE:
472 /* generate cuts for violated constraints (at sol) that are active or still violated at the projection, since
473 * a suboptimal solution or numerical issues could give a solution of the projection problem where constraints
474 * are not active; if the solution of the projection problem is in the interior of the region, we do nothing
475 */
476
477 /* get solution: build SCIP_SOL out of nlpi sol */
478 SCIP_CALL( SCIPnlpiGetSolution(sepadata->nlpi, sepadata->nlpiprob, &nlpisol, NULL, NULL, NULL, NULL) );
479 assert(nlpisol != NULL);
480
481 SCIP_CALL( SCIPcreateSol(scip, &projection, NULL) );
482 for( i = 0; i < nlpinvars; i++ )
483 {
484 SCIP_VAR* var;
485
486 var = sepadata->nlpivars[i];
487 assert(SCIPhashmapExists(sepadata->var2nlpiidx, (void*)var) );
488
489 SCIP_CALL( SCIPsetSolVal(scip, projection, var,
490 nlpisol[SCIPhashmapGetImageInt(sepadata->var2nlpiidx, (void *)var)]) );
491 }
492 SCIPdebug( SCIPprintSol(scip, projection, NULL, TRUE) );
493
494 /* check for active or violated constraints */
495 for( i = 0; i < sepadata->nnlrows; ++i )
496 {
497 SCIP_NLROW* nlrow;
498 CONVEXSIDE convexside;
499 SCIP_Real activity;
500
501 /* ignore constraints that are not violated by `sol` */
502 if( SCIPisFeasZero(scip, sepadata->constraintviolation[i]) )
503 continue;
504
505 convexside = sepadata->convexsides[i];
506 nlrow = sepadata->nlrows[i];
507 assert(nlrow != NULL);
508
509 /* check for currently active constraints at projected point */
510 SCIP_CALL( SCIPgetNlRowSolActivity(scip, nlrow, projection, &activity) );
511
512 SCIPdebugMsg(scip, "NlRow activity at nlpi solution: %g <= %g <= %g\n", SCIPnlrowGetLhs(nlrow), activity,
513 SCIPnlrowGetRhs(nlrow) );
514
515 /* if nlrow is active or violates the projection, build gradient cut at projection */
516 if( (convexside == RHS && SCIPisFeasGE(scip, activity, SCIPnlrowGetRhs(nlrow)))
517 || (convexside == LHS && SCIPisFeasLE(scip, activity, SCIPnlrowGetLhs(nlrow))) )
518 {
519 SCIP_ROW* row;
520
521 SCIP_CALL( generateCut(scip, sepa, sepadata->exprinterpreter, projection, nlrow, convexside, activity,
522 &row) );
523
524 SCIPdebugMsg(scip, "active or violated nlrow: (sols vio: %e)\n", sepadata->constraintviolation[i]);
525 SCIPdebug( SCIP_CALL( SCIPprintNlRow(scip, nlrow, NULL) ) );
526 SCIPdebugMsg(scip, "cut with efficacy %g generated\n", SCIPgetCutEfficacy(scip, sol, row));
527 SCIPdebug( SCIPprintRow(scip, row, NULL) );
528
529 /* add cut if it is efficacious for the point we want to separate (sol) */
530 if( SCIPisCutEfficacious(scip, sol, row) )
531 {
532 SCIP_Bool infeasible;
533
534 SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
535
536 if( infeasible )
537 {
538 *result = SCIP_CUTOFF;
539 SCIP_CALL( SCIPreleaseRow(scip, &row) );
540 break;
541 }
542 else
543 {
544 *result = SCIP_SEPARATED;
545 }
546 }
547
548 /* release the row */
549 SCIP_CALL( SCIPreleaseRow(scip, &row) );
550 }
551 }
552
553 #ifdef SCIP_DEBUG
554 {
555 SCIP_Real distance;
556
557 /* compute distance between LP sol and its projection (only makes sense when it is optimal) */
558 distance = 0.0;
559 for( i = 0; i < SCIPgetNNLPVars(scip); ++i )
560 {
561 SCIP_VAR* var;
562
563 var = SCIPgetNLPVars(scip)[i];
564 assert(var != NULL);
565
566 /* assert NLP solution is within the bounds of the variable (only make sense when sol is optimal) */
567 if( !SCIPisInfinity(scip, -SCIPvarGetLbLocal(var)) )
568 assert(SCIPisFeasLE(scip, SCIPvarGetLbLocal(var), SCIPvarGetNLPSol(var)));
569 if( !SCIPisInfinity(scip, SCIPvarGetUbLocal(var)) )
570 assert(SCIPisFeasLE(scip, SCIPvarGetNLPSol(var), SCIPvarGetUbLocal(var)));
571
572 /*SCIPdebugMsg(scip, "NLP sol (LP sol): %s = %f (%g)\n", SCIPvarGetName(var),
573 * SCIPvarGetNLPSol(var), SCIPgetSolVal(scip, sol, var));
574 */
575
576 distance += SQR( SCIPvarGetNLPSol(var) - SCIPgetSolVal(scip, sol, var) );
577 }
578
579 SCIPdebugMsg(scip, "NLP objval: %e, distance: %e\n", SCIPgetNLPObjval(scip), distance);
580 }
581 #endif
582
583 /* free solution */
584 SCIP_CALL( SCIPfreeSol(scip, &projection) );
585 break;
586
587 case SCIP_NLPSOLSTAT_GLOBINFEASIBLE:
588 case SCIP_NLPSOLSTAT_LOCINFEASIBLE:
589 /* fallthrough;
590 * @todo: write what it means to be locinfeasible and why it can't be used to cutoff the node */
591 case SCIP_NLPSOLSTAT_UNKNOWN:
592 /* unknown... assume numerical issues */
593 nlpunstable = TRUE;
594 break;
595
596 case SCIP_NLPSOLSTAT_UNBOUNDED:
597 default:
598 SCIPerrorMessage("Projection NLP is not unbounded by construction, should not get here!\n");
599 SCIPABORT();
600 nlpunstable = TRUE;
601 }
602
603 /* if nlp is detected to be unstable, don't try to separate again */
604 if( nlpunstable )
605 {
606 /* @todo: maybe change objective function to \sum [(x_i - x_i^*)/max(|x_i^*|, 1)]^2
607 * or some other scaling when unstable and try again.
608 * maybe free it here */
609 sepadata->skipsepa = TRUE;
610 }
611
612 /* reset objective */
613 BMSclearMemoryArray(linvals, nlpinvars);
614 SCIP_CALL( SCIPnlpiChgLinearCoefs(sepadata->nlpi, sepadata->nlpiprob, -1, nlpinvars, lininds, linvals) );
615
616 CLEANUP:
617 /* free memory */
618 SCIPfreeBufferArray(scip, &lininds);
619 SCIPfreeBufferArray(scip, &linvals);
620
621 return SCIP_OKAY;
622 }
623
624 /** computes the violation and maximum violation of the convex nlrows stored in sepadata wrt sol */
625 static
computeMaxViolation(SCIP * scip,SCIP_SEPADATA * sepadata,SCIP_SOL * sol,SCIP_Real * maxviolation)626 SCIP_RETCODE computeMaxViolation(
627 SCIP* scip, /**< SCIP data structure */
628 SCIP_SEPADATA* sepadata, /**< separator data */
629 SCIP_SOL* sol, /**< solution that should be separated */
630 SCIP_Real* maxviolation /**< buffer to store maximum violation */
631 )
632 {
633 SCIP_NLROW* nlrow;
634 int i;
635
636 assert(sepadata != NULL);
637 assert(sepadata->nnlrows > 0);
638 assert(sepadata->nlrows != NULL);
639 assert(sepadata->convexsides != NULL);
640 assert(sepadata->constraintviolation != NULL);
641
642 *maxviolation = 0.0;
643 for( i = 0; i < sepadata->nnlrows; i++ )
644 {
645 SCIP_Real activity;
646 SCIP_Real violation;
647
648 nlrow = sepadata->nlrows[i];
649
650 /* get activity of nlrow */
651 SCIP_CALL( SCIPgetNlRowSolActivity(scip, nlrow, sol, &activity) );
652
653 /* violation = max{activity - rhs, 0.0} when convex and max{lhs - activity, 0.0} when concave */
654 if( sepadata->convexsides[i] == RHS )
655 {
656 assert(SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_CONVEX);
657 assert(!SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrow)));
658
659 violation = activity - SCIPnlrowGetRhs(nlrow);
660 sepadata->constraintviolation[i] = MAX(violation, 0.0);
661 }
662 if( sepadata->convexsides[i] == LHS )
663 {
664 assert(SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_CONCAVE);
665 assert(!SCIPisInfinity(scip, -SCIPnlrowGetLhs(nlrow)));
666
667 violation = SCIPnlrowGetLhs(nlrow) - activity;
668 sepadata->constraintviolation[i] = MAX(violation, 0.0);
669 }
670
671 /* compute maximum */
672 if( *maxviolation < sepadata->constraintviolation[i] )
673 *maxviolation = sepadata->constraintviolation[i];
674 }
675
676 SCIPdebugMsg(scip, "Maximum violation %g\n", *maxviolation);
677
678 return SCIP_OKAY;
679 }
680
681
682 /** stores, from the constraints represented by nlrows, the nonlinear convex ones in sepadata */
683 static
storeNonlinearConvexNlrows(SCIP * scip,SCIP_SEPADATA * sepadata,SCIP_NLROW ** nlrows,int nnlrows)684 SCIP_RETCODE storeNonlinearConvexNlrows(
685 SCIP* scip, /**< SCIP data structure */
686 SCIP_SEPADATA* sepadata, /**< separator data */
687 SCIP_NLROW** nlrows, /**< nlrows from which to store convex ones */
688 int nnlrows /**< number of nlrows */
689 )
690 {
691 int i;
692
693 assert(scip != NULL);
694 assert(sepadata != NULL);
695
696 SCIPdebugMsg(scip, "storing convex nlrows\n");
697
698 sepadata->nlrowssize = nnlrows;
699 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(sepadata->nlrows), nnlrows) );
700 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(sepadata->convexsides), nnlrows) );
701 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(sepadata->constraintviolation), nnlrows) );
702
703 /* count the number of nonlinear convex rows and store them */
704 sepadata->nnlrows = 0;
705 for( i = 0; i < nnlrows; ++i )
706 {
707 SCIP_NLROW* nlrow;
708
709 nlrow = nlrows[i];
710 assert(nlrow != NULL);
711
712 /* linear case */
713 if( SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_LINEAR ||
714 (SCIPnlrowGetNQuadElems(nlrow) == 0 && SCIPnlrowGetExprtree(nlrow) == NULL) )
715 continue;
716
717 /* nonlinear case */
718 if( !SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrow)) && SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_CONVEX )
719 {
720 sepadata->convexsides[sepadata->nnlrows] = RHS;
721 sepadata->nlrows[sepadata->nnlrows] = nlrow;
722 ++(sepadata->nnlrows);
723 }
724 else if( !SCIPisInfinity(scip, -SCIPnlrowGetLhs(nlrow)) && SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_CONCAVE )
725 {
726 sepadata->convexsides[sepadata->nnlrows] = LHS;
727 sepadata->nlrows[sepadata->nnlrows] = nlrow;
728 ++(sepadata->nnlrows);
729 }
730 }
731
732 return SCIP_OKAY;
733 }
734
735
736 /*
737 * Callback methods of separator
738 */
739
740
741 /** destructor of separator to free user data (called when SCIP is exiting) */
742 static
SCIP_DECL_SEPAFREE(sepaFreeConvexproj)743 SCIP_DECL_SEPAFREE(sepaFreeConvexproj)
744 { /*lint --e{715}*/
745 SCIP_SEPADATA* sepadata;
746
747 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
748
749 /* free separator data */
750 sepadata = SCIPsepaGetData(sepa);
751 assert(sepadata != NULL);
752
753 SCIP_CALL( sepadataClear(scip, sepadata) );
754
755 SCIPfreeBlockMemory(scip, &sepadata);
756
757 SCIPsepaSetData(sepa, NULL);
758
759 return SCIP_OKAY;
760 }
761
762 /** solving process deinitialization method of separator (called before branch and bound process data is freed) */
763 static
SCIP_DECL_SEPAEXITSOL(sepaExitsolConvexproj)764 SCIP_DECL_SEPAEXITSOL(sepaExitsolConvexproj)
765 { /*lint --e{715}*/
766 SCIP_SEPADATA* sepadata;
767
768 assert(sepa != NULL);
769
770 sepadata = SCIPsepaGetData(sepa);
771
772 assert(sepadata != NULL);
773
774 SCIP_CALL( sepadataClear(scip, sepadata) );
775
776 return SCIP_OKAY;
777 }
778
779
780 /** LP solution separation method of separator */
781 static
SCIP_DECL_SEPAEXECLP(sepaExeclpConvexproj)782 SCIP_DECL_SEPAEXECLP(sepaExeclpConvexproj)
783 { /*lint --e{715}*/
784 SCIP_Real maxviolation;
785 SCIP_SOL* lpsol;
786 SCIP_SEPADATA* sepadata;
787
788 *result = SCIP_DIDNOTRUN;
789
790 sepadata = SCIPsepaGetData(sepa);
791 assert(sepadata != NULL);
792
793 /* do not run if there is no interesting convex relaxation (with at least one nonlinear convex constraint),
794 * or if we have found it to be numerically unstable
795 * @todo: should it be with at least 2 nonlinear convex constraints?
796 */
797 if( sepadata->skipsepa )
798 {
799 SCIPdebugMsg(scip, "not running because convex relaxation is uninteresting or numerically unstable\n");
800 return SCIP_OKAY;
801 }
802
803 /* the separator needs an NLP solver */
804 if( SCIPgetNNlpis(scip) == 0 )
805 return SCIP_OKAY;
806
807 /* only call separator up to a maximum depth */
808 if( sepadata->maxdepth >= 0 && SCIPgetDepth(scip) > sepadata->maxdepth )
809 return SCIP_OKAY;
810
811 /* only call separator, if we are not close to terminating */
812 if( SCIPisStopped(scip) )
813 return SCIP_OKAY;
814
815 /* do not run if SCIP does not have constructed an NLP */
816 if( !SCIPisNLPConstructed(scip) )
817 {
818 SCIPdebugMsg(scip, "NLP not constructed, skipping convex projection separator\n");
819 return SCIP_OKAY;
820 }
821
822 /* recompute convex NLP relaxation if the variable set changed and we are still at the root node
823 * @todo: does it make sense to do this??? */
824 if( sepadata->nlpiprob != NULL && SCIPgetNVars(scip) != sepadata->nlpinvars && SCIPgetDepth(scip) == 0 )
825 {
826 SCIP_CALL( sepadataClear(scip, sepadata) );
827 assert(sepadata->nlpiprob == NULL);
828 }
829
830 /* create or update convex NLP relaxation */
831 if( sepadata->nlpiprob == NULL )
832 {
833 /* store convex nonlinear constraints */
834 SCIP_CALL( storeNonlinearConvexNlrows(scip, sepadata, SCIPgetNLPNlRows(scip), SCIPgetNNLPNlRows(scip)) );
835
836 /* check that convex NLP relaxation is interesting (more than one nonlinear constraint) */
837 if( sepadata->nnlrows < 1 )
838 {
839 SCIPdebugMsg(scip, "convex relaxation uninteresting, don't run\n");
840 sepadata->skipsepa = TRUE;
841 return SCIP_OKAY;
842 }
843
844 /* create the expression interpreter */
845 SCIP_CALL( SCIPexprintCreate(SCIPblkmem(scip), &sepadata->exprinterpreter) );
846
847 sepadata->nlpinvars = SCIPgetNVars(scip);
848 sepadata->nlpi = SCIPgetNlpis(scip)[0];
849 assert(sepadata->nlpi != NULL);
850
851 SCIP_CALL( SCIPnlpiCreateProblem(sepadata->nlpi, &sepadata->nlpiprob, "convexproj-nlp") );
852 SCIP_CALL( SCIPhashmapCreate(&sepadata->var2nlpiidx, SCIPblkmem(scip), sepadata->nlpinvars) );
853 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &sepadata->nlpivars, SCIPgetVars(scip), sepadata->nlpinvars) ); /*lint !e666*/
854
855 SCIP_CALL( SCIPcreateNlpiProb(scip, sepadata->nlpi, SCIPgetNLPNlRows(scip), SCIPgetNNLPNlRows(scip),
856 sepadata->nlpiprob, sepadata->var2nlpiidx, NULL, NULL, SCIPgetCutoffbound(scip), FALSE, TRUE) );
857
858 /* add rows of the LP */
859 if( SCIPgetDepth(scip) == 0 )
860 {
861 SCIP_CALL( SCIPaddNlpiProbRows(scip, sepadata->nlpi, sepadata->nlpiprob, sepadata->var2nlpiidx,
862 SCIPgetLPRows(scip), SCIPgetNLPRows(scip)) );
863 }
864
865 /* set quadratic part of objective function */
866 SCIP_CALL( setQuadraticObj(scip, sepadata) );
867 }
868 else
869 {
870 SCIP_CALL( SCIPupdateNlpiProb(scip, sepadata->nlpi, sepadata->nlpiprob, sepadata->var2nlpiidx,
871 sepadata->nlpivars, sepadata->nlpinvars, SCIPgetCutoffbound(scip)) );
872 }
873
874 /* assert that the lp solution satisfies the cutoff bound; if this fails then we shouldn't have a cutoff bound in the
875 * nlpi, since then the projection could be in the interior of the actual convex relaxation */
876 assert(SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL ||
877 SCIPisLE(scip, SCIPgetLPObjval(scip), SCIPgetCutoffbound(scip)));
878
879 /* get current sol: LP or pseudo solution if LP sol is not available */
880 SCIP_CALL( SCIPcreateCurrentSol(scip, &lpsol, NULL) );
881
882 /* do not run if current solution's violation is small */
883 SCIP_CALL( computeMaxViolation(scip, sepadata, lpsol, &maxviolation) );
884 if( maxviolation < VIOLATIONFAC * SCIPfeastol(scip) )
885 {
886 SCIPdebugMsg(scip, "solution doesn't violate constraints enough, do not separate\n");
887 SCIP_CALL( SCIPfreeSol(scip, &lpsol) );
888 return SCIP_OKAY;
889 }
890
891 /* run the separator */
892 *result = SCIP_DIDNOTFIND;
893
894 /* separateCuts computes the projection and then gradient cuts on each constraint that was originally violated */
895 SCIP_CALL( separateCuts(scip, sepa, lpsol, result) );
896
897 /* free memory */
898 SCIP_CALL( SCIPfreeSol(scip, &lpsol) );
899
900 return SCIP_OKAY;
901 }
902
903
904 /*
905 * separator specific interface methods
906 */
907
908 /** creates the convexproj separator and includes it in SCIP */
SCIPincludeSepaConvexproj(SCIP * scip)909 SCIP_RETCODE SCIPincludeSepaConvexproj(
910 SCIP* scip /**< SCIP data structure */
911 )
912 {
913 SCIP_SEPADATA* sepadata;
914 SCIP_SEPA* sepa;
915
916 /* create convexproj separator data */
917 SCIP_CALL( SCIPallocBlockMemory(scip, &sepadata) );
918
919 /* this sets all data in sepadata to 0 */
920 BMSclearMemory(sepadata);
921
922 /* include separator */
923 SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST,
924 SEPA_USESSUBSCIP, SEPA_DELAY,
925 sepaExeclpConvexproj, NULL,
926 sepadata) );
927 assert(sepa != NULL);
928
929 /* set non fundamental callbacks via setter functions */
930 SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeConvexproj) );
931 SCIP_CALL( SCIPsetSepaExitsol(scip, sepa, sepaExitsolConvexproj) );
932
933 /* add convexproj separator parameters */
934 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxdepth",
935 "maximal depth at which the separator is applied (-1: unlimited)",
936 &sepadata->maxdepth, FALSE, DEFAULT_MAXDEPTH, -1, INT_MAX, NULL, NULL) );
937
938 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nlpiterlimit",
939 "iteration limit of NLP solver; 0 for no limit",
940 &sepadata->nlpiterlimit, TRUE, DEFAULT_NLPITERLIM, 0, INT_MAX, NULL, NULL) );
941
942 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/nlptimelimit",
943 "time limit of NLP solver; 0.0 for no limit",
944 &sepadata->nlptimelimit, TRUE, DEFAULT_NLPTIMELIMIT, 0.0, SCIP_REAL_MAX, NULL, NULL) );
945
946 return SCIP_OKAY;
947 }
948