[ANT-992] Add ramping for thermal clusters

src/libs/antares/study/parts/thermal/cluster.cpp

@@ -137,7 +137,11 @@ Data::ThermalCluster::ThermalCluster(Area* parent) :
     forcedLaw(thermalLawUniform),
     plannedLaw(thermalLawUniform),
     PthetaInf(HOURS_PER_YEAR, 0),
-    costsTimeSeries(1, CostsTimeSeries())
+    costsTimeSeries(1, CostsTimeSeries()),
+    maxUpwardPowerRampingRate(0.),
+    maxDownwardPowerRampingRate(0.),
+    powerIncreaseCost(0.),
+    powerDecreaseCost(0.)
 {
     // assert
     assert(parent and "A parent for a thermal dispatchable cluster can not be null");
@@ -182,6 +186,12 @@ void Data::ThermalCluster::copyFrom(const ThermalCluster& cluster)
     minUpTime = cluster.minUpTime;
     minDownTime = cluster.minDownTime;
 
+    // ramping
+    maxUpwardPowerRampingRate = cluster.maxUpwardPowerRampingRate;
+    maxDownwardPowerRampingRate = cluster.maxDownwardPowerRampingRate;
+    powerDecreaseCost = cluster.powerDecreaseCost;
+    powerIncreaseCost = cluster.powerIncreaseCost;
+
     // spinning
     spinning = cluster.spinning;
 
@@ -479,6 +489,11 @@ void Data::ThermalCluster::reset()
     marketBidCost = 0.;
     variableomcost = 0.;
     costsTimeSeries.resize(1, CostsTimeSeries());
+    powerIncreaseCost = 0;
+    powerDecreaseCost = 0;
+
+    maxUpwardPowerRampingRate = 0;
+    maxDownwardPowerRampingRate = 0;
 
     // modulation
     modulation.resize(thermalModulationMax, HOURS_PER_YEAR);
@@ -608,6 +623,36 @@ bool Data::ThermalCluster::integrityCheck()
             ret = false;
     }*/
 
+    // ramping
+    if (maxUpwardPowerRampingRate <= 0)
+    {
+        logs.error() << "Thermal cluster: " << parentArea->name << '/' << pName
+                     << ": The maximum upward power ramping rate must greater than zero.";
+        maxUpwardPowerRampingRate = 1.;
+        ret = false;
+    }
+    if (maxDownwardPowerRampingRate <= 0)
+    {
+        logs.error() << "Thermal cluster: " << parentArea->name << '/' << pName
+                     << ": The maximum downward power ramping rate must greater than zero.";
+        maxDownwardPowerRampingRate = 1.;
+        ret = false;
+    }
+    if (powerIncreaseCost < 0)
+    {
+        logs.error() << "Thermal cluster: " << parentArea->name << '/' << pName
+                     << ": The ramping power increase cost must be positive or null.";
+        powerIncreaseCost = 0.;
+        ret = false;
+    }
+    if (powerDecreaseCost < 0)
+    {
+        logs.error() << "Thermal cluster: " << parentArea->name << '/' << pName
+                     << ": The ramping power decrease cost must be positive or null.";
+        powerDecreaseCost = 0.;
+        ret = false;
+    }
+
     return ret;
 }
 

src/libs/antares/study/parts/thermal/cluster.h

@@ -340,6 +340,16 @@ class ThermalCluster final : public Cluster, public std::enable_shared_from_this
     double marketBidCost = 0;
     //! Variable O&M cost (euros/MWh)
     double variableomcost = 0;
