Fix thermal constraints in MILP [ANT-2580]

docs/user-guide/solver/static-modeler/04-parameters.md

@@ -433,6 +433,17 @@ _**This section is under construction**_
 
 > _**Note:**_ You can find more information on this parameter [here](../03-appendix.md#details-on-the-include-exportmps-parameter).
 
+---
+#### include-export-solutions
+- **Expected value:** `true` or `false`
+- **Required:** no
+- **Default value:** `false`
+- **Usage:** set to `true` to activate writing the raw optimization results, that is
+  - values and reduced costs for each variable
+  - marginal cost for each constraint
+
+This is an advanced option intended to help developers and advanced users better understand their simulation results.
+
 ---
 #### include-split-exported-mps
 [//]: # (TODO: document this parameter, seems to belong to another category)

src/libs/antares/study/include/antares/study/parameters.h

@@ -397,6 +397,7 @@ class Parameters final
         //! Enum to define unfeasible problem behavior \see UnfeasibleProblemBehavior
         UnfeasibleProblemBehavior unfeasibleProblemBehavior;
 
+        bool exportSolutions;
     } include;
 
     struct Compatibility

src/libs/antares/study/parameters.cpp

@@ -368,6 +368,7 @@ void Parameters::reset()
 
     include.exportMPS = mpsExportStatus::NO_EXPORT;
     include.exportStructure = false;
+    include.exportSolutions = false;
     namedProblems = false;
 
     include.unfeasibleProblemBehavior = UnfeasibleProblemBehavior::ERROR_MPS;
@@ -721,6 +722,11 @@ static bool SGDIntLoadFamily_Optimization(Parameters& d,
         return true;
     }
 
+    if (key == "include-export-solutions")
+    {
+        return value.to<bool>(d.include.exportSolutions);
+    }
+
     if (key == "include-exportstructure")
     {
         return value.to<bool>(d.include.exportStructure);
@@ -1769,6 +1775,10 @@ void Parameters::prepareForSimulation(const StudyLoadOptions& options)
     {
         logs.info() << "  :: ignoring hurdle costs";
     }
+    if (!include.exportSolutions)
+    {
+        logs.info() << "  :: ignoring solution export";
+    }
 
     logs.info() << "  :: solver " << options.optOptions.ortoolsSolver
                 << " is used for problem resolution";
@@ -1888,7 +1898,9 @@ void Parameters::saveToINI(IniFile& ini) const
         section->add("include-primaryreserve", include.reserve.primary);
 
         section->add("include-exportmps", mpsExportStatusToString(include.exportMPS));
+
         section->add("include-exportstructure", include.exportStructure);
+        section->add("include-export-solutions", include.exportSolutions);
 
         // Unfeasible problem behavior
         section->add("include-unfeasible-problem-behavior",

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
 

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

src/solver/optimisation/opt_optimisation_lineaire.cpp

@@ -62,6 +62,20 @@ void OPT_EcrireResultatFonctionObjectiveAuFormatTXT(
     writer.addEntryFromBuffer(filename, buffer);
 }
 
+void OPT_WriteSolution(const PROBLEME_ANTARES_A_RESOUDRE& pb,
+                       const OptPeriodStringGenerator& optPeriodStringGenerator,
+                       int optimizationNumber,
+                       Solver::IResultWriter& writer)
+{
+    Yuni::Clob buffer;
+    auto filename = createSolutionFilename(optPeriodStringGenerator, optimizationNumber);
+    for (int var = 0; var < pb.NombreDeVariables; var++)
+    {
+        buffer.appendFormat("%s\t%11.10e\n", pb.NomDesVariables[var].c_str(), pb.X[var]);
+    }
+    writer.addEntryFromBuffer(filename, buffer);
+}
+
 namespace
 {
 void notifyProblemHebdo(const PROBLEME_HEBDO* problemeHebdo,
@@ -141,6 +155,13 @@ bool runWeeklyOptimization(const OptimizationOptions& options,
                                                            optimizationNumber,
                                                            writer);
         }
+        if (problemeHebdo->exportSolutions)
+        {
+            OPT_WriteSolution(*problemeHebdo->ProblemeAResoudre,
+                              *optPeriodStringGenerator,
+                              optimizationNumber,
+                              writer);
+        }
     }
     return true;
 }

src/solver/simulation/include/antares/solver/simulation/sim_structure_probleme_economique.h

@@ -532,6 +532,7 @@ struct PROBLEME_HEBDO
     bool exportMPSOnError = false;
     bool ExportStructure = false;
     bool NamedProblems = false;
+    bool exportSolutions = false;
 
     uint32_t HeureDansLAnnee = 0;
     bool LeProblemeADejaEteInstancie = false;

src/solver/simulation/sim_calcul_economique.cpp

@@ -122,6 +122,7 @@ void SIM_InitialisationProblemeHebdo(Data::Study& study,
     problem.NombreDeContraintesCouplantes = activeConstraints.size();
 
     problem.ExportMPS = study.parameters.include.exportMPS;
+    problem.exportSolutions = study.parameters.include.exportSolutions;
     problem.ExportStructure = study.parameters.include.exportStructure;
     problem.NamedProblems = study.parameters.namedProblems;
     problem.exportMPSOnError = Data::exportMPS(parameters.include.unfeasibleProblemBehavior);

src/solver/utils/filename.cpp

@@ -65,3 +65,9 @@ std::string createMPSfilename(const OptPeriodStringGenerator& optPeriodStringGen
 {
     return createOptimizationFilename("problem", optPeriodStringGenerator, optNumber, "mps");
 }
+
+std::string createSolutionFilename(const OptPeriodStringGenerator& optPeriodStringGenerator,
+                                   const unsigned int optNumber)
+{
+    return createOptimizationFilename("solution", optPeriodStringGenerator, optNumber, "csv");
+}

src/solver/utils/include/antares/solver/utils/filename.h

@@ -32,5 +32,9 @@ std::shared_ptr<OptPeriodStringGenerator> createOptPeriodAsString(bool isOptimiz
 
 std::string createCriterionFilename(const OptPeriodStringGenerator& optPeriodStringGenerator,
                                     const unsigned int optNumber);
+
 std::string createMPSfilename(const OptPeriodStringGenerator& optPeriodStringGenerator,
                               const unsigned int optNumber);
+
+std::string createSolutionFilename(const OptPeriodStringGenerator& optPeriodStringGenerator,
+                                   const unsigned int optNumber);

src/solver/variable/include/antares/solver/variable/state.h

@@ -121,14 +121,6 @@ class State
       const std::array<uint, HOURS_PER_YEAR>& ON_min,
       const std::array<uint, HOURS_PER_YEAR>& ON_max) const;
 
-    /*!
-    ** \brief Smooth the thermal units run after resolutions
-    ** using heuristics
-    **
-    ** \param areaWideIndex Index of the thermal cluster for the current area
-    */
-    void yearEndSmoothDispatchedUnitsCount(const unsigned int areaWideIndex, uint numSpace);
-
 public:
     /*!
     ** \brief Reset internal data

src/solver/variable/state.cpp

@@ -316,25 +316,24 @@ void State::yearEndBuildFromThermalClusterIndex(const uint clusterAreaWideIndex)
                            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:
@@ -344,13 +343,13 @@ void State::yearEndBuildFromThermalClusterIndex(const uint clusterAreaWideIndex)
         }
         }
 
-        ON_max[h] = static_cast<uint>(std::ceil(thermalClusterAvailableProduction
-                                                / currentCluster->nominalCapacityWithSpinning));
+        ON_max[h] = static_cast<uint>(Utils::ceilDiv(thermalClusterAvailableProduction,
+                                                     currentCluster->nominalCapacityWithSpinning));
 
         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;