Fix milp bug (#2489)

Problème constaté par Amine sur le MILP avec Xpress. L'erreur est du à un mauvais ceiling sur des valeurs égales à n+10^-10 qui deviennent n+1 avec le ceil au lieu de n. --------- Co-authored-by: Juliette-Gerbaux <[email protected]> Co-authored-by: payetvin <[email protected]> Co-authored-by: Vincent Payet <[email protected]> Co-authored-by: Florian Omnès <[email protected]>
AntaresSimulatorTeam · Jan 23, 2025 · 99a50c0 · 99a50c0
1 parent d13d286
commit 99a50c0
Show file tree

Hide file tree

Showing 3 changed files with 25 additions and 10 deletions.
diff --git a/src/libs/antares/utils/include/antares/utils/utils.h b/src/libs/antares/utils/include/antares/utils/utils.h
@@ -53,8 +53,12 @@ std::vector<std::pair<std::string, std::string>> splitStringIntoPairs(const std:
 
 namespace Utils
 {
+
 bool isZero(double d);
 double round(double d, unsigned precision);
+double ceilDiv(double numerator, double denominator);
+double floorDiv(double numerator, double denominator);
+
 } // namespace Utils
 } // namespace Antares
 

diff --git a/src/libs/antares/utils/utils.cpp b/src/libs/antares/utils/utils.cpp
@@ -160,5 +160,17 @@ double round(double d, unsigned precision)
     return std::round(d * factor) / factor;
 }
 
+static constexpr double largeValue = 1000000;
+
+double ceilDiv(double numerator, double denominator)
+{
+    return std::ceil(std::round(numerator / denominator * largeValue) / largeValue);
+}
+
+double floorDiv(double numerator, double denominator)
+{
+    return std::floor(std::round(numerator / denominator * largeValue) / largeValue);
+}
+
 } // namespace Utils
 } // namespace Antares
diff --git a/src/solver/variable/state.cpp b/src/solver/variable/state.cpp
@@ -316,25 +316,24 @@ void State::yearEndBuildFromThermalClusterIndex(const uint clusterEnabledIndex)
                            static_cast<uint>(
                              std::ceil(thermalClusterAvailableProduction
                                        / currentCluster->nominalCapacityWithSpinning))),
-                  static_cast<uint>(std::ceil(thermalClusterProduction
-                                              / currentCluster->nominalCapacityWithSpinning)));
+                  static_cast<uint>(Utils::ceilDiv(thermalClusterProduction,
+                                                   currentCluster->nominalCapacityWithSpinning)));
             }
             else
             {
-                ON_min[h] = static_cast<uint>(std::ceil(
-                  thermalClusterProduction / currentCluster->nominalCapacityWithSpinning));
+                ON_min[h] = static_cast<uint>(
+                  Utils::ceilDiv(thermalClusterProduction,
+                                 currentCluster->nominalCapacityWithSpinning));
             }
             break;
         }
         case Antares::Data::UnitCommitmentMode::ucMILP:
         case Antares::Data::UnitCommitmentMode::ucHeuristicAccurate:
         {
             ON_min[h] = std::max(
-              static_cast<uint>(
-                std::ceil(thermalClusterProduction / currentCluster->nominalCapacityWithSpinning)),
-              thermalClusterDispatchedUnitsCountForYear[h]); // eq. to thermalClusterON for
-            // that hour
-
+              static_cast<uint>(Utils::ceilDiv(thermalClusterProduction,
+                                               currentCluster->nominalCapacityWithSpinning)),
+              thermalClusterDispatchedUnitsCountForYear[h]); // eq to thermalClusterON for that hour
             break;
         }
         case Antares::Data::UnitCommitmentMode::ucUnknown:
@@ -350,7 +349,7 @@ void State::yearEndBuildFromThermalClusterIndex(const uint clusterEnabledIndex)
         if (currentCluster->minStablePower > 0.)
         {
             maxUnitNeeded = static_cast<uint>(
-              std::floor(thermalClusterProduction / currentCluster->minStablePower));
+              Utils::floorDiv(thermalClusterProduction, currentCluster->minStablePower));
             if (ON_max[h] > maxUnitNeeded)
             {
                 ON_max[h] = maxUnitNeeded;