+    //! Cost of power increase (euros/MW)
+    double powerIncreaseCost = 0;
+    //! Cost of power decrease (euros/MW)
+    double powerDecreaseCost = 0;
+
+    //! Maximum hourly upward power ramping rate (MW/hour)
+    double maxUpwardPowerRampingRate = 0;
+    //! Maximum hourly downward power ramping rate (MW/hour)
+    double maxDownwardPowerRampingRate = 0;
+
     //@}
 
     /*!

src/libs/antares/study/parts/thermal/cluster_list.cpp

@@ -222,6 +222,15 @@ static bool ThermalClusterLoadFromProperty(ThermalCluster& cluster, const IniFil
     if (p->key == "startup-cost")
         return p->value.to<double>(cluster.startupCost);
 
+    if (p->key == "power-increase-cost")
+        return p->value.to<double>(cluster.powerIncreaseCost);
+    if (p->key == "power-decrease-cost")
+        return p->value.to<double>(cluster.powerDecreaseCost);
+    if (p->key == "max-upward-power-ramping-rate")
+        return p->value.to<double>(cluster.maxUpwardPowerRampingRate);
+    if (p->key == "max-downward-power-ramping-rate")
+        return p->value.to<double>(cluster.maxDownwardPowerRampingRate);
+
     if (p->key == "unitcount")
         return p->value.to<uint>(cluster.unitCount);
     if (p->key == "volatility.planned")
@@ -405,6 +414,15 @@ bool ThermalClusterList::saveToFolder(const AnyString& folder) const
             if (!Math::Zero(c.variableomcost))
                 s->add("variableomcost", Math::Round(c.variableomcost,3));
 
+            // ramping
+            if (not Math::Zero(c.powerIncreaseCost))
+                s->add("power-increase-cost", Math::Round(c.powerIncreaseCost, 3));
+            if (not Math::Zero(c.powerDecreaseCost))
+                s->add("power-decrease-cost", Math::Round(c.powerDecreaseCost, 3));
+            if (not Math::Zero(c.maxUpwardPowerRampingRate))
+                s->add("max-upward-power-ramping-rate", Math::Round(c.maxUpwardPowerRampingRate, 3));
+            if (not Math::Zero(c.maxDownwardPowerRampingRate))
+                s->add("max-downward-power-ramping-rate", Math::Round(c.maxDownwardPowerRampingRate, 3));
 
             //pollutant factor
             for (auto const& [key, val] : Pollutant::namesToEnum)

src/solver/optimisation/CMakeLists.txt

@@ -38,6 +38,11 @@ set(RTESOLVER_OPT
 		opt_decompte_variables_et_contraintes_couts_demarrage.cpp
 		opt_init_minmax_groupes_couts_demarrage.cpp
 		opt_nombre_min_groupes_demarres_couts_demarrage.cpp
+		opt_construction_contraintes_rampes_thermiques.cpp
+		opt_decompte_variables_et_contraintes_rampes_thermiques.cpp
+		opt_construction_variables_rampes_thermiques.cpp
+		opt_gestion_des_bornes_rampes_thermiques.cpp
+		opt_gestion_des_couts_rampes_thermiques.cpp
 		opt_export_structure.h
 		opt_export_structure.cpp
 		base_weekly_optimization.h
@@ -119,7 +124,12 @@ set(RTESOLVER_OPT
 		constraints/NbDispUnitsMinBoundSinceMinUpTime.cpp
 		constraints/MinDownTime.h
 		constraints/MinDownTime.cpp
-
+		constraints/RampingIncreaseRate.h
+		constraints/RampingIncreaseRate.cpp
+		constraints/RampingDecreaseRate.h
+		constraints/RampingDecreaseRate.cpp
+		constraints/PowerOutputVariation.h
+		constraints/PowerOutputVariation.cpp
 	)
 
 

src/solver/optimisation/constraints/ConstraintBuilder.cpp

@@ -58,6 +58,25 @@ ConstraintBuilder& ConstraintBuilder::DispatchableProduction(unsigned int index,
                                                              int delta)
 {
     AddVariable(GetVariableManager(offset, delta).DispatchableProduction(index), coeff);
+    // logs.info() << "dispatchable production idx " << GetVariableManager(offset, delta).DispatchableProduction(index);
+    return *this;
+}
+
+ConstraintBuilder& ConstraintBuilder::ProductionDecreaseAboveMin(unsigned int index,
+                                                                 double coeff,
+                                                                 int offset,
+                                                                 int delta)
+{
+    AddVariable(GetVariableManager(offset, delta).ProductionDecreaseAboveMin(index), coeff);
+    return *this;
+}
+
+ConstraintBuilder& ConstraintBuilder::ProductionIncreaseAboveMin(unsigned int index,
+                                                                 double coeff,
+                                                                 int offset,
+                                                                 int delta)
+{
+    AddVariable(GetVariableManager(offset, delta).ProductionIncreaseAboveMin(index), coeff);
     return *this;
 }
 

src/solver/optimisation/constraints/ConstraintBuilder.h

@@ -31,6 +31,16 @@ class VariableManager
         return nativeOptimVar.NumeroDeVariableDuPalierThermique[index];
     }
 
+    int ProductionIncreaseAboveMin(unsigned int index) const
+    {
+        return nativeOptimVar.powerRampingIncreaseIndex[index];
+    }
+
+    int ProductionDecreaseAboveMin(unsigned int index) const
+    {
+        return nativeOptimVar.powerRampingDecreaseIndex[index];
+    }
+
     int NumberOfDispatchableUnits(unsigned int index) const
     {
         return nativeOptimVar.NumeroDeVariableDuNombreDeGroupesEnMarcheDuPalierThermique[index];
@@ -183,6 +193,16 @@ class ConstraintBuilder
                                               int offset = 0,
                                               int delta = 0);
 
+   ConstraintBuilder& ProductionDecreaseAboveMin(unsigned int index,
+                                                 double coeff,
+                                                 int offset = 0,
+                                                 int delta = 0);
+
+    ConstraintBuilder& ProductionIncreaseAboveMin(unsigned int index,
+                                                  double coeff,
+                                                  int offset = 0,
+                                                  int delta = 0);
+
     ConstraintBuilder& NumberOfDispatchableUnits(unsigned int index,
                                                  double coeff,
                                                  int offset = 0,

src/solver/optimisation/constraints/PowerOutputVariation.cpp

@@ -0,0 +1,47 @@
+#include "PowerOutputVariation.h"
+
+void PowerOutputVariation::add(int pays, int cluster, int clusterIndex, int pdt, bool Simulation)
+{
+    if (!Simulation)
+    {
+        const PALIERS_THERMIQUES& PaliersThermiquesDuPays = problemeHebdo->PaliersThermiquesDuPays[pays];
+        double maxUpwardPowerRampingRate = PaliersThermiquesDuPays.maxUpwardPowerRampingRate[clusterIndex];
+        double pminDUnGroupeDuPalierThermique = PaliersThermiquesDuPays.pminDUnGroupeDuPalierThermique[clusterIndex];
+        // constraint : P(t) - P(t-1) - l * M^+(t) - P^+ + P^- = 0
+        if (pdt > 0)
+        {
+            builder.updateHourWithinWeek(pdt)
+              .DispatchableProduction(cluster, 1.0)
+              .DispatchableProduction(cluster, -1.0, -1, problemeHebdo->NombreDePasDeTempsPourUneOptimisation)
+              .NumberStartingDispatchableUnits(cluster, -pminDUnGroupeDuPalierThermique)
+              .ProductionIncreaseAboveMin(cluster, -1.0)
+              .ProductionDecreaseAboveMin(cluster, 1.0)
+              .equalTo();
+        }
+        else
+        {
+            builder.updateHourWithinWeek(pdt)
+              .DispatchableProduction(cluster, 1.0)
+              .NumberStartingDispatchableUnits(cluster, -pminDUnGroupeDuPalierThermique)
+              .ProductionIncreaseAboveMin(cluster, -1.0)
+              .ProductionDecreaseAboveMin(cluster, 1.0)
+              .equalTo();
+        }
+        if (builder.NumberOfVariables() > 0)
+        {
+            ConstraintNamer namer(problemeHebdo->ProblemeAResoudre->NomDesContraintes);
+
+            namer.UpdateTimeStep(problemeHebdo->weekInTheYear * 168 + pdt);
+            namer.UpdateArea(problemeHebdo->NomsDesPays[pays]);
+
+            namer.ProductionOutputVariation(problemeHebdo->ProblemeAResoudre->NombreDeContraintes, PaliersThermiquesDuPays.NomsDesPaliersThermiques[clusterIndex]);
+        }
+        builder.build();
+    }
+    else
+    {
+        int add = (pdt == 0) ? 4 : 5;
+        problemeHebdo->NbTermesContraintesPourLesRampes += add;
+        problemeHebdo->ProblemeAResoudre->NombreDeContraintes++;
+    }
+}

src/solver/optimisation/constraints/PowerOutputVariation.h

@@ -0,0 +1,20 @@
+#pragma once
+#include "ConstraintBuilder.h"
+
+/*!
+ * represent 'RampingIncreaseRate' Constraint type
+ */
+class PowerOutputVariation : private ConstraintFactory
+{
+public:
+    using ConstraintFactory::ConstraintFactory;
+
+    /*!
+     * @brief Add variables to the constraint and update constraints Matrix
+     * @param pays : area
+     * @param cluster : global index of the cluster
+     * @param pdt : timestep
+     * @param Simulation : ---
+     */
+    void add(int pays, int cluster, int clusterIndex, int pdt, bool Simulation);
+};

