add IEEE G1 model and expand flexibility of dynamic Generator model t…

src/fileInput/gridDynReadDYR.cpp

@@ -5,6 +5,7 @@
  */
 
 #include "ReaderInfo.h"
+#include "core/CoreExceptions.h"
 #include "core/CoreObject.h"
 #include "core/ObjectFactory.hpp"
 #include "core/coreDefinitions.hpp"
@@ -16,8 +17,12 @@
 #include "griddyn/Generator.h"
 #include "griddyn/Governor.h"
 #include "griddyn/GridBus.h"
+#include "griddyn/generators/DynamicGenerator.h"
+#include "griddyn/governors/GovernorIeeeG1.h"
+#include <cmath>
 #include <fstream>
 #include <iostream>
+#include <memory>
 #include <string>
 
 namespace griddyn {
@@ -27,6 +32,7 @@ namespace {
     void loadESST3A(CoreObject* parentObject, stringVec& tokens);
     void loadEXST1(CoreObject* parentObject, stringVec& tokens);
     void loadTGOV1(CoreObject* parentObject, stringVec& tokens);
+    void loadIEEEG1(CoreObject* parentObject, stringVec& tokens);
     void loadEXDC2(CoreObject* parentObject, stringVec& tokens);
     void loadSEXS(CoreObject* parentObject, stringVec& tokens);
 }  // namespace
@@ -79,6 +85,8 @@ void loadDyr(CoreObject* parentObject,
             loadEXDC2(parentObject, lineTokens);
         } else if (type == "'TGOV1'") {
             loadTGOV1(parentObject, lineTokens);
+        } else if (type == "'IEEEG1'") {
+            loadIEEEG1(parentObject, lineTokens);
         } else if (type == "'SEXS'") {
             loadSEXS(parentObject, lineTokens);
         } else {
@@ -291,6 +299,85 @@ namespace {
 
         gen->add(governorModel);
     }
+
+    void loadIEEEG1(CoreObject* parentObject, stringVec& tokens)
+    {
+        if (tokens.size() != 25U) {
+            throw InvalidParameterValue("IEEEG1 DYR record must contain 25 fields");
+        }
+        const int busId = std::stoi(tokens[0]);
+        const auto* bus = static_cast<GridBus*>(parentObject->findByUserID("bus", busId));
+        const int genId = std::stoi(tokens[2]);
+        if ((bus == nullptr) || (genId <= 0)) {
+            throw InvalidParameterValue("IEEEG1 primary generator identity");
+        }
+        auto* primary = dynamic_cast<DynamicGenerator*>(bus->getGen(genId - 1));
+        if (primary == nullptr) {
+            throw InvalidParameterValue("IEEEG1 requires a dynamic primary generator");
+        }
+
+        const auto params = gmlc::utilities::str2vector(tokens, kNullVal);
+        const int secondBusId = static_cast<int>(params[3]);
+        DynamicGenerator* secondary = nullptr;
+        if (secondBusId != 0) {
+            const auto* secondBus =
+                dynamic_cast<GridBus*>(parentObject->findByUserID("bus", secondBusId));
+            const int secondGenId = static_cast<int>(params[4]);
+            if ((secondBus == nullptr) || (secondGenId <= 0)) {
+                throw InvalidParameterValue("IEEEG1 secondary generator identity");
+            }
+            secondary = dynamic_cast<DynamicGenerator*>(secondBus->getGen(secondGenId - 1));
+            if ((secondary == nullptr) || (secondary == primary)) {
+                throw InvalidParameterValue(
+                    "IEEEG1 requires a distinct dynamic secondary generator");
+            }
+        } else {
+            constexpr double zeroTolerance = 1e-12;
+            if ((std::abs(params[15]) > zeroTolerance) || (std::abs(params[18]) > zeroTolerance) ||
+                (std::abs(params[21]) > zeroTolerance) || (std::abs(params[24]) > zeroTolerance)) {
+                throw InvalidParameterValue(
+                    "single-generator IEEEG1 requires K2, K4, K6, and K8 to be zero");
+            }
+        }
+
+        auto cof = CoreObjectFactory::instance();
+        std::unique_ptr<governors::GovernorIeeeG1> governor(
+            dynamic_cast<governors::GovernorIeeeG1*>(cof->createObject("governor", "ieeeg1")));
+        if (governor == nullptr) {
+            throw InvalidParameterValue("IEEEG1 factory registration");
+        }
+
+        // Exact frozen ANDES psse-dyr.yaml order after BUS and ID:
+        // BUS2, ID2, K, T1, T2, T3, UO, UC, PMAX, PMIN,
+        // T4, K1, K2, T5, K3, K4, T6, K5, K6, T7, K7, K8.
+        governor->set("k", params[5]);
+        governor->set("t1", params[6]);
+        governor->set("t2", params[7]);
+        governor->set("t3", params[8]);
+        governor->set("uo", params[9]);
+        governor->set("uc", params[10]);
+        governor->set("pmax", params[11]);
+        governor->set("pmin", params[12]);
+        governor->set("t4", params[13]);
+        governor->set("k1", params[14]);
+        governor->set("k2", params[15]);
+        governor->set("t5", params[16]);
+        governor->set("k3", params[17]);
+        governor->set("k4", params[18]);
+        governor->set("t6", params[19]);
+        governor->set("k5", params[20]);
+        governor->set("k6", params[21]);
+        governor->set("t7", params[22]);
+        governor->set("k7", params[23]);
+        governor->set("k8", params[24]);
+
+        auto* governorPointer = governor.release();
+        primary->add(governorPointer);
+        if (secondary != nullptr) {
+            secondary->setMechanicalPowerSource(governorPointer,
+                                                governors::GovernorIeeeG1::lpOutput);
+        }
+    }
 }  // namespace
 
 }  // namespace griddyn

