Refactor constraint building (#1607)

* start

* add FictitiousLoad

* shortTermStorageLevel

* flowDissociation

* bindingConstraintHour

* ghost

* this is not only about constraint builder

* pay attention to the time step

* specific issue

* bug free

* smooth green?

* M&M

* Max Pumping ok?

* Area Hydro Level (_-_)

* finally

* basic

* tiny

* clean code

* add {PMin,PMax}DispatchableGeneration

* clean code

* ConsistenceNODU ok

* NbUnitsOutageLessThanNbUnitsStop

* NbDispUnitsMinBoundSinceMinUpTime

* MinDownTime

* setup Access Control

* remove comments

* remove comments

* treat major code smells

* treat major code smells 2

* more on major code smells

* update after review

* about objects lifetime

* (does not compile) cut into small pieces

* Revert "(does not compile) cut into small pieces"

This reverts commit 5a047c0.

* aesthetic

* remove comment

* smart move

* **** move

* shallow

* customize error exception

* test new variable construction design

* refactor

* update build for NODU

* clean deprecated code

* update build for NumberStoppingDispatchableUnits

* on NumberStartingDispatchableUnits

* on NumberBreakingDownDispatchableUnits

* new build NTCDirect

* missing file  NTCDirect

* Interco Direct/Indirect Cost

* ShortTermStorageInjection

* ShortTermStorageWithdrawal

* ShortTermStorageLevel

* HydProd

* HydProdDown

* HydProdUp

* Pumping

* HydroLevel

* Overflow

* FinalStorage

* fix

* LayerStorage

* Psitive/Negative UnsuppliedEnergy

* remove old code

* remove code

* fix

* move

* clarity matters thanks to @sylvlecl

* There is no need for repetition

* it makes sense:  Constraint->ConstraintFactory

* update

* remove redundant check

* rename things

* unused

* debug

* fiiiiiiiix

* add explicit keyword

* inspect unique_ptr

* resolve conflicts

* quick fix

* suggestion by @flomnes

Co-authored-by: Florian Omnès <[email protected]>

* add docs

* docs for new constraint classes

* update

* complete new classes documentation

* update

* update

* small changes

* small changes 2

* remove comment

* default value for constraint op

* clean code

* smell code

* smell code

* fix

* pretty

* typo

---------

Co-authored-by: Florian Omnès <[email protected]>

src/solver/optimisation/CMakeLists.txt

@@ -70,6 +70,56 @@ set(RTESOLVER_OPT
 		opt_period_string_generator_base.h
 		opt_rename_problem.h
 		opt_rename_problem.cpp
+
+		constraints/ConstraintBuilder.cpp
+		constraints/ConstraintBuilder.h
+		constraints/AreaBalance.h
+		constraints/AreaBalance.cpp
+		constraints/FictitiousLoad.h
+		constraints/FictitiousLoad.cpp
+		constraints/ShortTermStorageLevel.h
+		constraints/ShortTermStorageLevel.cpp
+		constraints/FlowDissociation.h
+		constraints/FlowDissociation.cpp
+		constraints/BindingConstraintHour.h
+		constraints/BindingConstraintHour.cpp
+		constraints/BindingConstraintDay.h
+		constraints/BindingConstraintDay.cpp
+		constraints/BindingConstraintWeek.h
+		constraints/BindingConstraintWeek.cpp
+		constraints/HydroPower.h
+		constraints/HydroPower.cpp
+		constraints/HydroPowerSmoothingUsingVariationSum.h
+		constraints/HydroPowerSmoothingUsingVariationSum.cpp
+		constraints/HydroPowerSmoothingUsingVariationMaxDown.h
+		constraints/HydroPowerSmoothingUsingVariationMaxDown.cpp
+		constraints/HydroPowerSmoothingUsingVariationMaxUp.h
+		constraints/HydroPowerSmoothingUsingVariationMaxUp.cpp
+		constraints/MinHydroPower.h
+		constraints/MinHydroPower.cpp
+		constraints/MaxHydroPower.h
+		constraints/MaxHydroPower.cpp
+		constraints/MaxPumping.h
+		constraints/MaxPumping.cpp
+		constraints/AreaHydroLevel.h
+		constraints/AreaHydroLevel.cpp
+		constraints/FinalStockEquivalent.h
+		constraints/FinalStockEquivalent.cpp
+		constraints/FinalStockExpression.h
+		constraints/FinalStockExpression.cpp
+		constraints/PMaxDispatchableGeneration.h
+		constraints/PMaxDispatchableGeneration.cpp
+		constraints/PMinDispatchableGeneration.h
+		constraints/PMinDispatchableGeneration.cpp
+		constraints/ConsistenceNumberOfDispatchableUnits.h
+		constraints/ConsistenceNumberOfDispatchableUnits.cpp
+		constraints/NbUnitsOutageLessThanNbUnitsStop.h
+		constraints/NbUnitsOutageLessThanNbUnitsStop.cpp
+		constraints/NbDispUnitsMinBoundSinceMinUpTime.h
+		constraints/NbDispUnitsMinBoundSinceMinUpTime.cpp
+		constraints/MinDownTime.h
+		constraints/MinDownTime.cpp
+
 	)
 
 

src/solver/optimisation/constraints/AreaBalance.cpp

@@ -0,0 +1,52 @@
+#include "AreaBalance.h"
+
+static void shortTermStorageBalance(const ::ShortTermStorage::AREA_INPUT& shortTermStorageInput,
+                                    ConstraintBuilder& constraintBuilder)
+{
+    for (const auto& storage : shortTermStorageInput)
+    {
+        unsigned index = storage.clusterGlobalIndex;
+        constraintBuilder.ShortTermStorageInjection(index, 1.0)
+        .ShortTermStorageWithdrawal(index, -1.0);
+    }
+}
+
+void AreaBalance::add(int pdt, int pays)
+{
+    CORRESPONDANCES_DES_CONTRAINTES& CorrespondanceCntNativesCntOptim
+      = problemeHebdo->CorrespondanceCntNativesCntOptim[pdt];
+    CorrespondanceCntNativesCntOptim.NumeroDeContrainteDesBilansPays[pays]
+      = problemeHebdo->ProblemeAResoudre->NombreDeContraintes;
+
+        ConstraintNamer namer(problemeHebdo->ProblemeAResoudre->NomDesContraintes);
+        namer.UpdateTimeStep(problemeHebdo->weekInTheYear * 168 + pdt);
+        namer.UpdateArea(problemeHebdo->NomsDesPays[pays]);
+        namer.AreaBalance(problemeHebdo->ProblemeAResoudre->NombreDeContraintes);
+
+        builder.updateHourWithinWeek(pdt);
+
+        int interco = problemeHebdo->IndexDebutIntercoOrigine[pays];
+        while (interco >= 0)
+        {
+        builder.NTCDirect(interco, 1.0);
+        interco = problemeHebdo->IndexSuivantIntercoOrigine[interco];
+        }
+
+    interco = problemeHebdo->IndexDebutIntercoExtremite[pays];
+    while (interco >= 0)
+    {
+        builder.NTCDirect(interco, -1.0);
+        interco = problemeHebdo->IndexSuivantIntercoExtremite[interco];
+    }
+
+    exportPaliers(*problemeHebdo, builder, pays);
+    builder.HydProd(pays, -1.0)
+      .Pumping(pays, 1.0)
+      .PositiveUnsuppliedEnergy(pays, -1.0)
+      .NegativeUnsuppliedEnergy(pays, 1.0);
+
+    shortTermStorageBalance(problemeHebdo->ShortTermStorage[pays], builder);
+
+    builder.equalTo();
+    builder.build();
+}

src/solver/optimisation/constraints/AreaBalance.h

@@ -0,0 +1,18 @@
+#pragma once
+#include "ConstraintBuilder.h"
+
+/*!
+ * represent 'Area Balance' constraint type
+ */
+class AreaBalance : private ConstraintFactory
+{
+public:
+    using ConstraintFactory::ConstraintFactory;
+
+    /*!
+     * @brief Add variables to the constraint and update constraints Matrix
+     * @param pdt : timestep
+     * @param pays : area
+     */
+    void add(int pdt, int pays);
+};

src/solver/optimisation/constraints/AreaHydroLevel.cpp

@@ -0,0 +1,32 @@
+#include "AreaHydroLevel.h"
+
+void AreaHydroLevel::add(int pays, int pdt)
+{
+    CORRESPONDANCES_DES_CONTRAINTES& CorrespondanceCntNativesCntOptim
+      = problemeHebdo->CorrespondanceCntNativesCntOptim[pdt];
+    CorrespondanceCntNativesCntOptim.NumeroDeContrainteDesNiveauxPays[pays]
+      = problemeHebdo->ProblemeAResoudre->NombreDeContraintes;
+    if (problemeHebdo->CaracteristiquesHydrauliques[pays].SuiviNiveauHoraire)
+    {
+        builder.updateHourWithinWeek(pdt).HydroLevel(pays, 1.0);
+        if (pdt > 0)
+        {
+            builder.updateHourWithinWeek(pdt - 1).HydroLevel(pays, -1.0);
+        }
+        ConstraintNamer namer(problemeHebdo->ProblemeAResoudre->NomDesContraintes);
+
+        namer.UpdateArea(problemeHebdo->NomsDesPays[pays]);
+        namer.UpdateTimeStep(problemeHebdo->weekInTheYear * 168 + pdt);
+        namer.AreaHydroLevel(problemeHebdo->ProblemeAResoudre->NombreDeContraintes);
+        CorrespondanceCntNativesCntOptim.NumeroDeContrainteDesNiveauxPays[pays]
+          = problemeHebdo->ProblemeAResoudre->NombreDeContraintes;
+        builder.updateHourWithinWeek(pdt)
+          .HydProd(pays, 1.0)
+          .Pumping(pays, -problemeHebdo->CaracteristiquesHydrauliques[pays].PumpingRatio)
+          .Overflow(pays, 1.)
+          .equalTo()
+          .build();
+    }
+    else
+        CorrespondanceCntNativesCntOptim.NumeroDeContrainteDesNiveauxPays[pays] = -1;
+}

src/solver/optimisation/constraints/AreaHydroLevel.h

@@ -0,0 +1,18 @@
+#pragma once
+#include "ConstraintBuilder.h"
+
+/*!
+ * represent 'Area Hydraulic Level' constraint type
+ */
+class AreaHydroLevel : private ConstraintFactory
+{
+    public:
+    using ConstraintFactory::ConstraintFactory;
+
+    /*!
+     * @brief Add variables to the constraint and update constraints Matrix
+     * @param pdt : timestep
+     * @param pays : area
+     */
+    void add(int pays, int pdt);
+};

src/solver/optimisation/constraints/BindingConstraintDay.cpp

@@ -0,0 +1,69 @@
+#include "BindingConstraintDay.h"
+
+void BindingConstraintDay::add(int cntCouplante)
+{
+    const CONTRAINTES_COUPLANTES& MatriceDesContraintesCouplantes
+      = problemeHebdo->MatriceDesContraintesCouplantes[cntCouplante];
+    if (MatriceDesContraintesCouplantes.TypeDeContrainteCouplante != CONTRAINTE_JOURNALIERE)
+        return;
+
+    const int nbInterco
+      = MatriceDesContraintesCouplantes.NombreDInterconnexionsDansLaContrainteCouplante;
+    const int nbClusters
+      = MatriceDesContraintesCouplantes.NombreDePaliersDispatchDansLaContrainteCouplante;
+
+    const int NombreDePasDeTempsPourUneOptimisation
+      = problemeHebdo->NombreDePasDeTempsPourUneOptimisation;
+    const int NombreDePasDeTempsDUneJournee = problemeHebdo->NombreDePasDeTempsDUneJournee;
+    int pdtDebut = 0;
+    while (pdtDebut < NombreDePasDeTempsPourUneOptimisation)
+    {
+        int jour = problemeHebdo->NumeroDeJourDuPasDeTemps[pdtDebut];
+        CORRESPONDANCES_DES_CONTRAINTES_JOURNALIERES& CorrespondanceCntNativesCntOptimJournalieres
+          = problemeHebdo->CorrespondanceCntNativesCntOptimJournalieres[jour];
+
+        for (int index = 0; index < nbInterco; index++)
+        {
+            int interco = MatriceDesContraintesCouplantes.NumeroDeLInterconnexion[index];
+            double poids = MatriceDesContraintesCouplantes.PoidsDeLInterconnexion[index];
+            int offset = MatriceDesContraintesCouplantes.OffsetTemporelSurLInterco[index];
+
+            for (int pdt = pdtDebut; pdt < pdtDebut + NombreDePasDeTempsDUneJournee; pdt++)
+            {
+                builder.updateHourWithinWeek(pdt).NTCDirect(
+                  interco, poids, offset, problemeHebdo->NombreDePasDeTemps);
+            }
+        }
+
+        for (int index = 0; index < nbClusters; index++)
+        {
+            int pays = MatriceDesContraintesCouplantes.PaysDuPalierDispatch[index];
+            const PALIERS_THERMIQUES& PaliersThermiquesDuPays
+              = problemeHebdo->PaliersThermiquesDuPays[pays];
+            const int palier
+              = PaliersThermiquesDuPays.NumeroDuPalierDansLEnsembleDesPaliersThermiques
+                  [MatriceDesContraintesCouplantes.NumeroDuPalierDispatch[index]];
+            double poids = MatriceDesContraintesCouplantes.PoidsDuPalierDispatch[index];
+            int offset = MatriceDesContraintesCouplantes.OffsetTemporelSurLePalierDispatch[index];
+
+            for (int pdt = pdtDebut; pdt < pdtDebut + NombreDePasDeTempsDUneJournee; pdt++)
+            {
+                builder.updateHourWithinWeek(pdt).DispatchableProduction(
+                  palier, poids, offset, problemeHebdo->NombreDePasDeTemps);
+            }
+        }
+
+        CorrespondanceCntNativesCntOptimJournalieres
+          .NumeroDeContrainteDesContraintesCouplantes[cntCouplante]
+          = problemeHebdo->ProblemeAResoudre->NombreDeContraintes;
+
+        char op = MatriceDesContraintesCouplantes.SensDeLaContrainteCouplante;
+        builder.SetOperator(op);
+        ConstraintNamer namer(problemeHebdo->ProblemeAResoudre->NomDesContraintes);
+        namer.UpdateTimeStep(jour);
+        namer.BindingConstraintDay(problemeHebdo->ProblemeAResoudre->NombreDeContraintes,
+                                   MatriceDesContraintesCouplantes.NomDeLaContrainteCouplante);
+        builder.build();
+        pdtDebut += problemeHebdo->NombreDePasDeTempsDUneJournee;
+    }
+}

src/solver/optimisation/constraints/BindingConstraintDay.h

@@ -0,0 +1,17 @@
+#pragma once
+#include "ConstraintBuilder.h"
+
+/*!
+ * represent 'Daily Binding Constraint' type
+ */
+class BindingConstraintDay : private ConstraintFactory
+{
+public:
+    using ConstraintFactory::ConstraintFactory;
+
+    /*!
+     * @brief Add variables to the constraint and update constraints Matrix
+     * @param cntCouplante : the binding constraint number
+     */
+    void add(int cntCouplante);
+};

src/solver/optimisation/constraints/BindingConstraintHour.cpp

@@ -0,0 +1,53 @@
+#include "BindingConstraintHour.h"
+
+void BindingConstraintHour::add(int pdt, int cntCouplante)
+{
+    const CONTRAINTES_COUPLANTES& MatriceDesContraintesCouplantes
+      = problemeHebdo->MatriceDesContraintesCouplantes[cntCouplante];
+    CORRESPONDANCES_DES_CONTRAINTES& CorrespondanceCntNativesCntOptim
+      = problemeHebdo->CorrespondanceCntNativesCntOptim[pdt];
+    CorrespondanceCntNativesCntOptim.NumeroDeContrainteDesContraintesCouplantes[cntCouplante]
+      = problemeHebdo->ProblemeAResoudre->NombreDeContraintes;
+
+    if (MatriceDesContraintesCouplantes.TypeDeContrainteCouplante != CONTRAINTE_HORAIRE)
+        return;
+
+    builder.updateHourWithinWeek(pdt);
+    // Links
+    const int nbInterco
+      = MatriceDesContraintesCouplantes.NombreDInterconnexionsDansLaContrainteCouplante;
+    for (int index = 0; index < nbInterco; index++)
+    {
+        const int interco = MatriceDesContraintesCouplantes.NumeroDeLInterconnexion[index];
+        const double poids = MatriceDesContraintesCouplantes.PoidsDeLInterconnexion[index];
+        const int offset = MatriceDesContraintesCouplantes.OffsetTemporelSurLInterco[index];
+
+        builder.updateHourWithinWeek(pdt).NTCDirect(
+          interco, poids, offset, problemeHebdo->NombreDePasDeTemps);
+    }
+
+    // Thermal clusters
+    const int nbClusters
+      = MatriceDesContraintesCouplantes.NombreDePaliersDispatchDansLaContrainteCouplante;
+    for (int index = 0; index < nbClusters; index++)
+    {
+        const int pays = MatriceDesContraintesCouplantes.PaysDuPalierDispatch[index];
+        const PALIERS_THERMIQUES& PaliersThermiquesDuPays
+          = problemeHebdo->PaliersThermiquesDuPays[pays];
+        const int palier = PaliersThermiquesDuPays.NumeroDuPalierDansLEnsembleDesPaliersThermiques
+                             [MatriceDesContraintesCouplantes.NumeroDuPalierDispatch[index]];
+        const double poids = MatriceDesContraintesCouplantes.PoidsDuPalierDispatch[index];
+        const int offset = MatriceDesContraintesCouplantes.OffsetTemporelSurLePalierDispatch[index];
+
+        builder.updateHourWithinWeek(pdt).DispatchableProduction(
+          palier, poids, offset, problemeHebdo->NombreDePasDeTemps);
+    }
+
+    char op = MatriceDesContraintesCouplantes.SensDeLaContrainteCouplante;
+    builder.SetOperator(op);
+    ConstraintNamer namer(problemeHebdo->ProblemeAResoudre->NomDesContraintes);
+    namer.UpdateTimeStep(problemeHebdo->weekInTheYear * 168 + pdt);
+    namer.BindingConstraintHour(problemeHebdo->ProblemeAResoudre->NombreDeContraintes,
+                                MatriceDesContraintesCouplantes.NomDeLaContrainteCouplante);
+    builder.build();
+}

src/solver/optimisation/constraints/BindingConstraintHour.h

@@ -0,0 +1,18 @@
+#pragma once
+#include "ConstraintBuilder.h"
+
+/*!
+ * represent 'Hourly Binding Constraint' type
+ */
+class BindingConstraintHour : private ConstraintFactory
+{
+    public:
+    using ConstraintFactory::ConstraintFactory;
+
+    /*!
+     * @brief Add variables to the constraint and update constraints Matrix
+     * @param pdt : timestep
+     * @param cntCouplante : the binding constraint number
+     */
+    void add(int pdt, int cntCouplante);
+};