src/solver/optimisation/constraints/RampingDecreaseRate.cpp

@@ -0,0 +1,48 @@
+#include "RampingDecreaseRate.h"
+
+void RampingDecreaseRate::add(int pays, int cluster, int clusterIndex, int pdt, bool Simulation)
+{
+    if (!Simulation)
+    {
+        const PALIERS_THERMIQUES& PaliersThermiquesDuPays = problemeHebdo->PaliersThermiquesDuPays[pays];
+        double maxDownwardPowerRampingRate = PaliersThermiquesDuPays.maxDownwardPowerRampingRate[clusterIndex];
+        double pmaxDUnGroupeDuPalierThermique = PaliersThermiquesDuPays.PmaxDUnGroupeDuPalierThermique[clusterIndex];
+        // constraint : P(t) - P(t-1) + R^- * M(t) + u * M^-(t) + u * M^--(t) > 0
+        if (pdt > 0)
+        {
+            builder.updateHourWithinWeek(pdt)
+              .DispatchableProduction(cluster, 1.0)
+              .DispatchableProduction(cluster, -1.0, -1, problemeHebdo->NombreDePasDeTempsPourUneOptimisation)
+              .NumberOfDispatchableUnits(cluster, maxDownwardPowerRampingRate)
+              .NumberStoppingDispatchableUnits(cluster, pmaxDUnGroupeDuPalierThermique)
+              .NumberBreakingDownDispatchableUnits(cluster, pmaxDUnGroupeDuPalierThermique)
+              .greaterThan();
+        }
+        else
+        {
+            builder.updateHourWithinWeek(pdt)
+              .DispatchableProduction(cluster, 1.0)
+              .NumberOfDispatchableUnits(cluster, maxDownwardPowerRampingRate)
+              .NumberStoppingDispatchableUnits(cluster, pmaxDUnGroupeDuPalierThermique)
+              .NumberBreakingDownDispatchableUnits(cluster, pmaxDUnGroupeDuPalierThermique)
+              .greaterThan();
+        }
+        if (builder.NumberOfVariables() > 0)
+        {
+            ConstraintNamer namer(problemeHebdo->ProblemeAResoudre->NomDesContraintes);
+
+            namer.UpdateTimeStep(problemeHebdo->weekInTheYear * 168 + pdt);
+            namer.UpdateArea(problemeHebdo->NomsDesPays[pays]);
+
+            namer.RampingDecreaseRate(problemeHebdo->ProblemeAResoudre->NombreDeContraintes, 
+                PaliersThermiquesDuPays.NomsDesPaliersThermiques[clusterIndex]);
+        }
+        builder.build();
+    }
+    else
+    {
+        int add = (pdt == 0) ? 4 : 5;
+        problemeHebdo->NbTermesContraintesPourLesRampes += add;
+        problemeHebdo->ProblemeAResoudre->NombreDeContraintes++;
+    }
+}