src/griddyn/CMakeLists.txt

@@ -107,6 +107,7 @@ set(gov_headers
     Governor.h
     governors/GovernorTypes.h
     governors/GovernorHydro.h
+    governors/GovernorIeeeG1.h
     governors/GovernorIeeeSimple.h
     governors/GovernorReheat.h
     governors/GovernorSteamNR.h
@@ -118,6 +119,7 @@ set(gov_sources
     ${gov_headers}
     governors/Governor.cpp
     governors/GovernorHydro.cpp
+    governors/GovernorIeeeG1.cpp
     governors/GovernorIeeeSimple.cpp
     governors/GovernorReheat.cpp
     governors/GovernorSteamNR.cpp

src/griddyn/GridSubModel.cpp

@@ -65,6 +65,11 @@ void GridSubModel::dynInitializeB(const IOdata& inputs,
     }
 }
 
+bool GridSubModel::setOutputInitializationTarget(index_t /*outputIndex*/, double /*target*/)
+{
+    return false;
+}
+
 double GridSubModel::get(std::string_view param, units::unit unitType) const
 {
     auto fptr = getObjectFunction(this, param);

src/griddyn/GridSubModel.h

@@ -36,6 +36,14 @@ class GridSubModel: public GridComponent {
                                 const IOdata& desiredOutput,
                                 IOdata& fieldSet) override final;
 
+    /** Supply an initialization target for one output of a shared submodel.
+     *
+     * Most submodels do not need this hook and return false. Multi-output
+     * controllers such as IEEEG1 use it to collect the initialized mechanical
+     * powers of generators which consume outputs owned by another generator.
+     */
+    virtual bool setOutputInitializationTarget(index_t outputIndex, double target);
+
     virtual double get(std::string_view param,
                        units::unit unitType = units::defunit) const override;
 };

src/griddyn/generators/DynamicGenerator.cpp

@@ -74,6 +74,11 @@ CoreObject* DynamicGenerator::clone(CoreObject* obj) const
     if (gen == nullptr) {
         return obj;
     }
+    gen->mechanicalPowerSourceExplicit = mechanicalPowerSourceExplicit;
+    gen->mechanicalPowerOutput = mechanicalPowerOutput;
+    gen->mechanicalPowerSourceName = mechanicalPowerSourceName;
+    gen->mechanicalPowerSource =
+        (mechanicalPowerSourceExplicit && (mechanicalPowerSource == gov)) ? gen->gov : nullptr;
     return gen;
 }
 namespace {
@@ -224,6 +229,7 @@ void DynamicGenerator::dynObjectInitializeA(CoreTime time0, std::uint32_t flags)
     if (genModel == nullptr) {
         add(new GenModel());
     }
+    resolveMechanicalPowerSource();
     if (gov != nullptr) {
         if (!genModel->checkFlag(GenModel::GenModelFlags::INTERNAL_FREQUENCY_CALCULATION)) {
             opFlags.set(USES_BUS_FREQUENCY);
@@ -278,6 +284,10 @@ void DynamicGenerator::dynObjectInitializeB(const IOdata& inputs,
     //  genModel->guessState (prevTime, m_state.data (), m_dstate_dt.data (), cLocalbSolverMode);
 
     Pset = m_Pmech / scale;
+    if (mechanicalPowerSourceExplicit && (mechanicalPowerSource != nullptr) &&
+        (mechanicalPowerSource != gov)) {
+        mechanicalPowerSource->setOutputInitializationTarget(mechanicalPowerOutput, m_Pmech);
+    }
     if (isoc != nullptr) {
         Pset -= isoc->getOutput();
     }
@@ -426,7 +436,12 @@ void DynamicGenerator::add(GridSubModel* obj)
             obj->set("xs", m_Xs);
         }
     } else if (dynamic_cast<Governor*>(obj) != nullptr) {
+        const bool reconnectExplicitSource =
+            mechanicalPowerSourceExplicit && (mechanicalPowerSource == gov);
         gov = static_cast<Governor*>(replaceModel(obj, gov, GOVERNOR_LOC));
+        if (reconnectExplicitSource) {
+            setMechanicalPowerSource(gov, mechanicalPowerOutput);
+        }
         // mesh up the Pmax and Pmin giving priority to the new gov
         const double govpmax = gov->get("pmax");
         const double govpmin = gov->get("pmin");
@@ -462,6 +477,85 @@ void DynamicGenerator::add(GridSubModel* obj)
     }
 }
 
+void DynamicGenerator::setMechanicalPowerSource(GridSubModel* source, index_t outputIndex)
+{
+    if (source == nullptr) {
+        clearMechanicalPowerSource();
+        return;
+    }
+    if ((outputIndex < 0) || (outputIndex >= source->numOutputs())) {
+        throw InvalidParameterValue("mechanical power output");
+    }
+
+    mechanicalPowerSource = source;
+    mechanicalPowerOutput = outputIndex;
+    mechanicalPowerSourceExplicit = true;
+    mechanicalPowerSourceName =
+        (source->getParent() != nullptr) ? fullObjectName(source) : source->getName();
+    subInputs.seqID = 0;
+    subInputLocs.seqID = 0;
+}
+
+void DynamicGenerator::setMechanicalPowerSource(std::string_view sourceName, index_t outputIndex)
+{
+    if ((sourceName.empty()) || (sourceName == "default") || (sourceName == "local")) {
+        clearMechanicalPowerSource();
+        return;
+    }
+    if (outputIndex < 0) {
+        throw InvalidParameterValue("mechanical power output");
+    }
+
+    mechanicalPowerSource = nullptr;
+    mechanicalPowerOutput = outputIndex;
+    mechanicalPowerSourceExplicit = true;
+    mechanicalPowerSourceName = sourceName;
+    subInputs.seqID = 0;
+    subInputLocs.seqID = 0;
+}
+
+void DynamicGenerator::clearMechanicalPowerSource()
+{
+    mechanicalPowerSource = nullptr;
+    mechanicalPowerOutput = 0;
+    mechanicalPowerSourceName.clear();
+    mechanicalPowerSourceExplicit = false;
+    subInputs.seqID = 0;
+    subInputLocs.seqID = 0;
+}
+
+GridSubModel* DynamicGenerator::getMechanicalPowerSource() const
+{
+    return mechanicalPowerSourceExplicit ? mechanicalPowerSource : gov;
+}
+
+index_t DynamicGenerator::getMechanicalPowerOutput() const
+{
+    return mechanicalPowerSourceExplicit ? mechanicalPowerOutput : 0;
+}
+
+bool DynamicGenerator::hasExplicitMechanicalPowerSource() const
+{
+    return mechanicalPowerSourceExplicit;
+}
+
+void DynamicGenerator::resolveMechanicalPowerSource()
+{
+    if (!mechanicalPowerSourceExplicit || (mechanicalPowerSource != nullptr)) {
+        return;
+    }
+
+    auto* source =
+        dynamic_cast<GridSubModel*>(locateObject(mechanicalPowerSourceName, getRoot(), false));
+    if (source == nullptr) {
+        throw InvalidParameterValue("mechanical power source '" + mechanicalPowerSourceName + "'");
+    }
+    if (mechanicalPowerOutput >= source->numOutputs()) {
+        throw InvalidParameterValue("mechanical power output");
+    }
+    mechanicalPowerSource = source;
+}
+
 GridSubModel* DynamicGenerator::replaceModel(GridSubModel* newObject,
                                              GridSubModel* oldObject,
                                              index_t newIndex)
@@ -485,6 +579,9 @@ void DynamicGenerator::set(std::string_view param, std::string_view val)
             throw(InvalidParameterValue(val));
         }
         buildDynModel(dmodel);
