Merge pull request #178 from coin-or/primalcut

Improve primal reduction cut heuristic
coin-or · Jan 3, 2025 · a81275b · a81275b
2 parents 750b43e + bb6f7f2
commit a81275b
Show file tree

Hide file tree

Showing 5 changed files with 89 additions and 29 deletions.
diff --git a/src/Enums.h b/src/Enums.h
@@ -372,6 +372,12 @@ enum class ES_PrimalNLPStrategy
     IterationOrTimeAndAllFeasibleSolutions
 };
 
+enum class ES_ReductionCutStrategy
+{
+    Fraction,
+    GoldenRatio
+};
+
 enum class ES_ReformulationBinaryMonomials
 {
     None,

diff --git a/src/Report.cpp b/src/Report.cpp
@@ -103,7 +103,7 @@ void Report::outputIterationDetail(int iterationNumber, std::string iterationDes
 
         if(lineType == E_IterationLineType::DualRepair)
         {
-            combDualCuts = fmt::format("Repairs: {:<4d}", dualCutsAdded);
+            combDualCuts = fmt::format("Relaxed constraints: {:<4d}", dualCutsAdded);
         }
         else if(lineType == E_IterationLineType::DualReductionCut)
         {
@@ -154,6 +154,11 @@ void Report::outputIterationDetail(int iterationNumber, std::string iterationDes
             env->output->outputInfo(fmt::format("{:>6d}: {:<10s}{:^10.2f}{:^13s}{:>27s}{:>19s}{:<32s}", iterationNumber,
                 iterationDesc, totalTime, combDualCuts, "", "", ""));
             env->output->outputDebug("");
+
+            env->output->outputInfo(fmt::format("{} {} {} / {}", std::string(27, ' '),
+                "Total primal improvements after repair / reduction cut:  ",
+                env->solutionStatistics.numberOfPrimalImprovementsAfterInfeasibilityRepair,
+                env->solutionStatistics.numberOfPrimalImprovementsAfterReductionCut));
         }
         else
         {

diff --git a/src/Solver.cpp b/src/Solver.cpp
@@ -896,15 +896,28 @@ void Solver::initializeSettings()
         "These settings control the added dual reduction cuts from the primal solution that will try to force a better "
         "primal solution. This functionality is only used if SHOT cannot deduce that the problem is nonconvex .");
 
-    env->settings->createSetting("ReductionCut.MaxIterations", "Dual", 5,
+    env->settings->createSetting(
+        "ReductionCut.Use", "Dual", true, "Enable the dual reduction cut strategy for nonconvex problems");
+
+    VectorString enumReductionCutStrategy;
+    enumReductionCutStrategy.push_back("Fraction");
+    enumReductionCutStrategy.push_back("GoldenRatio");
+
+    ES_ReductionCutStrategy reductionCutStrategy;
+
+    reductionCutStrategy = ES_ReductionCutStrategy::Fraction;
+
+    env->settings->createSetting("ReductionCut.Strategy", "Dual", static_cast<int>(reductionCutStrategy),
+        "The reduction cut strategy to use", enumReductionCutStrategy,
+        static_cast<int>(ES_ReductionCutStrategy::Fraction));
+    enumMIPSolver.clear();
+
+    env->settings->createSetting("ReductionCut.MaxIterations", "Dual", 20,
         "Max number of primal cut reduction without primal improvement", 0, SHOT_INT_MAX);
 
     env->settings->createSetting(
         "ReductionCut.ReductionFactor", "Dual", 0.001, "The factor used to reduce the cutoff value", 0, 1.0);
 
-    env->settings->createSetting(
-        "ReductionCut.Use", "Dual", true, "Enable the dual reduction cut strategy for nonconvex problems");
-
     // Dual strategy settings: Relaxation strategies
 
     env->settings->createSettingGroup("Dual", "Relaxation", "Relaxation strategies",
@@ -1141,9 +1154,6 @@ void Solver::initializeSettings()
         "Which formulation to use in eigenvalue decomposition", enumEigenValueStrategy, 0);
     enumEigenValueStrategy.clear();
 
-    env->settings->createSetting("Reformulation.Quadratics.EigenValueDecomposition.Method", "Model", false,
-        "Whether to use the eigen value decomposition of convex quadratic functions");
-
     env->settings->createSetting("Reformulation.Quadratics.EigenValueDecomposition.Tolerance", "Model", 1e-6,
         "Variables with eigenvalues smaller than this value will be ignored", 0.0, SHOT_DBL_MAX);
 

diff --git a/src/Tasks/TaskAddPrimalReductionCut.cpp b/src/Tasks/TaskAddPrimalReductionCut.cpp
@@ -72,42 +72,80 @@ void TaskAddPrimalReductionCut::run()
         return;
     }
 
-    double relativeGap = env->results->getRelativeCurrentObjectiveGap();
+    double cutOffToUse;
 
-    // Different logic if gap is large
-    if(relativeGap <= 1.0 && relativeGap > 0.1)
+    if(env->settings->getSetting<int>("ReductionCut.Strategy", "Dual") == (int)ES_ReductionCutStrategy::Fraction)
     {
-        double factor
-            = ((double)env->solutionStatistics.numberOfPrimalReductionCutsUpdatesWithoutEffect) / (maxIterations + 1.0);
+        double relativeGap = env->results->getRelativeCurrentObjectiveGap();
 
-        env->dualSolver->cutOffToUse
-            = factor * env->results->currentDualBound + (1 - factor) * env->results->currentPrimalBound;
-
-        if(env->reformulatedProblem->objectiveFunction->properties.isMinimize)
+        // Different logic if gap is large
+        if(relativeGap <= 1.0 && relativeGap > 0.1)
         {
-            env->results->currentDualBound = SHOT_DBL_MIN;
+            double factor = ((double)env->solutionStatistics.numberOfPrimalReductionCutsUpdatesWithoutEffect)
+                / (maxIterations + 1.0);
+
+            cutOffToUse = factor * env->results->currentDualBound + (1 - factor) * env->results->currentPrimalBound;
         }
         else
         {
-            env->results->currentDualBound = SHOT_DBL_MAX;
+            double reductionFactor = env->settings->getSetting<double>("ReductionCut.ReductionFactor", "Dual");
+
+            if(env->reformulatedProblem->objectiveFunction->properties.isMinimize)
+            {
+                cutOffToUse = env->dualSolver->cutOffToUse - reductionFactor * std::abs(env->dualSolver->cutOffToUse);
+            }
+            else
+            {
+                cutOffToUse = env->dualSolver->cutOffToUse + reductionFactor * std::abs(env->dualSolver->cutOffToUse);
+            }
         }
     }
-    else
+    else if(env->settings->getSetting<int>("ReductionCut.Strategy", "Dual")
+        == (int)ES_ReductionCutStrategy::GoldenRatio)
     {
-        double reductionFactor = env->settings->getSetting<double>("ReductionCut.ReductionFactor", "Dual");
+        double factor = 0.618;
 
-        if(env->reformulatedProblem->objectiveFunction->properties.isMinimize)
+        // If first cut iteration after PB update
+        if(env->solutionStatistics.numberOfPrimalReductionCutsUpdatesWithoutEffect == 0)
         {
-            env->dualSolver->cutOffToUse
-                = env->dualSolver->cutOffToUse - reductionFactor * std::abs(env->dualSolver->cutOffToUse);
-            env->results->currentDualBound = SHOT_DBL_MIN;
+            if(env->reformulatedProblem->objectiveFunction->properties.isMinimize)
+            {
+                currentLowerBoundForReductionCut = std::max(SHOT_DBL_MIN, env->results->globalDualBound);
+            }
+            else
+            {
+                currentLowerBoundForReductionCut = std::min(SHOT_DBL_MAX, env->results->globalDualBound);
+            }
         }
-        else
+
+        cutOffToUse = (1 - factor) * currentLowerBoundForReductionCut + factor * env->results->currentPrimalBound;
+
+        env->output->outputDebug(fmt::format(
+            "        {} {}", "Cut off difference:", std::abs(currentLowerBoundForReductionCut - cutOffToUse)));
+
+        if(std::abs(currentLowerBoundForReductionCut - cutOffToUse) < 1e-6)
         {
-            env->dualSolver->cutOffToUse
-                = env->dualSolver->cutOffToUse + reductionFactor * std::abs(env->dualSolver->cutOffToUse);
-            env->results->currentDualBound = SHOT_DBL_MAX;
+            env->tasks->setNextTask(taskIDIfFalse);
+            return;
         }
+
+        currentLowerBoundForReductionCut = cutOffToUse;
+    }
+    else
+    {
+        env->tasks->setNextTask(taskIDIfFalse);
+        return;
+    }
+
+    env->dualSolver->cutOffToUse = cutOffToUse;
+
+    if(env->reformulatedProblem->objectiveFunction->properties.isMinimize)
+    {
+        env->results->currentDualBound = SHOT_DBL_MIN;
+    }
+    else
+    {
+        env->results->currentDualBound = SHOT_DBL_MAX;
     }
 
     std::stringstream tmpType;

diff --git a/src/Tasks/TaskAddPrimalReductionCut.h b/src/Tasks/TaskAddPrimalReductionCut.h
@@ -27,5 +27,6 @@ class TaskAddPrimalReductionCut : public TaskBase
     std::string taskIDIfTrue;
     std::string taskIDIfFalse;
     int totalReductionCutUpdates = 0;
+    double currentLowerBoundForReductionCut = SHOT_DBL_INF;
 };
 } // namespace SHOT