1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 2 /* */ 3 /* This file is part of the program and library */ 4 /* SCIP --- Solving Constraint Integer Programs */ 5 /* */ 6 /* Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB) */ 7 /* */ 8 /* Licensed under the Apache License, Version 2.0 (the "License"); */ 9 /* you may not use this file except in compliance with the License. */ 10 /* You may obtain a copy of the License at */ 11 /* */ 12 /* http://www.apache.org/licenses/LICENSE-2.0 */ 13 /* */ 14 /* Unless required by applicable law or agreed to in writing, software */ 15 /* distributed under the License is distributed on an "AS IS" BASIS, */ 16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */ 17 /* See the License for the specific language governing permissions and */ 18 /* limitations under the License. */ 19 /* */ 20 /* You should have received a copy of the Apache-2.0 license */ 21 /* along with SCIP; see the file LICENSE. If not visit scipopt.org. */ 22 /* */ 23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 24 25 /**@file sepa_lagromory.c 26 * @ingroup DEFPLUGINS_SEPA 27 * @brief Lagromory separator 28 * @author Suresh Bolusani 29 * @author Mark Ruben Turner 30 * @author Mathieu Besançon 31 * 32 * This separator is based on the following article that discusses Lagromory separation using the relax-and-cut 33 * framework. Multiple enhancements have been implemented on top of the basic algorithm described in the article. 34 * 35 * Fischetti M. and Salvagnin D. (2011).@n 36 * A relax-and-cut framework for Gomory mixed-integer cuts.@n 37 * Mathematical Programming Computation, 3, 79-102. 38 * 39 * Consider the following linear relaxation at a node: 40 * 41 * \f[ 42 * \begin{array}{rrl} 43 * \min & c^T x &\\ 44 * & x & \in P, 45 * \end{array} 46 * \f] 47 * 48 * where \f$P\f$ is the feasible region of the relaxation. Let the following be the cuts generated so far in the current 49 * separation round. 50 * 51 * \f[ 52 * {\alpha^i}^T x \leq \alpha^i_0, i = 1, 2, \hdots, M 53 * \f] 54 * 55 * Then, the following is the Lagrangian dual problem considered in the relax-and-cut framework used in the separator. 56 * 57 * \f[ 58 * z_D := \max\limits_{u \geq 0} \left\{L(u) := \min \left\{c^T x + \sum\limits_{i = 1}^{M} \left(u_i 59 * \left({\alpha^i}^T x - \alpha^i_0\right) \right) \mid x \in P\right\} \right\}, 60 * \f] 61 * 62 * where \f$u\f$ are the Lagrangian multipliers (referred to as \a dualvector in this separator) used for penalizing the 63 * violation of the generated cuts, and \f$z_D\f$ is the optimal objective value (which is approximated via \a ubparam in this separator). 64 * Then, the following are the steps of the relax-and-cut algorithm implemented in this separator. 65 * 66 * \begin{itemize} 67 * \item Generate an initial pool of cuts to build the initial Lagrangian dual problem. 68 * \item Select initial values for Lagrangian multipliers \f$u^0\f$ (e.g., all zeroes vector). 69 * \item In the outer main loop \f$i\f$ of the algorithm: 70 * \begin{enumerate} 71 * \item Solve the Lagrangian dual problem until certain termination criterion is met. This results in an inner 72 * subgradient loop, whose iteration \f$j\f$ is described below. 73 * \begin{enumerate} 74 * \item Fix \f$u^j\f$, and solve the LP corresponding to the Lagrangian dual with fixed multipliers. 75 * Gather its optimal simplex tableau and optimal objective value (i.e., the Lagrangian value) 76 * \f$L(u^j)\f$. 77 * \item Update \f$u^j\f$ to \f$u^{j+1}\f$ as follows. 78 * \f[ 79 * u^{j+1} = \left(u^j + \lambda_j s^k\right)_+, 80 * \f] 81 * where \f$\lambda_j\f$ is the step length: 82 * \f[ 83 * \lambda_j = \frac{\mu_j (UB - L(u^j))}{\|s^j\|^2_2}, 84 * \f] 85 * where \f$mu_j\f$ is a factor (i.e., \a muparam) such that \f$0 < \mu_j \leq 2\f$, UB is \p ubparam, 86 * and \f$s^j\f$ is the subgradient vector defined as: 87 * \f[ 88 * s^j_k = \left({\alpha^k}^T x - \alpha^k_0\right), k = 1, 2, \hdots, M. 89 * \f] 90 * The factor \f$mu_j\f$ is updated as below. 91 * \f[ 92 * mu_j = \begin{cases} 93 * \p mubacktrackfactor * mu_j & \text{if } L(u^j) < bestLB - \delta\\ 94 * \begin{cases} 95 * \p muslab1factor * mu_j & \text{if } bestLB - avgLB < \p deltaslab1ub * delta\\ 96 * \p muslab2factor * mu_j & \text{if } \p deltaslab1ub * \delta \leq bestLB - avgLB < \p deltaslab2ub * delta\\ 97 * \p muslab3factor * mu_j & \text{otherwise} 98 * \end{cases} & \text{otherwise}, 99 * \end{cases} 100 * \f] 101 * where \f$bestLB\f$ and \f$avgLB\f$ are best and average Lagrangian values found so far, and 102 * \f$\delta = UB - bestLB\f$. 103 * \item Stabilize \f$u^{j+1}\f$ by projecting onto a norm ball followed by taking a convex combination 104 * with a core vector of Lagrangian multipliers. 105 * \item Generate GMI cuts based on the optimal simplex tableau. 106 * \item Relax the newly generated cuts by penalizing and adding them to the objective function. 107 * \item Go to the next iteration \f$j+1\f$. 108 * \end{enumerate} 109 * \item Gather all the generated cuts and build an LP by adding all these cuts to the node relaxation. 110 * \item Solve this LP to obtain its optimal primal and dual solutions. 111 * \item If this primal solution is MIP primal feasible, then add this solution to the solution pool, add all 112 * the generated cuts to the cutpool or sepastore as needed, and exit the separator. 113 * \item Otherwise, update the Lagrangian multipliers based on this optimal dual solution, and go to the next 114 * iteration \f$i+1\f$. 115 * \end{enumerate} 116 * \end{itemize} 117 * 118 * @todo store all LP sols in a data structure, and send them to fix-and-propagate at the end. 119 * 120 * @todo test heuristics such as feasibility pump with multiple input solutions. 121 * 122 * @todo find dual degenerate problems and test the separator on these problems. 123 * 124 * @todo identify instance classes where these cuts work better. 125 * 126 * @todo add termination criteria based on failed efforts. 127 * 128 * @todo for warm starting, if there are additional rows/cols, set their basis status to non-basic and then set WS info. 129 * 130 * @todo filter cuts using multiple explored LP solutions. 131 * 132 * @todo for bases on optimal face only, aggregate to get a new basis and separate it. 133 * 134 * @todo generate other separators in addition to GMI cuts (0-1/2) 135 * 136 * @todo: convert iters from int to SCIP_Longint 137 */ 138 139 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 140 141 #include "scip/heur_trysol.h" 142 #include "scip/sepa_lagromory.h" 143 144 #define SEPA_NAME "lagromory" 145 #define SEPA_DESC "separator for Lagromory cuts for MIP relaxations" 146 #define SEPA_PRIORITY -8000 147 #define SEPA_FREQ -1 148 #define SEPA_MAXBOUNDDIST 1.0 149 #define SEPA_USESSUBSCIP FALSE /**< does the separator use a secondary SCIP instance? */ 150 #define SEPA_DELAY FALSE /**< should separation method be delayed, if other separators found cuts? */ 151 152 /* generic parameters */ 153 #define DEFAULT_AWAY 0.01 /**< minimal integrality violation of a basis variable to try separation */ 154 #define DEFAULT_DELAYEDCUTS FALSE /**< should cuts be added to the delayed cut pool? */ 155 #define DEFAULT_SEPARATEROWS TRUE /**< separate rows with integral slack? */ 156 #define DEFAULT_SORTCUTOFFSOL TRUE /**< sort fractional integer columns based on fractionality? */ 157 #define DEFAULT_SIDETYPEBASIS TRUE /**< choose side types of row (lhs/rhs) based on basis information? */ 158 #define DEFAULT_DYNAMICCUTS TRUE /**< should generated cuts be removed from the LP if they are no longer tight? */ 159 #define DEFAULT_MAKEINTEGRAL FALSE /**< try to scale all cuts to integral coefficients? */ 160 #define DEFAULT_FORCECUTS FALSE /**< force cuts to be added to the LP? */ 161 #define DEFAULT_ALLOWLOCAL FALSE /**< should locally valid cuts be generated? */ 162 163 /* parameters related to the separator's termination check */ 164 #define DEFAULT_MAXROUNDSROOT 1 /**< maximal number of separation rounds in the root node (-1: unlimited) */ 165 #define DEFAULT_MAXROUNDS 1 /**< maximal number of separation rounds per node (-1: unlimited) */ 166 #define DEFAULT_DUALDEGENERACYRATETHRESHOLD 0.5 /**< minimum dual degeneracy rate for separator execution */ 167 #define DEFAULT_VARCONSRATIOTHRESHOLD 1.0 /**< minimum variable-constraint ratio on optimal face for separator execution */ 168 #define DEFAULT_MINRESTART 1 /**< minimum restart round for separator execution (0: from beginning of the 169 instance solving, >= n with n >= 1: from restart round n) */ 170 #define DEFAULT_PERLPMAXCUTSROOT 50 /**< maximal number of cuts separated per Lagromory LP in the root node */ 171 #define DEFAULT_PERLPMAXCUTS 10 /**< maximal number of cuts separated per Lagromory LP in the non-root node */ 172 #define DEFAULT_PERROUNDLPITERLIMITFACTOR -1.0 /**< factor w.r.t. root node LP iterations for maximal separating LP iterations 173 per separation round (negative for no limit) */ 174 #define DEFAULT_ROOTLPITERLIMITFACTOR -1.0 /**< factor w.r.t. root node LP iterations for maximal separating LP iterations 175 in the root node (negative for no limit) */ 176 #define DEFAULT_TOTALLPITERLIMITFACTOR -1.0 /**< factor w.r.t. root node LP iterations for maximal separating LP iterations 177 in the tree (negative for no limit) */ 178 #define DEFAULT_PERROUNDMAXLPITERS 50000 /**< maximal number of separating LP iterations per separation round (-1: unlimited) */ 179 #define DEFAULT_PERROUNDCUTSFACTORROOT 1.0 /**< factor w.r.t. number of integer columns for number of cuts separated per 180 separation round in root node */ 181 #define DEFAULT_PERROUNDCUTSFACTOR 0.5 /**< factor w.r.t. number of integer columns for number of cuts separated per 182 separation round at a non-root node */ 183 #define DEFAULT_TOTALCUTSFACTOR 50.0 /**< factor w.r.t. number of integer columns for total number of cuts separated */ 184 #define DEFAULT_MAXMAINITERS 4 /**< maximal number of main loop iterations of the relax-and-cut algorithm */ 185 #define DEFAULT_MAXSUBGRADIENTITERS 6 /**< maximal number of subgradient loop iterations of the relax-and-cut algorithm */ 186 187 /* parameters related to the relax-and-cut algorithm */ 188 #define DEFAULT_MUPARAMCONST TRUE /**< is the mu parameter (factor for step length) constant? */ 189 #define DEFAULT_MUPARAMINIT 0.01 /**< initial value of the mu parameter (factor for step length) */ 190 #define DEFAULT_MUPARAMLB 0.0 /**< lower bound for the mu parameter (factor for step length) */ 191 #define DEFAULT_MUPARAMUB 2.0 /**< upper bound for the mu parameter (factor for step length) */ 192 #define DEFAULT_MUBACKTRACKFACTOR 0.5 /**< factor of mu while backtracking the mu parameter (factor for step length) - 193 see updateMuSteplengthParam() */ 194 #define DEFAULT_MUSLAB1FACTOR 10.0 /**< factor of mu parameter (factor for step length) for larger increment - see 195 updateMuSteplengthParam() */ 196 #define DEFAULT_MUSLAB2FACTOR 2.0 /**< factor of mu parameter (factor for step length) for smaller increment - see 197 updateMuSteplengthParam() */ 198 #define DEFAULT_MUSLAB3FACTOR 0.5 /**< factor of mu parameter (factor for step length) for reduction - see 199 updateMuSteplengthParam() */ 200 #define DEFAULT_DELTASLAB1UB 0.001 /**< factor of delta deciding larger increment of mu parameter (factor for step 201 length) - see updateMuSteplengthParam() */ 202 #define DEFAULT_DELTASLAB2UB 0.01 /**< factor of delta deciding smaller increment of mu parameter (factor for step 203 length) - see updateMuSteplengthParam() */ 204 #define DEFAULT_UBPARAMPOSFACTOR 2.0 /**< factor for positive upper bound used as an estimate for the optimal 205 Lagrangian dual value */ 206 #define DEFAULT_UBPARAMNEGFACTOR 0.5 /**< factor for negative upper bound used as an estimate for the optimal 207 Lagrangian dual value */ 208 #define DEFAULT_MAXLAGRANGIANVALSFORAVG 2 /**< maximal number of iterations for rolling average of Lagrangian value */ 209 #define DEFAULT_MAXCONSECITERSFORMUUPDATE 10 /**< consecutive number of iterations used to determine if mu needs to be backtracked */ 210 #define DEFAULT_PERROOTLPITERFACTOR 0.2 /**< factor w.r.t. root node LP iterations for iteration limit of each separating 211 LP (negative for no limit) */ 212 #define DEFAULT_PERLPITERFACTOR 0.1 /**< factor w.r.t. non-root node LP iterations for iteration limit of each 213 separating LP (negative for no limit) */ 214 #define DEFAULT_CUTGENFREQ 1 /**< frequency of subgradient iterations for generating cuts */ 215 #define DEFAULT_CUTADDFREQ 1 /**< frequency of subgradient iterations for adding cuts to objective function */ 216 #define DEFAULT_CUTSFILTERFACTOR 1.0 /**< fraction of generated cuts per explored basis to accept from separator */ 217 #define DEFAULT_OPTIMALFACEPRIORITY 2 /**< priority of the optimal face for separator execution (0: low priority, 1: 218 medium priority, 2: high priority) */ 219 #define DEFAULT_AGGREGATECUTS TRUE /**< aggregate all generated cuts using the Lagrangian multipliers? */ 220 /* parameters for stabilization of the Lagrangian multipliers */ 221 #define DEFAULT_PROJECTIONTYPE 2 /**< the ball into which the Lagrangian multipliers are projected for 222 stabilization (0: no projection, 1: L1-norm ball projection, 2: L2-norm ball 223 projection, 3: L_inf-norm ball projection) */ 224 #define DEFAULT_STABILITYCENTERTYPE 1 /**< type of stability center for taking weighted average of Lagrangian multipliers for 225 stabilization (0: no weighted stabilization, 1: best Lagrangian multipliers) */ 226 #define DEFAULT_RADIUSINIT 0.5 /**< initial radius of the ball used in stabilization of Lagrangian multipliers */ 227 #define DEFAULT_RADIUSMAX 20.0 /**< maximum radius of the ball used in stabilization of Lagrangian multipliers */ 228 #define DEFAULT_RADIUSMIN 1e-6 /**< minimum radius of the ball used in stabilization of Lagrangian multipliers */ 229 #define DEFAULT_CONST 2.0 /**< a constant for stablity center based stabilization of Lagrangian multipliers */ 230 #define DEFAULT_RADIUSUPDATEWEIGHT 0.98 /**< multiplier to evaluate cut violation score used for updating ball radius */ 231 232 /* macros that are used directly */ 233 #define RANDSEED 42 /**< random seed */ 234 #define MAKECONTINTEGRAL FALSE /**< convert continuous variable to integral variables in SCIPmakeRowIntegral()? */ 235 #define POSTPROCESS TRUE /**< apply postprocessing after MIR calculation? - see SCIPcalcMIR() */ 236 #define BOUNDSWITCH 0.9999 /**< threshold for bound switching - see SCIPcalcMIR() */ 237 #define USEVBDS TRUE /**< use variable bounds? - see SCIPcalcMIR() */ 238 #define FIXINTEGRALRHS FALSE /**< try to generate an integral rhs? - see SCIPcalcMIR() */ 239 #define MAXAGGRLEN(ncols) (0.1*(ncols)+1000) /**< maximal length of base inequality */ 240 241 /* 242 * Data structures 243 */ 244 245 /** separator data */ 246 struct SCIP_SepaData 247 { 248 /* generic variables */ 249 SCIP_Real away; /**< minimal integrality violation of a basis variable to try separation */ 250 SCIP_Bool delayedcuts; /**< should cuts be added to the delayed cut pool? */ 251 SCIP_Bool separaterows; /**< separate rows with integral slack? */ 252 SCIP_Bool sortcutoffsol; /**< sort fractional integer columns based on fractionality? */ 253 SCIP_Bool sidetypebasis; /**< choose side types of row (lhs/rhs) based on basis information? */ 254 SCIP_Bool dynamiccuts; /**< should generated cuts be removed from the LP if they are no longer tight? */ 255 SCIP_Bool makeintegral; /**< try to scale all cuts to integral coefficients? */ 256 SCIP_Bool forcecuts; /**< force cuts to be added to the LP? */ 257 SCIP_Bool allowlocal; /**< should locally valid cuts be generated? */ 258 SCIP_RANDNUMGEN* randnumgen; /**< random number generator */ 259 SCIP_HEUR* heurtrysol; /**< a pointer to the trysol heuristic, if available */ 260 261 /* used to define separating LPs */ 262 SCIP_LPI* lpiwithsoftcuts; /**< pointer to the lpi interface of Lagrangian dual with fixed multipliers */ 263 SCIP_LPI* lpiwithhardcuts; /**< pointer to the lpi interface of LP with all generated cuts */ 264 int nrowsinhardcutslp; /**< nrows of \a lpiwithhardcuts */ 265 int nrunsinsoftcutslp; /**< number of branch-and-bound runs on current instance */ 266 267 /* used for termination checks */ 268 SCIP_Longint ncalls; /**< number of calls to the separator */ 269 int maxroundsroot; /**< maximal number of separation rounds in the root node (-1: unlimited) */ 270 int maxrounds; /**< maximal number of separation rounds per node (-1: unlimited) */ 271 SCIP_Real dualdegeneracyratethreshold; /**< minimum dual degeneracy rate for separator execution */ 272 SCIP_Real varconsratiothreshold; /**< minimum variable-constraint ratio on optimal face for separator execution */ 273 int minrestart; /**< minimum restart round for separator execution (0: from beginning of the instance 274 solving, >= n with n >= 1: from restart round n) */ 275 int nmaxcutsperlproot; /**< maximal number of cuts separated per Lagromory LP in the root node */ 276 int nmaxcutsperlp; /**< maximal number of cuts separated per Lagromory LP in the non-root node */ 277 SCIP_Real perroundlpiterlimitfactor; /**< factor w.r.t. root node LP iterations for maximal separating LP iterations per 278 separation round (negative for no limit) */ 279 SCIP_Real rootlpiterlimitfactor; /**< factor w.r.t. root node LP iterations for maximal separating LP iterations in the 280 root node (negative for no limit) */ 281 SCIP_Real totallpiterlimitfactor; /**< factor w.r.t. root node LP iterations for maximal separating LP iterations in the 282 tree (negative for no limit) */ 283 int perroundnmaxlpiters; /**< maximal number of separating LP iterations per separation round (-1: unlimited) */ 284 SCIP_Real perroundcutsfactorroot; /**< factor w.r.t. number of integer columns for number of cuts separated per 285 separation round in root node */ 286 SCIP_Real perroundcutsfactor; /**< factor w.r.t. number of integer columns for number of cuts separated per 287 separation round at a non-root node */ 288 SCIP_Real totalcutsfactor; /**< factor w.r.t. number of integer columns for total number of cuts separated */ 289 int nmaxmainiters; /**< maximal number of main loop iterations of the relax-and-cut algorithm */ 290 int nmaxsubgradientiters; /**< maximal number of subgradient loop iterations of the relax-and-cut algorithm */ 291 int nmaxperroundlpiters; /**< maximal number of separating LP iterations per separation round */ 292 int nmaxrootlpiters; /**< maximal number of separating LP iterations in the root node */ 293 int nrootlpiters; /**< number of separating LP iterations in the root node */ 294 int nmaxtotallpiters; /**< maximal number of separating LP iterations in the tree */ 295 int ntotallpiters; /**< number of separating LP iterations in the tree */ 296 int nmaxperroundcutsroot; /**< maximal number of cuts separated per separation round in root node */ 297 int nmaxperroundcuts; /**< maximal number of cuts separated per separation round */ 298 int nmaxtotalcuts; /**< maximal number of cuts separated in the tree */ 299 int ntotalcuts; /**< number of cuts separated in the tree */ 300 301 /* used for the relax-and-cut algorithm */ 302 SCIP_Bool muparamconst; /**< is the mu parameter (factor for step length) constant? */ 303 SCIP_Real muparaminit; /**< initial value of the mu parameter (factor for step length) */ 304 SCIP_Real muparamlb; /**< lower bound for the mu parameter (factor for step length) */ 305 SCIP_Real muparamub; /**< upper bound for the mu parameter (factor for step length) */ 306 SCIP_Real mubacktrackfactor; /**< factor of mu while backtracking the mu parameter (factor for step length) - see 307 updateMuSteplengthParam() */ 308 SCIP_Real muslab1factor; /**< factor of mu parameter (factor for step length) for larger increment - see 309 updateMuSteplengthParam() */ 310 SCIP_Real muslab2factor; /**< factor of mu parameter (factor for step length) for smaller increment - see 311 updateMuSteplengthParam() */ 312 SCIP_Real muslab3factor; /**< factor of mu parameter (factor for step length) for reduction - see updateMuSteplengthParam() */ 313 SCIP_Real deltaslab1ub; /**< factor of delta deciding larger increment of mu parameter (factor for step 314 length) - see updateMuSteplengthParam() */ 315 SCIP_Real deltaslab2ub; /**< factor of delta deciding smaller increment of mu parameter (factor for step 316 length) - see updateMuSteplengthParam() */ 317 SCIP_Real ubparamposfactor; /**< factor for positive upper bound used as an estimate for the optimal Lagrangian 318 dual value */ 319 SCIP_Real ubparamnegfactor; /**< factor for negative upper bound used as an estimate for the optimal Lagrangian 320 dual value */ 321 int nmaxlagrangianvalsforavg; /**< maximal number of iterations for rolling average of Lagrangian value */ 322 int nmaxconsecitersformuupdate; /**< consecutive number of iterations used to determine if mu needs to be backtracked */ 323 SCIP_Real perrootlpiterfactor; /**< factor w.r.t. root node LP iterations for iteration limit of each separating LP 324 (negative for no limit) */ 325 SCIP_Real perlpiterfactor; /**< factor w.r.t. non-root node LP iterations for iteration limit of each separating 326 LP (negative for no limit) */ 327 int cutgenfreq; /**< frequency of subgradient iterations for generating cuts */ 328 int cutaddfreq; /**< frequency of subgradient iterations for adding cuts to objective function */ 329 SCIP_Real cutsfilterfactor; /**< fraction of generated cuts per explored basis to accept from separator */ 330 int optimalfacepriority; /**< priority of the optimal face for separator execution (0: low priority, 1: medium 331 priority, 2: high priority) */ 332 SCIP_Bool aggregatecuts; /**< aggregate all generated cuts using the Lagrangian multipliers? */ 333 334 /* for stabilization of Lagrangian multipliers */ 335 int projectiontype; /**< the ball into which the Lagrangian multipliers are projected for stabilization 336 (0: no projection, 1: L1-norm ball projection, 2: L2-norm ball projection, 3: 337 L_inf-norm ball projection) */ 338 int stabilitycentertype; /**< type of stability center for taking weighted average of Lagrangian multipliers for 339 stabilization (0: no weighted stabilization, 1: best Lagrangian multipliers) */ 340 SCIP_Real radiusinit; /**< initial radius of the ball used in stabilization of Lagrangian multipliers */ 341 SCIP_Real radiusmax; /**< maximum radius of the ball used in stabilization of Lagrangian multipliers */ 342 SCIP_Real radiusmin; /**< minimum radius of the ball used in stabilization of Lagrangian multipliers */ 343 SCIP_Real constant; /**< a constant for stablity center based stabilization of Lagrangian multipliers */ 344 SCIP_Real radiusupdateweight; /**< multiplier to evaluate cut violation score used for updating ball radius */ 345 }; 346 347 348 /* 349 * Local methods 350 */ 351 352 /** start the diving mode for solving LPs corresponding to the Lagrangian dual with fixed multipliers in the subgradient 353 * loop of the separator, and update some sepadata values */ 354 static 355 SCIP_RETCODE createLPWithSoftCuts( 356 SCIP* scip, /**< SCIP data structure */ 357 SCIP_SEPADATA* sepadata /**< separator data structure */ 358 ) 359 { 360 SCIP_ROW** rows; 361 SCIP_COL** cols; 362 int nruns; 363 int runnum; 364 int nrows; 365 int ncols; 366 unsigned int nintcols; 367 SCIP_Longint nrootlpiters; 368 369 assert(scip != NULL); 370 assert(sepadata != NULL); 371 372 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) ); 373 SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) ); 374 runnum = SCIPgetNRuns(scip); /* current run number of SCIP (starts from 1) indicating restarts */ 375 nruns = sepadata->nrunsinsoftcutslp; /* previous run number of SCIP in which the diving LP was created */ 376 nintcols = 0; 377 378 /* start diving mode so that the diving LP can be used for all subgradient iterations */ 379 SCIP_CALL( SCIPstartDive(scip) ); 380 381 /* store LPI pointer to be able to use LPI functions directly (e.g., setting time limit) */ 382 SCIP_CALL( SCIPgetLPI(scip, &(sepadata->lpiwithsoftcuts)) ); 383 384 /* if called for the first time in a restart (including the initial run), set certain sepadata values */ 385 if( nruns != runnum ) 386 { 387 /* get number of LP iterations of root node's first LP solving */ 388 nrootlpiters = SCIPgetNRootFirstLPIterations(scip); 389 390 /* calculate maximum number of LP iterations allowed for all separation calls in the root node */ 391 if( (sepadata->rootlpiterlimitfactor >= 0.0) && !SCIPisInfinity(scip, sepadata->rootlpiterlimitfactor) ) 392 { 393 sepadata->nmaxrootlpiters = (int)(sepadata->rootlpiterlimitfactor * nrootlpiters); 394 } 395 else 396 { 397 sepadata->nmaxrootlpiters = -1; /* no finite limit */ 398 } 399 400 /* calculate maximum number of LP iterations allowed for all separation calls in the entire tree */ 401 if( (sepadata->totallpiterlimitfactor >= 0.0) && !SCIPisInfinity(scip, sepadata->totallpiterlimitfactor) ) 402 { 403 sepadata->nmaxtotallpiters = (int)(sepadata->totallpiterlimitfactor * nrootlpiters); 404 } 405 else 406 { 407 sepadata->nmaxtotallpiters = -1; /* no finite limit */ 408 } 409 410 /* calculate maximum number of LP iterations allowed per separation call */ 411 if( (sepadata->perroundlpiterlimitfactor >= 0.0) && !SCIPisInfinity(scip, sepadata->perroundlpiterlimitfactor) ) 412 { 413 sepadata->nmaxperroundlpiters = (int)(sepadata->perroundlpiterlimitfactor * nrootlpiters); 414 } 415 else 416 { 417 sepadata->nmaxperroundlpiters = -1; /* no finite limit */ 418 } 419 420 /* update maximum number of LP iterations allowed per separation call using absolute limits */ 421 if( sepadata->perroundnmaxlpiters > 0 ) 422 { 423 sepadata->nmaxperroundlpiters = ((sepadata->nmaxperroundlpiters >= 0) ? MIN(sepadata->nmaxperroundlpiters, 424 sepadata->perroundnmaxlpiters) : sepadata->perroundnmaxlpiters); 425 } 426 427 /* set maximum number of cuts allowed to generate per round in root and non-root nodes as well as the total tree */ 428 for( int i = 0; i < ncols; ++i ) 429 { 430 nintcols += SCIPcolIsIntegral(cols[i]); 431 } 432 sepadata->nmaxperroundcutsroot = (int)(sepadata->perroundcutsfactorroot * nintcols); 433 sepadata->nmaxperroundcuts = (int)(sepadata->perroundcutsfactor * nintcols); 434 435 if( sepadata->ncalls == 0 ) 436 { 437 sepadata->nmaxtotalcuts = (int)(sepadata->totalcutsfactor * nintcols); 438 sepadata->ntotalcuts = 0; 439 } 440 441 /* update the run number of solving to represent the restart number in which the above limits were set */ 442 sepadata->nrunsinsoftcutslp = runnum; 443 } 444 445 return SCIP_OKAY; 446 } 447 448 /** end the diving mode that was used for solving LPs corresponding to the Lagrangian dual with fixed multipliers */ 449 static 450 SCIP_RETCODE deleteLPWithSoftCuts( 451 SCIP* scip, /**< SCIP data structure */ 452 SCIP_SEPADATA* sepadata /**< separator data structure */ 453 ) 454 { 455 assert(scip != NULL); 456 assert(sepadata != NULL); 457 458 SCIP_CALL( SCIPendDive(scip) ); 459 460 return SCIP_OKAY; 461 } 462 463 /** set up LP interface to solve LPs in the (outer) main loop of the relax-and-cut algorithm; these LPs are built by 464 * adding all the generated cuts to the node relaxation */ 465 /* @todo add lpi iters to global statistics */ 466 static 467 SCIP_RETCODE createLPWithHardCuts( 468 SCIP* scip, /**< SCIP data structure */ 469 SCIP_SEPADATA* sepadata, /**< separator data structure */ 470 SCIP_ROW** cuts, /**< generated cuts to be added to the LP */ 471 int ncuts /**< number of generated cuts to be added to the LP */ 472 ) 473 { 474 SCIP_ROW** rows; 475 SCIP_COL** cols; 476 SCIP_Real* collb; 477 SCIP_Real* colub; 478 SCIP_Real* colobj; 479 SCIP_Real* rowlhs; 480 SCIP_Real* rowrhs; 481 SCIP_Real rowconst; 482 SCIP_Real* rowvals; 483 SCIP_Real* rowval; 484 SCIP_COL** rowcols; 485 SCIP_LPI* wslpi; 486 SCIP_LPISTATE* lpistate; 487 BMS_BLKMEM* blkmem; 488 SCIP_Real pinf; 489 SCIP_Real ninf; 490 int nrows; 491 int ncols; 492 int nrownonz; 493 int collppos; 494 int* rowcolinds; 495 int* rowbegs; 496 497 assert(scip != NULL); 498 assert(sepadata != NULL); 499 500 SCIP_LPI** lpi = &(sepadata->lpiwithhardcuts); 501 blkmem = SCIPblkmem(scip); 502 503 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) ); 504 SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) ); 505 506 /* if this function is called for the first time in this separation round, then create an LPI and add cols & rows */ 507 if( ncuts == 0 ) 508 { 509 if( *lpi != NULL ) 510 { 511 SCIP_CALL( SCIPlpiFree(lpi) ); 512 *lpi = NULL; 513 } 514 assert(*lpi == NULL); 515 516 /* create an LPI with appropriate objective sense */ 517 if( SCIPgetObjsense(scip) == SCIP_OBJSENSE_MAXIMIZE ) 518 { 519 SCIP_CALL( SCIPlpiCreate(lpi, SCIPgetMessagehdlr(scip), "node LP with generated cuts", SCIP_OBJSEN_MAXIMIZE) ); 520 } 521 else 522 { 523 SCIP_CALL( SCIPlpiCreate(lpi, SCIPgetMessagehdlr(scip), "node LP with generated cuts", SCIP_OBJSEN_MINIMIZE) ); 524 } 525 526 /* add cols to the LP interface */ 527 SCIP_CALL( SCIPallocBufferArray(scip, &colobj, ncols) ); 528 SCIP_CALL( SCIPallocBufferArray(scip, &collb, ncols) ); 529 SCIP_CALL( SCIPallocBufferArray(scip, &colub, ncols) ); 530 /* gather required column information */ 531 for( int i = 0; i < ncols; ++i ) 532 { 533 colobj[i] = SCIPcolGetObj(cols[i]); 534 collb[i] = SCIPcolGetLb(cols[i]); 535 colub[i] = SCIPcolGetUb(cols[i]); 536 } 537 /* add cols */ 538 SCIP_CALL( SCIPlpiAddCols(*lpi, ncols, colobj, collb, colub, NULL, 0, NULL, NULL, NULL) ); 539 SCIPfreeBufferArray(scip, &colub); 540 SCIPfreeBufferArray(scip, &collb); 541 SCIPfreeBufferArray(scip, &colobj); 542 543 /* add rows to the LP interface */ 544 /* find number of nonzeroes in rows */ 545 nrownonz = 0; 546 for( int i = 0; i < nrows; ++i ) 547 { 548 assert(!(SCIPisInfinity(scip, -SCIProwGetLhs(rows[i])) && SCIPisInfinity(scip, SCIProwGetRhs(rows[i])))); 549 nrownonz += SCIProwGetNLPNonz(rows[i]); 550 } 551 SCIP_CALL( SCIPallocBufferArray(scip, &rowcolinds, nrownonz) ); 552 SCIP_CALL( SCIPallocBufferArray(scip, &rowvals, nrownonz) ); 553 SCIP_CALL( SCIPallocBufferArray(scip, &rowbegs, nrows + 1) ); 554 SCIP_CALL( SCIPallocBufferArray(scip, &rowlhs, nrows) ); 555 SCIP_CALL( SCIPallocBufferArray(scip, &rowrhs, nrows) ); 556 /* gather required row information */ 557 rowbegs[0] = 0; 558 pinf = SCIPlpiInfinity(*lpi); 559 ninf = -SCIPlpiInfinity(*lpi); 560 for( int i = 0; i < nrows; ++i ) 561 { 562 nrownonz = SCIProwGetNLPNonz(rows[i]); 563 assert(nrownonz <= ncols); 564 rowval = SCIProwGetVals(rows[i]); 565 rowcols = SCIProwGetCols(rows[i]); 566 567 rowbegs[i + 1] = rowbegs[i] + nrownonz; 568 rowconst = SCIProwGetConstant(rows[i]); 569 rowlhs[i] = SCIPisInfinity(scip, -SCIProwGetLhs(rows[i])) ? ninf : SCIProwGetLhs(rows[i]) - rowconst; 570 rowrhs[i] = SCIPisInfinity(scip, SCIProwGetRhs(rows[i])) ? pinf : SCIProwGetRhs(rows[i]) - rowconst; 571 572 for( int j = 0; j < nrownonz; ++j ) 573 { 574 collppos = SCIPcolGetLPPos(rowcols[j]); 575 assert(collppos >= 0); 576 assert(collppos <= ncols); 577 578 rowcolinds[rowbegs[i] + j] = collppos; 579 rowvals[rowbegs[i] + j] = rowval[j]; 580 } 581 } 582 /* add rows */ 583 SCIP_CALL( SCIPlpiAddRows(*lpi, nrows, rowlhs, rowrhs, NULL, rowbegs[nrows], rowbegs, rowcolinds, rowvals) ); 584 585 /* get warm starting info */ 586 SCIP_CALL( SCIPgetLPI(scip, &wslpi) ); 587 SCIP_CALL( SCIPlpiGetState(wslpi, blkmem, &lpistate) ); 588 } 589 /* if there are any cuts, then add the cuts that were not added earlier to the LPI */ 590 else 591 { 592 assert(nrows + ncuts >= sepadata->nrowsinhardcutslp); 593 594 /* get warm starting info */ 595 wslpi = *lpi; 596 SCIP_CALL( SCIPlpiGetState(wslpi, blkmem, &lpistate) ); 597 598 /* find number of nonzeros in cuts and allocate memory */ 599 nrownonz = 0; 600 pinf = SCIPlpiInfinity(*lpi); 601 ninf = -SCIPlpiInfinity(*lpi); 602 for( int i = sepadata->nrowsinhardcutslp - nrows; i < ncuts; ++i ) 603 { 604 assert(!(SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i])) && SCIPisInfinity(scip, SCIProwGetRhs(cuts[i])))); 605 assert(SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i]))); 606 assert(!SCIPisInfinity(scip, SCIProwGetRhs(cuts[i]))); 607 nrownonz += SCIProwGetNNonz(cuts[i]); 608 } 609 SCIP_CALL( SCIPallocBufferArray(scip, &rowcolinds, nrownonz) ); 610 SCIP_CALL( SCIPallocBufferArray(scip, &rowvals, nrownonz) ); 611 SCIP_CALL( SCIPallocBufferArray(scip, &rowbegs, (ncuts - sepadata->nrowsinhardcutslp + nrows) + 1) ); 612 SCIP_CALL( SCIPallocBufferArray(scip, &rowlhs, (ncuts - sepadata->nrowsinhardcutslp + nrows)) ); 613 SCIP_CALL( SCIPallocBufferArray(scip, &rowrhs, (ncuts - sepadata->nrowsinhardcutslp + nrows)) ); 614 615 /* gather required cut information */ 616 rowbegs[0] = 0; 617 for( int i = sepadata->nrowsinhardcutslp - nrows; i < ncuts; ++i ) 618 { 619 nrownonz = SCIProwGetNNonz(cuts[i]); 620 assert(nrownonz <= ncols); 621 rowval = SCIProwGetVals(cuts[i]); 622 rowcols = SCIProwGetCols(cuts[i]); 623 624 rowbegs[i - sepadata->nrowsinhardcutslp + nrows + 1] = rowbegs[i - sepadata->nrowsinhardcutslp + nrows] + 625 nrownonz; 626 rowconst = SCIProwGetConstant(cuts[i]); 627 rowlhs[i - sepadata->nrowsinhardcutslp + nrows] = SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i])) ? ninf : 628 SCIProwGetLhs(cuts[i]) - rowconst; 629 rowrhs[i - sepadata->nrowsinhardcutslp + nrows] = SCIPisInfinity(scip, SCIProwGetRhs(cuts[i])) ? pinf : 630 SCIProwGetRhs(cuts[i]) - rowconst; 631 632 for( int j = 0; j < nrownonz; ++j ) 633 { 634 collppos = SCIPcolGetLPPos(rowcols[j]); 635 assert(collppos >= 0); 636 assert(collppos <= ncols); 637 638 rowcolinds[rowbegs[i - sepadata->nrowsinhardcutslp + nrows] + j] = collppos; 639 rowvals[rowbegs[i - sepadata->nrowsinhardcutslp + nrows] + j] = rowval[j]; 640 } 641 } 642 643 /* add cuts */ 644 SCIP_CALL( SCIPlpiAddRows(*lpi, (ncuts - sepadata->nrowsinhardcutslp + nrows), rowlhs, rowrhs, NULL, 645 rowbegs[(ncuts - sepadata->nrowsinhardcutslp + nrows)], rowbegs, rowcolinds, rowvals) ); 646 } 647 648 /* set warm starting basis */ 649 SCIP_CALL( SCIPlpiSetState(*lpi, blkmem, lpistate) ); 650 651 /* reset remaining sepadata values */ 652 sepadata->nrowsinhardcutslp = nrows + ncuts; 653 654 /* free memory */ 655 SCIP_CALL( SCIPlpiFreeState(*lpi, blkmem, &lpistate) ); 656 SCIPfreeBufferArray(scip, &rowrhs); 657 SCIPfreeBufferArray(scip, &rowlhs); 658 SCIPfreeBufferArray(scip, &rowbegs); 659 SCIPfreeBufferArray(scip, &rowvals); 660 SCIPfreeBufferArray(scip, &rowcolinds); 661 662 return SCIP_OKAY; 663 } 664 665 /** free separator data */ 666 static 667 SCIP_RETCODE sepadataFree( 668 SCIP* scip, /**< SCIP data structure */ 669 SCIP_SEPADATA** sepadata /**< separator data structure */ 670 ) 671 { 672 assert(scip != NULL); 673 assert(sepadata != NULL); 674 assert(*sepadata != NULL); 675 676 if( (*sepadata)->lpiwithhardcuts != NULL ) 677 { 678 SCIP_CALL( SCIPlpiFree(&((*sepadata)->lpiwithhardcuts)) ); 679 } 680 681 (*sepadata)->nrowsinhardcutslp = 0; 682 (*sepadata)->nrunsinsoftcutslp = 0; 683 (*sepadata)->ncalls = 0; 684 (*sepadata)->nmaxperroundlpiters = 0; 685 (*sepadata)->nmaxrootlpiters = 0; 686 (*sepadata)->nrootlpiters = 0; 687 (*sepadata)->nmaxtotallpiters = 0; 688 (*sepadata)->ntotallpiters = 0; 689 (*sepadata)->nmaxperroundcutsroot = 0; 690 (*sepadata)->nmaxperroundcuts = 0; 691 (*sepadata)->nmaxtotalcuts = 0; 692 (*sepadata)->ntotalcuts = 0; 693 694 SCIPfreeBlockMemory(scip, sepadata); 695 696 return SCIP_OKAY; 697 } 698 699 /** update mu parameter which is used as a factor in the step length calculation; refer to the top of the file for a 700 * description of the formula. 701 */ 702 /* @todo some adaptive strategy like constant after certain changes? */ 703 static 704 SCIP_RETCODE updateMuSteplengthParam( 705 SCIP* scip, /**< SCIP data structure */ 706 SCIP_SEPADATA* sepadata, /**< separator data structure */ 707 int subgradientiternum, /**< subgradient iteration number */ 708 SCIP_Real ubparam, /**< estimate of the optimal Lagrangian dual value */ 709 SCIP_Real* lagrangianvals, /**< vector of Lagrangian values found so far */ 710 SCIP_Real bestlagrangianval, /**< best Lagrangian value found so far */ 711 SCIP_Real avglagrangianval, /**< rolling average of the Lagrangian values found so far */ 712 SCIP_Real* muparam, /**< mu parameter to be updated */ 713 SCIP_Bool* backtrack /**< whether mu parameter has been backtracked */ 714 ) 715 { 716 SCIP_Real delta; 717 SCIP_Real deltaslab1ub; 718 SCIP_Real deltaslab2ub; 719 SCIP_Real muslab1factor; 720 SCIP_Real muslab2factor; 721 SCIP_Real muslab3factor; 722 int maxiters; 723 int i; 724 725 *backtrack = FALSE; 726 727 /* update the mu parameter only if it is not set to be a constant value */ 728 if( !sepadata->muparamconst ) 729 { 730 delta = ubparam - bestlagrangianval; 731 deltaslab1ub = MIN(sepadata->deltaslab1ub, sepadata->deltaslab2ub); 732 deltaslab2ub = MAX(sepadata->deltaslab1ub, sepadata->deltaslab2ub); 733 /* ensure that the ordering of different user input parameters is as expected */ 734 if( SCIPisPositive(scip, sepadata->muslab1factor - sepadata->muslab2factor) ) 735 { 736 if( SCIPisPositive(scip, sepadata->muslab2factor - sepadata->muslab3factor) ) 737 { 738 muslab1factor = sepadata->muslab1factor; 739 muslab2factor = sepadata->muslab2factor; 740 muslab3factor = sepadata->muslab3factor; 741 } 742 else 743 { 744 if( SCIPisPositive(scip, sepadata->muslab1factor - sepadata->muslab3factor) ) 745 { 746 muslab1factor = sepadata->muslab1factor; 747 muslab2factor = sepadata->muslab3factor; 748 muslab3factor = sepadata->muslab2factor; 749 } 750 else 751 { 752 muslab1factor = sepadata->muslab3factor; 753 muslab2factor = sepadata->muslab1factor; 754 muslab3factor = sepadata->muslab2factor; 755 } 756 } 757 } 758 else 759 { 760 if( SCIPisPositive(scip, sepadata->muslab1factor - sepadata->muslab3factor) ) 761 { 762 muslab1factor = sepadata->muslab2factor; 763 muslab2factor = sepadata->muslab1factor; 764 muslab3factor = sepadata->muslab3factor; 765 } 766 else 767 { 768 if( SCIPisPositive(scip, sepadata->muslab2factor - sepadata->muslab3factor) ) 769 { 770 muslab1factor = sepadata->muslab2factor; 771 muslab2factor = sepadata->muslab3factor; 772 muslab3factor = sepadata->muslab1factor; 773 } 774 else 775 { 776 muslab1factor = sepadata->muslab3factor; 777 muslab2factor = sepadata->muslab2factor; 778 muslab3factor = sepadata->muslab1factor; 779 } 780 } 781 } 782 783 maxiters = MIN(sepadata->nmaxconsecitersformuupdate, sepadata->nmaxlagrangianvalsforavg); 784 i = -1; 785 786 /* if certain number of iterations are done, then check for a possibility of backtracking and apply accordingly */ 787 if( subgradientiternum >= maxiters ) 788 { 789 for( i = subgradientiternum - maxiters; i < subgradientiternum; i++ ) 790 { 791 if( SCIPisGE(scip, lagrangianvals[i], bestlagrangianval - delta) ) 792 break; 793 } 794 795 if( i == subgradientiternum ) 796 { 797 *muparam *= sepadata->mubacktrackfactor; 798 *backtrack = TRUE; 799 } 800 } 801 802 /* update mu parameter based on the different between best and average Lagrangian values */ 803 if( (subgradientiternum < maxiters) || (i >= 0 && i < subgradientiternum) ) 804 { 805 if( bestlagrangianval - avglagrangianval < deltaslab1ub * delta ) 806 *muparam *= muslab1factor; 807 else if( bestlagrangianval - avglagrangianval < deltaslab2ub * delta ) 808 *muparam *= muslab2factor; 809 else 810 *muparam *= muslab3factor; 811 } 812 813 /* reset the mu parameter to within its bounds */ 814 *muparam = MAX(*muparam, sepadata->muparamlb); 815 *muparam = MIN(*muparam, sepadata->muparamub); 816 } 817 818 return SCIP_OKAY; 819 } 820 821 /** update subgradient, i.e., residuals of generated cuts */ 822 /* @note: assumed that \f$i^{th}\f$ cut is of the form \f${\alpha^i}^T x \leq {\alpha^i_0}\f$ */ 823 static 824 void updateSubgradient( 825 SCIP* scip, /**< SCIP data structure */ 826 SCIP_SOL* sol, /**< LP solution used in updating subgradient vector */ 827 SCIP_ROW** cuts, /**< cuts generated so far */ 828 int ncuts, /**< number of cuts generated so far */ 829 SCIP_Real* subgradient, /**< vector of subgradients to be updated */ 830 SCIP_Real* dualvector, /**< Lagrangian multipliers */ 831 SCIP_Bool* subgradientzero, /**< whether the subgradient vector is all zero */ 832 int* ncutviols, /**< number of violations of generated cuts */ 833 SCIP_Real* maxcutviol, /**< maximum violation of generated cuts */ 834 int* nnzsubgradientdualprod, /**< number of nonzero products of subgradient vector and Lagrangian multipliers (i.e., 835 number of complementarity slackness violations) */ 836 SCIP_Real* maxnzsubgradientdualprod /**< maximum value of nonzero products of subgradient vector and Lagrangian multipliers 837 (i.e., maximum value of complementarity slackness violations) */ 838 ) 839 { 840 int nzerosubgradient; 841 SCIP_Real term; 842 843 assert(subgradientzero != NULL); 844 assert(ncutviols != NULL); 845 assert(maxcutviol != NULL); 846 assert(nnzsubgradientdualprod != NULL); 847 assert(maxnzsubgradientdualprod != NULL); 848 849 *ncutviols = 0; 850 *maxcutviol = 0.0; 851 *nnzsubgradientdualprod = 0; 852 *maxnzsubgradientdualprod = 0.0; 853 nzerosubgradient = 0; 854 *subgradientzero = FALSE; 855 856 /* for each cut, calculate the residual along with various violation metrics */ 857 for( int i = 0; i < ncuts; i++ ) 858 { 859 assert(SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i]))); 860 assert(!SCIPisInfinity(scip, SCIProwGetRhs(cuts[i]))); 861 subgradient[i] = SCIPgetRowSolActivity(scip, cuts[i], sol) + SCIProwGetConstant(cuts[i]) - SCIProwGetRhs(cuts[i]); 862 if( SCIPisFeasZero(scip, subgradient[i]) ) 863 { 864 subgradient[i] = 0.0; 865 nzerosubgradient++; 866 } 867 else 868 { 869 /* check for cut violation */ 870 if( SCIPisFeasPositive(scip, subgradient[i]) ) 871 { 872 (*ncutviols)++; 873 *maxcutviol = MAX(*maxcutviol, subgradient[i]); 874 } 875 876 /* check for violation of complementarity slackness associated with the cut */ 877 if( !SCIPisZero(scip, subgradient[i] * dualvector[i]) ) 878 { 879 (*nnzsubgradientdualprod)++; 880 term = REALABS(subgradient[i] * dualvector[i]); 881 *maxnzsubgradientdualprod = MAX(*maxnzsubgradientdualprod, term); 882 } 883 } 884 } 885 886 /* indicator for all zero subgradient vector */ 887 if( nzerosubgradient == ncuts ) 888 { 889 *subgradientzero = TRUE; 890 } 891 } 892 893 /** update Lagrangian value, i.e., optimal value of the Lagrangian dual with fixed multipliers */ 894 static 895 void updateLagrangianValue( 896 SCIP* scip, /**< SCIP data structure */ 897 SCIP_Real objval, /**< objective value of the Lagrangian dual with fixed multipliers */ 898 SCIP_Real* dualvector, /**< Lagrangian multipliers */ 899 SCIP_ROW** cuts, /**< cuts generated so far */ 900 int ncuts, /**< number of cuts generated so far */ 901 SCIP_Real* lagrangianval /**< Lagrangian value to be updated */ 902 ) 903 { 904 *lagrangianval = objval; 905 906 for( int i = 0; i < ncuts; i++ ) 907 { 908 assert(SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i]))); 909 assert(!SCIPisInfinity(scip, SCIProwGetRhs(cuts[i]))); 910 *lagrangianval += dualvector[i] * (SCIProwGetConstant(cuts[i]) - SCIProwGetRhs(cuts[i])); 911 } 912 } 913 914 /** update step length based on various input arguments; refer to the top of the file for an expression */ 915 static 916 void updateStepLength( 917 SCIP* scip, /**< SCIP data structure */ 918 SCIP_Real muparam, /**< mu parameter used as a factor for step length */ 919 SCIP_Real ubparam, /**< estimate of the optimal Lagrangian dual value */ 920 SCIP_Real lagrangianval, /**< Lagrangian value */ 921 SCIP_Real* subgradient, /**< subgradient vector */ 922 int ncuts, /**< number of cuts generated so far */ 923 SCIP_Real* steplength /**< step length to be updated */ 924 ) 925 { 926 SCIP_Real normsquared = 0.0; 927 928 for( int i = 0; i < ncuts; i++ ) 929 { 930 normsquared += SQR(subgradient[i]); 931 } 932 933 if( !SCIPisFeasZero(scip, normsquared) ) 934 { 935 *steplength = (muparam * (ubparam - lagrangianval))/(normsquared); /*lint !e795*/ 936 } 937 } 938 939 /** update the ball radius (based on various violation metrics) that is used for stabilization of Lagrangian multipliers */ 940 static 941 void updateBallRadius( 942 SCIP* scip, /**< SCIP data structure */ 943 SCIP_SEPADATA* sepadata, /**< separator data structure */ 944 SCIP_Real maxviolscore, /**< weighted average of maximum value of generated cut violations and maximum value of 945 complementarity slackness violations, in the current iteration */ 946 SCIP_Real maxviolscoreold, /**< weighted average of maximum value of generated cut violations and maximum value of 947 complementarity slackness violations, in the previous iteration */ 948 SCIP_Real nviolscore, /**< weighted average of number of generated cut violations and number of complementarity 949 slackness violations, in the current iteration */ 950 SCIP_Real nviolscoreold, /**< weighted average of number of generated cut violations and number of complementarity 951 slackness violations, in the previous iteration */ 952 int nlpiters, /**< number of LP iterations taken for solving the Lagrangian dual with fixed multipliers in 953 current iteration */ 954 SCIP_Real* ballradius /**< norm ball radius to be updated */ 955 ) 956 { 957 SCIP_Bool maxviolscoreimproved; 958 SCIP_Bool nviolscoreimproved; 959 960 assert(ballradius != NULL); 961 962 maxviolscoreimproved = !SCIPisNegative(scip, maxviolscoreold - maxviolscore); 963 nviolscoreimproved = !SCIPisNegative(scip, nviolscoreold - nviolscore); 964 965 if( maxviolscoreimproved && nviolscoreimproved ) 966 { 967 /* both the maximum violation and number of violations scores have become better, so, increase the radius */ 968 if( sepadata->optimalfacepriority <= 1 ) 969 { 970 *ballradius *= 2.0; 971 *ballradius = MIN(*ballradius, sepadata->radiusmax); 972 } 973 else 974 { 975 *ballradius *= 1.5; 976 *ballradius = MIN(*ballradius, sepadata->radiusmax/2.0); 977 } 978 } 979 else if( !maxviolscoreimproved && !nviolscoreimproved ) 980 { 981 /* both the maximum violation and number of violations scores have become worse, so, decrease the radius */ 982 *ballradius *= 0.5; 983 *ballradius = MAX(*ballradius, sepadata->radiusmin); 984 } 985 else if( nlpiters == 0 ) 986 { 987 /* only one among the maximum violation and number of violations scores has become better, and the LP basis did 988 * not change (i.e., nlpters = 0), so, increase the radius slightly */ 989 if( sepadata->optimalfacepriority <= 1 ) 990 { 991 *ballradius *= 1.5; 992 *ballradius = MIN(*ballradius, sepadata->radiusmax); 993 } 994 else 995 { 996 *ballradius *= 1.2; 997 *ballradius = MIN(*ballradius, sepadata->radiusmax/2.0); 998 } 999 } 1000 } 1001 1002 /** projection of Lagrangian multipliers onto L1-norm ball. This algorithm is based on the following article. 1003 * 1004 * Condat L. (2016).@n 1005 * Fast projection onto the simplex and the \f$l_1\f$ ball.@n 1006 * Mathematical Programming, 158, 1-2, 575-585. 1007 * 1008 */ 1009 static 1010 SCIP_RETCODE l1BallProjection( 1011 SCIP* scip, /**< SCIP data structure */ 1012 SCIP_Real* dualvector, /**< Lagrangian multipliers to be projected onto L1-norm vall */ 1013 int dualvectorlen, /**< length of the Lagrangian multipliers vector */ 1014 SCIP_Real radius /**< radius of the L1-norm ball */ 1015 ) 1016 { 1017 SCIP_Real* temp1vals; 1018 SCIP_Real* temp2vals; 1019 SCIP_Real pivotparam; 1020 SCIP_Real val; 1021 SCIP_Real term; 1022 SCIP_Bool temp1changed; 1023 int ntemp1removed; 1024 int* temp1inds; 1025 int* temp2inds; 1026 int temp1len; 1027 int temp2len; 1028 1029 assert(!SCIPisNegative(scip, radius)); 1030 val = REALABS(dualvector[0]); 1031 /* calculate the L1-norm of the Lagrangian multipliers */ 1032 for( int i = 1; i < dualvectorlen; i++ ) 1033 { 1034 val += REALABS(dualvector[i]); 1035 } 1036 1037 /* if the vector of Lagrangian multipliers lies outside the L1-norm ball, then do the projection */ 1038 if( SCIPisGT(scip, val, radius) ) 1039 { 1040 SCIP_CALL( SCIPallocCleanBufferArray(scip, &temp1vals, dualvectorlen) ); 1041 SCIP_CALL( SCIPallocCleanBufferArray(scip, &temp2vals, dualvectorlen) ); 1042 SCIP_CALL( SCIPallocBufferArray(scip, &temp1inds, dualvectorlen) ); 1043 SCIP_CALL( SCIPallocBufferArray(scip, &temp2inds, dualvectorlen) ); 1044 for( int i = 0; i < dualvectorlen; i++ ) 1045 { 1046 temp1inds[i] = -1; 1047 temp2inds[i] = -1; 1048 } 1049 temp2len = 0; 1050 1051 temp1vals[0] = REALABS(dualvector[0]); 1052 temp1inds[0] = 0; 1053 temp1len = 1; 1054 pivotparam = REALABS(dualvector[0]) - radius; 1055 1056 for( int i = 1; i < dualvectorlen; i++ ) 1057 { 1058 if( SCIPisGT(scip, REALABS(dualvector[i]), pivotparam) ) 1059 { 1060 pivotparam += ((REALABS(dualvector[i]) - pivotparam) / (temp1len + 1)); 1061 if( SCIPisGT(scip, pivotparam, REALABS(dualvector[i]) - radius) ) 1062 { 1063 temp1vals[temp1len] = REALABS(dualvector[i]); 1064 temp1inds[temp1len] = i; 1065 temp1len++; 1066 } 1067 else 1068 { 1069 for( int j = 0; j < temp1len; j++ ) 1070 { 1071 temp2vals[temp2len + j] = temp1vals[j]; 1072 temp2inds[temp2len + j] = temp1inds[j]; 1073 } 1074 temp2len += temp1len; 1075 temp1vals[0] = REALABS(dualvector[i]); 1076 temp1inds[0] = i; 1077 temp1len = 1; 1078 pivotparam = REALABS(dualvector[i]) - radius; 1079 } 1080 } 1081 } 1082 1083 for( int i = 0; i < temp2len; i++ ) 1084 { 1085 if( SCIPisGT(scip, temp2vals[i], pivotparam) ) 1086 { 1087 temp1vals[temp1len] = temp2vals[i]; 1088 temp1inds[temp1len] = temp2inds[i]; 1089 temp1len++; 1090 pivotparam += ((temp2vals[i] - pivotparam) / temp1len); 1091 } 1092 } 1093 1094 temp1changed = TRUE; 1095 ntemp1removed = 0; 1096 while( temp1changed ) 1097 { 1098 temp1changed = FALSE; 1099 1100 for( int i = 0; i < temp1len; i++ ) 1101 { 1102 if( (temp1inds[i] >= 0) && SCIPisLE(scip, temp1vals[i], pivotparam) ) 1103 { 1104 temp1inds[i] = -1; 1105 temp1changed = TRUE; 1106 ntemp1removed++; 1107 assert(temp1len - ntemp1removed > 0); 1108 pivotparam += ((pivotparam - temp1vals[i]) / (temp1len - ntemp1removed)); 1109 } 1110 } 1111 } 1112 1113 for( int i = 0; i < dualvectorlen; i++ ) 1114 { 1115 term = REALABS(dualvector[i]); 1116 val = MAX(term - pivotparam, 0.0); 1117 1118 if( SCIPisPositive(scip, dualvector[i]) ) 1119 { 1120 dualvector[i] = val; 1121 } 1122 else if( SCIPisNegative(scip, dualvector[i]) ) 1123 { 1124 dualvector[i] = -val; 1125 } 1126 } 1127 1128 /* free memory */ 1129 for( int i = 0; i < dualvectorlen; i++ ) 1130 { 1131 temp2vals[i] = 0.0; 1132 temp1vals[i] = 0.0; 1133 } 1134 SCIPfreeBufferArray(scip, &temp2inds); 1135 SCIPfreeBufferArray(scip, &temp1inds); 1136 SCIPfreeCleanBufferArray(scip, &temp2vals); 1137 SCIPfreeCleanBufferArray(scip, &temp1vals); 1138 } 1139 1140 return SCIP_OKAY; 1141 } 1142 1143 /** projection of Lagrangian multipliers onto L2-norm ball */ 1144 static 1145 void l2BallProjection( 1146 SCIP* scip, /**< SCIP data structure */ 1147 SCIP_Real* dualvector, /**< Lagrangian multipliers to be projected onto L2-norm vall */ 1148 int dualvectorlen, /**< length of the Lagrangian multipliers vector */ 1149 SCIP_Real radius /**< radius of the L2-norm ball */ 1150 ) 1151 { 1152 SCIP_Real l2norm; 1153 SCIP_Real factor; 1154 1155 assert(!SCIPisNegative(scip, radius)); 1156 1157 l2norm = 0.0; 1158 /* calculate the L2-norm of the Lagrangian multipliers */ 1159 for( int i = 0; i < dualvectorlen; i++ ) 1160 { 1161 l2norm += SQR(dualvector[i]); 1162 } 1163 l2norm = SQRT(l2norm); 1164 factor = radius/(1.0 + l2norm); 1165 1166 /* if the vector of Lagrangian multipliers is outside the L2-norm ball, then do the projection */ 1167 if( SCIPisGT(scip, l2norm, radius) && SCIPisLT(scip, factor, 1.0) ) 1168 { 1169 for( int i = 0; i < dualvectorlen; i++ ) 1170 { 1171 dualvector[i] *= factor; 1172 } 1173 } 1174 } 1175 1176 /** projection of Lagrangian multipliers onto L_infinity-norm ball */ 1177 static 1178 void linfBallProjection( 1179 SCIP* scip, /**< SCIP data structure */ 1180 SCIP_Real* dualvector, /**< Lagrangian multipliers to be projected onto L_infinity-norm vall */ 1181 int dualvectorlen, /**< length of the Lagrangian multipliers vector */ 1182 SCIP_Real radius /**< radius of the L_infinity-norm ball */ 1183 ) 1184 { 1185 assert(!SCIPisNegative(scip, radius)); 1186 1187 /* if the vector of Lagrangian multipliers is outside the L_infinity-norm ball, then do the projection */ 1188 for( int i = 0; i < dualvectorlen; i++ ) 1189 { 1190 if( SCIPisLT(scip, dualvector[i], -radius) ) 1191 { 1192 dualvector[i] = -radius; 1193 } 1194 else if( SCIPisGT(scip, dualvector[i], radius) ) 1195 { 1196 dualvector[i] = radius; 1197 } 1198 } 1199 } 1200 1201 /** weighted Lagrangian multipliers based on a given vector as stability center */ 1202 /* @todo calculate weight outside this function and pass it (so that this function becomes generic and independent of 1203 * the terminology related to best Lagrangian multipliers) 1204 */ 1205 static 1206 SCIP_RETCODE weightedDualVector( 1207 SCIP* scip, /**< SCIP data structure */ 1208 SCIP_SEPADATA* sepadata, /**< separator data structure */ 1209 SCIP_Real* dualvector, /**< Lagrangian multipliers */ 1210 int dualvectorlen, /**< length of the Lagrangian multipliers vector */ 1211 SCIP_Real* stabilitycenter, /**< stability center (i.e., core vector of Lagrangian multipliers) */ 1212 int stabilitycenterlen, /**< length of the stability center */ 1213 int nbestdualupdates, /**< number of best Lagrangian values found so far */ 1214 int totaliternum /**< total number of iterations of the relax-and-cut algorithm performed so far */ 1215 ) 1216 { 1217 SCIP_Real constant; 1218 SCIP_Real weight; 1219 SCIP_Real alpha; 1220 1221 constant = MAX(2.0, sepadata->constant); 1222 /* weight factor from the literature on Dantzig-Wolfe decomposition stabilization schemes */ 1223 weight = MIN(constant, (totaliternum + 1 + nbestdualupdates) / 2.0); 1224 alpha = 1.0 / weight; 1225 1226 assert(dualvectorlen >= stabilitycenterlen); 1227 1228 /* weighted Lagrangian multipliers */ 1229 for( int i = 0; i < stabilitycenterlen; i++ ) 1230 { 1231 dualvector[i] = alpha * dualvector[i] + (1 - alpha) * stabilitycenter[i]; 1232 } 1233 for( int i = stabilitycenterlen; i < dualvectorlen; i++ ) 1234 { 1235 dualvector[i] = alpha * dualvector[i]; 1236 } 1237 1238 return SCIP_OKAY; 1239 } 1240 1241 /** stabilize Lagrangian multipliers */ 1242 static 1243 SCIP_RETCODE stabilizeDualVector( 1244 SCIP* scip, /**< SCIP data structure */ 1245 SCIP_SEPADATA* sepadata, /**< separator data structure */ 1246 SCIP_Real* dualvector, /**< Lagrangian multipliers */ 1247 int dualvectorlen, /**< length of the Lagrangian multipliers vector */ 1248 SCIP_Real* bestdualvector, /**< best Lagrangian multipliers found so far */ 1249 int bestdualvectorlen, /**< length of the best Lagrangian multipliers vector */ 1250 int nbestdualupdates, /**< number of best Lagrangian values found so far */ 1251 int subgradientiternum, /**< iteration number of the subgradient algorithm */ 1252 int totaliternum, /**< total number of iterations of the relax-and-cut algorithm performed so far */ 1253 SCIP_Real maxviolscore, /**< weighted average of maximum value of generated cut violations and maximum value of 1254 complementarity slackness violations, in the current iteration */ 1255 SCIP_Real maxviolscoreold, /**< weighted average of maximum value of generated cut violations and maximum value of 1256 complementarity slackness violations, in the previous iteration */ 1257 SCIP_Real nviolscore, /**< weighted average of number of generated cut violations and number of complementarity 1258 slackness violations, in the current iteration */ 1259 SCIP_Real nviolscoreold, /**< weighted average of number of generated cut violations and number of complementarity 1260 slackness violations, in the previous iteration */ 1261 int nlpiters, /**< number of LP iterations taken for solving the Lagrangian dual with fixed multipliers in 1262 current iteration */ 1263 SCIP_Real* ballradius /**< norm ball radius */ 1264 ) 1265 { 1266 if( sepadata->projectiontype > 0 ) 1267 { 1268 if( subgradientiternum >= 1 ) 1269 { 1270 /* update the ball radius */ 1271 updateBallRadius(scip, sepadata, maxviolscore, maxviolscoreold, nviolscore, nviolscoreold, nlpiters, 1272 ballradius); 1273 } 1274 1275 if( sepadata->projectiontype == 1 ) 1276 { 1277 /* projection of Lagrangian multipliers onto L1-norm ball */ 1278 SCIP_CALL( l1BallProjection(scip, dualvector, dualvectorlen, *ballradius) ); 1279 } 1280 else if( sepadata->projectiontype == 2 ) 1281 { 1282 /* projection of Lagrangian multipliers onto L2-norm ball */ 1283 l2BallProjection(scip, dualvector, dualvectorlen, *ballradius); 1284 } 1285 else if( sepadata->projectiontype == 3 ) 1286 { 1287 /* projection of Lagrangian multipliers onto L_inf-norm ball */ 1288 linfBallProjection(scip, dualvector, dualvectorlen, *ballradius); 1289 } 1290 } 1291 1292 if( sepadata->stabilitycentertype == 1 ) 1293 { 1294 /* weighted Lagrangian multipliers based on best Langrangian multipliers as stability center */ 1295 SCIP_CALL( weightedDualVector(scip, sepadata, dualvector, dualvectorlen, bestdualvector, 1296 bestdualvectorlen, nbestdualupdates, totaliternum) ); 1297 } 1298 1299 return SCIP_OKAY; 1300 } 1301 1302 /** update Lagrangian multipliers */ 1303 static 1304 SCIP_RETCODE updateDualVector( 1305 SCIP* scip, /**< SCIP data structure */ 1306 SCIP_SEPADATA* sepadata, /**< separator data structure */ 1307 SCIP_Real* dualvector1, /**< Lagrangian multipliers vector to be updated */ 1308 SCIP_Real* dualvector2, /**< Lagrangian multipliers vector used for backtracking */ 1309 int dualvector2len, /**< length of the Lagrangian multipliers vector used for backtracking */ 1310 int ndualvector2updates,/**< number of best Lagrangian values found so far */ 1311 int subgradientiternum, /**< iteration number of the subgradient algorithm */ 1312 int totaliternum, /**< total number of iterations of the relax-and-cut algorithm performed so far */ 1313 SCIP_Real steplength, /**< step length used for updating Lagrangian multipliers */ 1314 SCIP_Real* subgradient, /**< subgradient vector */ 1315 int ncuts, /**< number of generated cuts so far */ 1316 SCIP_Bool backtrack, /**< whether the Lagrangian multipliers need to be backtracked */ 1317 SCIP_Real maxviolscore, /**< weighted average of maximum value of generated cut violations and maximum value of 1318 complementarity slackness violations, in the current iteration */ 1319 SCIP_Real maxviolscoreold, /**< weighted average of maximum value of generated cut violations and maximum value of 1320 complementarity slackness violations, in the previous iteration */ 1321 SCIP_Real nviolscore, /**< weighted average of number of generated cut violations and number of complementarity 1322 slackness violations, in the current iteration */ 1323 SCIP_Real nviolscoreold, /**< weighted average of number of generated cut violations and number of complementarity 1324 slackness violations, in the previous iteration */ 1325 int nlpiters, /**< number of LP iterations taken for solving the Lagrangian dual with fixed multipliers in 1326 current iteration */ 1327 SCIP_Bool* dualvecsdiffer, /**< whether the updated Lagrangian multipliers differ from the old one */ 1328 SCIP_Real* ballradius /**< norm ball radius */ 1329 ) 1330 { 1331 SCIP_Real* dualvector1copy; 1332 1333 assert(dualvector2len <= ncuts); 1334 assert(dualvecsdiffer != NULL); 1335 1336 *dualvecsdiffer = FALSE; 1337 /* @todo do allocation and free operations outside only once instead of every time this function is called? */ 1338 /* copy of the Lagrangian multipliers to be used to check if the updated vector is different than this */ 1339 SCIP_CALL( SCIPallocCleanBufferArray(scip, &dualvector1copy, ncuts) ); 1340 for( int i = 0; i < ncuts; i++ ) 1341 { 1342 dualvector1copy[i] = dualvector1[i]; 1343 } 1344 1345 /* if backtracking was not identified at the time of the mu parameter update, then update the Lagrangian multipliers 1346 * based on the given subgradient vector 1347 */ 1348 if( !backtrack ) 1349 { 1350 assert((subgradient != NULL) || (ncuts == 0)); 1351 assert(subgradientiternum >= 0); 1352 1353 /* update Lagrangian multipliers */ 1354 for( int i = 0; i < ncuts; i++ ) 1355 { 1356 dualvector1[i] += steplength * subgradient[i]; 1357 } 1358 1359 /* projection onto non-negative orthant */ 1360 for( int i = 0; i < ncuts; i++ ) 1361 { 1362 dualvector1[i] = MAX(dualvector1[i], 0.0); 1363 } 1364 1365 /* stabilization of Lagrangian multipliers */ 1366 SCIP_CALL( stabilizeDualVector(scip, sepadata, dualvector1, ncuts, dualvector2, dualvector2len, 1367 ndualvector2updates, subgradientiternum, totaliternum, maxviolscore, maxviolscoreold, nviolscore, 1368 nviolscoreold, nlpiters, ballradius) ); 1369 1370 /* projection onto non-negative orthant again in case stabilization changed some components negative*/ 1371 for( int i = 0; i < ncuts; i++ ) 1372 { 1373 dualvector1[i] = MAX(dualvector1[i], 0.0); 1374 } 1375 } 1376 /* if backtracking was identified at the time of the mu parameter update, then backtrack the Lagrangian multipliers 1377 * based on the given backtracking multipliers 1378 */ 1379 else 1380 { 1381 for( int i = 0; i < dualvector2len; i++ ) 1382 { 1383 dualvector1[i] = dualvector2[i]; 1384 } 1385 1386 for( int i = dualvector2len; i < ncuts; i++ ) 1387 { 1388 dualvector1[i] = 0.0; 1389 } 1390 } 1391 1392 /* identify if the vector of Lagrangian multipliers is indeed different after updating */ 1393 for( int i = 0; i < ncuts; i++ ) 1394 { 1395 if( !SCIPisEQ(scip, dualvector1[i], dualvector1copy[i]) ) 1396 { 1397 *dualvecsdiffer = TRUE; 1398 break; 1399 } 1400 } 1401 1402 /* free memory */ 1403 for( int i = 0; i < ncuts; i++ ) 1404 { 1405 dualvector1copy[i] = 0.0; 1406 } 1407 SCIPfreeCleanBufferArray(scip, &dualvector1copy); 1408 1409 return SCIP_OKAY; 1410 } 1411 1412 /** check different termination criteria */ 1413 /* @note: the criterion based on objvecsdiffer assumes deterministic solving process (i.e., we would get same LP solution 1414 * for "Lagrangian dual with fixed Lagrangian multipliers" when the objective vector remains the same across iterations). 1415 */ 1416 /* @todo nlpssolved criterion? */ 1417 static 1418 SCIP_RETCODE checkLagrangianDualTermination( 1419 SCIP_SEPADATA* sepadata, /**< separator data structure */ 1420 int nnewaddedsoftcuts, /**< number of cuts that were recently penalized and added to the Lagrangian dual's 1421 objective function */ 1422 int nyettoaddsoftcuts, /**< number of cuts that are yet to be penalized and added to the Lagrangian dual's 1423 objective function */ 1424 SCIP_Bool objvecsdiffer, /**< whether the Lagrangian dual's objective function has changed */ 1425 int ngeneratedcurrroundcuts, /**< number of cuts generated in the current separation round */ 1426 int nmaxgeneratedperroundcuts,/**< maximal number of cuts allowed to generate per separation round */ 1427 int ncurrroundlpiters, /**< number of separating LP iterations in the current separation round */ 1428 int depth, /**< depth of the current node */ 1429 SCIP_Bool* terminate /**< whether to terminate the subgradient algorithm loop */ 1430 ) 1431 { 1432 *terminate = FALSE; 1433 1434 /* check if no new cuts were added to the Lagrangian dual, no cuts are remaining to be added, and the objective 1435 * function of the Lagrangian dual with fixed multipliers had not changed from the previous iteration 1436 */ 1437 if( (nnewaddedsoftcuts == 0) && (nyettoaddsoftcuts == 0) && !objvecsdiffer ) 1438 *terminate = TRUE; 1439 1440 /* check if allowed number of cuts in this separation round have already been generated */ 1441 if( ngeneratedcurrroundcuts >= nmaxgeneratedperroundcuts ) 1442 *terminate = TRUE; 1443 1444 /* check if allowed number of cuts in the tree have already been generated */ 1445 if( sepadata->ntotalcuts >= sepadata->nmaxtotalcuts ) 1446 *terminate = TRUE; 1447 1448 /* check if allowed number of simplex iterations in this separation round have already been used up */ 1449 if( (sepadata->nmaxperroundlpiters >= 0) && (ncurrroundlpiters >= sepadata->nmaxperroundlpiters) ) 1450 *terminate = TRUE; 1451 1452 /* check if allowed number of simplex iterations in the root node have already been used up */ 1453 if( (depth == 0) && (sepadata->nmaxrootlpiters >= 0) && (sepadata->nrootlpiters >= sepadata->nmaxrootlpiters) ) 1454 *terminate = TRUE; 1455 1456 /* check if allowed number of simplex iterations in the tree have already been used up */ 1457 if( (sepadata->nmaxtotallpiters >= 0) && (sepadata->ntotallpiters >= sepadata->nmaxtotallpiters) ) 1458 *terminate = TRUE; 1459 1460 return SCIP_OKAY; 1461 } 1462 1463 /** solve the LP corresponding to the Lagrangian dual with fixed Lagrangian multipliers */ 1464 static 1465 SCIP_RETCODE solveLagromoryLP( 1466 SCIP* scip, /**< SCIP data structure */ 1467 SCIP_SEPADATA* sepadata, /**< separator data structure */ 1468 int depth, /**< depth of the current node in the tree */ 1469 SCIP_Real origobjoffset, /**< objective offset in the current node's relaxation */ 1470 SCIP_Bool* solfound, /**< whether an LP optimal solution has been found */ 1471 SCIP_SOL* sol, /**< data structure to store LP optimal solution, if found */ 1472 SCIP_Real* solvals, /**< values of the LP optimal solution, if found */ 1473 SCIP_Real* objval, /**< optimal objective value of the LP optimal solution, if found */ 1474 int* ncurrroundlpiters /**< number of LP iterations taken for solving Lagrangian dual problems with fixed multipliers 1475 in the current separator round */ 1476 ) 1477 { 1478 SCIP_Real timelimit; 1479 SCIP_COL** cols; 1480 SCIP_COL* col; 1481 SCIP_VAR* var; 1482 SCIP_LPI* lpi; 1483 SCIP_Bool lperror; 1484 SCIP_Bool cutoff; 1485 SCIP_LPSOLSTAT stat; 1486 SCIP_Longint ntotallpiters; 1487 SCIP_Longint nlpiters; 1488 int ncols; 1489 int iterlimit; 1490 1491 assert(solfound != NULL); 1492 assert(sol != NULL); 1493 assert(solvals != NULL); 1494 assert(ncurrroundlpiters != NULL); 1495 assert(objval != NULL); 1496 1497 *solfound = FALSE; 1498 lperror = FALSE; 1499 cutoff = FALSE; 1500 iterlimit = -1; 1501 lpi = sepadata->lpiwithsoftcuts; 1502 1503 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) ); 1504 1505 /* set time limit */ 1506 SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) ); 1507 if( !SCIPisInfinity(scip, timelimit) ) 1508 { 1509 timelimit -= SCIPgetSolvingTime(scip); 1510 if( timelimit <= 0.0 ) 1511 { 1512 SCIPdebugMsg(scip, "skip Lagromory cut generation since no time left\n"); 1513 goto TERMINATE; 1514 } 1515 /* @note: the following direct LPI call is being used because of the lack of an equivalent function call in 1516 * scip_lp.c (lpSetRealpar exists in lp.c though) 1517 */ 1518 SCIP_CALL( SCIPlpiSetRealpar(lpi, SCIP_LPPAR_LPTILIM, timelimit) ); 1519 } 1520 1521 /* find iteration limit */ 1522 if( (depth == 0) && 1523 (sepadata->perrootlpiterfactor >= 0.0 && !SCIPisInfinity(scip, sepadata->perrootlpiterfactor)) ) 1524 { 1525 iterlimit = (int)(sepadata->perrootlpiterfactor * SCIPgetNRootFirstLPIterations(scip)); 1526 } 1527 else if( (depth > 0) && 1528 (sepadata->perlpiterfactor >= 0.0 && !SCIPisInfinity(scip, sepadata->perlpiterfactor)) ) 1529 { 1530 iterlimit = (int)(sepadata->perlpiterfactor * SCIPgetNNodeInitLPIterations(scip)); 1531 } 1532 if( sepadata->nmaxperroundlpiters >= 0 ) 1533 { 1534 if( sepadata->nmaxperroundlpiters - *ncurrroundlpiters >= 0 ) 1535 { 1536 if( iterlimit >= 0 ) 1537 { 1538 iterlimit = MIN(iterlimit, sepadata->nmaxperroundlpiters - *ncurrroundlpiters); 1539 } 1540 else 1541 { 1542 iterlimit = sepadata->nmaxperroundlpiters - *ncurrroundlpiters; 1543 } 1544 } 1545 else 1546 { 1547 iterlimit = 0; 1548 } 1549 } 1550 /* @todo impose a finite iteration limit only when the dualvector changes from zero to non-zero for the first time because 1551 * many simplex pivots are performed in this case even with warm starting (compared to the case when the 1552 * dualvector changes from non-zero to non-zero). 1553 */ 1554 1555 /* solve the LP with an iteration limit and get number of simplex iterations taken */ 1556 ntotallpiters = SCIPgetNLPIterations(scip); 1557 1558 SCIP_CALL( SCIPsolveDiveLP(scip, iterlimit, &lperror, &cutoff) ); 1559 1560 nlpiters = SCIPgetNLPIterations(scip) - ntotallpiters; 1561 1562 /* get the solution and objective value if optimal */ 1563 stat = SCIPgetLPSolstat(scip); 1564 /* @todo is there any way to accept terminations due to iterlimit and timelimit as well? It is not possible 1565 * currently because primal sol is not saved in these cases. 1566 */ 1567 /* @note: ideally, only primal feasibility is sufficient. But, there is no such option with SCIPgetLPSolstat. */ 1568 if( stat == SCIP_LPSOLSTAT_OPTIMAL ) 1569 { 1570 if( SCIPisLPSolBasic(scip) ) 1571 { 1572 *solfound = TRUE; 1573 1574 /* update sol */ 1575 for( int i = 0; i < ncols; ++i ) 1576 { 1577 col = cols[i]; 1578 assert(col != NULL); 1579 1580 var = SCIPcolGetVar(col); 1581 1582 solvals[i] = SCIPcolGetPrimsol(col); 1583 SCIP_CALL( SCIPsetSolVal(scip, sol, var, solvals[i]) ); 1584 } 1585 1586 *objval = SCIPgetLPObjval(scip); 1587 *objval = *objval + origobjoffset; 1588 } 1589 } 1590 1591 /* update some statistics */ 1592 if( depth == 0 ) 1593 { 1594 sepadata->nrootlpiters += (int)nlpiters; 1595 } 1596 sepadata->ntotallpiters += (int)nlpiters; 1597 *ncurrroundlpiters += (int)nlpiters; 1598 1599 TERMINATE: 1600 return SCIP_OKAY; 1601 } 1602 1603 /** solve the LP corresponding to the node relaxation upon adding all the generated cuts */ 1604 static 1605 SCIP_RETCODE solveLPWithHardCuts( 1606 SCIP* scip, /**< SCIP data structure */ 1607 SCIP_SEPADATA* sepadata, /**< separator data structure */ 1608 SCIP_Bool* solfound, /**< whether an LP optimal solution has been found */ 1609 SCIP_SOL* sol, /**< data structure to store LP optimal solution, if found */ 1610 SCIP_Real* solvals /**< values of the LP optimal solution, if found */ 1611 ) 1612 { 1613 SCIP_Real timelimit; 1614 SCIP_COL** cols; 1615 SCIP_COL* col; 1616 SCIP_VAR* var; 1617 int ncols; 1618 1619 assert(solfound != NULL); 1620 assert(sol != NULL); 1621 assert(solvals != NULL); 1622 1623 *solfound = FALSE; 1624 1625 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) ); 1626 1627 /* set time limit */ 1628 SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) ); 1629 if( !SCIPisInfinity(scip, timelimit) ) 1630 { 1631 timelimit -= SCIPgetSolvingTime(scip); 1632 if( timelimit <= 0.0 ) 1633 { 1634 SCIPdebugMsg(scip, "skip Lagromory cut generation since no time left\n"); 1635 goto TERMINATE; 1636 } 1637 SCIP_CALL( SCIPlpiSetRealpar(sepadata->lpiwithhardcuts, SCIP_LPPAR_LPTILIM, timelimit) ); 1638 } 1639 1640 /* solve the LP */ 1641 SCIPlpiSolvePrimal(sepadata->lpiwithhardcuts); /*lint !e534*/ 1642 1643 /* get the solution if primal feasible */ 1644 if( SCIPlpiIsPrimalFeasible(sepadata->lpiwithhardcuts) ) 1645 { 1646 *solfound = TRUE; 1647 SCIP_CALL( SCIPlpiGetSol(sepadata->lpiwithhardcuts, NULL, solvals, NULL, NULL, NULL) ); 1648 1649 /* update sol */ 1650 for( int i = 0; i < ncols; ++i ) 1651 { 1652 col = cols[i]; 1653 assert(col != NULL); 1654 1655 var = SCIPcolGetVar(col); 1656 1657 SCIP_CALL( SCIPsetSolVal(scip, sol, var, solvals[i]) ); 1658 } 1659 } 1660 1661 TERMINATE: 1662 return SCIP_OKAY; 1663 } 1664 1665 /** construct a cut based on the input cut coefficients, sides, etc */ 1666 static 1667 SCIP_RETCODE constructCutRow( 1668 SCIP* scip, /**< SCIP data structure */ 1669 SCIP_SEPA* sepa, /**< pointer to the separator */ 1670 SCIP_SEPADATA* sepadata, /**< separator data structure */ 1671 int mainiternum, /**< iteration number of the outer loop of the relax-and-cut algorithm */ 1672 int subgradientiternum, /**< iteration number of the subgradient algorithm */ 1673 int cutnnz, /**< number of nonzeros in cut */ 1674 int* cutinds, /**< column indices in cut */ 1675 SCIP_Real* cutcoefs, /**< cut cofficients */ 1676 SCIP_Real cutefficacy, /**< cut efficacy */ 1677 SCIP_Real cutrhs, /**< RHS of cut */ 1678 SCIP_Bool cutislocal, /**< whether cut is local */ 1679 int cutrank, /**< rank of cut */ 1680 SCIP_ROW** generatedcuts, /**< array of generated cuts */ 1681 SCIP_Real* generatedcutefficacies, /**< array of generated cut efficacies w.r.t. the respective LP bases used for cut 1682 generations */ 1683 int ngeneratedcurrroundcuts, /**< number of cuts generated until the previous basis in the current separation round */ 1684 int* ngeneratednewcuts, /**< number of new cuts generated using the current basis */ 1685 SCIP_Bool* cutoff /**< should the current node be cutoff? */ 1686 ) 1687 { 1688 SCIP_COL** cols; 1689 SCIP_Real minactivity; 1690 SCIP_Real maxactivity; 1691 SCIP_Real lhs; 1692 SCIP_Real rhs; 1693 1694 assert(scip != NULL); 1695 assert(mainiternum >= 0); 1696 assert(ngeneratednewcuts != NULL); 1697 assert(cutoff != NULL); 1698 1699 *cutoff = FALSE; 1700 1701 if( cutnnz == 0 && SCIPisFeasNegative(scip, cutrhs) ) /*lint !e644*/ 1702 { 1703 SCIPdebugMsg(scip, " -> Lagromory cut detected node infeasibility with cut 0 <= %g.\n", cutrhs); 1704 *cutoff = TRUE; 1705 return SCIP_OKAY; 1706 } 1707 1708 /* only take efficient cuts */ 1709 if( SCIPisEfficacious(scip, cutefficacy) ) 1710 { 1711 SCIP_ROW* cut; 1712 SCIP_VAR* var; 1713 char cutname[SCIP_MAXSTRLEN]; 1714 int v; 1715 1716 /* construct cut name */ 1717 if( subgradientiternum >= 0 ) 1718 { 1719 (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "%s_%" SCIP_LONGINT_FORMAT "_%d" "_%d" "_%d", SCIPsepaGetName(sepa), 1720 sepadata->ncalls, mainiternum, subgradientiternum, ngeneratedcurrroundcuts + *ngeneratednewcuts); 1721 } 1722 else 1723 { 1724 (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "%s_%" SCIP_LONGINT_FORMAT "_%d" "_%d", SCIPsepaGetName(sepa), 1725 sepadata->ncalls, mainiternum, ngeneratedcurrroundcuts + *ngeneratednewcuts); 1726 } 1727 1728 /* create empty cut */ 1729 SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &cut, sepa, cutname, -SCIPinfinity(scip), cutrhs, cutislocal, FALSE, 1730 sepadata->dynamiccuts) ); 1731 1732 /* set cut rank */ 1733 SCIProwChgRank(cut, cutrank); /*lint !e644*/ 1734 1735 /* cache the row extension and only flush them if the cut gets added */ 1736 SCIP_CALL( SCIPcacheRowExtensions(scip, cut) ); 1737 1738 cols = SCIPgetLPCols(scip); 1739 1740 /* collect all non-zero coefficients */ 1741 for( v = 0; v < cutnnz; ++v ) 1742 { 1743 var = SCIPcolGetVar(cols[cutinds[v]]); 1744 SCIP_CALL( SCIPaddVarToRow(scip, cut, var, cutcoefs[v]) ); 1745 } 1746 1747 /* flush all changes before adding the cut */ 1748 SCIP_CALL( SCIPflushRowExtensions(scip, cut) ); 1749 1750 if( SCIProwGetNNonz(cut) == 0 ) 1751 { 1752 assert( SCIPisFeasNegative(scip, cutrhs) ); 1753 SCIPdebugMsg(scip, " -> Lagromory cut detected node infeasibility with cut 0 <= %g.\n", cutrhs); 1754 *cutoff = TRUE; 1755 return SCIP_OKAY; 1756 } 1757 else 1758 { 1759 /* gathering lhs and rhs again in case the separator is extended later to add other cuts/constraints that may 1760 * have non-inf lhs or inf rhs */ 1761 lhs = SCIProwGetLhs(cut); 1762 rhs = SCIProwGetRhs(cut); 1763 assert(SCIPisInfinity(scip, -lhs)); 1764 assert(!SCIPisInfinity(scip, rhs)); 1765 1766 SCIPdebugMsg(scip, " -> %s cut <%s>: rhs=%f, eff=%f\n", "lagromory", cutname, cutrhs, cutefficacy); 1767 1768 /* check if the cut leads to infeasibility (i.e., *cutoff = TRUE) */ 1769 { 1770 /* modifiable cuts cannot be declared infeasible, since we don't know all coefficients */ 1771 if( SCIProwIsModifiable(cut) ) 1772 *cutoff = FALSE; 1773 1774 /* check for activity infeasibility */ 1775 minactivity = SCIPgetRowMinActivity(scip, cut); 1776 maxactivity = SCIPgetRowMaxActivity(scip, cut); 1777 1778 if( (!SCIPisInfinity(scip, rhs) && SCIPisFeasPositive(scip, minactivity - rhs)) || 1779 (!SCIPisInfinity(scip, -lhs) && SCIPisFeasNegative(scip, maxactivity - lhs)) ) 1780 { 1781 SCIPdebugMsg(scip, "cut <%s> is infeasible (sides=[%g,%g], act=[%g,%g])\n", 1782 SCIProwGetName(cut), lhs, rhs, minactivity, maxactivity); 1783 *cutoff = TRUE; 1784 } 1785 } 1786 1787 /* store the newly generated cut in an array and update some statistics */ 1788 generatedcuts[ngeneratedcurrroundcuts + *ngeneratednewcuts] = cut; 1789 generatedcutefficacies[ngeneratedcurrroundcuts + *ngeneratednewcuts] = cutefficacy; 1790 ++(*ngeneratednewcuts); 1791 } 1792 } 1793 1794 return SCIP_OKAY; 1795 } 1796 1797 /** aggregated generated cuts based on the best Lagrangian multipliers */ 1798 static 1799 SCIP_RETCODE aggregateGeneratedCuts( 1800 SCIP* scip, /**< SCIP data structure */ 1801 SCIP_SEPA* sepa, /**< pointer to the separator */ 1802 SCIP_SEPADATA* sepadata, /**< separator data structure */ 1803 SCIP_ROW** generatedcuts, /**< cuts generated in the current separation round */ 1804 SCIP_Real* bestdualvector, /**< best Lagrangian multipliers vector */ 1805 int bestdualvectorlen, /**< length of the best Lagrangian multipliers vector */ 1806 SCIP_ROW** aggrcuts, /**< aggregated cuts generated so far in the current separation round */ 1807 int* naggrcuts, /**< number of aggregated cuts generated so far in the current separation round */ 1808 SCIP_Bool* cutoff /**< should the current node be cutoff? */ 1809 ) 1810 { 1811 SCIP_Real* cutvals; /**< cut cofficients */ 1812 SCIP_COL** cutcols; 1813 SCIP_COL** cols; 1814 SCIP_VAR* var; 1815 SCIP_ROW* cut; 1816 SCIP_ROW* aggrcut; 1817 SCIP_Real minactivity; 1818 SCIP_Real maxactivity; 1819 SCIP_Real cutlhs; 1820 SCIP_Real cutrhs; 1821 SCIP_Real cutconst; 1822 SCIP_Real aggrcutlhs; 1823 SCIP_Real aggrcutrhs; 1824 SCIP_Real aggrcutconst; 1825 SCIP_Real* aggrcutvals; 1826 SCIP_Real* aggrcutcoefs; 1827 SCIP_Real multiplier; 1828 SCIP_Bool aggrcutislocal; 1829 SCIP_Bool aggrindicator; 1830 SCIP_Real* tmpcutvals; 1831 SCIP_Real QUAD(quadterm); 1832 SCIP_Real QUAD(tmpcutconst); 1833 SCIP_Real QUAD(tmpcutrhs); 1834 SCIP_Real QUAD(quadprod); 1835 char aggrcutname[SCIP_MAXSTRLEN]; 1836 int cutnnz; /**< number of nonzeros in cut */ 1837 int aggrcutnnz; 1838 int* aggrcutinds; /**< column indices in cut */ 1839 int aggrcutrank; /**< rank of cut */ 1840 int cutrank; 1841 int ncols; 1842 int collppos; 1843 int nlocalcuts; 1844 1845 assert(scip != NULL); 1846 assert(naggrcuts != NULL); 1847 assert(cutoff != NULL); 1848 1849 *cutoff = FALSE; 1850 aggrcutlhs = -SCIPinfinity(scip); 1851 aggrcutnnz = 0; 1852 aggrcutrank = -1; 1853 nlocalcuts = 0; 1854 aggrindicator = FALSE; 1855 1856 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) ); 1857 1858 /* allocate memory */ 1859 SCIP_CALL( SCIPallocCleanBufferArray(scip, &tmpcutvals, QUAD_ARRAY_SIZE(ncols)) ); 1860 QUAD_ASSIGN(tmpcutconst, 0.0); 1861 QUAD_ASSIGN(tmpcutrhs, 0.0); 1862 SCIP_CALL( SCIPallocCleanBufferArray(scip, &aggrcutvals, ncols) ); 1863 SCIP_CALL( SCIPallocBufferArray(scip, &aggrcutcoefs, ncols) ); 1864 SCIP_CALL( SCIPallocBufferArray(scip, &aggrcutinds, ncols) ); 1865 1866 /* aggregate cuts based on the input Lagrangian multipliers */ 1867 for( int i = 0; i < bestdualvectorlen; i++ ) 1868 { 1869 multiplier = bestdualvector[i]; 1870 if( SCIPisGE(scip, multiplier, 1e-4) ) 1871 { 1872 cut = generatedcuts[i]; 1873 cutnnz = SCIProwGetNNonz(cut); 1874 cutlhs = SCIProwGetLhs(cut); 1875 cutrhs = SCIProwGetRhs(cut); 1876 assert(SCIPisInfinity(scip, -cutlhs)); 1877 assert(!SCIPisInfinity(scip, cutrhs)); 1878 cutvals = SCIProwGetVals(cut); 1879 cutcols = SCIProwGetCols(cut); 1880 cutconst = SCIProwGetConstant(cut); 1881 1882 for( int j = 0; j < cutnnz; j++ ) 1883 { 1884 collppos = SCIPcolGetLPPos(cutcols[j]); 1885 assert(collppos >= 0); 1886 assert(collppos <= ncols); 1887 1888 QUAD_ARRAY_LOAD(quadterm, tmpcutvals, collppos); 1889 SCIPquadprecProdDD(quadprod, multiplier, cutvals[j]); 1890 SCIPquadprecSumQQ(quadterm, quadterm, quadprod); 1891 QUAD_ARRAY_STORE(tmpcutvals, collppos, quadterm); 1892 } 1893 1894 SCIPquadprecProdDD(quadprod, multiplier, cutconst); 1895 SCIPquadprecSumQQ(tmpcutconst, tmpcutconst, quadprod); 1896 SCIPquadprecProdDD(quadprod, multiplier, cutrhs); 1897 SCIPquadprecSumQQ(tmpcutrhs, tmpcutrhs, quadprod); 1898 1899 cutrank = SCIProwGetRank(cut); 1900 aggrcutrank = MAX(aggrcutrank, cutrank); 1901 nlocalcuts += (SCIProwIsLocal(cut) ? 1 : 0); 1902 aggrindicator = TRUE; 1903 } 1904 } 1905 /* if at least one cut is local, then the aggregated cut is local */ 1906 aggrcutislocal = (nlocalcuts > 0 ? TRUE : FALSE); 1907 1908 /* if the aggregation was successful, then create an empty row and build a cut */ 1909 if( aggrindicator ) 1910 { 1911 aggrcutconst = QUAD_TO_DBL(tmpcutconst); 1912 aggrcutrhs = QUAD_TO_DBL(tmpcutrhs); 1913 1914 /* build sparse representation of the aggregated cut */ 1915 for( int i = 0; i < ncols; i++ ) 1916 { 1917 QUAD_ARRAY_LOAD(quadterm, tmpcutvals, i); 1918 aggrcutvals[i] = QUAD_TO_DBL(quadterm); 1919 if( !SCIPisZero(scip, aggrcutvals[i]) ) 1920 { 1921 aggrcutcoefs[aggrcutnnz] = aggrcutvals[i]; 1922 aggrcutinds[aggrcutnnz] = i; 1923 aggrcutnnz++; 1924 } 1925 } 1926 1927 if( aggrcutnnz == 0 && SCIPisFeasNegative(scip, aggrcutrhs) ) /*lint !e644*/ 1928 { 1929 SCIPdebugMsg(scip, " -> Lagromory cut detected node infeasibility with cut 0 <= %g.\n", aggrcutrhs); 1930 *cutoff = TRUE; 1931 goto TERMINATE; 1932 } 1933 1934 /* construct aggregated cut name */ 1935 (void) SCIPsnprintf(aggrcutname, SCIP_MAXSTRLEN, "%s_%" SCIP_LONGINT_FORMAT "_aggregated", SCIPsepaGetName(sepa), 1936 sepadata->ncalls); 1937 1938 /* create empty cut */ 1939 SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &aggrcut, sepa, aggrcutname, aggrcutlhs - aggrcutconst, aggrcutrhs - 1940 aggrcutconst, aggrcutislocal, FALSE, sepadata->dynamiccuts) ); 1941 1942 /* set cut rank */ 1943 SCIProwChgRank(aggrcut, aggrcutrank); /*lint !e644*/ 1944 1945 /* cache the row extension and only flush them if the cut gets added */ 1946 SCIP_CALL( SCIPcacheRowExtensions(scip, aggrcut) ); 1947 1948 /* collect all non-zero coefficients */ 1949 for( int i = 0; i < aggrcutnnz; i++ ) 1950 { 1951 var = SCIPcolGetVar(cols[aggrcutinds[i]]); 1952 SCIP_CALL( SCIPaddVarToRow(scip, aggrcut, var, aggrcutcoefs[i]) ); 1953 } 1954 1955 /* flush all changes before adding the cut */ 1956 SCIP_CALL( SCIPflushRowExtensions(scip, aggrcut) ); 1957 1958 if( SCIProwGetNNonz(aggrcut) == 0 ) 1959 { 1960 assert( SCIPisFeasNegative(scip, aggrcutrhs) ); 1961 SCIPdebugMsg(scip, " -> Lagromory cut detected node infeasibility with cut 0 <= %g.\n", aggrcutrhs); 1962 *cutoff = TRUE; 1963 goto TERMINATE; 1964 } 1965 else 1966 { 1967 /* gathering lhs and rhs again in case the separator is extended later to add other cuts/constraints that may 1968 * have non-inf lhs or inf rhs */ 1969 cutlhs = SCIProwGetLhs(aggrcut); 1970 cutrhs = SCIProwGetRhs(aggrcut); 1971 assert(SCIPisInfinity(scip, -cutlhs)); 1972 assert(!SCIPisInfinity(scip, cutrhs)); 1973 1974 SCIPdebugMsg(scip, " -> %s cut <%s>: rhs=%f\n", "lagromory", aggrcutname, aggrcutrhs); 1975 1976 /* check if the cut leads to infeasibility (i.e., *cutoff = TRUE) */ 1977 { 1978 /* modifiable cuts cannot be declared infeasible, since we don't know all coefficients */ 1979 if( SCIProwIsModifiable(aggrcut) ) 1980 *cutoff = FALSE; 1981 1982 /* check for activity infeasibility */ 1983 minactivity = SCIPgetRowMinActivity(scip, aggrcut); 1984 maxactivity = SCIPgetRowMaxActivity(scip, aggrcut); 1985 1986 if( (!SCIPisInfinity(scip, cutrhs) && SCIPisFeasPositive(scip, minactivity - cutrhs)) || 1987 (!SCIPisInfinity(scip, -cutlhs) && SCIPisFeasNegative(scip, maxactivity - cutlhs)) ) 1988 { 1989 SCIPdebugMsg(scip, "cut <%s> is infeasible (sides=[%g,%g], act=[%g,%g])\n", 1990 SCIProwGetName(aggrcut), cutlhs, cutrhs, minactivity, maxactivity); 1991 *cutoff = TRUE; 1992 } 1993 } 1994 1995 /* add the aggregated cut to a separate data structure */ 1996 aggrcuts[*naggrcuts] = aggrcut; 1997 (*naggrcuts)++; 1998 } 1999 2000 QUAD_ASSIGN(quadterm, 0.0); 2001 for( int i = 0; i < ncols; i++ ) 2002 { 2003 aggrcutvals[i] = 0.0; 2004 QUAD_ARRAY_STORE(tmpcutvals, i, quadterm); 2005 } 2006 } 2007 2008 TERMINATE: 2009 /* free memory */ 2010 SCIPfreeBufferArray(scip, &aggrcutinds); 2011 SCIPfreeBufferArray(scip, &aggrcutcoefs); 2012 SCIPfreeCleanBufferArray(scip, &aggrcutvals); 2013 SCIPfreeCleanBufferArray(scip, &tmpcutvals); 2014 2015 return SCIP_OKAY; 2016 } 2017 2018 /** main method: LP solution separation method of separator */ 2019 static 2020 SCIP_RETCODE generateGMICuts( 2021 SCIP* scip, /**< SCIP data structure */ 2022 SCIP_SEPA* sepa, /**< pointer to the separator */ 2023 SCIP_SEPADATA* sepadata, /**< separator data structure */ 2024 int mainiternum, /**< iteration number of the outer loop of the relax-and-cut algorithm */ 2025 int subgradientiternum, /**< iteration number of the subgradient algorithm */ 2026 SCIP_SOL* sol, /**< LP solution to be used for cut generation */ 2027 SCIP_Real* solvals, /**< values of the LP solution to be used for cut generation */ 2028 int nmaxgeneratedperroundcuts,/**< maximal number of cuts allowed to generate per separation round */ 2029 SCIP_Bool allowlocal, /**< should locally valid cuts be generated? */ 2030 SCIP_ROW** generatedcurrroundcuts, /**< cuts generated in the current separation round */ 2031 SCIP_Real* generatedcutefficacies, /**< array of generated cut efficacies w.r.t. the respective LP bases used for cut 2032 generations */ 2033 int ngeneratedcurrroundcuts, /**< number of cuts generated until the previous basis in the current separation round */ 2034 int* ngeneratednewcuts, /**< number of new cuts generated using the current basis */ 2035 int depth, /**< depth of the current node in the tree */ 2036 SCIP_Bool* cutoff /**< should the current node be cutoff? */ 2037 ) 2038 { 2039 SCIP_Real minfrac; 2040 SCIP_Real maxfrac; 2041 SCIP_Real frac; 2042 SCIP_Real* basisfrac; 2043 SCIP_Real* binvrow; 2044 SCIP_AGGRROW* aggrrow; 2045 SCIP_Bool success; 2046 SCIP_Real* cutcoefs; 2047 SCIP_Real cutrhs; 2048 SCIP_Real cutefficacy; 2049 SCIP_Bool cutislocal; 2050 SCIP_ROW** rows; 2051 SCIP_COL** cols; 2052 SCIP_VAR* var; 2053 SCIP_ROW* row; 2054 int cutrank; 2055 int cutnnz; 2056 int* cutinds; 2057 int* basisind; 2058 int* inds; 2059 int nrows; 2060 int ncols; 2061 int* basisperm; 2062 int k; 2063 int c; 2064 int ninds; 2065 int nmaxcutsperlp; 2066 2067 assert(ngeneratednewcuts != NULL); 2068 2069 minfrac = sepadata->away; 2070 maxfrac = 1.0 - sepadata->away; 2071 SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) ); 2072 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) ); 2073 *ngeneratednewcuts = 0; 2074 nmaxcutsperlp = ((depth == 0) ? sepadata->nmaxcutsperlproot : sepadata->nmaxcutsperlp); 2075 2076 /* allocate memory */ 2077 SCIP_CALL( SCIPallocBufferArray(scip, &basisperm, nrows) ); 2078 SCIP_CALL( SCIPallocBufferArray(scip, &basisfrac, nrows) ); 2079 SCIP_CALL( SCIPallocBufferArray(scip, &binvrow, nrows) ); 2080 SCIP_CALL( SCIPallocBufferArray(scip, &inds, nrows) ); 2081 SCIP_CALL( SCIPallocBufferArray(scip, &basisind, nrows) ); 2082 SCIP_CALL( SCIPallocBufferArray(scip, &cutcoefs, ncols) ); 2083 SCIP_CALL( SCIPallocBufferArray(scip, &cutinds, ncols) ); 2084 SCIP_CALL( SCIPaggrRowCreate(scip, &aggrrow) ); 2085 2086 SCIP_CALL( SCIPgetLPBasisInd(scip, basisind) ); 2087 2088 /* check if the rows of the simplex tableau are suitable for cut generation and build an array of fractions */ 2089 for( int i = 0; i < nrows; ++i ) 2090 { 2091 frac = 0.0; 2092 2093 c = basisind[i]; 2094 2095 basisperm[i] = i; 2096 2097 /* if the simplex tableau row corresponds to an LP column */ 2098 if( c >= 0 ) 2099 { 2100 assert(c < ncols); 2101 2102 var = SCIPcolGetVar(cols[c]); 2103 /* if the column is non-continuous one */ 2104 if( SCIPvarGetType(var) != SCIP_VARTYPE_CONTINUOUS ) 2105 { 2106 frac = SCIPfeasFrac(scip, solvals[c]); 2107 frac = MIN(frac, 1.0 - frac); 2108 } 2109 } 2110 /* if the simplex tableau row corresponds to an LP row and separation on rows is allowed */ 2111 else if( sepadata->separaterows ) 2112 { 2113 assert(0 <= -c-1); 2114 assert(-c-1 < nrows); 2115 2116 row = rows[-c-1]; 2117 /* if the row is suitable for cut generation */ 2118 if( SCIProwIsIntegral(row) && !SCIProwIsModifiable(row) ) 2119 { 2120 frac = SCIPfeasFrac(scip, SCIPgetRowActivity(scip, row)); 2121 frac = MIN(frac, 1.0 - frac); 2122 } 2123 } 2124 2125 if( frac >= minfrac ) 2126 { 2127 /* slightly change fractionality to have random order for equal fractions */ 2128 basisfrac[i] = frac + SCIPrandomGetReal(sepadata->randnumgen, -1e-6, 1e-6); 2129 } 2130 else 2131 { 2132 basisfrac[i] = 0.0; 2133 } 2134 } 2135 2136 /* if there is a need to sort the fractionalities */ 2137 if( sepadata->sortcutoffsol ) 2138 { 2139 /* sort basis indices by fractionality */ 2140 SCIPsortDownRealInt(basisfrac, basisperm, nrows); 2141 } 2142 2143 /* for all basic columns belonging to integer variables, try to generate a GMI cut */ 2144 for( int i = 0; i < nrows && !SCIPisStopped(scip) && !*cutoff; ++i ) 2145 { 2146 if( (ngeneratedcurrroundcuts + *ngeneratednewcuts >= nmaxgeneratedperroundcuts) || 2147 (sepadata->ntotalcuts + *ngeneratednewcuts >= sepadata->nmaxtotalcuts) || 2148 (*ngeneratednewcuts >= nmaxcutsperlp) ) 2149 break; 2150 2151 ninds = -1; 2152 cutefficacy = 0.0; 2153 2154 /* either break the loop or proceed to the next iteration if the fractionality is zero */ 2155 if( basisfrac[i] == 0.0 ) 2156 { 2157 if( sepadata->sortcutoffsol ) 2158 break; 2159 else 2160 continue; 2161 } 2162 2163 k = basisperm[i]; 2164 2165 /* get the row of B^-1 for this basic integer variable with fractional solution value and call aggregate function */ 2166 SCIP_CALL( SCIPgetLPBInvRow(scip, k, binvrow, inds, &ninds) ); 2167 2168 SCIP_CALL( SCIPaggrRowSumRows(scip, aggrrow, binvrow, inds, ninds, sepadata->sidetypebasis, allowlocal, 2, 2169 (int) MAXAGGRLEN(ncols), &success) ); 2170 2171 if( !success ) 2172 continue; 2173 2174 /* @todo Currently we are using the SCIPcalcMIR() function to compute the coefficients of the Gomory 2175 * cut. Alternatively, we could use the direct version (see thesis of Achterberg formula (8.4)) which 2176 * leads to cut a of the form \sum a_i x_i \geq 1. Rumor has it that these cuts are better. 2177 */ 2178 2179 /* try to create GMI cut out of the aggregation row */ 2180 SCIP_CALL( SCIPcalcMIR(scip, sol, POSTPROCESS, BOUNDSWITCH, USEVBDS, allowlocal, FIXINTEGRALRHS, NULL, 2181 NULL, minfrac, maxfrac, 1.0, aggrrow, cutcoefs, &cutrhs, cutinds, &cutnnz, &cutefficacy, &cutrank, 2182 &cutislocal, &success) ); 2183 2184 if( success ) 2185 { 2186 assert(allowlocal || !cutislocal); /*lint !e644*/ 2187 2188 SCIP_CALL( constructCutRow(scip, sepa, sepadata, mainiternum, subgradientiternum, cutnnz, cutinds, cutcoefs, 2189 cutefficacy, cutrhs, cutislocal, cutrank, generatedcurrroundcuts, generatedcutefficacies, 2190 ngeneratedcurrroundcuts, ngeneratednewcuts, cutoff)); 2191 } 2192 } 2193 2194 /* free temporary memory */ 2195 SCIPfreeBufferArray(scip, &cutinds); 2196 SCIPfreeBufferArray(scip, &cutcoefs); 2197 SCIPfreeBufferArray(scip, &basisind); 2198 SCIPfreeBufferArray(scip, &inds); 2199 SCIPfreeBufferArray(scip, &binvrow); 2200 SCIPfreeBufferArray(scip, &basisfrac); 2201 SCIPfreeBufferArray(scip, &basisperm); 2202 SCIPaggrRowFree(scip, &aggrrow); 2203 2204 return SCIP_OKAY; 2205 } 2206 2207 /** update objective vector w.r.t. the fixed Lagrangian multipliers */ 2208 static 2209 SCIP_RETCODE updateObjectiveVector( 2210 SCIP* scip, /**< SCIP data structure */ 2211 SCIP_Real* dualvector, /**< Lagrangian multipliers vector */ 2212 SCIP_ROW** cuts, /**< cuts generated so far in the current separation round */ 2213 int ncuts, /**< number of cuts generated so far in the current separation round */ 2214 SCIP_Real* origobjcoefs, /**< original objective function coefficients of the node linear relaxation */ 2215 SCIP_Bool* objvecsdiffer /**< whether the updated objective function coefficients differ from the old ones */ 2216 ) 2217 { 2218 SCIP_Real* objvals; 2219 SCIP_Real* prod; 2220 SCIP_Real* oldobjvals; 2221 SCIP_Real* cutvals; 2222 SCIP_COL** cutcols; 2223 SCIP_COL** cols; 2224 SCIP_VAR* var; 2225 int cutnnz; 2226 int collppos; 2227 int ncols; 2228 2229 assert(objvecsdiffer != NULL); 2230 assert(ncuts > 0); 2231 2232 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) ); 2233 2234 /* allocate memory */ 2235 SCIP_CALL( SCIPallocBufferArray(scip, &objvals, ncols) ); 2236 SCIP_CALL( SCIPallocBufferArray(scip, &oldobjvals, ncols) ); 2237 SCIP_CALL( SCIPallocCleanBufferArray(scip, &prod, ncols) ); 2238 *objvecsdiffer = FALSE; 2239 2240 /* find the product of Lagrangian multipliers and cut coefficients */ 2241 for( int i = 0; i < ncuts; i++ ) 2242 { 2243 if( !SCIPisZero(scip, dualvector[i]) ) 2244 { 2245 assert(!(SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i])) && SCIPisInfinity(scip, SCIProwGetRhs(cuts[i])))); 2246 assert(SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i]))); 2247 assert(!SCIPisInfinity(scip, SCIProwGetRhs(cuts[i]))); 2248 2249 cutnnz = SCIProwGetNNonz(cuts[i]); 2250 assert(cutnnz <= ncols); 2251 cutvals = SCIProwGetVals(cuts[i]); 2252 cutcols = SCIProwGetCols(cuts[i]); 2253 2254 for( int j = 0; j < cutnnz; ++j ) 2255 { 2256 collppos = SCIPcolGetLPPos(cutcols[j]); 2257 assert(collppos >= 0); 2258 assert(collppos <= ncols); 2259 2260 prod[collppos] += dualvector[i] * cutvals[j]; 2261 } 2262 } 2263 } 2264 2265 /* change objective coefficients */ 2266 for( int i = 0; i < ncols; i++ ) 2267 { 2268 var = SCIPcolGetVar(cols[i]); 2269 oldobjvals[i] = SCIPgetVarObjDive(scip, var); 2270 objvals[i] = origobjcoefs[i] + prod[i]; 2271 2272 SCIP_CALL( SCIPchgVarObjDive(scip, var, objvals[i]) ); 2273 2274 /* identify if the updated objective vector is indeed different from the previous one */ 2275 if( !(*objvecsdiffer) && !SCIPisEQ(scip, oldobjvals[i], objvals[i]) ) 2276 *objvecsdiffer = TRUE; 2277 } 2278 2279 for( int i = 0; i < ncols; i++) 2280 { 2281 prod[i] = 0.0; 2282 } 2283 2284 /* free memory */ 2285 SCIPfreeCleanBufferArray(scip, &prod); 2286 SCIPfreeBufferArray(scip, &oldobjvals); 2287 SCIPfreeBufferArray(scip, &objvals); 2288 2289 return SCIP_OKAY; 2290 } 2291 2292 /** add GMI cuts to the objective function of the Lagrangian dual problem by introducing new Lagrangian multipliers */ 2293 static 2294 SCIP_RETCODE addGMICutsAsSoftConss( 2295 SCIP_Real* dualvector, /**< Lagrangian multipliers vector */ 2296 int ngeneratedcuts, /**< number of cuts generated so far in the current separation round */ 2297 int* naddedcuts, /**< number of cuts added so far in the current separation round to the Lagrangian dual problem 2298 upon penalization */ 2299 int* nnewaddedsoftcuts /**< number of cuts added newly to the Lagrangian dual problem upon penalization */ 2300 ) 2301 { 2302 assert(*naddedcuts <= ngeneratedcuts); 2303 2304 /* set the initial penalty of the newly penalized cuts as zero */ 2305 for( int i = *naddedcuts; i < ngeneratedcuts; i++ ) 2306 dualvector[i] = 0.0; 2307 2308 *nnewaddedsoftcuts = ngeneratedcuts - *naddedcuts; 2309 *naddedcuts = ngeneratedcuts; 2310 2311 return SCIP_OKAY; 2312 } 2313 2314 /** solve the Lagrangian dual problem */ 2315 static 2316 SCIP_RETCODE solveLagrangianDual( 2317 SCIP* scip, /**< SCIP data structure */ 2318 SCIP_SEPA* sepa, /**< pointer to the separator */ 2319 SCIP_SEPADATA* sepadata, /**< separator data structure */ 2320 SCIP_SOL* sol, /**< data structure to store an LP solution upon solving a Lagrangian dual problem with 2321 fixed Lagrangian multipliers */ 2322 SCIP_Real* solvals, /**< values of the LP solution obtained upon solving a Lagrangian dual problem with fixed 2323 Lagrangian multipliers */ 2324 int mainiternum, /**< iteration number of the outer loop of the relax-and-cut algorithm */ 2325 SCIP_Real ubparam, /**< estimate of the optimal Lagrangian dual value */ 2326 int depth, /**< depth of the current node in the tree */ 2327 SCIP_Bool allowlocal, /**< should locally valid cuts be generated? */ 2328 int nmaxgeneratedperroundcuts,/**< maximal number of cuts allowed to generate per separation round */ 2329 SCIP_Real* origobjcoefs, /**< original objective function coefficients of the node linear relaxation */ 2330 SCIP_Real origobjoffset, /**< original objective function offset of the node linear relaxation */ 2331 SCIP_Real* dualvector, /**< Lagrangian multipliers vector */ 2332 int* nsoftcuts, /**< number of generated cuts that were penalized and added to the Lagrangian dual problem */ 2333 SCIP_ROW** generatedcurrroundcuts, /**< cuts generated in the current separation round */ 2334 SCIP_Real* generatedcutefficacies, /**< array of generated cut efficacies w.r.t. the respective LP bases used for cut 2335 generations */ 2336 int* ngeneratedcutsperiter, /**< number of cuts generated per subgradient iteration in the current separation round */ 2337 int* ngeneratedcurrroundcuts, /**< number of cuts generated so far in the current separation round */ 2338 int* ncurrroundlpiters, /**< number of LP iterations taken for solving Lagrangian dual problems with fixed 2339 multipliers in the current separator round */ 2340 SCIP_Bool* cutoff, /**< should the current node be cutoff? */ 2341 SCIP_Real* bestlagrangianval, /**< best Lagrangian value found so far */ 2342 SCIP_Real* bestdualvector, /**< Lagrangian multipliers corresponding to the best Lagrangian value found so far */ 2343 int* bestdualvectorlen, /**< length of the Lagrangian multipliers corresponding to the best Lagrangian value 2344 found so far */ 2345 int* nbestdualupdates, /**< number of best Lagrangian values found so far */ 2346 int* totaliternum /**< total number of iterations of the relax-and-cut algorithm performed so far */ 2347 ) 2348 { 2349 SCIP_Real* subgradient; 2350 SCIP_Real muparam; 2351 SCIP_Real steplength; 2352 SCIP_Real objval; 2353 SCIP_Real lagrangianval; 2354 SCIP_Real* lagrangianvals; 2355 SCIP_Real avglagrangianval; 2356 SCIP_Real maxsoftcutviol; 2357 SCIP_Real maxnzsubgradientdualprod; 2358 SCIP_Real maxviolscore; 2359 SCIP_Real maxviolscoreold; 2360 SCIP_Real nviolscore; 2361 SCIP_Real nviolscoreold; 2362 SCIP_Real scoreweight; 2363 SCIP_Real ballradius; 2364 SCIP_Bool solfound; 2365 SCIP_Bool backtrack; 2366 SCIP_Bool terminate; 2367 SCIP_Bool subgradientzero; 2368 SCIP_Bool objvecsdiffer; 2369 SCIP_Bool dualvecsdiffer; 2370 SCIP_Bool solvelp; 2371 int ncurrroundlpiterslast; 2372 int nlpiters; 2373 int ngeneratednewcuts; 2374 int nnewaddedsoftcuts; 2375 int nsoftcutviols; 2376 int nnzsubgradientdualprod; 2377 2378 SCIP_CALL( SCIPallocBufferArray(scip, &lagrangianvals, sepadata->nmaxsubgradientiters) ); 2379 SCIP_CALL( SCIPallocCleanBufferArray(scip, &subgradient, nmaxgeneratedperroundcuts) ); 2380 2381 muparam = sepadata->muparaminit; 2382 steplength = 0.0; 2383 objval = 0.0; 2384 avglagrangianval = 0.0; 2385 maxsoftcutviol = 0.0; 2386 maxnzsubgradientdualprod = 0.0; 2387 maxviolscore = 0.0; 2388 nviolscore = 0.0; 2389 scoreweight = 1.0; 2390 ballradius = sepadata->radiusinit; 2391 ngeneratednewcuts = 0; 2392 nsoftcutviols = 0; 2393 nnzsubgradientdualprod = 0; 2394 terminate = FALSE; 2395 subgradientzero = FALSE; 2396 objvecsdiffer = FALSE; 2397 dualvecsdiffer = FALSE; 2398 solvelp = TRUE; 2399 2400 /* update objective vector based on input Lagrangian multipliers */ 2401 if( *nsoftcuts > 0 ) 2402 { 2403 SCIP_CALL( updateObjectiveVector(scip, dualvector, generatedcurrroundcuts, *nsoftcuts, origobjcoefs, &objvecsdiffer) ); 2404 } 2405 2406 /* termination check */ 2407 SCIP_CALL( checkLagrangianDualTermination(sepadata, -1, -1, FALSE, *ngeneratedcurrroundcuts, 2408 nmaxgeneratedperroundcuts, *ncurrroundlpiters, depth, &terminate) ); 2409 2410 /* the subgradient algorithm loop */ 2411 for( int i = 0; i < sepadata->nmaxsubgradientiters && !SCIPisStopped(scip) && !terminate; i++ ) 2412 { 2413 solfound = FALSE; 2414 subgradientzero = FALSE; 2415 objvecsdiffer = FALSE; 2416 dualvecsdiffer = FALSE; 2417 nnewaddedsoftcuts = 0; 2418 scoreweight *= sepadata->radiusupdateweight; 2419 2420 ncurrroundlpiterslast = *ncurrroundlpiters; 2421 if( solvelp ) 2422 { 2423 /* solve Lagrangian dual for fixed Lagrangian multipliers */ 2424 SCIP_CALL( solveLagromoryLP(scip, sepadata, depth, origobjoffset, &solfound, sol, solvals, &objval, 2425 ncurrroundlpiters) ); 2426 } 2427 nlpiters = *ncurrroundlpiters - ncurrroundlpiterslast; 2428 2429 /* if an optimal solution has been found, then generate cuts and do other operations */ 2430 if( solfound ) 2431 { 2432 /* generate GMI cuts if a new basis solution is found */ 2433 if( (nlpiters >= 1) && (i % sepadata->cutgenfreq == 0) ) 2434 { 2435 ngeneratednewcuts = 0; 2436 SCIP_CALL( generateGMICuts(scip, sepa, sepadata, mainiternum, i, sol, solvals, 2437 nmaxgeneratedperroundcuts, allowlocal, generatedcurrroundcuts, generatedcutefficacies, 2438 *ngeneratedcurrroundcuts, &ngeneratednewcuts, depth, cutoff)); 2439 sepadata->ntotalcuts += ngeneratednewcuts; 2440 *ngeneratedcurrroundcuts += ngeneratednewcuts; 2441 ngeneratedcutsperiter[mainiternum * sepadata->nmaxsubgradientiters + i + 1] = ngeneratednewcuts; 2442 } 2443 2444 /* update subgradient, i.e., find the residuals of the penalized cuts, and determine various violations */ 2445 updateSubgradient(scip, sol, generatedcurrroundcuts, *nsoftcuts, subgradient, dualvector, &subgradientzero, 2446 &nsoftcutviols, &maxsoftcutviol, &nnzsubgradientdualprod, &maxnzsubgradientdualprod); 2447 2448 /* calculate Lagrangian value for the fixed Lagrangian multipliers, and update best and avg values */ 2449 updateLagrangianValue(scip, objval, dualvector, generatedcurrroundcuts, *nsoftcuts, &lagrangianval); 2450 if( SCIPisPositive(scip, lagrangianval - *bestlagrangianval) ) 2451 { 2452 *bestlagrangianval = lagrangianval; 2453 for( int j = 0; j < *nsoftcuts; j++ ) 2454 { 2455 bestdualvector[j] = dualvector[j]; 2456 } 2457 *bestdualvectorlen = *nsoftcuts; 2458 (*nbestdualupdates)++; 2459 } 2460 lagrangianvals[i] = lagrangianval; 2461 if( i < sepadata->nmaxlagrangianvalsforavg ) 2462 { 2463 avglagrangianval = (avglagrangianval * i + lagrangianval)/(i+1); 2464 } 2465 else 2466 { 2467 avglagrangianval = (avglagrangianval * sepadata->nmaxlagrangianvalsforavg - 2468 lagrangianvals[i - sepadata->nmaxlagrangianvalsforavg] + 2469 lagrangianval)/(sepadata->nmaxlagrangianvalsforavg); 2470 } 2471 2472 /* if the subgradient vector is non-zero, then update the mu parameter and the Lagrangian multipliers */ 2473 if( !subgradientzero ) 2474 { 2475 /* update mu param */ 2476 SCIP_CALL( updateMuSteplengthParam(scip, sepadata, i, ubparam, lagrangianvals, *bestlagrangianval, avglagrangianval, 2477 &muparam, &backtrack) ); 2478 2479 /* update step length */ 2480 updateStepLength(scip, muparam, ubparam, lagrangianval, subgradient, *nsoftcuts, &steplength); 2481 2482 /* update scores to determine how to update the stabilization ball radius */ 2483 maxviolscoreold = maxviolscore; 2484 nviolscoreold = nviolscore; 2485 maxviolscore = (1.0 - scoreweight) * maxsoftcutviol + scoreweight * maxnzsubgradientdualprod; 2486 nviolscore = (1.0 - scoreweight) * nsoftcutviols + scoreweight * nnzsubgradientdualprod; 2487 2488 /* update Lagrangian multipliers */ 2489 SCIP_CALL( updateDualVector(scip, sepadata, dualvector, bestdualvector, *bestdualvectorlen, 2490 *nbestdualupdates, i, *totaliternum, steplength, subgradient, *nsoftcuts, backtrack, maxviolscore, 2491 maxviolscoreold, nviolscore, nviolscoreold, nlpiters, &dualvecsdiffer, &ballradius) ); 2492 2493 /* update objective vector based on updated Lagrangian multipliers */ 2494 if( dualvecsdiffer ) 2495 { 2496 SCIP_CALL( updateObjectiveVector(scip, dualvector, generatedcurrroundcuts, *nsoftcuts, origobjcoefs, &objvecsdiffer) ); 2497 } 2498 } 2499 /* if the subgradient vector if zero, then simply mark that the Lagrangian multipliers and the objective 2500 * function of the Lagrangian dual did not change */ 2501 else 2502 { 2503 dualvecsdiffer = FALSE; 2504 objvecsdiffer = FALSE; 2505 } 2506 2507 2508 /* add generated GMI cuts to the objective function of the Lagrangian dual problem by introducing new 2509 * Lagrangian multipliers */ 2510 if( (i % sepadata->cutaddfreq == 0) || (!dualvecsdiffer && !objvecsdiffer && 2511 (*ngeneratedcurrroundcuts - *nsoftcuts > 0)) ) 2512 { 2513 SCIP_CALL( addGMICutsAsSoftConss(dualvector, *ngeneratedcurrroundcuts, nsoftcuts, &nnewaddedsoftcuts) ); 2514 } 2515 } 2516 else 2517 { 2518 /* add any remaining generated GMI cuts to the objective function of the Lagrangian dual problem by introducing 2519 * new Lagrangian multipliers */ 2520 if( (*ngeneratedcurrroundcuts - *nsoftcuts) > 0 ) 2521 { 2522 SCIP_CALL( addGMICutsAsSoftConss(dualvector, *ngeneratedcurrroundcuts, nsoftcuts, &nnewaddedsoftcuts) ); 2523 } 2524 2525 solvelp = FALSE; 2526 } 2527 2528 /* termination check */ 2529 SCIP_CALL( checkLagrangianDualTermination(sepadata, nnewaddedsoftcuts, *ngeneratedcurrroundcuts - *nsoftcuts, 2530 objvecsdiffer, *ngeneratedcurrroundcuts, nmaxgeneratedperroundcuts, *ncurrroundlpiters, depth, 2531 &terminate) ); 2532 2533 (*totaliternum)++; 2534 } 2535 2536 /* add any remaining generated GMI cuts to the objective function of the Lagrangian dual problem by introducing new 2537 * Lagrangian multipliers */ 2538 if( (*ngeneratedcurrroundcuts - *nsoftcuts) > 0 ) 2539 { 2540 SCIP_CALL( addGMICutsAsSoftConss(dualvector, *ngeneratedcurrroundcuts, nsoftcuts, &nnewaddedsoftcuts) ); 2541 } 2542 2543 /* free memory */ 2544 for( int i = 0; i < nmaxgeneratedperroundcuts; i++) 2545 { 2546 subgradient[i] = 0.0; 2547 } 2548 SCIPfreeCleanBufferArray(scip, &subgradient); 2549 SCIPfreeBufferArray(scip, &lagrangianvals); 2550 2551 return SCIP_OKAY; 2552 } 2553 2554 /** generates initial cut pool before solving the Lagrangian dual */ 2555 static 2556 SCIP_RETCODE generateInitCutPool( 2557 SCIP* scip, /**< SCIP data structure */ 2558 SCIP_SEPA* sepa, /**< separator */ 2559 SCIP_SEPADATA* sepadata, /**< separator data structure */ 2560 int mainiternum, /**< iteration number of the outer loop of the relax-and-cut algorithm */ 2561 SCIP_SOL* sol, /**< LP solution to be used for cut generation */ 2562 SCIP_Real* solvals, /**< values of the LP solution to be used for cut generation */ 2563 int nmaxgeneratedperroundcuts,/**< maximal number of cuts allowed to generate per separation round */ 2564 SCIP_Bool allowlocal, /**< should locally valid cuts be generated? */ 2565 SCIP_ROW** generatedcurrroundcuts, /**< cuts generated in the current separation round */ 2566 SCIP_Real* generatedcutefficacies, /**< array of generated cut efficacies w.r.t. the respective LP bases used for cut 2567 generations */ 2568 int* ngeneratedcutsperiter, /**< number of cuts generated per subgradient iteration in the current separation round */ 2569 int* ngeneratedcurrroundcuts, /**< number of cuts generated so far in the current separation round */ 2570 int depth, /**< depth of the current node in the tree */ 2571 SCIP_Bool* cutoff /**< should the current node be cutoff? */ 2572 ) 2573 { 2574 int ngeneratednewcuts; 2575 2576 ngeneratednewcuts = 0; 2577 2578 /* generate initial set of cuts */ 2579 SCIP_CALL( generateGMICuts(scip, sepa, sepadata, mainiternum, -1, sol, solvals, nmaxgeneratedperroundcuts, 2580 allowlocal, generatedcurrroundcuts, generatedcutefficacies, *ngeneratedcurrroundcuts, &ngeneratednewcuts, 2581 depth, cutoff) ); 2582 2583 /* update certain statistics */ 2584 sepadata->ntotalcuts += ngeneratednewcuts; 2585 *ngeneratedcurrroundcuts += ngeneratednewcuts; 2586 ngeneratedcutsperiter[sepadata->nmaxsubgradientiters * mainiternum] = ngeneratednewcuts; 2587 2588 return SCIP_OKAY; 2589 } 2590 2591 /** add cuts to SCIP */ 2592 static 2593 SCIP_RETCODE addCuts( 2594 SCIP* scip, /**< SCIP data structure */ 2595 SCIP_SEPADATA* sepadata, /**< separator data structure */ 2596 SCIP_ROW** cuts, /**< cuts generated so far in the current separation round */ 2597 int ncuts, /**< number of cuts generated so far in the current separation round */ 2598 SCIP_Longint maxdnom, /**< maximum denominator in the rational representation of cuts */ 2599 SCIP_Real maxscale, /**< maximal scale factor to scale the cuts to integral values */ 2600 int* naddedcuts, /**< number of cuts added to either global cutpool or sepastore */ 2601 SCIP_Bool* cutoff /**< should the current node be cutoff? */ 2602 ) 2603 { 2604 SCIP_ROW* cut; 2605 SCIP_Bool madeintegral; 2606 2607 assert(scip != NULL); 2608 assert(sepadata != NULL); 2609 assert(naddedcuts != NULL); 2610 assert(cutoff != NULL); 2611 2612 *cutoff = FALSE; 2613 madeintegral = FALSE; 2614 2615 for( int i = 0; i < ncuts && !*cutoff; i++ ) 2616 { 2617 cut = cuts[i]; 2618 2619 if( SCIProwGetNNonz(cut) == 1 ) 2620 { 2621 /* Add the bound change as cut to avoid that the LP gets modified. This would mean that the LP is not flushed 2622 * and the method SCIPgetLPBInvRow() fails; SCIP internally will apply this bound change automatically. */ 2623 SCIP_CALL( SCIPaddRow(scip, cut, TRUE, cutoff) ); 2624 ++(*naddedcuts); 2625 } 2626 else 2627 { 2628 if( sepadata->makeintegral && SCIPgetRowNumIntCols(scip, cut) == SCIProwGetNNonz(cut) ) 2629 { 2630 /* try to scale the cut to integral values */ 2631 SCIP_CALL( SCIPmakeRowIntegral(scip, cut, -SCIPepsilon(scip), SCIPsumepsilon(scip), 2632 maxdnom, maxscale, MAKECONTINTEGRAL, &madeintegral) ); 2633 2634 /* if RHS = plus infinity (due to scaling), the cut is useless, so we are not adding it */ 2635 if( madeintegral && SCIPisInfinity(scip, SCIProwGetRhs(cut)) ) 2636 return SCIP_OKAY; 2637 } 2638 2639 if( SCIPisCutNew(scip, cut) ) 2640 { 2641 /* add global cuts which are not implicit bound changes to the cut pool */ 2642 if( !SCIProwIsLocal(cut) ) 2643 { 2644 if( sepadata->delayedcuts ) 2645 { 2646 SCIP_CALL( SCIPaddDelayedPoolCut(scip, cut) ); 2647 } 2648 else 2649 { 2650 SCIP_CALL( SCIPaddPoolCut(scip, cut) ); 2651 } 2652 } 2653 else 2654 { 2655 /* local cuts we add to the sepastore */ 2656 SCIP_CALL( SCIPaddRow(scip, cut, sepadata->forcecuts, cutoff) ); 2657 } 2658 ++(*naddedcuts); 2659 } 2660 } 2661 } 2662 2663 return SCIP_OKAY; 2664 } 2665 2666 /** check different termination criteria */ 2667 /* @todo nlpssolved criterion? */ 2668 static 2669 SCIP_RETCODE checkMainLoopTermination( 2670 SCIP_SEPADATA* sepadata, /**< separator data structure */ 2671 SCIP_Bool cutoff, /**< should the current node be cutoff? */ 2672 SCIP_Bool dualvecsdiffer, /**< whether the updated Lagrangian multipliers differ from the old one */ 2673 int ngeneratedcurrroundcuts, /**< number of cuts generated in the current separation round */ 2674 int nsoftcuts, /**< number of generated cuts that were penalized and added to the Lagrangian dual problem */ 2675 int nmaxgeneratedperroundcuts,/**< maximal number of cuts allowed to generate per separation round */ 2676 int ncurrroundlpiters, /**< number of LP iterations taken for solving Lagrangian dual problems with fixed 2677 multipliers in the current separator round */ 2678 int depth, /**< depth of the current node in the tree */ 2679 SCIP_Bool* terminate /**< whether to terminate the relax-and-cut algorithm */ 2680 ) 2681 { 2682 *terminate = FALSE; 2683 2684 /* check if the node has been identified to be cutoff */ 2685 if( cutoff ) 2686 *terminate = TRUE; 2687 2688 /* check if the Lagrangian multipliers do not differ from the previous iteration and no new cuts exist for penalizing */ 2689 if( !dualvecsdiffer && (ngeneratedcurrroundcuts == nsoftcuts) ) 2690 *terminate = TRUE; 2691 2692 /* check if allowed number of cuts in this separation round have already been generated */ 2693 if( ngeneratedcurrroundcuts >= nmaxgeneratedperroundcuts ) 2694 *terminate = TRUE; 2695 2696 /* check if allowed number of cuts in the tree have already been generated */ 2697 if( sepadata->ntotalcuts >= sepadata->nmaxtotalcuts ) 2698 *terminate = TRUE; 2699 2700 /* check if allowed number of simplex iterations in this separation round have already been used up */ 2701 if( (sepadata->nmaxperroundlpiters >= 0) && (ncurrroundlpiters >= sepadata->nmaxperroundlpiters) ) 2702 *terminate = TRUE; 2703 2704 /* check if allowed number of simplex iterations in the root node have already been used up */ 2705 if( (depth == 0) && (sepadata->nmaxrootlpiters >= 0) && (sepadata->nrootlpiters >= sepadata->nmaxrootlpiters) ) 2706 *terminate = TRUE; 2707 2708 /* check if allowed number of simplex iterations in the tree have already been used up */ 2709 if( (sepadata->nmaxtotallpiters >= 0) && (sepadata->ntotallpiters >= sepadata->nmaxtotallpiters) ) 2710 *terminate = TRUE; 2711 2712 return SCIP_OKAY; 2713 } 2714 2715 /** Searches and tries to add Lagromory cuts */ 2716 static 2717 SCIP_RETCODE separateCuts( 2718 SCIP* scip, /**< SCIP data structure */ 2719 SCIP_SEPA* sepa, /**< separator */ 2720 SCIP_SEPADATA* sepadata, /**< separator data structure */ 2721 SCIP_Real ubparam, /**< estimate of the optimal Lagrangian dual value */ 2722 int depth, /**< depth of the current node in the tree */ 2723 SCIP_Bool allowlocal, /**< should locally valid cuts be generated? */ 2724 SCIP_RESULT* result /**< final result of the separation round */ 2725 ) 2726 { 2727 SCIP_ROW** generatedcurrroundcuts; 2728 SCIP_ROW** aggregatedcurrroundcuts; 2729 SCIP_Real* generatedcutefficacies; 2730 SCIP_Bool solfound; 2731 SCIP_SOL* softcutslpsol; 2732 SCIP_Real* softcutslpsolvals; 2733 SCIP_SOL* hardcutslpsol; 2734 SCIP_Real* hardcutslpsolvals; 2735 SCIP_Real* dualsol; 2736 SCIP_Real* dualvector; 2737 SCIP_Real* bestdualvector; 2738 SCIP_Real bestlagrangianval; 2739 SCIP_Real* origobjcoefs; 2740 SCIP_Real origobjoffset; 2741 SCIP_Real objval; 2742 SCIP_Real maxscale; 2743 SCIP_Longint maxdnom; 2744 SCIP_Bool cutoff; 2745 SCIP_Bool cutoff2; 2746 SCIP_Bool dualvecsdiffer; 2747 SCIP_Bool terminate; 2748 SCIP_COL** cols; 2749 2750 int* ngeneratedcutsperiter; 2751 int bestdualvectorlen; 2752 int nbestdualupdates; 2753 int totaliternum; 2754 int* cutindsperm; 2755 int nprocessedcuts; 2756 int ncols; 2757 int nrows; 2758 int ngeneratedcurrroundcuts; 2759 int nselectedcurrroundcuts; 2760 int nselectedcuts; 2761 int naddedcurrroundcuts; 2762 int naggregatedcurrroundcuts; 2763 int nmaxgeneratedperroundcuts; 2764 int ncurrroundlpiters; 2765 int nsoftcuts; 2766 int nsoftcutsold; 2767 int maxdepth; 2768 2769 assert(*result == SCIP_DIDNOTRUN); 2770 assert(sepadata != NULL); 2771 sepadata->ncalls = SCIPsepaGetNCalls(sepa); 2772 objval = SCIPinfinity(scip); 2773 bestlagrangianval = -SCIPinfinity(scip); 2774 bestdualvectorlen = 0; 2775 nbestdualupdates = 0; 2776 totaliternum = 0; 2777 ngeneratedcurrroundcuts = 0; 2778 naddedcurrroundcuts = 0; 2779 naggregatedcurrroundcuts = 0; 2780 ncurrroundlpiters = 0; 2781 nsoftcuts = 0; 2782 solfound = FALSE; 2783 cutoff = FALSE; 2784 cutoff2 = FALSE; 2785 dualvecsdiffer = FALSE; 2786 terminate = FALSE; 2787 2788 SCIPdebugMsg(scip, "Separating cuts...\n"); 2789 2790 /* initialize the LP that will be used to solve the Lagrangian dual with fixed Lagrangian multipliers */ 2791 SCIP_CALL( createLPWithSoftCuts(scip, sepadata) ); 2792 /* create the LP that represents the node relaxation including all the generated cuts in this separator */ 2793 SCIP_CALL( createLPWithHardCuts(scip, sepadata, NULL, 0) ); 2794 2795 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) ); 2796 nrows = SCIPgetNLPRows(scip); 2797 2798 /* get the maximal number of cuts allowed in a separation round */ 2799 nmaxgeneratedperroundcuts = ((depth == 0) ? sepadata->nmaxperroundcutsroot : sepadata->nmaxperroundcuts); 2800 2801 /* allocate memory */ 2802 SCIP_CALL( SCIPallocBufferArray(scip, &generatedcurrroundcuts, nmaxgeneratedperroundcuts) ); 2803 SCIP_CALL( SCIPallocBufferArray(scip, &aggregatedcurrroundcuts, nmaxgeneratedperroundcuts) ); 2804 SCIP_CALL( SCIPallocBufferArray(scip, &generatedcutefficacies, nmaxgeneratedperroundcuts) ); 2805 SCIP_CALL( SCIPallocBufferArray(scip, &cutindsperm, nmaxgeneratedperroundcuts) ); 2806 SCIP_CALL( SCIPallocCleanBufferArray(scip, &ngeneratedcutsperiter, sepadata->nmaxmainiters * 2807 (sepadata->nmaxsubgradientiters + 1)) ); 2808 SCIP_CALL( SCIPallocCleanBufferArray(scip, &dualsol, nrows + nmaxgeneratedperroundcuts) ); 2809 SCIP_CALL( SCIPallocCleanBufferArray(scip, &dualvector, nmaxgeneratedperroundcuts) ); 2810 SCIP_CALL( SCIPallocCleanBufferArray(scip, &bestdualvector, nmaxgeneratedperroundcuts) ); 2811 SCIP_CALL( SCIPallocBufferArray(scip, &softcutslpsolvals, ncols) ); 2812 SCIP_CALL( SCIPcreateSol(scip, &softcutslpsol, NULL) ); 2813 SCIP_CALL( SCIPallocBufferArray(scip, &hardcutslpsolvals, ncols) ); 2814 SCIP_CALL( SCIPcreateSol(scip, &hardcutslpsol, NULL) ); 2815 SCIP_CALL( SCIPallocBufferArray(scip, &origobjcoefs, ncols) ); 2816 2817 /* store current objective function */ 2818 for( int i = 0; i < ncols; i++ ) 2819 { 2820 origobjcoefs[i] = SCIPcolGetObj(cols[i]); 2821 } 2822 origobjoffset = SCIPgetTransObjoffset(scip); 2823 2824 /* solve node LP relaxation to have an initial simplex tableau */ 2825 SCIP_CALL( solveLagromoryLP(scip, sepadata, depth, origobjoffset, &solfound, softcutslpsol, softcutslpsolvals, &objval, 2826 &ncurrroundlpiters)); 2827 2828 /* generate initial cut pool */ 2829 SCIP_CALL( generateInitCutPool(scip, sepa, sepadata, 0, softcutslpsol, softcutslpsolvals, nmaxgeneratedperroundcuts, allowlocal, 2830 generatedcurrroundcuts, generatedcutefficacies, ngeneratedcutsperiter, &ngeneratedcurrroundcuts, depth, 2831 &cutoff) ); 2832 2833 /* termination check */ 2834 SCIP_CALL( checkMainLoopTermination(sepadata, cutoff, TRUE, ngeneratedcurrroundcuts, nsoftcuts, 2835 nmaxgeneratedperroundcuts, ncurrroundlpiters, depth, &terminate) ); 2836 2837 /* compute cuts for each integer col with fractional val */ 2838 for( int i = 0; i < sepadata->nmaxmainiters && !SCIPisStopped(scip) && !terminate; ++i ) 2839 { 2840 nsoftcutsold = nsoftcuts; 2841 nsoftcuts = ngeneratedcurrroundcuts; 2842 dualvecsdiffer = FALSE; 2843 2844 /* solve Lagrangain dual */ 2845 SCIP_CALL( solveLagrangianDual(scip, sepa, sepadata, softcutslpsol, softcutslpsolvals, i, ubparam, depth, allowlocal, 2846 nmaxgeneratedperroundcuts, origobjcoefs, origobjoffset, dualvector, &nsoftcuts, generatedcurrroundcuts, 2847 generatedcutefficacies, ngeneratedcutsperiter, &ngeneratedcurrroundcuts, &ncurrroundlpiters, &cutoff, 2848 &bestlagrangianval, bestdualvector, &bestdualvectorlen, &nbestdualupdates, &totaliternum) ); 2849 2850 /* @todo filter cuts before adding them to the new LP that was created based on the node relaxation? */ 2851 2852 /* update the LP representing the node relaxation by adding newly generated cuts */ 2853 if( !cutoff && (ngeneratedcurrroundcuts - nsoftcutsold > 0) ) 2854 { 2855 SCIP_CALL( createLPWithHardCuts(scip, sepadata, generatedcurrroundcuts, ngeneratedcurrroundcuts) ); 2856 2857 /* solve the LP and get relevant information */ 2858 SCIP_CALL( solveLPWithHardCuts(scip, sepadata, &solfound, hardcutslpsol, hardcutslpsolvals) ); 2859 2860 /* if primal solution is found, then pass it to trysol heuristic */ 2861 /* @note if trysol heuristic is was not present, then the solution found above, which can potentially be a MIP 2862 * primal feasible solution, will go to waste. 2863 */ 2864 if( solfound && sepadata->heurtrysol != NULL ) 2865 { 2866 SCIP_CALL( SCIPheurPassSolTrySol(scip, sepadata->heurtrysol, hardcutslpsol) ); 2867 } 2868 2869 /* if dual feasible, then fetch dual solution and reset Lagrangian multipliers based on it. otherwise, retain the 2870 * Lagrangian multipliers and simply initialize the new multipliers to zeroes. */ 2871 if( SCIPlpiIsDualFeasible(sepadata->lpiwithhardcuts) ) 2872 { 2873 SCIP_CALL( SCIPlpiGetSol(sepadata->lpiwithhardcuts, NULL, NULL, dualsol, NULL, NULL) ); 2874 SCIP_CALL( updateDualVector(scip, sepadata, dualvector, &(dualsol[nrows]), 2875 ngeneratedcurrroundcuts, 0, -1, -1, 0.0, NULL, ngeneratedcurrroundcuts, TRUE, 0.0, 0.0, 0.0, 0.0, -1, 2876 &dualvecsdiffer, NULL) ); 2877 } 2878 else 2879 { 2880 SCIP_CALL( updateDualVector(scip, sepadata, dualvector, dualvector, nsoftcuts, 0, -1, -1, 0.0, NULL, 2881 ngeneratedcurrroundcuts, TRUE, 0.0, 0.0, 0.0, 0.0, -1, &dualvecsdiffer, NULL) ); 2882 } 2883 } 2884 2885 /* termination check */ 2886 SCIP_CALL( checkMainLoopTermination(sepadata, cutoff, dualvecsdiffer, ngeneratedcurrroundcuts, nsoftcuts, 2887 nmaxgeneratedperroundcuts, ncurrroundlpiters, depth, &terminate) ); 2888 } 2889 2890 /* set the maximal denominator in rational representation of gomory cut and the maximal scale factor to 2891 * scale resulting cut to integral values to avoid numerical instabilities 2892 */ 2893 /**@todo find better but still stable gomory cut settings: look at dcmulti, gesa3, khb0525, misc06, p2756 */ 2894 /* @note: above todo was copied from sepa_gomory.c. So, if gomory code is changed, same changes can be done here. */ 2895 maxdepth = SCIPgetMaxDepth(scip); 2896 if( depth == 0 ) 2897 { 2898 maxdnom = 1000; 2899 maxscale = 1000.0; 2900 } 2901 else if( depth <= maxdepth/4 ) 2902 { 2903 maxdnom = 1000; 2904 maxscale = 1000.0; 2905 } 2906 else if( depth <= maxdepth/2 ) 2907 { 2908 maxdnom = 100; 2909 maxscale = 100.0; 2910 } 2911 else 2912 { 2913 maxdnom = 10; 2914 maxscale = 10.0; 2915 } 2916 2917 /* filter cuts by calling cut selection algorithms and add cuts accordingly */ 2918 /* @todo an idea is to filter cuts after every main iter */ 2919 /* @todo we can remove !cutoff criterion by adding forcedcuts */ 2920 if( !cutoff && (ngeneratedcurrroundcuts > 0) ) 2921 { 2922 if( SCIPisGE(scip, sepadata->cutsfilterfactor, 1.0) ) 2923 { 2924 nselectedcurrroundcuts = ngeneratedcurrroundcuts; 2925 SCIP_CALL( addCuts(scip, sepadata, generatedcurrroundcuts, nselectedcurrroundcuts, maxdnom, maxscale, 2926 &naddedcurrroundcuts, &cutoff2) ); 2927 cutoff = cutoff2; 2928 } 2929 else if( SCIPisPositive(scip, sepadata->cutsfilterfactor) ) 2930 { 2931 nprocessedcuts = 0; 2932 for( int i = 0; i < sepadata->nmaxmainiters * (sepadata->nmaxsubgradientiters + 1); i++ ) 2933 { 2934 if( ngeneratedcutsperiter[i] != 0 ) 2935 { 2936 for( int j = 0; j < ngeneratedcutsperiter[i]; j++ ) 2937 cutindsperm[j] = j + nprocessedcuts; 2938 2939 /* sort cut efficacies by fractionality */ 2940 SCIPsortDownRealInt(&generatedcutefficacies[nprocessedcuts], cutindsperm, ngeneratedcutsperiter[i]); 2941 2942 nselectedcuts = (int)SCIPceil(scip, sepadata->cutsfilterfactor * ngeneratedcutsperiter[i]); 2943 2944 SCIP_CALL( addCuts(scip, sepadata, &generatedcurrroundcuts[nprocessedcuts], nselectedcuts, maxdnom, 2945 maxscale, &naddedcurrroundcuts, &cutoff2) ); 2946 cutoff = cutoff2; 2947 2948 nprocessedcuts += ngeneratedcutsperiter[i]; 2949 } 2950 } 2951 } 2952 } 2953 else if( ngeneratedcurrroundcuts > 0 ) 2954 { 2955 nselectedcurrroundcuts = ngeneratedcurrroundcuts; 2956 SCIP_CALL( addCuts(scip, sepadata, generatedcurrroundcuts, nselectedcurrroundcuts, maxdnom, maxscale, 2957 &naddedcurrroundcuts, &cutoff2) ); 2958 } 2959 2960 /* add an aggregated cut based on best Lagrangian multipliers */ 2961 if( sepadata->aggregatecuts && (ngeneratedcurrroundcuts > 0) && (bestdualvectorlen > 0) ) 2962 { 2963 assert(bestdualvectorlen <= ngeneratedcurrroundcuts); 2964 SCIP_CALL( aggregateGeneratedCuts(scip, sepa, sepadata, generatedcurrroundcuts, bestdualvector, bestdualvectorlen, 2965 aggregatedcurrroundcuts, &naggregatedcurrroundcuts, &cutoff2) ); 2966 cutoff = (!cutoff ? cutoff2 : cutoff); 2967 if( naggregatedcurrroundcuts > 0 ) 2968 { 2969 SCIP_CALL( addCuts(scip, sepadata, aggregatedcurrroundcuts, naggregatedcurrroundcuts, maxdnom, maxscale, 2970 &naddedcurrroundcuts, &cutoff2) ); 2971 cutoff = (!cutoff ? cutoff2 : cutoff); 2972 } 2973 } 2974 2975 if( cutoff ) 2976 { 2977 *result = SCIP_CUTOFF; 2978 } 2979 else if( naddedcurrroundcuts > 0 ) 2980 { 2981 *result = SCIP_SEPARATED; 2982 } 2983 else 2984 { 2985 *result = SCIP_DIDNOTFIND; 2986 } 2987 2988 /* delete the LP representing the Lagrangian dual problem with fixed Lagrangian multipliers */ 2989 SCIP_CALL( deleteLPWithSoftCuts(scip, sepadata) ); 2990 2991 /* release the rows in aggregatedcurrroundcuts */ 2992 for( int i = 0; i < naggregatedcurrroundcuts; ++i ) 2993 { 2994 SCIP_CALL( SCIPreleaseRow(scip, &(aggregatedcurrroundcuts[i])) ); 2995 } 2996 /* release the rows in generatedcurrroundcuts */ 2997 for( int i = 0; i < ngeneratedcurrroundcuts; ++i ) 2998 { 2999 SCIP_CALL( SCIPreleaseRow(scip, &(generatedcurrroundcuts[i])) ); 3000 } 3001 3002 for( int i = 0; i < sepadata->nmaxmainiters * (sepadata->nmaxsubgradientiters + 1); i++ ) 3003 { 3004 ngeneratedcutsperiter[i] = 0; 3005 } 3006 for( int i = 0; i < nrows; i++ ) 3007 { 3008 dualsol[i] = 0.0; 3009 } 3010 for( int i = 0; i < nmaxgeneratedperroundcuts; i++ ) 3011 { 3012 dualsol[nrows + i] = 0.0; 3013 dualvector[i] = 0.0; 3014 bestdualvector[i] = 0.0; 3015 } 3016 /* free memory */ 3017 SCIPfreeBufferArray(scip, &origobjcoefs); 3018 SCIPfreeBufferArray(scip, &hardcutslpsolvals); 3019 SCIP_CALL( SCIPfreeSol(scip, &hardcutslpsol) ); 3020 SCIPfreeBufferArray(scip, &softcutslpsolvals); 3021 SCIP_CALL( SCIPfreeSol(scip, &softcutslpsol) ); 3022 SCIPfreeCleanBufferArray(scip, &ngeneratedcutsperiter); 3023 SCIPfreeBufferArray(scip, &cutindsperm); 3024 SCIPfreeBufferArray(scip, &generatedcutefficacies); 3025 SCIPfreeBufferArray(scip, &aggregatedcurrroundcuts); 3026 SCIPfreeBufferArray(scip, &generatedcurrroundcuts); 3027 SCIPfreeCleanBufferArray(scip, &dualvector); 3028 SCIPfreeCleanBufferArray(scip, &bestdualvector); 3029 SCIPfreeCleanBufferArray(scip, &dualsol); 3030 3031 return SCIP_OKAY; 3032 } 3033 3034 /** creates separator data */ 3035 static 3036 SCIP_RETCODE sepadataCreate( 3037 SCIP* scip, /**< SCIP data structure */ 3038 SCIP_SEPADATA** sepadata /**< separator data structure */ 3039 ) 3040 { 3041 assert(scip != NULL); 3042 assert(sepadata != NULL); 3043 3044 SCIP_CALL( SCIPallocBlockMemory(scip, sepadata) ); 3045 BMSclearMemory(*sepadata); 3046 3047 return SCIP_OKAY; 3048 } 3049 3050 3051 /* 3052 * Callback methods of separator 3053 */ 3054 3055 /** copy method for separator plugins (called when SCIP copies plugins) */ 3056 static 3057 SCIP_DECL_SEPACOPY(sepaCopyLagromory) 3058 { /*lint --e{715}*/ 3059 assert(scip != NULL); 3060 assert(sepa != NULL); 3061 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0); 3062 3063 /* call inclusion method of constraint handler */ 3064 SCIP_CALL( SCIPincludeSepaLagromory(scip) ); 3065 3066 return SCIP_OKAY; 3067 } 3068 3069 3070 /** destructor of separator to free user data (called when SCIP is exiting) */ 3071 static 3072 SCIP_DECL_SEPAFREE(sepaFreeLagromory) 3073 { 3074 SCIP_SEPADATA* sepadata; 3075 3076 sepadata = SCIPsepaGetData(sepa); 3077 assert(sepadata != NULL); 3078 3079 SCIP_CALL( sepadataFree(scip, &sepadata) ); 3080 SCIPsepaSetData(sepa, NULL); 3081 3082 return SCIP_OKAY; 3083 } 3084 3085 3086 /** initialization method of separator (called after problem was transformed) */ 3087 static 3088 SCIP_DECL_SEPAINIT(sepaInitLagromory) 3089 { 3090 SCIP_SEPADATA* sepadata; 3091 3092 sepadata = SCIPsepaGetData(sepa); 3093 assert(scip != NULL); 3094 assert(sepadata != NULL); 3095 3096 /* create and initialize random number generator */ 3097 SCIP_CALL( SCIPcreateRandom(scip, &(sepadata->randnumgen), RANDSEED, TRUE) ); 3098 3099 /* find trysol heuristic */ 3100 if ( sepadata->heurtrysol == NULL ) 3101 { 3102 sepadata->heurtrysol = SCIPfindHeur(scip, "trysol"); 3103 } 3104 3105 return SCIP_OKAY; 3106 } 3107 3108 /** deinitialization method of separator (called before transformed problem is freed) */ 3109 static 3110 SCIP_DECL_SEPAEXIT(sepaExitLagromory) 3111 { /*lint --e{715}*/ 3112 SCIP_SEPADATA* sepadata; 3113 3114 sepadata = SCIPsepaGetData(sepa); 3115 assert(sepadata != NULL); 3116 3117 SCIPfreeRandom(scip, &(sepadata->randnumgen)); 3118 3119 return SCIP_OKAY; 3120 } 3121 3122 /** LP solution separation method of separator */ 3123 static 3124 SCIP_DECL_SEPAEXECLP(sepaExeclpLagromory) 3125 { 3126 /*lint --e{715}*/ 3127 3128 SCIP_SEPADATA* sepadata; 3129 SCIP_SOL* bestsol; 3130 SCIP_COL** cols; 3131 SCIP_COL* col; 3132 SCIP_VAR* var; 3133 SCIP_Real ubparam; 3134 SCIP_Real dualdegeneracyrate; 3135 SCIP_Real varconsratio; 3136 SCIP_Real threshold1; 3137 SCIP_Real threshold2; 3138 int nrows; 3139 int ncols; 3140 int ncalls; 3141 int runnum; 3142 3143 assert(sepa != NULL); 3144 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0); 3145 sepadata = SCIPsepaGetData(sepa); 3146 assert(sepadata != NULL); 3147 3148 assert(result != NULL); 3149 *result = SCIP_DIDNOTRUN; 3150 3151 assert(scip != NULL); 3152 3153 ncalls = SCIPsepaGetNCallsAtNode(sepa); 3154 runnum = SCIPgetNRuns(scip); 3155 dualdegeneracyrate = 0.0; 3156 varconsratio = 0.0; 3157 3158 /* only generate Lagromory cuts if we are not close to terminating */ 3159 if( SCIPisStopped(scip) ) 3160 return SCIP_OKAY; 3161 3162 /* only call the separator starting sepadata->minrestart runs */ 3163 if( (sepadata->minrestart >= 1) && (runnum < sepadata->minrestart + 1) ) 3164 return SCIP_OKAY; 3165 3166 /* only call the separator a given number of times at each node */ 3167 if( (depth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot) 3168 || (depth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) ) 3169 return SCIP_OKAY; 3170 3171 /* only generate cuts if an optimal LP solution is at hand */ 3172 if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL ) 3173 return SCIP_OKAY; 3174 3175 /* only generate cuts if the LP solution is basic */ 3176 if( ! SCIPisLPSolBasic(scip) ) 3177 return SCIP_OKAY; 3178 3179 /* get LP data */ 3180 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) ); 3181 assert(cols != NULL); 3182 3183 nrows = SCIPgetNLPRows(scip); 3184 3185 /* return if LP has no columns or no rows */ 3186 if( ncols == 0 || nrows == 0 ) 3187 return SCIP_OKAY; 3188 3189 /* return if dual degeneracy metrics are below threshold values */ 3190 threshold1 = sepadata->dualdegeneracyratethreshold; 3191 threshold2 = sepadata->varconsratiothreshold; 3192 SCIP_CALL( SCIPgetLPDualDegeneracy(scip, &dualdegeneracyrate, &varconsratio) ); 3193 if( (dualdegeneracyrate < threshold1) && (varconsratio < threshold2) ) 3194 return SCIP_OKAY; 3195 3196 bestsol = SCIPgetBestSol(scip); 3197 ubparam = 0.0; 3198 3199 /* if a MIP primal solution exists, use it to estimate the optimal value of the Lagrangian dual problem */ 3200 if( bestsol != NULL ) 3201 { 3202 for( int i = 0; i < ncols; ++i ) 3203 { 3204 col = cols[i]; 3205 assert(col != NULL); 3206 3207 var = SCIPcolGetVar(col); 3208 3209 ubparam += SCIPgetSolVal(scip, bestsol, var) * SCIPcolGetObj(col); 3210 } 3211 3212 ubparam += SCIPgetTransObjoffset(scip); 3213 } 3214 /* if a MIP primal solution does not exist, then use the node relaxation's LP solutin to estimate the optimal value 3215 * of the Lagrangian dual problem 3216 */ 3217 else 3218 { 3219 for( int i = 0; i < ncols; ++i ) 3220 { 3221 col = cols[i]; 3222 assert(col != NULL); 3223 3224 ubparam += SCIPcolGetPrimsol(col) * SCIPcolGetObj(col); 3225 } 3226 3227 ubparam += SCIPgetTransObjoffset(scip); 3228 ubparam *= SCIPisPositive(scip, ubparam) ? sepadata->ubparamposfactor : sepadata->ubparamnegfactor; 3229 } 3230 3231 /* the main separation function of the separator */ 3232 SCIP_CALL( separateCuts(scip, sepa, sepadata, ubparam, depth, (allowlocal && sepadata->allowlocal), result) ); 3233 3234 return SCIP_OKAY; 3235 } 3236 3237 3238 /* 3239 * separator specific interface methods 3240 */ 3241 3242 /** creates the Lagromory separator and includes it in SCIP */ 3243 SCIP_RETCODE SCIPincludeSepaLagromory( 3244 SCIP* scip /**< SCIP data structure */ 3245 ) 3246 { 3247 SCIP_SEPADATA* sepadata; 3248 SCIP_SEPA* sepa; 3249 3250 /* create separator data */ 3251 SCIP_CALL( sepadataCreate(scip, &sepadata) ); 3252 3253 sepadata->heurtrysol = NULL; 3254 3255 /* include separator */ 3256 SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST, 3257 SEPA_USESSUBSCIP, SEPA_DELAY, sepaExeclpLagromory, NULL, sepadata) ); 3258 3259 assert(sepa != NULL); 3260 3261 /* set non fundamental callbacks via setter functions */ 3262 SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyLagromory) ); 3263 SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeLagromory) ); 3264 SCIP_CALL( SCIPsetSepaInit(scip, sepa, sepaInitLagromory) ); 3265 SCIP_CALL( SCIPsetSepaExit(scip, sepa, sepaExitLagromory) ); 3266 3267 /* add separator parameters */ 3268 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/away", "minimal integrality violation of a basis " 3269 "variable to try separation", &sepadata->away, FALSE, DEFAULT_AWAY, 0.0, 1.0, NULL, NULL) ); 3270 3271 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/rootlpiterlimitfactor", "factor w.r.t. root node LP " 3272 "iterations for maximal separating LP iterations in the root node (negative for no limit)", 3273 &sepadata->rootlpiterlimitfactor, TRUE, DEFAULT_ROOTLPITERLIMITFACTOR, -1.0, SCIP_REAL_MAX, NULL, NULL) ); 3274 3275 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/totallpiterlimitfactor", "factor w.r.t. root node LP " 3276 "iterations for maximal separating LP iterations in the tree (negative for no limit)", 3277 &sepadata->totallpiterlimitfactor, TRUE, DEFAULT_TOTALLPITERLIMITFACTOR, -1.0, SCIP_REAL_MAX, NULL, NULL) ); 3278 3279 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/perroundlpiterlimitfactor", "factor w.r.t. root node LP " 3280 "iterations for maximal separating LP iterations per separation round (negative for no limit)", 3281 &sepadata->perroundlpiterlimitfactor, TRUE, DEFAULT_PERROUNDLPITERLIMITFACTOR, -1.0, SCIP_REAL_MAX, NULL, 3282 NULL) ); 3283 3284 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/perroundcutsfactorroot", "factor w.r.t. number of integer " 3285 "columns for number of cuts separated per separation round in root node", &sepadata->perroundcutsfactorroot, 3286 TRUE, DEFAULT_PERROUNDCUTSFACTORROOT, 0.0, SCIP_REAL_MAX, NULL, NULL) ); 3287 3288 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/perroundcutsfactor", "factor w.r.t. number of integer " 3289 "columns for number of cuts separated per separation round at a non-root node", 3290 &sepadata->perroundcutsfactor, TRUE, DEFAULT_PERROUNDCUTSFACTOR, 0.0, SCIP_REAL_MAX, NULL, NULL) ); 3291 3292 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/totalcutsfactor", "factor w.r.t. number of integer " 3293 "columns for total number of cuts separated", &sepadata->totalcutsfactor, TRUE, DEFAULT_TOTALCUTSFACTOR, 0.0, 3294 SCIP_REAL_MAX, NULL, NULL) ); 3295 3296 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muparaminit", "initial value of the mu parameter (factor " 3297 "for step length)", &sepadata->muparaminit, TRUE, DEFAULT_MUPARAMINIT, 0.0, 100.0, NULL, NULL) ); 3298 3299 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muparamlb", "lower bound of the mu parameter (factor for " 3300 "step length)", &sepadata->muparamlb, TRUE, DEFAULT_MUPARAMLB, 0.0, 1.0, NULL, NULL) ); 3301 3302 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muparamub", "upper bound of the mu parameter (factor for " 3303 "step length)", &sepadata->muparamub, TRUE, DEFAULT_MUPARAMUB, 1.0, 10.0, NULL, NULL) ); 3304 3305 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/mubacktrackfactor", "factor of mu while backtracking the " 3306 "mu parameter (factor for step length)", &sepadata->mubacktrackfactor, TRUE, DEFAULT_MUBACKTRACKFACTOR, 3307 0.0, 1.0, NULL, NULL) ); 3308 3309 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muslab1factor", "factor of mu parameter (factor for step " 3310 "length) for larger increment" , &sepadata->muslab1factor, TRUE, DEFAULT_MUSLAB1FACTOR, 0.0, SCIP_REAL_MAX, NULL, 3311 NULL)); 3312 3313 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muslab2factor", "factor of mu parameter (factor for step " 3314 "length) for smaller increment", &sepadata->muslab2factor, TRUE, DEFAULT_MUSLAB2FACTOR, 0.0, SCIP_REAL_MAX, NULL, 3315 NULL) ); 3316 3317 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muslab3factor", "factor of mu parameter (factor for step " 3318 "length) for reduction", &sepadata->muslab3factor, TRUE, DEFAULT_MUSLAB3FACTOR, 0.0, SCIP_REAL_MAX, NULL, NULL) ); 3319 3320 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/deltaslab1ub", "factor of delta deciding larger increment " 3321 "of mu parameter (factor for step length)", &sepadata->deltaslab1ub, TRUE, DEFAULT_DELTASLAB1UB, 0.0, 1.0, 3322 NULL, NULL) ); 3323 3324 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/deltaslab2ub", "factor of delta deciding smaller " 3325 "increment of mu parameter (factor for step length)", &sepadata->deltaslab2ub, TRUE, DEFAULT_DELTASLAB2UB, 3326 0.0, 1.0, NULL, NULL) ); 3327 3328 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/ubparamposfactor", "factor for positive upper bound used " 3329 "as an estimate for the optimal Lagrangian dual value", &sepadata->ubparamposfactor, TRUE, DEFAULT_UBPARAMPOSFACTOR, 3330 1.0, 100.0, NULL, NULL) ); 3331 3332 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/ubparamnegfactor", "factor for negative upper bound used " 3333 "as an estimate for the optimal Lagrangian dual value", &sepadata->ubparamnegfactor, TRUE, DEFAULT_UBPARAMNEGFACTOR, 3334 0.0, 1.0, NULL, NULL) ); 3335 3336 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/perrootlpiterfactor", "factor w.r.t. root node LP " 3337 "iterations for iteration limit of each separating LP (negative for no limit)", 3338 &sepadata->perrootlpiterfactor, TRUE, DEFAULT_PERROOTLPITERFACTOR, -1.0, SCIP_REAL_MAX, NULL, NULL) ); 3339 3340 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/perlpiterfactor", "factor w.r.t. non-root node LP " 3341 "iterations for iteration limit of each separating LP (negative for no limit)", &sepadata->perlpiterfactor, TRUE, 3342 DEFAULT_PERLPITERFACTOR, -1.0, SCIP_REAL_MAX, NULL, NULL) ); 3343 3344 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/cutsfilterfactor", "fraction of generated cuts per " 3345 "explored basis to accept from separator", &sepadata->cutsfilterfactor, TRUE, DEFAULT_CUTSFILTERFACTOR, 0.0, 3346 1.0, NULL, NULL) ); 3347 3348 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/radiusinit", "initial radius of the ball used in " 3349 "stabilization of Lagrangian multipliers", &sepadata->radiusinit, TRUE, DEFAULT_RADIUSINIT, 0.0, 1.0, NULL, 3350 NULL) ); 3351 3352 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/radiusmax", "maximum radius of the ball used in " 3353 "stabilization of Lagrangian multipliers", &sepadata->radiusmax, TRUE, DEFAULT_RADIUSMAX, 0.0, SCIP_REAL_MAX, 3354 NULL, NULL) ); 3355 3356 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/radiusmin", "minimum radius of the ball used in " 3357 "stabilization of Lagrangian multipliers", &sepadata->radiusmin, TRUE, DEFAULT_RADIUSMIN, 0.0, SCIP_REAL_MAX, 3358 NULL, NULL) ); 3359 3360 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/constant", "a constant for stablity center based " 3361 "stabilization of Lagrangian multipliers", &sepadata->constant, TRUE, DEFAULT_CONST, 2.0, SCIP_REAL_MAX, 3362 NULL, NULL) ); 3363 3364 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/radiusupdateweight", "multiplier to evaluate cut " 3365 "violation score used for updating ball radius", &sepadata->radiusupdateweight, TRUE, 3366 DEFAULT_RADIUSUPDATEWEIGHT, 0.0, 1.0, NULL, NULL) ); 3367 3368 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/dualdegeneracyratethreshold", "minimum dual degeneracy " 3369 "rate for separator execution", &sepadata->dualdegeneracyratethreshold, FALSE, 3370 DEFAULT_DUALDEGENERACYRATETHRESHOLD, 0.0, 1.0, NULL, NULL) ); 3371 3372 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/varconsratiothreshold", "minimum variable-constraint " 3373 "ratio on optimal face for separator execution", &sepadata->varconsratiothreshold, FALSE, 3374 DEFAULT_VARCONSRATIOTHRESHOLD, 1.0, SCIP_REAL_MAX, NULL, NULL) ); 3375 3376 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/muparamconst", "is the mu parameter (factor for step " 3377 "length) constant?" , &sepadata->muparamconst, TRUE, DEFAULT_MUPARAMCONST, NULL, NULL) ); 3378 3379 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/separaterows", "separate rows with integral slack?", 3380 &sepadata->separaterows, TRUE, DEFAULT_SEPARATEROWS, NULL, NULL) ); 3381 3382 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/sortcutoffsol", "sort fractional integer columns" 3383 "based on fractionality?" , &sepadata->sortcutoffsol, TRUE, DEFAULT_SORTCUTOFFSOL, NULL, NULL) ); 3384 3385 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/sidetypebasis", "choose side types of row (lhs/rhs) " 3386 "based on basis information?", &sepadata->sidetypebasis, TRUE, DEFAULT_SIDETYPEBASIS, NULL, NULL) ); 3387 3388 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/dynamiccuts", "should generated cuts be removed from " 3389 "LP if they are no longer tight?", &sepadata->dynamiccuts, FALSE, DEFAULT_DYNAMICCUTS, NULL, NULL) ); 3390 3391 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/makeintegral", "try to scale all cuts to integral " 3392 "coefficients?", &sepadata->makeintegral, TRUE, DEFAULT_MAKEINTEGRAL, NULL, NULL) ); 3393 3394 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/forcecuts", "force cuts to be added to the LP?", 3395 &sepadata->forcecuts, TRUE, DEFAULT_FORCECUTS, NULL, NULL) ); 3396 3397 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/delayedcuts", "should cuts be added to the delayed cut " 3398 "pool", &sepadata->delayedcuts, TRUE, DEFAULT_DELAYEDCUTS, NULL, NULL) ); 3399 3400 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/allowlocal", "should locally valid cuts be generated?", 3401 &sepadata->allowlocal, TRUE, DEFAULT_ALLOWLOCAL, NULL, NULL) ); 3402 3403 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/aggregatecuts", "aggregate all generated cuts using the " 3404 "Lagrangian multipliers?", &sepadata->aggregatecuts, TRUE, DEFAULT_AGGREGATECUTS, NULL, NULL) ); 3405 3406 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxrounds", "maximal number of separation rounds per node " 3407 "(-1: unlimited)", &sepadata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) ); 3408 3409 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxroundsroot", "maximal number of separation rounds in " 3410 "the root node (-1: unlimited)", &sepadata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, 3411 NULL) ); 3412 3413 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/perroundnmaxlpiters", "maximal number of separating LP " 3414 "iterations per separation round (-1: unlimited)", &sepadata->perroundnmaxlpiters, FALSE, 3415 DEFAULT_PERROUNDMAXLPITERS, -1, INT_MAX, NULL, NULL) ); 3416 3417 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxcutsperlp", "maximal number of cuts separated per " 3418 "Lagromory LP in the non-root node", &sepadata->nmaxcutsperlp, FALSE, DEFAULT_PERLPMAXCUTS, 0, INT_MAX, NULL, 3419 NULL)); 3420 3421 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxcutsperlproot", "maximal number of cuts separated per " 3422 "Lagromory LP in the root node", &sepadata->nmaxcutsperlproot, FALSE, DEFAULT_PERLPMAXCUTSROOT, 0, INT_MAX, 3423 NULL, NULL)); 3424 3425 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxmainiters", "maximal number of main loop iterations of " 3426 "the relax-and-cut algorithm", &sepadata->nmaxmainiters, TRUE, DEFAULT_MAXMAINITERS, 0, INT_MAX, NULL, NULL) 3427 ); 3428 3429 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxsubgradientiters", "maximal number of subgradient loop " 3430 "iterations of the relax-and-cut algorithm", &sepadata->nmaxsubgradientiters, TRUE, 3431 DEFAULT_MAXSUBGRADIENTITERS, 0, INT_MAX, NULL, NULL) ); 3432 3433 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/cutgenfreq", "frequency of subgradient iterations for " 3434 "generating cuts", &sepadata->cutgenfreq, TRUE, DEFAULT_CUTGENFREQ, 0, INT_MAX, NULL, NULL) ); 3435 3436 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/cutaddfreq", "frequency of subgradient iterations for " 3437 "adding cuts to objective function", &sepadata->cutaddfreq, TRUE, DEFAULT_CUTADDFREQ, 0, INT_MAX, NULL, NULL) ); 3438 3439 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxlagrangianvalsforavg", "maximal number of iterations " 3440 "for rolling average of Lagrangian value", &sepadata->nmaxlagrangianvalsforavg, TRUE, 3441 DEFAULT_MAXLAGRANGIANVALSFORAVG, 0, INT_MAX, NULL, NULL) ); 3442 3443 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxconsecitersformuupdate", "consecutive number of " 3444 "iterations used to determine if mu needs to be backtracked", &sepadata->nmaxconsecitersformuupdate, TRUE, 3445 DEFAULT_MAXCONSECITERSFORMUUPDATE, 0, INT_MAX, NULL, NULL) ); 3446 3447 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/projectiontype", "the ball into which the Lagrangian multipliers " 3448 "are projected for stabilization (0: no projection, 1: L1-norm ball projection, 2: L2-norm ball projection, 3: " 3449 "L_inf-norm ball projection)", &sepadata->projectiontype, TRUE, DEFAULT_PROJECTIONTYPE, 0, 3, NULL, NULL)); 3450 3451 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/stabilitycentertype", "type of stability center for " 3452 "taking weighted average of Lagrangian multipliers for stabilization (0: no weighted stabilization, 1: best " 3453 "Lagrangian multipliers)", &sepadata->stabilitycentertype, TRUE, DEFAULT_STABILITYCENTERTYPE, 0, 1, NULL, 3454 NULL) ); 3455 3456 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/optimalfacepriority", "priority of the optimal face for " 3457 "separator execution (0: low priority, 1: medium priority, 2: high priority)", 3458 &sepadata->optimalfacepriority, TRUE, DEFAULT_OPTIMALFACEPRIORITY, 0, 2, NULL, NULL) ); 3459 3460 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/minrestart", "minimum restart round for separator " 3461 "execution (0: from beginning of the instance solving, >= n with n >= 1: from restart round n)", 3462 &sepadata->minrestart, TRUE, DEFAULT_MINRESTART, 0, INT_MAX, NULL, NULL) ); 3463 3464 return SCIP_OKAY; 3465 } 3466