Add golden ratio search

Author: andreaslundell

src/Enums.h

@@ -372,6 +372,12 @@ enum class ES_PrimalNLPStrategy
     IterationOrTimeAndAllFeasibleSolutions
 };
 
+enum class ES_ReductionCutStrategy
+{
+    Fraction,
+    GoldenSection
+};
+
 enum class ES_ReformulationBinaryMonomials
 {
     None,

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.Use", "Dual", true, "Enable the dual reduction cut strategy for nonconvex problems");
+
+    VectorString enumReductionCutStrategy;
+    enumReductionCutStrategy.push_back("Fraction");
+    enumReductionCutStrategy.push_back("GoldenSection");
+
+    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", 5,
         "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",

src/Tasks/TaskAddPrimalReductionCut.cpp

@@ -72,42 +72,78 @@ 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::GoldenSection)
     {
-        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;
+
+        std::cout << "        Cut off diff: " << std::abs(currentLowerBoundForReductionCut - cutOffToUse) << std::endl;
+        if(std::abs(currentLowerBoundForReductionCut - cutOffToUse) < 1e-4)
         {
-            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;
+        std::cout << "        Cut off to use: " << cutOffToUse << std::endl;
+        std::cout << "        Primal improv. after repair: "
+                  << env->solutionStatistics.numberOfPrimalImprovementsAfterInfeasibilityRepair << std::endl;
+        std::cout << "        Primal improv. after red. cut: "
+                  << env->solutionStatistics.numberOfPrimalImprovementsAfterReductionCut << std::endl;
+    }
+
+    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;

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