+    } else if ((param == "mechanical_power_source") || (param == "mechanicalpowersource") ||
+               (param == "pmech_source") || (param == "pmechsource")) {
+        setMechanicalPowerSource(val, mechanicalPowerOutput);
     } else {
         try {
             Generator::set(param, val);
@@ -519,7 +616,10 @@ void DynamicGenerator::timestep(CoreTime time, const IOdata& inputs, const Solve
 
         if ((gov != nullptr) && (gov->isEnabled())) {
             gov->timestep(time, {omega, Pset / scale}, sMode);
-            m_Pmech = gov->getOutput();
+        }
+        auto* pmechSource = getMechanicalPowerSource();
+        if ((pmechSource != nullptr) && (pmechSource->isEnabled())) {
+            m_Pmech = pmechSource->getOutput(getMechanicalPowerOutput());
         }
 
         if ((pss != nullptr) && (pss->isEnabled())) {
@@ -670,6 +770,19 @@ void DynamicGenerator::set(std::string_view param, double val, unit unitType)
         }
     } else if (param == "eft") {
         m_Eft = val;
+    } else if ((param == "mechanical_power_output") || (param == "mechanicalpoweroutput") ||
+               (param == "pmech_output") || (param == "pmechoutput")) {
+        const auto outputIndex = static_cast<index_t>(val);
+        if ((val < 0.0) || (static_cast<double>(outputIndex) != val)) {
+            throw InvalidParameterValue("mechanical power output");
+        }
+        if ((mechanicalPowerSource != nullptr) &&
+            (outputIndex >= mechanicalPowerSource->numOutputs())) {
+            throw InvalidParameterValue("mechanical power output");
+        }
+        mechanicalPowerOutput = outputIndex;
+        subInputs.seqID = 0;
+        subInputLocs.seqID = 0;
     } else if (param == "vref") {
         if (ext != nullptr) {
             ext->set(param, val, unitType);
@@ -1128,8 +1241,11 @@ void DynamicGenerator::generateSubModelInputs(const IOdata& inputs,
     subInputs.inputs[GOVERNOR_LOC][govpSetInLocation] = pcontrol * scale;
 
     double pmech = pcontrol * scale;
-    if ((gov != nullptr) && (gov->isEnabled())) {
-        pmech = gov->getOutput(subInputs.inputs[GOVERNOR_LOC], stateDataValue, sMode, 0);
+    auto* pmechSource = getMechanicalPowerSource();
+    if ((pmechSource != nullptr) && (pmechSource->isEnabled())) {
+        const auto& sourceInputs = (pmechSource == gov) ? subInputs.inputs[GOVERNOR_LOC] : noInputs;
+        pmech =
+            pmechSource->getOutput(sourceInputs, stateDataValue, sMode, getMechanicalPowerOutput());
     }
     if (std::abs(pmech) > 1e25) {
         pmech = 0.0;
@@ -1202,8 +1318,12 @@ void DynamicGenerator::generateSubModelInputLocs(const IOlocs& inputLocs,
             subInputLocs.inputLocs[GOVERNOR_LOC][govOmegaInLocation] = floc;
         }
         subInputLocs.inputLocs[GOVERNOR_LOC][govpSetInLocation] = pSetLocation(sMode);
+    }
+
+    auto* pmechSource = getMechanicalPowerSource();
+    if ((pmechSource != nullptr) && (pmechSource->isEnabled())) {
         subInputLocs.inputLocs[GEN_MODEL_LOC][genModelPmechInLocation] =
-            gov->getOutputLoc(sMode, 0);
+            pmechSource->getOutputLoc(sMode, getMechanicalPowerOutput());
     } else {
         subInputLocs.inputLocs[GEN_MODEL_LOC][genModelPmechInLocation] = pSetLocation(sMode);
     }