Verify dual bound

coin-or · Oct 9, 2024 · fe62762 · fe62762
1 parent 2a3f4f6
commit fe62762
Show file tree

Hide file tree

Showing 17 changed files with 553 additions and 213 deletions.
diff --git a/src/DualSolver.cpp b/src/DualSolver.cpp
@@ -69,8 +69,11 @@ void DualSolver::checkDualSolutionCandidates()
 
         if(updateDual)
         {
-            // New dual solution
-            env->results->setDualBound(C.objValue);
+            if(C.sourceType == E_DualSolutionSource::ConvexBounding)
+                env->results->setDualBound(C.objValue, true); // Force valid dual bound update
+            else
+                env->results->setDualBound(C.objValue);
+
             currDualBound = C.objValue;
 
             if(env->results->getNumberOfIterations() > 0)
@@ -83,7 +86,8 @@ void DualSolver::checkDualSolutionCandidates()
 
             if(C.sourceType == E_DualSolutionSource::MIPSolutionOptimal
                 || C.sourceType == E_DualSolutionSource::LPSolution
-                || C.sourceType == E_DualSolutionSource::MIPSolverBound)
+                || C.sourceType == E_DualSolutionSource::MIPSolverBound
+                || C.sourceType == E_DualSolutionSource::ConvexBounding)
             {
                 env->results->addDualSolution(C);
             }
@@ -104,11 +108,21 @@ void DualSolver::checkDualSolutionCandidates()
             case E_DualSolutionSource::MIPSolverBound:
                 sourceDesc = "MIP solver bound";
                 break;
+            case E_DualSolutionSource::ConvexBounding:
+                sourceDesc = "Convex MIP bounding";
+                break;
             default:
                 break;
             }
 
-            env->output->outputDebug(fmt::format("        New dual bound {}, source: {}", C.objValue, sourceDesc));
+            if(C.sourceType == E_DualSolutionSource::ConvexBounding)
+            {
+                env->output->outputInfo(fmt::format("        New dual bound {}, source: {}", C.objValue, sourceDesc));
+            }
+            else
+            {
+                env->output->outputDebug(fmt::format("        New dual bound {}, source: {}", C.objValue, sourceDesc));
+            }
         }
     }
 
@@ -232,6 +246,11 @@ void DualSolver::addGeneratedHyperplane(const Hyperplane& hyperplane)
     currentIteration->totNumHyperplanes++;
     env->solutionStatistics.iterationLastDualCutAdded = currentIteration->iterationNumber;
 
+    if(hyperplane.isSourceConvex)
+        env->solutionStatistics.numberOfHyperplanesWithConvexSource++;
+    else
+        env->solutionStatistics.numberOfHyperplanesWithNonconvexSource++;
+
     env->output->outputTrace("        Hyperplane generated from: " + source);
 }
 

diff --git a/src/Enums.h b/src/Enums.h
@@ -46,7 +46,8 @@ enum class E_DualSolutionSource
     LPSolution,
     MIPSolutionOptimal,
     ObjectiveConstraint,
-    MIPSolverBound
+    MIPSolverBound,
+    ConvexBounding
 };
 
 enum class E_EventType
@@ -152,7 +153,8 @@ enum class E_PrimalSolutionSource
     MIPSolutionPool,
     LPFixedIntegers,
     MIPCallback,
-    InteriorPointSearch
+    InteriorPointSearch,
+    ConvexBounding
 };
 
 enum class E_ProblemConvexity