src/solver/optimisation/constraints/RampingDecreaseRate.h

@@ -0,0 +1,20 @@
+#pragma once
+#include "ConstraintBuilder.h"
+
+/*!
+ * represent 'RampingIncreaseRate' Constraint type
+ */
+class RampingDecreaseRate : private ConstraintFactory
+{
+public:
+    using ConstraintFactory::ConstraintFactory;
+
+    /*!
+     * @brief Add variables to the constraint and update constraints Matrix
+     * @param pays : area
+     * @param cluster : global index of the cluster
+     * @param pdt : timestep
+     * @param Simulation : ---
+     */
+    void add(int pays, int cluster, int clusterIndex, int pdt, bool Simulation);
+};

src/solver/optimisation/constraints/RampingIncreaseRate.cpp

@@ -0,0 +1,47 @@
+#include "RampingIncreaseRate.h"
+
+void RampingIncreaseRate::add(int pays, int cluster, int clusterIndex, int pdt, bool Simulation)
+{
+    if (!Simulation)
+    {
+        const PALIERS_THERMIQUES& PaliersThermiquesDuPays = problemeHebdo->PaliersThermiquesDuPays[pays];
+        double maxUpwardPowerRampingRate = PaliersThermiquesDuPays.maxUpwardPowerRampingRate[clusterIndex];
+        double pminDUnGroupeDuPalierThermique = PaliersThermiquesDuPays.pminDUnGroupeDuPalierThermique[clusterIndex];
+        // constraint : P(t) - P(t-1) - R^+ * M(t) - l * M^+(t) < 0
+        // logs.info() << "add up ramping constraint for cluster " << cluster << " pdt " << pdt;
+
+        if (pdt > 0)
+        {
+            builder.updateHourWithinWeek(pdt)
+              .DispatchableProduction(cluster, 1.0)
+              .DispatchableProduction(
+                cluster, -1.0, -1, problemeHebdo->NombreDePasDeTempsPourUneOptimisation)
+              .NumberOfDispatchableUnits(cluster, -maxUpwardPowerRampingRate)
+              .NumberStartingDispatchableUnits(cluster, -pminDUnGroupeDuPalierThermique)
+              .lessThan();
+        }
+        else
+        {
+            builder.updateHourWithinWeek(pdt)
+              .DispatchableProduction(cluster, 1.0)
+              .NumberOfDispatchableUnits(cluster, -maxUpwardPowerRampingRate)
+              .NumberStartingDispatchableUnits(cluster, -pminDUnGroupeDuPalierThermique)
+              .lessThan();
+        }
+        if (builder.NumberOfVariables() > 0)
+        {
+            ConstraintNamer namer(problemeHebdo->ProblemeAResoudre->NomDesContraintes);
+            namer.UpdateTimeStep(problemeHebdo->weekInTheYear * 168 + pdt);
+            namer.UpdateArea(problemeHebdo->NomsDesPays[pays]);
+            namer.RampingIncreaseRate(problemeHebdo->ProblemeAResoudre->NombreDeContraintes, 
+                PaliersThermiquesDuPays.NomsDesPaliersThermiques[clusterIndex]);
+        }
+        builder.build();
+    }
+    else
+    {
+        int add = (pdt == 0) ? 3 : 4;
+        problemeHebdo->NbTermesContraintesPourLesRampes += add;
+        problemeHebdo->ProblemeAResoudre->NombreDeContraintes++;
+    }
+}