From 99a50c00c6de694746b689ed2a01867a19743f79 Mon Sep 17 00:00:00 2001
From: Juliette-Gerbaux <130555142+Juliette-Gerbaux@users.noreply.github.com>
Date: Thu, 23 Jan 2025 14:11:20 +0100
Subject: [PATCH] Fix milp bug (#2489)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

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 <juliette.gerbaux_externe@rte-france.com>
Co-authored-by: payetvin <113102157+payetvin@users.noreply.github.com>
Co-authored-by: Vincent Payet <vincent.payet@rte-france.com>
Co-authored-by: Florian Omnès <florian.omnes@rte-france.com>
---
 .../utils/include/antares/utils/utils.h       |  4 ++++
 src/libs/antares/utils/utils.cpp              | 12 ++++++++++++
 src/solver/variable/state.cpp                 | 19 +++++++++----------
 3 files changed, 25 insertions(+), 10 deletions(-)

diff --git a/src/libs/antares/utils/include/antares/utils/utils.h b/src/libs/antares/utils/include/antares/utils/utils.h
index 47695188a7..4627783680 100644
--- 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
index c2c1286021..f7b2388eb9 100644
--- 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
index 0f96aef811..7774e9ac89 100644
--- a/src/solver/variable/state.cpp
+++ b/src/solver/variable/state.cpp
@@ -316,13 +316,14 @@ 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;
         }
@@ -330,11 +331,9 @@ void State::yearEndBuildFromThermalClusterIndex(const uint clusterEnabledIndex)
         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;