diff --git a/src/MIPSolver/MIPSolverCplex.cpp b/src/MIPSolver/MIPSolverCplex.cpp
@@ -1139,9 +1139,7 @@ void MIPSolverCplex::setTimeLimit(double seconds)
 {
     try
     {
-        if(seconds > 1e+75)
-        {
-        }
+        if(seconds > 1e+75) { }
         else if(seconds > 0)
         {
             cplexInstance.setParam(IloCplex::Param::TimeLimit, seconds);

diff --git a/src/PrimalSolver.cpp b/src/PrimalSolver.cpp
@@ -116,7 +116,7 @@ bool PrimalSolver::checkPrimalSolutionPoint(PrimalSolution primalSol)
         sourceDesc = "NLP fixed";
         break;
     case E_PrimalSolutionSource::MIPSolutionPool:
-        sourceDesc = "MILP sol. pool";
+        sourceDesc = "MIP sol. pool";
         break;
     case E_PrimalSolutionSource::LPFixedIntegers:
         sourceDesc = "LP fixed";
@@ -127,6 +127,9 @@ bool PrimalSolver::checkPrimalSolutionPoint(PrimalSolution primalSol)
     case E_PrimalSolutionSource::InteriorPointSearch:
         sourceDesc = "Interior point search";
         break;
+    case E_PrimalSolutionSource::ConvexBounding:
+        sourceDesc = "Convex MIP bounding";
+        break;
     default:
         sourceDesc = "other";
         break;

diff --git a/src/Results.cpp b/src/Results.cpp
@@ -121,14 +121,16 @@ void Results::addPrimalSolution(PrimalSolution solution)
     if(env->problem->objectiveFunction->properties.isMinimize)
     {
         std::sort(this->primalSolutions.begin(), this->primalSolutions.end(),
-            [](const PrimalSolution& firstSolution, const PrimalSolution& secondSolution)
-            { return (firstSolution.objValue < secondSolution.objValue); });
+            [](const PrimalSolution& firstSolution, const PrimalSolution& secondSolution) {
+                return (firstSolution.objValue < secondSolution.objValue);
+            });
     }
     else
     {
         std::sort(this->primalSolutions.begin(), this->primalSolutions.end(),
-            [](const PrimalSolution& firstSolution, const PrimalSolution& secondSolution)
-            { return (firstSolution.objValue > secondSolution.objValue); });
+            [](const PrimalSolution& firstSolution, const PrimalSolution& secondSolution) {
+                return (firstSolution.objValue > secondSolution.objValue);
+            });
     }
 
     env->solutionStatistics.numberOfFoundPrimalSolutions++;
@@ -1135,7 +1137,7 @@ double Results::getCurrentDualBound() { return (this->currentDualBound); }
 
 double Results::getGlobalDualBound() { return (globalDualBound); }
 
-void Results::setDualBound(double value)
+void Results::setDualBound(double value, bool forceGlobal)
 {
     double primalBound = this->getPrimalBound();
 
@@ -1152,8 +1154,12 @@ void Results::setDualBound(double value)
 
     this->currentDualBound = value;
 
-    if(this->solutionIsGlobal)
+    if(this->solutionIsGlobal || forceGlobal)
+    {
+        env->output->outputInfo(fmt::format(" New dual bound: {} ", this->currentDualBound));
         this->globalDualBound = value;
+        this->solutionIsGlobal = true;
+    }
 
     env->solutionStatistics.numberOfIterationsWithDualStagnation = 0;
 

diff --git a/src/Results.h b/src/Results.h
@@ -46,7 +46,7 @@ class DllExport Results
     void addDualSolution(DualSolution solution);
     double getCurrentDualBound();
     double getGlobalDualBound();
-    void setDualBound(double value);
+    void setDualBound(double value, bool forceGlobal = false);
 
     // For minimization problems, the lower bound is the dual while the upper bound is the primal objective value
     // for maximization problems, the lower bound is the primal while the upper bound is the dual objective value

diff --git a/src/SolutionStrategy/SolutionStrategyMIQCQP.cpp b/src/SolutionStrategy/SolutionStrategyMIQCQP.cpp
@@ -63,6 +63,8 @@
 #include "../Settings.h"
 #include "../Timing.h"
 
+#include "../Model/Problem.h"
+
 namespace SHOT
 {
 

diff --git a/src/SolutionStrategy/SolutionStrategyMultiTree.cpp b/src/SolutionStrategy/SolutionStrategyMultiTree.cpp
@@ -66,6 +66,8 @@
 
 #include "../Tasks/TaskAddIntegerCuts.h"
 
+#include "../Tasks/TaskPerformDualBounding.h"
+
 #include "../Output.h"
 #include "../Model/Problem.h"
 #include "../Model/ObjectiveFunction.h"
@@ -137,6 +139,12 @@ SolutionStrategyMultiTree::SolutionStrategyMultiTree(EnvironmentPtr envPtr)
     auto tSolveIteration = std::make_shared<TaskSolveIteration>(env);
     env->tasks->addTask(tSolveIteration, "SolveIter");
 
+    // if(env->reformulatedProblem->properties.convexity != E_ProblemConvexity::Convex)
+    //    {
+    auto tPerformDualBounding = std::make_shared<TaskPerformDualBounding>(env);
+    env->tasks->addTask(tPerformDualBounding, "DualBounding");
+    //}
+
     auto tSelectPrimSolPool = std::make_shared<TaskSelectPrimalCandidatesFromSolutionPool>(env);
     env->tasks->addTask(tSelectPrimSolPool, "SelectPrimSolPool");
     std::dynamic_pointer_cast<TaskSequential>(tFinalizeSolution)->addTask(tSelectPrimSolPool);

diff --git a/src/Solver.cpp b/src/Solver.cpp
@@ -1939,6 +1939,9 @@ void Solver::setConvexityBasedSettings()
 
             env->settings->updateSetting("BoundTightening.FeasibilityBased.TimeLimit", "Model", 5.0);
 
+            // Need to save these to perform dual bound updates
+            // env->settings->updateSetting("HyperplaneCuts.SaveHyperplanePoints", "Dual", true);
+
 #ifdef HAS_CPLEX
 
             if(static_cast<ES_MIPSolver>(env->settings->getSetting<int>("MIP.Solver", "Dual")) == ES_MIPSolver::Cplex)

diff --git a/src/Structs.h b/src/Structs.h
@@ -179,23 +179,17 @@ struct Hyperplane
     int sourceConstraintIndex; // -1 if objective function
     VectorDouble generatedPoint;
     double objectiveFunctionValue; // Used for the objective cuts only
-    E_HyperplaneSource source;
+    E_HyperplaneSource source = E_HyperplaneSource::None;
     bool isObjectiveHyperplane = false;
     bool isSourceConvex = false;
     double pointHash;
 };
 
-struct GeneratedHyperplane
+struct GeneratedHyperplane : Hyperplane
 {
-    NumericConstraintPtr sourceConstraint;
-    int sourceConstraintIndex; // -1 if objective function
-    VectorDouble generatedPoint;
-    E_HyperplaneSource source = E_HyperplaneSource::None;
     bool isLazy = false;
     bool isRemoved = false;
-    bool isSourceConvex = false;
     int iterationGenerated = -1;
-    double pointHash;
 };
 
 struct IntegerCut
@@ -231,6 +225,8 @@ struct SolutionStatistics
     int numberOfConstraintsRemovedInPresolve = 0;
     int numberOfVariableBoundsTightenedInPresolve = 0;
 
+    int numberOfHyperplanesWithConvexSource = 0;
+    int numberOfHyperplanesWithNonconvexSource = 0;
     int numberOfIntegerCuts = 0;
 
     int numberOfIterationsWithDualStagnation = 0;