diff --git a/docs/developer-guide/andes-compatibility.md b/docs/developer-guide/andes-compatibility.md index 6466511b..7fb4f58a 100644 --- a/docs/developer-guide/andes-compatibility.md +++ b/docs/developer-guide/andes-compatibility.md @@ -80,7 +80,8 @@ adapters still have attachment, parameter, equation, and validation gaps. | `ESAC6A`, `SCRX`, `EXPIC1` | `SEXS` | `ExciterSEXS` | **Planned compatibility approximations.** ANDES marks these conversions as TODO/approximate and currently discards most source parameters. Reproduce this only as an explicit compatibility mode and emit a diagnostic; do not describe it as exact model support. | | `TGOV1` | `TGOV1` | `GovernorTgov1` | **Implemented.** Matches the ANDES v2.0.0 equations and PSS/e schema `R, T1, VMAX, VMIN, T2, T3, Dt`; DYR attachment, initialization, limiter, analytic-Jacobian, and isolated speed-step trajectory regressions are covered. | | `HYGOV` | `HYGOV` | `GovernorHydro` is a candidate | **Planned.** Compare equations, water-column dynamics, gate/rate limits, and parameter units before mapping it. | -| `IEESGO`, `IEEEG1` | Same named ANDES models | `GovernorReheat` and steam-governor classes are only candidates | **Planned.** No exact equivalence has been established; audit block diagrams before choosing reuse versus new models. | +| `IEESGO` | `IEESGO` | `GovernorReheat` and steam-governor classes are only candidates | **Planned.** No exact equivalence has been established; audit block diagrams before choosing reuse versus a new model. | +| `IEEEG1` | `IEEEG1` | `GovernorIeeeG1` | **Implemented.** Frozen-ANDES equations and DYR order, initialization and perturbed equations, valve rate/position limits, analytic Jacobian, one- and two-machine attachment, and mixed synchronous-machine coupling are covered. Native ANDES import, alphanumeric DYR IDs, unequal machine bases, and a captured disturbed trajectory remain open. | | `GAST` | `GAST` | None exact | **No direct analogue.** Add a gas-turbine governor implementation and DYR adapter. | | `GGOV1` | `TGOV1` | `GovernorTgov1` | **Planned compatibility approximation.** ANDES currently retains only `R` when converting this record. If reproduced, emit a diagnostic and keep exact GGOV1 support as a separate task. | | `IEEEST`, `ST2CUT` | Same named ANDES models | Base `Stabilizer` only | **No direct analogue.** The base class is not a usable PSS implementation; add the two models and their generator/exciter signal connections. | @@ -175,7 +176,7 @@ seconds, and continuing through 2.0 seconds. The five controller chains are: | 2 | Robust RAW/DYR identity and base handling, IEEE 14 power-flow parity, and a controller-free GENROU trajectory | Planned | PR 1 | | 3 | Generator/controller signal plumbing, PSS-to-exciter routing, and validated `TGOV1` | Implemented | PR 2 | | 4 | Complete `ESST3A` and `EXST1` models plus DYR adapters | Models/adapters implemented; trajectories open | PR 3 | -| 5 | Complete `IEEEG1`, `ST2CUT`, and `IEEEST` models plus DYR adapters | Planned | PR 4 | +| 5 | Complete `IEEEG1`, `ST2CUT`, and `IEEEST` models plus DYR adapters | `IEEEG1` implemented; stabilizers planned | PR 4 | | 6 | DYR `Toggle`, complete IEEE 14 initialization/equilibrium, and the two-second trajectory regression | Planned | PR 5 | PRs must merge in this order. Each PR must pass without relying on production @@ -437,6 +438,34 @@ Required models: - `IEEEST`, used at bus 3: implement selectable input, second-order filters, two lead-lag stages, gain, washout/lag, output limiting, and voltage gating. +#### IEEEG1 implementation and parameter contract + +`GovernorIeeeG1` follows frozen ANDES v2.0.0 +`andes/models/governor/ieeeg1.py`. Its DYR adapter uses the exact frozen +`psse-dyr.yaml` order `BUS, ID, BUS2, ID2, K, T1, T2, T3, UO, UC, PMAX, +PMIN, T4, K1, K2, T5, K3, K4, T6, K5, K6, T7, K7, K8`. The eight power +coefficients are normalized by their sum, as in ANDES. `T1=0` and zero +turbine-stage time constants are exact algebraic bypasses and do not allocate +unnecessary differential states. + +The valve rate is limited to `[UC, UO]`; valve-position anti-windup holds at +`[PMIN, PMAX]` only while requested motion points farther into the active +limit. No equation discrepancy from frozen ANDES is known. ANDES marks +governor scheduling as unsupported, so GridDyn likewise holds the initialized +reference and does not claim a dynamic `paux`/dispatch input. + +With `BUS2=0`, the adapter uses the high-pressure output on the primary +generator and requires `K2`, `K4`, `K6`, and `K8` to be zero. With a second +machine, the governor remains owned and evaluated once by the primary +generator, while the generic indexed mechanical-power connection routes its +low-pressure output to the secondary. The secondary may use a different +synchronous-generator model class and may retain an unrelated local governor; +the explicit mechanical source selects which output drives its machine. +`GenModel4` plus `GenModel6` is covered. Exact two-machine operation currently +requires equal machine MBASE values; base conversion for unequal machines, +alphanumeric DYR machine IDs, native ANDES input, and a captured disturbed +trajectory remain open. + Do not map `IEEEG1` to `GovernorReheat` or either stabilizer to the base `Stabilizer` unless an equation-by-equation audit proves exact equivalence. The actual IEEE 14 `ST2CUT` records use rotor-speed input, while its `IEEEST` @@ -633,29 +662,30 @@ can be split and checked off as soon as models in the same family diverge. A model is not compatible merely because a similarly named GridDyn model exists: it also needs native-input mapping, initialization, and a trajectory test. -| ANDES model(s) | GridDyn mapping / next action | Status | -| --------------------------------------------------------------------------------------------------------------------- | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ------------------ | -| `GENCLS` | `GenModelClassical`; parameter and trajectory comparison required. | Partial | -| `GENROU` | `GenModelGENROU` equations, DYR mapping, and five-machine initialization match ANDES references; fix exact machine identity/base handling, native ANDES import, and trajectory parity. | Partial | -| `TGOV1` | `GovernorTgov1`; ANDES-equation, DYR-order, limiter, initialization, Jacobian, and isolated trajectory regressions are complete. | Implemented | -| `EXDC2` | `ExciterDC2A` candidate; add omitted transducer/switch/saturation behavior and prove equation equivalence. | Partial | -| `ZIP`, `FLoad` | GridDyn static/dynamic load models; identify exact parameter and frequency-response equivalence. | Partial | -| `Motor3`, `Motor5` | `MotorLoad3`, `MotorLoad5`; parameter mapping and trajectory comparisons required. | Partial | -| `ACE`, `ACEc`, `COI` | Map area-control and center-of-inertia services. | Untriaged | -| `PLBVFU1`, `IEEEVC` | No direct voltage-compensator/playback analogue; define exciter input and playback interfaces. | No direct analogue | -| `TG2`, `TGOV1DB`, `TGOV1N`, `TGOV1NDB`, `IEEEG1`, `IEESGO`, `GAST`, `HYGOV`, `HYGOVDB`, `HYGOV4` | Existing hydro/reheat/steam classes are candidates only; each needs an equation audit or a new model. | Planned | -| `SHAFT5` | Multi-mass shaft model; map states and mechanical interfaces. | Untriaged | -| `ESDC1A`, `ESDC2A`, `SEXS`, `IEEET1`, `IEEET3` | Existing DC/IEEE/SEXS exciters are candidates; complete model-specific audits and DYR/native mappings. | Partial | -| `ESST3A` | `ExciterESST3A`; exact GENROU path plus documented reduced-order synchronous-machine signal approximations. | Implemented | -| `EXST1` | `ExciterEXST1`; positive-time-constant equations and DYR mapping, with the documented corrected regulator-output limiter; zero-time-constant blocks, native import, and a captured trajectory remain. | Implemented | -| `IEEEX1`, `EXAC1`, `EXAC2`, `EXAC4`, `ESST4B`, `AC8B`, `ESAC1A`, `ESST1A`, `ESAC5A` | No exact named GridDyn implementations; implement individually rather than mapping by family name. | No direct analogue | -| `IEEEST`, `ST2CUT` | Base `Stabilizer` is not a functional equivalent; implement both models and signal connections. | No direct analogue | -| `BusFreq`, `BusROCOF`, `PMU`, `PLL1`, `PLL2`, `FreqDiv` | Measurement and frequency-estimation models. | Untriaged | -| `REGCA1`, `REGCP1`, `REECA1`, `REECA1E`, `REECA1G`, `REECB1`, `REPCA1`, `REGCV1`, `REGCV2`, `REGF1`, `REGF2`, `REGF3` | Generic `GenModelInverter` is insufficient; add composable generator, electrical, plant, and frequency controls. | No direct analogue | -| `WTDTA1`, `WTDS`, `WTARA1`, `WTPTA1`, `WTTQA1`, `WTARV1` | Add wind-turbine drive-train, aerodynamic, pitch, torque, and renewable-voltage submodels and interfaces. | No direct analogue | -| `PVD1`, `ESD1`, `EV1`, `EV2`, `DGPRCT1`, `DGPRCTExt` | Distributed energy-resource and protection models. | Untriaged | -| `Fault`, `Alter`, `TimeSeries`, `Toggle` | Event/action semantics and time-series input mapping. | Untriaged | -| `Summary`, `Output` | Reporting configuration; define output-channel mapping after model compatibility. | Untriaged | +| ANDES model(s) | GridDyn mapping / next action | Status | +| --------------------------------------------------------------------------------------------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ------------------ | +| `GENCLS` | `GenModelClassical`; parameter and trajectory comparison required. | Partial | +| `GENROU` | `GenModelGENROU` equations, DYR mapping, and five-machine initialization match ANDES references; fix exact machine identity/base handling, native ANDES import, and trajectory parity. | Partial | +| `TGOV1` | `GovernorTgov1`; ANDES-equation, DYR-order, limiter, initialization, Jacobian, and isolated trajectory regressions are complete. | Implemented | +| `EXDC2` | `ExciterDC2A` candidate; add omitted transducer/switch/saturation behavior and prove equation equivalence. | Partial | +| `ZIP`, `FLoad` | GridDyn static/dynamic load models; identify exact parameter and frequency-response equivalence. | Partial | +| `Motor3`, `Motor5` | `MotorLoad3`, `MotorLoad5`; parameter mapping and trajectory comparisons required. | Partial | +| `ACE`, `ACEc`, `COI` | Map area-control and center-of-inertia services. | Untriaged | +| `PLBVFU1`, `IEEEVC` | No direct voltage-compensator/playback analogue; define exciter input and playback interfaces. | No direct analogue | +| `IEEEG1` | `GovernorIeeeG1`; frozen equations, limits, DYR mapping, one-/two-machine connections, initialization, and Jacobian checks are complete; unequal MBASE, native import, and disturbed trajectory remain. | Implemented | +| `TG2`, `TGOV1DB`, `TGOV1N`, `TGOV1NDB`, `IEESGO`, `GAST`, `HYGOV`, `HYGOVDB`, `HYGOV4` | Existing hydro/reheat/steam classes are candidates only; each needs an equation audit or a new model. | Planned | +| `SHAFT5` | Multi-mass shaft model; map states and mechanical interfaces. | Untriaged | +| `ESDC1A`, `ESDC2A`, `SEXS`, `IEEET1`, `IEEET3` | Existing DC/IEEE/SEXS exciters are candidates; complete model-specific audits and DYR/native mappings. | Partial | +| `ESST3A` | `ExciterESST3A`; exact GENROU path plus documented reduced-order synchronous-machine signal approximations. | Implemented | +| `EXST1` | `ExciterEXST1`; positive-time-constant equations and DYR mapping, with the documented corrected regulator-output limiter; zero-time-constant blocks, native import, and a captured trajectory remain. | Implemented | +| `IEEEX1`, `EXAC1`, `EXAC2`, `EXAC4`, `ESST4B`, `AC8B`, `ESAC1A`, `ESST1A`, `ESAC5A` | No exact named GridDyn implementations; implement individually rather than mapping by family name. | No direct analogue | +| `IEEEST`, `ST2CUT` | Base `Stabilizer` is not a functional equivalent; implement both models and signal connections. | No direct analogue | +| `BusFreq`, `BusROCOF`, `PMU`, `PLL1`, `PLL2`, `FreqDiv` | Measurement and frequency-estimation models. | Untriaged | +| `REGCA1`, `REGCP1`, `REECA1`, `REECA1E`, `REECA1G`, `REECB1`, `REPCA1`, `REGCV1`, `REGCV2`, `REGF1`, `REGF2`, `REGF3` | Generic `GenModelInverter` is insufficient; add composable generator, electrical, plant, and frequency controls. | No direct analogue | +| `WTDTA1`, `WTDS`, `WTARA1`, `WTPTA1`, `WTTQA1`, `WTARV1` | Add wind-turbine drive-train, aerodynamic, pitch, torque, and renewable-voltage submodels and interfaces. | No direct analogue | +| `PVD1`, `ESD1`, `EV1`, `EV2`, `DGPRCT1`, `DGPRCTExt` | Distributed energy-resource and protection models. | Untriaged | +| `Fault`, `Alter`, `TimeSeries`, `Toggle` | Event/action semantics and time-series input mapping. | Untriaged | +| `Summary`, `Output` | Reporting configuration; define output-channel mapping after model compatibility. | Untriaged | ## Numerical-regression policy diff --git a/src/fileInput/gridDynReadDYR.cpp b/src/fileInput/gridDynReadDYR.cpp index 77020815..01283fce 100644 --- a/src/fileInput/gridDynReadDYR.cpp +++ b/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 #include #include +#include #include 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(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(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(params[3]); + DynamicGenerator* secondary = nullptr; + if (secondBusId != 0) { + const auto* secondBus = + dynamic_cast(parentObject->findByUserID("bus", secondBusId)); + const int secondGenId = static_cast(params[4]); + if ((secondBus == nullptr) || (secondGenId <= 0)) { + throw InvalidParameterValue("IEEEG1 secondary generator identity"); + } + secondary = dynamic_cast(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 governor( + dynamic_cast(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 diff --git a/src/griddyn/CMakeLists.txt b/src/griddyn/CMakeLists.txt index abc6a074..b04b0b18 100644 --- a/src/griddyn/CMakeLists.txt +++ b/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 diff --git a/src/griddyn/GridSubModel.cpp b/src/griddyn/GridSubModel.cpp index b4946113..8860e82e 100644 --- a/src/griddyn/GridSubModel.cpp +++ b/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); diff --git a/src/griddyn/GridSubModel.h b/src/griddyn/GridSubModel.h index d2f52971..9d099d24 100644 --- a/src/griddyn/GridSubModel.h +++ b/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; }; diff --git a/src/griddyn/generators/DynamicGenerator.cpp b/src/griddyn/generators/DynamicGenerator.cpp index 685098b7..4b9ab5c5 100644 --- a/src/griddyn/generators/DynamicGenerator.cpp +++ b/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(obj) != nullptr) { + const bool reconnectExplicitSource = + mechanicalPowerSourceExplicit && (mechanicalPowerSource == gov); gov = static_cast(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(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(val); + if ((val < 0.0) || (static_cast(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); } diff --git a/src/griddyn/generators/DynamicGenerator.h b/src/griddyn/generators/DynamicGenerator.h index b1e841dd..c4068cee 100644 --- a/src/griddyn/generators/DynamicGenerator.h +++ b/src/griddyn/generators/DynamicGenerator.h @@ -58,6 +58,10 @@ class DynamicGenerator: public Generator { Source* pSetControl = nullptr; //!< source for throttle control Source* vSetControl = nullptr; //!< source for voltage level control IsocController* isoc = nullptr; //!< pointer to a isochronous controller + GridSubModel* mechanicalPowerSource = nullptr; //!< optional non-owning Pmech source + index_t mechanicalPowerOutput = 0; //!< selected output on the Pmech source + std::string mechanicalPowerSourceName; //!< clone/configuration path for the Pmech source + bool mechanicalPowerSourceExplicit = false; //!< use the source instead of the local governor // const double *m_stateTemp = nullptr; //!< temporary state // vector(assumed not writable) const double *m_dstate_dt_Temp = nullptr; //!< // a temporary deriv vector; @@ -104,6 +108,17 @@ class DynamicGenerator: public Generator { @throw unrecognizedObjectError is object is not valid*/ virtual void add(GridSubModel* obj) override; + /** Connect this machine's single mechanical-power input to a submodel output. + The source remains owned and evaluated by its normal parent. */ + void setMechanicalPowerSource(GridSubModel* source, index_t outputIndex = 0); + /** Store a source path for resolution after the complete object tree is available. */ + void setMechanicalPowerSource(std::string_view sourceName, index_t outputIndex = 0); + /** Restore the conventional local-governor-output-zero connection. */ + void clearMechanicalPowerSource(); + GridSubModel* getMechanicalPowerSource() const; + index_t getMechanicalPowerOutput() const; + bool hasExplicitMechanicalPowerSource() const; + virtual void algebraicUpdate(const IOdata& inputs, const StateData& stateDataValue, double update[], @@ -213,6 +228,8 @@ class DynamicGenerator: public Generator { GridSubModel* replaceModel(GridSubModel* newObject, GridSubModel* oldObject, index_t newIndex); + void resolveMechanicalPowerSource(); + void buildDynModel(DynModel dynModel); }; diff --git a/src/griddyn/governors/Governor.cpp b/src/griddyn/governors/Governor.cpp index f8252a2a..6a1a2ad7 100644 --- a/src/griddyn/governors/Governor.cpp +++ b/src/griddyn/governors/Governor.cpp @@ -36,6 +36,9 @@ namespace governors { static ChildTypeFactory gFgov5("governor", std::to_array({"tgov1"})); + static ChildTypeFactory + gFgov6("governor", std::to_array({"ieeeg1"})); + } // namespace governors using units::convert; using units::puHz; diff --git a/src/griddyn/governors/GovernorIeeeG1.cpp b/src/griddyn/governors/GovernorIeeeG1.cpp new file mode 100644 index 00000000..e8f43c84 --- /dev/null +++ b/src/griddyn/governors/GovernorIeeeG1.cpp @@ -0,0 +1,599 @@ +/* + * Copyright (c) 2014-2026, Lawrence Livermore National Security + * See the top-level NOTICE for additional details. All rights reserved. + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include "GovernorIeeeG1.h" + +#include "core/CoreExceptions.h" +#include "core/CoreObjectTemplates.hpp" +#include "utilities/MatrixData.hpp" +#include +#include +#include +#include +#include + +namespace griddyn::governors { +// The equations below mirror the published block diagram. Extra precedence +// parentheses obscure the signal flow without changing the calculation. +// NOLINTBEGIN(readability-math-missing-parentheses) +namespace { + constexpr double initializationTolerance = 1e-7; +} + +GovernorIeeeG1::GovernorIeeeG1(const std::string& objName): Governor(objName) +{ + K = 20.0; + T1 = 1.0; + T2 = 1.0; + T3 = 0.1; + Pmax = 5.0; + Pmin = 0.0; + m_outputSize = 2; + opFlags.set(IGNORE_DEADBAND); + opFlags.set(IGNORE_FILTER); + opFlags.set(IGNORE_THROTTLE); + opFlags.set(USES_POWER_LIMITS); + opFlags.set(USES_RAMP_LIMITS); +} + +CoreObject* GovernorIeeeG1::clone(CoreObject* obj) const +{ + auto* governorClone = cloneBase(this, obj); + if (governorClone == nullptr) { + return obj; + } + governorClone->Uo = Uo; + governorClone->Uc = Uc; + governorClone->turbineTime = turbineTime; + governorClone->powerFraction = powerFraction; + governorClone->initializationTarget = initializationTarget; + governorClone->initializationTargetSet = initializationTargetSet; + return governorClone; +} + +void GovernorIeeeG1::dynObjectInitializeA(CoreTime time0, std::uint32_t /*flags*/) +{ + const double fractionSum = std::accumulate(powerFraction.begin(), powerFraction.end(), 0.0); + if ((Pmax < Pmin) || !std::isfinite(fractionSum) || (fractionSum <= 0.0)) { + throw InvalidParameterValue("IEEEG1 limits or power fractions"); + } + + index_t nextState = 0; + leadLagState = (T1 > 0.0) ? nextState++ : kInvalidLocation; + valveState = nextState++; + for (index_t stage = 0; std::cmp_less(stage, turbineTime.size()); ++stage) { + turbineState[stage] = (turbineTime[stage] > 0.0) ? nextState++ : kInvalidLocation; + } + + auto& local = offsets.local().local; + local.algSize = 2; + local.diffSize = nextState; + local.algRoots = 2; + local.diffRoots = 0; + local.jacSize = 48; + prevTime = time0; + initializationActive = false; + opFlags.reset(RATE_LIMITED); + opFlags.reset(RATE_LIMIT_HIGH); + opFlags.reset(POWER_LIMITED); + opFlags.reset(POWER_LIMIT_HIGH); +} + +void GovernorIeeeG1::dynObjectInitializeB(const IOdata& /*inputs*/, + const IOdata& desiredOutput, + IOdata& fieldSet) +{ + if (desiredOutput.empty() || !std::isfinite(desiredOutput[hpOutput])) { + throw InvalidParameterValue("IEEEG1 initial high-pressure output"); + } + initializationTarget[hpOutput] = desiredOutput[hpOutput]; + initializationTargetSet[hpOutput] = true; + initializationActive = true; + initializeStatesFromTargets(); + + fieldSet.resize(2); + fieldSet[govpSetInLocation] = Pset; +} + +bool GovernorIeeeG1::setOutputInitializationTarget(index_t outputIndex, double target) +{ + if ((outputIndex < 0) || (outputIndex >= 2) || !std::isfinite(target)) { + throw InvalidParameterValue("IEEEG1 output initialization target"); + } + initializationTarget[outputIndex] = target; + initializationTargetSet[outputIndex] = true; + if (initializationActive) { + initializeStatesFromTargets(); + } + return true; +} + +std::array GovernorIeeeG1::normalizedFractions() const +{ + const double sum = std::accumulate(powerFraction.begin(), powerFraction.end(), 0.0); + std::array normalized{}; + std::ranges::transform(powerFraction, normalized.begin(), [sum](double value) { + return value / sum; + }); + return normalized; +} + +void GovernorIeeeG1::initializeStatesFromTargets() +{ + const auto fraction = normalizedFractions(); + const double hpFraction = fraction[0] + fraction[2] + fraction[4] + fraction[6]; + const double lpFraction = fraction[1] + fraction[3] + fraction[5] + fraction[7]; + + bool totalSet = false; + double totalPower = 0.0; + const auto addTarget = [&](index_t outputIndex, double branchFraction) { + if (!initializationTargetSet[outputIndex]) { + return; + } + if (branchFraction <= initializationTolerance) { + if (std::abs(initializationTarget[outputIndex]) > initializationTolerance) { + throw InvalidParameterValue("IEEEG1 initial output has zero power fraction"); + } + return; + } + const double candidate = initializationTarget[outputIndex] / branchFraction; + if (totalSet && + (std::abs(candidate - totalPower) > + initializationTolerance * (1.0 + std::abs(totalPower)))) { + throw InvalidParameterValue("IEEEG1 initial HP/LP powers are not proportional"); + } + totalPower = candidate; + totalSet = true; + }; + addTarget(hpOutput, hpFraction); + addTarget(lpOutput, lpFraction); + if (!totalSet || (totalPower < Pmin - initializationTolerance) || + (totalPower > Pmax + initializationTolerance)) { + throw InvalidParameterValue("IEEEG1 initial valve position outside limits"); + } + + Pset = totalPower; + double* diffState = m_state.data() + offsets.getDiffOffset(cLocalSolverMode); + if (leadLagState != kInvalidLocation) { + diffState[leadLagState] = 0.0; + } + diffState[valveState] = totalPower; + for (const auto stageState : turbineState) { + if (stageState != kInvalidLocation) { + diffState[stageState] = totalPower; + } + } + + const auto output = powerOutputs(diffState); + double* algState = m_state.data() + offsets.getAlgOffset(cLocalSolverMode); + algState[hpOutput] = output[hpOutput]; + algState[lpOutput] = output[lpOutput]; +} + +double GovernorIeeeG1::leadLagOutput(const IOdata& inputs, const double diffState[]) const +{ + const double speedDeviation = 1.0 - inputs[govOmegaInLocation]; + if (leadLagState == kInvalidLocation) { + return K * speedDeviation; + } + const double state = diffState[leadLagState]; + return K * (state + (T2 / T1) * (speedDeviation - state)); +} + +double GovernorIeeeG1::unlimitedValveRate(const IOdata& inputs, const double diffState[]) const +{ + return (leadLagOutput(inputs, diffState) + Pset - diffState[valveState]) / T3; +} + +double GovernorIeeeG1::limitedValveRate(const IOdata& inputs, const double diffState[]) const +{ + return std::clamp(unlimitedValveRate(inputs, diffState), + static_cast(Uc), + static_cast(Uo)); +} + +std::array GovernorIeeeG1::turbineOutputs(const double diffState[]) const +{ + std::array output{}; + double input = diffState[valveState]; + for (index_t stage = 0; std::cmp_less(stage, output.size()); ++stage) { + output[stage] = + (turbineState[stage] != kInvalidLocation) ? diffState[turbineState[stage]] : input; + input = output[stage]; + } + return output; +} + +std::array GovernorIeeeG1::powerOutputs(const double diffState[]) const +{ + const auto fraction = normalizedFractions(); + const auto stage = turbineOutputs(diffState); + return {{fraction[0] * stage[0] + fraction[2] * stage[1] + fraction[4] * stage[2] + + fraction[6] * stage[3], + fraction[1] * stage[0] + fraction[3] * stage[1] + fraction[5] * stage[2] + + fraction[7] * stage[3]}}; +} + +int GovernorIeeeG1::rateLimitStatus(const IOdata& inputs, const double diffState[]) const +{ + const double rate = unlimitedValveRate(inputs, diffState); + if (rate >= Uo) { + return 1; + } + return (rate <= Uc) ? -1 : 0; +} + +int GovernorIeeeG1::valveLimitStatus(const IOdata& inputs, const double diffState[]) const +{ + const double rate = limitedValveRate(inputs, diffState); + if ((diffState[valveState] >= Pmax) && (rate >= 0.0)) { + return 1; + } + if ((diffState[valveState] <= Pmin) && (rate <= 0.0)) { + return -1; + } + return 0; +} + +void GovernorIeeeG1::residual(const IOdata& inputs, + const StateData& stateData, + double resid[], + const SolverMode& sMode) +{ + const auto locations = offsets.getLocations(stateData, resid, sMode, this); + if (hasAlgebraic(sMode)) { + const auto output = powerOutputs(locations.diffStateLoc); + locations.destLoc[hpOutput] = output[hpOutput] - locations.algStateLoc[hpOutput]; + locations.destLoc[lpOutput] = output[lpOutput] - locations.algStateLoc[lpOutput]; + } + if (!hasDifferential(sMode)) { + return; + } + derivative(inputs, stateData, resid, sMode); + for (index_t state = 0; state < locations.diffSize; ++state) { + locations.destDiffLoc[state] -= locations.dstateLoc[state]; + } +} + +void GovernorIeeeG1::derivative(const IOdata& inputs, + const StateData& stateData, + double deriv[], + const SolverMode& sMode) +{ + if (!hasDifferential(sMode)) { + return; + } + const auto locations = offsets.getLocations(stateData, deriv, sMode, this); + const double* state = locations.diffStateLoc; + double* stateDerivative = locations.destDiffLoc; + + if (leadLagState != kInvalidLocation) { + stateDerivative[leadLagState] = + (1.0 - inputs[govOmegaInLocation] - state[leadLagState]) / T1; + } + stateDerivative[valveState] = + (valveLimitStatus(inputs, state) == 0) ? limitedValveRate(inputs, state) : 0.0; + + double input = state[valveState]; + for (index_t stage = 0; std::cmp_less(stage, turbineState.size()); ++stage) { + if (turbineState[stage] != kInvalidLocation) { + stateDerivative[turbineState[stage]] = + (input - state[turbineState[stage]]) / turbineTime[stage]; + input = state[turbineState[stage]]; + } + } +} + +void GovernorIeeeG1::algebraicUpdate(const IOdata& /*inputs*/, + const StateData& stateData, + double update[], + const SolverMode& sMode, + double /*alpha*/) +{ + if (!hasAlgebraic(sMode)) { + return; + } + const auto locations = offsets.getLocations(stateData, update, sMode, this); + const auto output = powerOutputs(locations.diffStateLoc); + locations.destLoc[hpOutput] = output[hpOutput]; + locations.destLoc[lpOutput] = output[lpOutput]; +} + +void GovernorIeeeG1::jacobianElements(const IOdata& inputs, + const StateData& stateData, + MatrixData& matrixData, + const IOlocs& inputLocs, + const SolverMode& sMode) +{ + const auto locations = offsets.getLocations(stateData, sMode, this); + const index_t algOffset = locations.algOffset; + const index_t diffOffset = locations.diffOffset; + const double* state = locations.diffStateLoc; + + std::array stageSource{}; + index_t source = valveState; + for (index_t stage = 0; std::cmp_less(stage, stageSource.size()); ++stage) { + if (turbineState[stage] != kInvalidLocation) { + source = turbineState[stage]; + } + stageSource[stage] = source; + } + + if (hasAlgebraic(sMode)) { + const auto fraction = normalizedFractions(); + matrixData.assign(algOffset + hpOutput, algOffset + hpOutput, -1.0); + matrixData.assign(algOffset + lpOutput, algOffset + lpOutput, -1.0); + for (std::size_t stage = 0; stage < stageSource.size(); ++stage) { + matrixData.assign(algOffset + hpOutput, + diffOffset + stageSource[stage], + fraction[2 * stage]); + matrixData.assign(algOffset + lpOutput, + diffOffset + stageSource[stage], + fraction[2 * stage + 1]); + } + } + if (!hasDifferential(sMode)) { + return; + } + + if (leadLagState != kInvalidLocation) { + matrixData.assign(diffOffset + leadLagState, + diffOffset + leadLagState, + -1.0 / T1 - stateData.cj); + matrixData.assignCheckCol(diffOffset + leadLagState, + inputLocs[govOmegaInLocation], + -1.0 / T1); + } + + const index_t valveRow = diffOffset + valveState; + const bool valveLimited = valveLimitStatus(inputs, state) != 0; + const bool rateLimited = rateLimitStatus(inputs, state) != 0; + matrixData.assign(valveRow, valveRow, -stateData.cj); + if (!valveLimited && !rateLimited) { + matrixData.assign(valveRow, valveRow, -1.0 / T3); + if (leadLagState == kInvalidLocation) { + matrixData.assignCheckCol(valveRow, + inputLocs[govOmegaInLocation], + -static_cast(K) / T3); + } else { + const double leadRatio = T2 / T1; + matrixData.assign(valveRow, diffOffset + leadLagState, K * (1.0 - leadRatio) / T3); + matrixData.assignCheckCol(valveRow, inputLocs[govOmegaInLocation], -K * leadRatio / T3); + } + } + + source = valveState; + for (index_t stage = 0; stage < 4; ++stage) { + if (turbineState[stage] == kInvalidLocation) { + continue; + } + const index_t row = diffOffset + turbineState[stage]; + matrixData.assign(row, diffOffset + source, 1.0 / turbineTime[stage]); + matrixData.assign(row, row, -1.0 / turbineTime[stage] - stateData.cj); + source = turbineState[stage]; + } +} + +void GovernorIeeeG1::timestep(CoreTime time, const IOdata& inputs, const SolverMode& /*sMode*/) +{ + derivative(inputs, emptyStateData, m_dstate_dt.data(), cLocalSolverMode); + const double timeStep = time - prevTime; + const index_t diffOffset = offsets.getDiffOffset(cLocalSolverMode); + for (index_t state = 0; state < offsets.local().local.diffSize; ++state) { + m_state[diffOffset + state] += timeStep * m_dstate_dt[diffOffset + state]; + } + double* diffState = m_state.data() + diffOffset; + diffState[valveState] = + std::clamp(diffState[valveState], static_cast(Pmin), static_cast(Pmax)); + const auto output = powerOutputs(diffState); + const index_t algOffset = offsets.getAlgOffset(cLocalSolverMode); + m_state[algOffset + hpOutput] = output[hpOutput]; + m_state[algOffset + lpOutput] = output[lpOutput]; + updateLimitFlags(inputs, diffState); + prevTime = time; +} + +bool GovernorIeeeG1::updateLimitFlags(const IOdata& inputs, const double diffState[]) +{ + const int rateStatus = rateLimitStatus(inputs, diffState); + const int valveStatus = valveLimitStatus(inputs, diffState); + const bool changed = (opFlags[RATE_LIMITED] != (rateStatus != 0)) || + (opFlags[RATE_LIMIT_HIGH] != (rateStatus > 0)) || + (opFlags[POWER_LIMITED] != (valveStatus != 0)) || + (opFlags[POWER_LIMIT_HIGH] != (valveStatus > 0)); + opFlags.set(RATE_LIMITED, rateStatus != 0); + opFlags.set(RATE_LIMIT_HIGH, rateStatus > 0); + opFlags.set(POWER_LIMITED, valveStatus != 0); + opFlags.set(POWER_LIMIT_HIGH, valveStatus > 0); + return changed; +} + +void GovernorIeeeG1::rootTest(const IOdata& inputs, + const StateData& stateData, + double roots[], + const SolverMode& sMode) +{ + const auto locations = offsets.getLocations(stateData, sMode, this); + const double* state = locations.diffStateLoc; + const index_t rootOffset = offsets.getRootOffset(sMode); + const double unlimitedRate = unlimitedValveRate(inputs, state); + if (opFlags[RATE_LIMITED]) { + roots[rootOffset] = opFlags[RATE_LIMIT_HIGH] ? Uo - unlimitedRate : unlimitedRate - Uc; + } else { + roots[rootOffset] = std::min(Uo - unlimitedRate, unlimitedRate - Uc); + } + + const double limitedRate = limitedValveRate(inputs, state); + if (opFlags[POWER_LIMITED]) { + roots[rootOffset + 1] = opFlags[POWER_LIMIT_HIGH] ? limitedRate : -limitedRate; + } else { + roots[rootOffset + 1] = std::min(Pmax - state[valveState], state[valveState] - Pmin); + } +} + +void GovernorIeeeG1::rootTrigger(CoreTime /*time*/, + const IOdata& inputs, + const std::vector& rootMask, + const SolverMode& sMode) +{ + const index_t rootOffset = offsets.getRootOffset(sMode); + if ((rootMask[rootOffset] == 0) && (rootMask[rootOffset + 1] == 0)) { + return; + } + double* state = m_state.data() + offsets.getDiffOffset(cLocalSolverMode); + if (state[valveState] >= Pmax) { + state[valveState] = Pmax; + } else if (state[valveState] <= Pmin) { + state[valveState] = Pmin; + } + if (updateLimitFlags(inputs, state)) { + alert(this, JAC_COUNT_CHANGE); + } +} + +ChangeCode GovernorIeeeG1::rootCheck(const IOdata& inputs, + const StateData& /*stateData*/, + const SolverMode& /*sMode*/, + CheckLevel /*level*/) +{ + double* state = m_state.data() + offsets.getDiffOffset(cLocalSolverMode); + if (state[valveState] > Pmax) { + state[valveState] = Pmax; + } else if (state[valveState] < Pmin) { + state[valveState] = Pmin; + } + if (updateLimitFlags(inputs, state)) { + alert(this, JAC_COUNT_CHANGE); + return ChangeCode::JACOBIAN_CHANGE; + } + return ChangeCode::NO_CHANGE; +} + +void GovernorIeeeG1::set(std::string_view param, std::string_view val) +{ + Governor::set(param, val); +} + +void GovernorIeeeG1::set(std::string_view param, double val, units::unit unitType) +{ + if ((param == "k") || (param == "droop")) { + if (!std::isfinite(val) || (val <= 0.0)) { + throw InvalidParameterValue("IEEEG1 gain must be positive and finite"); + } + Governor::set(param, val, unitType); + } else if (param == "r") { + if (!std::isfinite(val) || (val <= 0.0)) { + throw InvalidParameterValue("IEEEG1 R must be positive and finite"); + } + Governor::set(param, val, unitType); + } else if ((param == "t1") || (param == "t2")) { + if (!std::isfinite(val) || (val < 0.0)) { + throw InvalidParameterValue("IEEEG1 T1/T2 must be nonnegative and finite"); + } + Governor::set(param, val, unitType); + } else if (param == "t3") { + if (!std::isfinite(val) || (val <= 0.0)) { + throw InvalidParameterValue("IEEEG1 T3 must be positive and finite"); + } + Governor::set(param, val, unitType); + } else if (param == "uo") { + if (!std::isfinite(val) || (val < 0.0)) { + throw InvalidParameterValue("IEEEG1 UO must be nonnegative and finite"); + } + Uo = val; + } else if (param == "uc") { + if (!std::isfinite(val) || (val > 0.0)) { + throw InvalidParameterValue("IEEEG1 UC must be nonpositive and finite"); + } + Uc = val; + } else if ((param == "pmax") || (param == "pmin")) { + if (!std::isfinite(val)) { + throw InvalidParameterValue("IEEEG1 valve limit must be finite"); + } + Governor::set(param, val, unitType); + } else if ((param.size() == 2) && (param[0] == 't') && (param[1] >= '4') && (param[1] <= '7')) { + if (!std::isfinite(val) || (val < 0.0)) { + throw InvalidParameterValue("IEEEG1 turbine time constant must be nonnegative"); + } + turbineTime[static_cast(param[1] - '4')] = val; + } else if ((param.size() == 2) && (param[0] == 'k') && (param[1] >= '1') && (param[1] <= '8')) { + if (!std::isfinite(val) || (val < 0.0)) { + throw InvalidParameterValue("IEEEG1 power fraction must be nonnegative"); + } + powerFraction[static_cast(param[1] - '1')] = val; + } else { + Governor::set(param, val, unitType); + } +} + +double GovernorIeeeG1::get(std::string_view param, units::unit unitType) const +{ + if (param == "uo") { + return Uo; + } + if (param == "uc") { + return Uc; + } + if ((param.size() == 2) && (param[0] == 't') && (param[1] >= '4') && (param[1] <= '7')) { + return turbineTime[static_cast(param[1] - '4')]; + } + if ((param.size() == 2) && (param[0] == 'k') && (param[1] >= '1') && (param[1] <= '8')) { + return powerFraction[static_cast(param[1] - '1')]; + } + if ((param.size() == 3) && (param[0] == 'k') && (param[1] >= '1') && (param[1] <= '8') && + (param[2] == 'n')) { + return normalizedFractions()[static_cast(param[1] - '1')]; + } + return Governor::get(param, unitType); +} + +stringVec GovernorIeeeG1::localStateNames() const +{ + stringVec names{"php", "plp"}; + if (leadLagState != kInvalidLocation) { + names.emplace_back("leadlag"); + } + names.emplace_back("valve"); + for (index_t stage = 0; stage < 4; ++stage) { + if (turbineState[stage] != kInvalidLocation) { + names.emplace_back("p" + std::to_string(stage + 4)); + } + } + return names; +} + +index_t GovernorIeeeG1::findIndex(std::string_view field, const SolverMode& sMode) const +{ + if ((field == "php") || (field == "hp") || (field == "pm") || (field == "pmech")) { + return offsets.getAlgOffset(sMode) + hpOutput; + } + if ((field == "plp") || (field == "lp")) { + return offsets.getAlgOffset(sMode) + lpOutput; + } + if ((field == "leadlag") || (field == "ll")) { + return (leadLagState != kInvalidLocation) ? offsets.getDiffOffset(sMode) + leadLagState : + kInvalidLocation; + } + if ((field == "valve") || (field == "pv")) { + return offsets.getDiffOffset(sMode) + valveState; + } + if ((field.size() == 2) && (field[0] == 'p') && (field[1] >= '4') && (field[1] <= '7')) { + const index_t stage = field[1] - '4'; + return (turbineState[stage] != kInvalidLocation) ? + offsets.getDiffOffset(sMode) + turbineState[stage] : + kInvalidLocation; + } + return kInvalidLocation; +} + +const std::vector& GovernorIeeeG1::outputNames() const +{ + static const std::vector names{{"php", "hp", "pmech", "power", "output"}, + {"plp", "lp"}}; + return names; +} +// NOLINTEND(readability-math-missing-parentheses) +} // namespace griddyn::governors diff --git a/src/griddyn/governors/GovernorIeeeG1.h b/src/griddyn/governors/GovernorIeeeG1.h new file mode 100644 index 00000000..9b662a05 --- /dev/null +++ b/src/griddyn/governors/GovernorIeeeG1.h @@ -0,0 +1,130 @@ +/* + * Copyright (c) 2014-2026, Lawrence Livermore National Security + * See the top-level NOTICE for additional details. All rights reserved. + * SPDX-License-Identifier: BSD-3-Clause + */ + +#pragma once + +#include "../Governor.h" +#include +#include +#include + +namespace griddyn::governors { +/** IEEE Type 1 speed-governing model (PSS/E IEEEG1). + * + * This class follows frozen ANDES v2.0.0 + * `andes/models/governor/ieeeg1.py`. It retains the published block structure + * while eliminating algebraic variables internal to ordinary lead-lag and lag + * blocks. A zero T1 or turbine-stage time constant is therefore an exact + * bypass and owns no differential state. + * + * With speed deviation @f$w_d=1-\omega@f$, lead-lag state @f$x_g@f$, + * valve position @f$p_v@f$, and turbine-stage outputs @f$p_4,\ldots,p_7@f$, + * + * @f[ + * T_1\dot x_g=w_d-x_g,\qquad + * y_g=K\left[x_g+\frac{T_2}{T_1}(w_d-x_g)\right], + * @f] + * @f[ + * v_s=\frac{y_g+P_0-p_v}{T_3},\quad + * \dot p_v=\operatorname{AW}_{[P_{MIN},P_{MAX}]} + * (\operatorname{limit}_{[U_C,U_O]}(v_s)), + * @f] + * @f[ + * T_i\dot p_i=p_{i-1}-p_i,\quad i=4,\ldots,7. + * @f] + * + * The two outputs are the normalized odd- and even-coefficient sums: + * @f$P_{HP}=K_{1n}p_4+K_{3n}p_5+K_{5n}p_6+K_{7n}p_7@f$ and + * @f$P_{LP}=K_{2n}p_4+K_{4n}p_5+K_{6n}p_6+K_{8n}p_7@f$. + */ +class GovernorIeeeG1: public Governor { + public: + enum IEEEG1Flags { + RATE_LIMITED = OBJECT_FLAG9, + RATE_LIMIT_HIGH = OBJECT_FLAG10, + }; + + static constexpr index_t hpOutput = 0; + static constexpr index_t lpOutput = 1; + + protected: + model_parameter Uo = 0.1; + model_parameter Uc = -0.1; + std::array turbineTime{{0.4, 8.0, 0.5, 0.05}}; + std::array powerFraction{{0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0}}; + + index_t leadLagState = kInvalidLocation; + index_t valveState = kInvalidLocation; + std::array turbineState{ + {kInvalidLocation, kInvalidLocation, kInvalidLocation, kInvalidLocation}}; + std::array initializationTarget{{0.0, 0.0}}; + std::array initializationTargetSet{{false, false}}; + bool initializationActive = false; + + public: + explicit GovernorIeeeG1(const std::string& objName = "govIeeeG1_#"); + CoreObject* clone(CoreObject* obj = nullptr) const override; + + void dynObjectInitializeA(CoreTime time0, std::uint32_t flags) override; + void dynObjectInitializeB(const IOdata& inputs, + const IOdata& desiredOutput, + IOdata& fieldSet) override; + bool setOutputInitializationTarget(index_t outputIndex, double target) override; + + void set(std::string_view param, std::string_view val) override; + void set(std::string_view param, double val, units::unit unitType = units::defunit) override; + double get(std::string_view param, units::unit unitType = units::defunit) const override; + + void residual(const IOdata& inputs, + const StateData& stateData, + double resid[], + const SolverMode& sMode) override; + void derivative(const IOdata& inputs, + const StateData& stateData, + double deriv[], + const SolverMode& sMode) override; + void algebraicUpdate(const IOdata& inputs, + const StateData& stateData, + double update[], + const SolverMode& sMode, + double alpha) override; + void jacobianElements(const IOdata& inputs, + const StateData& stateData, + MatrixData& matrixData, + const IOlocs& inputLocs, + const SolverMode& sMode) override; + void timestep(CoreTime time, const IOdata& inputs, const SolverMode& sMode) override; + + void rootTest(const IOdata& inputs, + const StateData& stateData, + double roots[], + const SolverMode& sMode) override; + void rootTrigger(CoreTime time, + const IOdata& inputs, + const std::vector& rootMask, + const SolverMode& sMode) override; + ChangeCode rootCheck(const IOdata& inputs, + const StateData& stateData, + const SolverMode& sMode, + CheckLevel level) override; + + stringVec localStateNames() const override; + index_t findIndex(std::string_view field, const SolverMode& sMode) const override; + const std::vector& outputNames() const override; + + private: + std::array normalizedFractions() const; + double leadLagOutput(const IOdata& inputs, const double diffState[]) const; + double unlimitedValveRate(const IOdata& inputs, const double diffState[]) const; + double limitedValveRate(const IOdata& inputs, const double diffState[]) const; + std::array turbineOutputs(const double diffState[]) const; + std::array powerOutputs(const double diffState[]) const; + int rateLimitStatus(const IOdata& inputs, const double diffState[]) const; + int valveLimitStatus(const IOdata& inputs, const double diffState[]) const; + bool updateLimitFlags(const IOdata& inputs, const double diffState[]); + void initializeStatesFromTargets(); +}; +} // namespace griddyn::governors diff --git a/src/griddyn/governors/GovernorTypes.h b/src/griddyn/governors/GovernorTypes.h index f1b03f8a..b31bf0fb 100644 --- a/src/griddyn/governors/GovernorTypes.h +++ b/src/griddyn/governors/GovernorTypes.h @@ -7,6 +7,7 @@ #pragma once #include "griddyn/governors/GovernorHydro.h" +#include "griddyn/governors/GovernorIeeeG1.h" #include "griddyn/governors/GovernorReheat.h" #include "griddyn/governors/GovernorSteamTCSR.h" #include "griddyn/governors/GovernorTgov1.h" diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 28487510..1d0a82c7 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -89,7 +89,9 @@ if(GRIDDYN_ENABLE_HELICS_EXECUTABLE) list(APPEND optional_component_gtest_sources ../src/helics/test/helicsTest.cpp) list(APPEND optional_component_gtest_sources exeTestHelper.cpp) endif() -if(GRIDDYN_ENABLE_FMI_EXPORT) +# The export test loads the generated FMU through the FMI import wrapper, so it requires both the +# exporter and the optional FMI import library. +if(GRIDDYN_ENABLE_FMI_EXPORT AND GRIDDYN_ENABLE_FMI) list(APPEND optional_component_gtest_sources ../src/fmi_export/testFMIExport.cpp) endif() diff --git a/test/andesTests/testAndesDyrReader.cpp b/test/andesTests/testAndesDyrReader.cpp index 2de45982..7a9e58ff 100644 --- a/test/andesTests/testAndesDyrReader.cpp +++ b/test/andesTests/testAndesDyrReader.cpp @@ -11,7 +11,9 @@ #include "griddyn/GridDynSimulation.h" #include "griddyn/exciters/ExciterESST3A.h" #include "griddyn/exciters/ExciterEXST1.h" +#include "griddyn/generators/DynamicGenerator.h" #include "griddyn/genmodels/GenModelGENROU.h" +#include "griddyn/governors/GovernorIeeeG1.h" #include "griddyn/governors/GovernorTgov1.h" #include #include @@ -20,6 +22,7 @@ #include #include #include +#include namespace { constexpr std::string_view andesTestDirectory{GRIDDYN_TEST_DIRECTORY "/andes_tests/"}; @@ -144,6 +147,60 @@ TEST(AndesDyrReaderTests, MapsTgov1ParametersInAndesDyrOrder) } } +TEST(AndesDyrReaderTests, MapsIeeeG1ParametersInFrozenAndesDyrOrder) +{ + auto simulation = std::make_unique(); + griddyn::loadFile(simulation.get(), makeAndesTestPath("ieee14.raw")); + griddyn::loadFile(simulation.get(), makeAndesTestPath("ieee14_genrou.dyr")); + griddyn::loadFile(simulation.get(), makeAndesTestPath("ieee14_ieeeg1.dyr")); + + auto* bus = dynamic_cast(simulation->findByUserID("bus", 2)); + ASSERT_NE(bus, nullptr); + auto* generator = dynamic_cast(bus->getGen(0)); + ASSERT_NE(generator, nullptr); + auto* governor = dynamic_cast(generator->find("governor")); + ASSERT_NE(governor, nullptr); + + const std::pair expected[]{ + {"k", 18.0}, {"t1", 0.12}, {"t2", 0.03}, {"t3", 0.22}, {"uo", 0.42}, + {"uc", -0.31}, {"pmax", 1.25}, {"pmin", 0.15}, {"t4", 0.14}, {"k1", 0.20}, + {"k2", 0.0}, {"t5", 0.25}, {"k3", 0.30}, {"k4", 0.0}, {"t6", 0.36}, + {"k5", 0.10}, {"k6", 0.0}, {"t7", 0.47}, {"k7", 0.40}, {"k8", 0.0}}; + for (const auto& parameter : expected) { + EXPECT_DOUBLE_EQ(governor->get(parameter.first), parameter.second) << parameter.first; + } + EXPECT_EQ(generator->getMechanicalPowerSource(), governor); + EXPECT_EQ(generator->getMechanicalPowerOutput(), griddyn::governors::GovernorIeeeG1::hpOutput); + + ASSERT_EQ(simulation->dynInitialize(), 0); + EXPECT_EQ(runResidualCheck(simulation, griddyn::cDaeSolverMode, false), 0); + EXPECT_EQ(runJacobianCheck(simulation, griddyn::cDaeSolverMode, false), 0); +} + +TEST(AndesDyrReaderTests, ConnectsIeeeG1LowPressureOutputToSecondGenerator) +{ + auto simulation = std::make_unique(); + griddyn::loadFile(simulation.get(), makeAndesTestPath("ieee14.raw")); + griddyn::loadFile(simulation.get(), makeAndesTestPath("ieee14_genrou.dyr")); + griddyn::loadFile(simulation.get(), makeAndesTestPath("ieee14_ieeeg1_cross.dyr")); + + auto* primaryBus = dynamic_cast(simulation->findByUserID("bus", 1)); + auto* secondaryBus = dynamic_cast(simulation->findByUserID("bus", 2)); + ASSERT_NE(primaryBus, nullptr); + ASSERT_NE(secondaryBus, nullptr); + auto* primary = dynamic_cast(primaryBus->getGen(0)); + auto* secondary = dynamic_cast(secondaryBus->getGen(0)); + ASSERT_NE(primary, nullptr); + ASSERT_NE(secondary, nullptr); + auto* governor = dynamic_cast(primary->find("governor")); + ASSERT_NE(governor, nullptr); + + EXPECT_EQ(primary->getMechanicalPowerSource(), governor); + EXPECT_EQ(primary->getMechanicalPowerOutput(), griddyn::governors::GovernorIeeeG1::hpOutput); + EXPECT_EQ(secondary->getMechanicalPowerSource(), governor); + EXPECT_EQ(secondary->getMechanicalPowerOutput(), griddyn::governors::GovernorIeeeG1::lpOutput); +} + TEST(AndesDyrReaderTests, LoadsEsst3aWithGenrouSignals) { auto simulation = std::make_unique(); diff --git a/test/componentTests/testGenerators.cpp b/test/componentTests/testGenerators.cpp index 8a433b55..a549c136 100644 --- a/test/componentTests/testGenerators.cpp +++ b/test/componentTests/testGenerators.cpp @@ -5,9 +5,13 @@ */ #include "../gtestHelper.h" +#include "core/CoreExceptions.h" #include "gmlc/utilities/TimeSeriesMulti.hpp" #include "griddyn/Generator.h" #include "griddyn/GridBus.h" +#include "griddyn/generators/DynamicGenerator.h" +#include "griddyn/genmodels/GenModel6.h" +#include "griddyn/governors/GovernorIeeeSimple.h" #include #include @@ -17,6 +21,61 @@ using namespace griddyn; class GeneratorTests: public GridDynSimulationTestFixture, public ::testing::Test {}; +TEST(DynamicGeneratorModelTests, MechanicalPowerSourceCanBeExternalAndIndexed) +{ + DynamicGenerator primary; + auto* governor = new governors::GovernorIeeeSimple(); + primary.add(governor); + + EXPECT_EQ(primary.getMechanicalPowerSource(), governor); + EXPECT_EQ(primary.getMechanicalPowerOutput(), 0); + EXPECT_FALSE(primary.hasExplicitMechanicalPowerSource()); + + auto* multiOutputSource = new genmodels::GenModel6(); + primary.add(multiOutputSource); + + DynamicGenerator secondary; + secondary.setMechanicalPowerSource(multiOutputSource, 1); + EXPECT_EQ(secondary.getMechanicalPowerSource(), multiOutputSource); + EXPECT_EQ(secondary.getMechanicalPowerOutput(), 1); + EXPECT_TRUE(secondary.hasExplicitMechanicalPowerSource()); + + EXPECT_THROW(secondary.setMechanicalPowerSource(multiOutputSource, 2), InvalidParameterValue); + + secondary.clearMechanicalPowerSource(); + EXPECT_EQ(secondary.getMechanicalPowerSource(), nullptr); + EXPECT_FALSE(secondary.hasExplicitMechanicalPowerSource()); +} + +TEST_F(GeneratorTests, MechanicalPowerSourcePathResolvesAfterTreeAssembly) +{ + const std::string fileName = + std::string(GRIDDYN_TEST_DIRECTORY "/governor_tests/test_gov_stability.xml"); + gds = readSimXMLFile(fileName); + + auto* bus1 = dynamic_cast(gds->find("bus1")); + auto* bus2 = dynamic_cast(gds->find("bus2")); + ASSERT_NE(bus1, nullptr); + ASSERT_NE(bus2, nullptr); + + auto* primary = dynamic_cast(bus1->getGen(0)); + auto* secondary = dynamic_cast(bus2->getGen(0)); + ASSERT_NE(primary, nullptr); + ASSERT_NE(secondary, nullptr); + + auto* governor = new governors::GovernorIeeeSimple("shared_governor"); + primary->add(governor); + secondary->add(new genmodels::GenModel6()); + ASSERT_NE(dynamic_cast(secondary->find("genmodel")), nullptr); + secondary->set("mechanical_power_output", 0.0); + secondary->set("mechanical_power_source", fullObjectName(governor)); + EXPECT_EQ(secondary->getMechanicalPowerSource(), nullptr); + + ASSERT_EQ(gds->dynInitialize(), 0); + EXPECT_EQ(secondary->getMechanicalPowerSource(), governor); + EXPECT_EQ(secondary->getMechanicalPowerOutput(), 0); +} + TEST_F(GeneratorTests, GenTestRemote) { std::string fileName = std::string(GEN_TEST_DIRECTORY "test_gen_remote.xml"); diff --git a/test/componentTests/testGovernors.cpp b/test/componentTests/testGovernors.cpp index 85985c14..a3233da3 100644 --- a/test/componentTests/testGovernors.cpp +++ b/test/componentTests/testGovernors.cpp @@ -7,11 +7,17 @@ #include "../gtestHelper.h" #include "core/ObjectFactory.hpp" #include "griddyn/Generator.h" +#include "griddyn/GridBus.h" +#include "griddyn/generators/DynamicGenerator.h" +#include "griddyn/genmodels/GenModel6.h" +#include "griddyn/governors/GovernorIeeeG1.h" #include "griddyn/governors/GovernorTgov1.h" #include "griddyn/simulation/Diagnostics.h" #include "utilities/MatrixDataSparse.hpp" +#include #include #include +#include #include #include #include @@ -34,8 +40,136 @@ void configureTgov1(governors::GovernorTgov1& governor) governor.set("t3", 2.1); governor.set("dt", 0.1); } + +void configureIeeeG1(governors::GovernorIeeeG1& governor) +{ + governor.set("k", 20.0); + governor.set("t1", 0.2); + governor.set("t2", 0.05); + governor.set("t3", 0.1); + governor.set("uo", 0.3); + governor.set("uc", -0.25); + governor.set("pmax", 1.2); + governor.set("pmin", 0.1); + governor.set("t4", 0.4); + governor.set("k1", 0.3); + governor.set("k2", 0.1); + governor.set("t5", 0.0); + governor.set("k3", 0.2); + governor.set("k4", 0.1); + governor.set("t6", 0.5); + governor.set("k5", 0.1); + governor.set("k6", 0.05); + governor.set("t7", 0.2); + governor.set("k7", 0.1); + governor.set("k8", 0.05); +} } // namespace +TEST(GovernorModelTests, IeeeG1MatchesAndesInitializationAndPerturbedEquations) +{ + governors::GovernorIeeeG1 governor; + configureIeeeG1(governor); + governor.dynInitializeA(0.0, 0); + governor.setOutputInitializationTarget(governors::GovernorIeeeG1::lpOutput, 0.3); + + IOdata inputs{1.0, 0.7}; + IOdata desiredOutput{0.7}; + IOdata fieldSet(2, 0.0); + governor.dynInitializeB(inputs, desiredOutput, fieldSet); + + // Frozen ANDES initializes the lead-lag state to zero and the valve and + // all non-bypassed turbine stages to tm0 + tm02. T5=0 is an exact bypass. + const auto& initialized = governor.getStates(); + const std::vector expectedInitial{0.7, 0.3, 0.0, 1.0, 1.0, 1.0, 1.0}; + ASSERT_EQ(initialized.size(), expectedInitial.size()); + for (std::size_t index = 0; index < expectedInitial.size(); ++index) { + EXPECT_NEAR(initialized[index], expectedInitial[index], 1e-14) << index; + } + EXPECT_DOUBLE_EQ(fieldSet[govpSetInLocation], 1.0); + + // Local order is [PHP, PLP, LL_x, valve, L4, L6, L7]. These references + // are a direct evaluation of frozen ANDES IEEEG1 away from equilibrium. + std::vector state{0.69, 0.31, 0.02, 0.90, 0.85, 0.80, 0.75}; + std::vector stateDerivative(state.size(), 0.0); + governor.setState(0.0, state.data(), stateDerivative.data(), cLocalSolverMode); + inputs[govOmegaInLocation] = 0.99; + + std::vector derivative(state.size(), 0.0); + governor.derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode); + EXPECT_NEAR(derivative[2], -0.05, 1e-14); + EXPECT_NEAR(derivative[3], 0.30, 1e-14); + EXPECT_NEAR(derivative[4], 0.125, 1e-14); + EXPECT_NEAR(derivative[5], 0.10, 1e-14); + EXPECT_NEAR(derivative[6], 0.25, 1e-14); + + std::vector residual(state.size(), 0.0); + governor.residual(inputs, emptyStateData, residual.data(), cLocalSolverMode); + EXPECT_NEAR(residual[0], -0.11, 1e-14); + EXPECT_NEAR(residual[1], -0.0625, 1e-14); + for (std::size_t index = 2; index < state.size(); ++index) { + EXPECT_NEAR(residual[index], derivative[index], 1e-14) << index; + } +} + +TEST(GovernorModelTests, IeeeG1RateAndAntiWindupLimitTransitions) +{ + governors::GovernorIeeeG1 governor; + configureIeeeG1(governor); + governor.dynInitializeA(0.0, 0); + governor.setRootOffset(0, cLocalSolverMode); + governor.setOutputInitializationTarget(governors::GovernorIeeeG1::lpOutput, 0.3); + + // A 10% underspeed drives the valve rate upward at Pmax, so the + // non-windup limiter must hold the valve state. + IOdata inputs{0.90, 0.7}; + IOdata fieldSet(2, 0.0); + governor.dynInitializeB(inputs, {0.7}, fieldSet); + + std::vector state{0.7, 0.3, 0.0, 1.2, 1.0, 1.0, 1.0}; + std::vector stateDerivative(state.size(), 0.0); + governor.setState(0.0, state.data(), stateDerivative.data(), cLocalSolverMode); + + std::array roots{}; + governor.rootTest(inputs, emptyStateData, roots.data(), cLocalSolverMode); + EXPECT_DOUBLE_EQ(roots[1], 0.0); + governor.rootTrigger(0.0, inputs, {1, 1}, cLocalSolverMode); + + std::vector derivative(state.size(), 0.0); + governor.derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode); + EXPECT_DOUBLE_EQ(derivative[3], 0.0); + + inputs[govOmegaInLocation] = 1.10; + governor.rootTest(inputs, emptyStateData, roots.data(), cLocalSolverMode); + EXPECT_LT(roots[1], 0.0); + governor.rootTrigger(0.0, inputs, {1, 1}, cLocalSolverMode); + governor.derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode); + EXPECT_DOUBLE_EQ(derivative[3], -0.25); +} + +TEST(GovernorModelTests, IeeeG1FactoryCloneAndParameterValidation) +{ + auto factory = CoreObjectFactory::instance(); + std::unique_ptr object(factory->createObject("governor", "ieeeg1")); + auto* governor = dynamic_cast(object.get()); + ASSERT_NE(governor, nullptr); + configureIeeeG1(*governor); + + std::unique_ptr clonedObject(governor->clone()); + auto* clone = dynamic_cast(clonedObject.get()); + ASSERT_NE(clone, nullptr); + EXPECT_DOUBLE_EQ(clone->get("uo"), 0.3); + EXPECT_DOUBLE_EQ(clone->get("t7"), 0.2); + EXPECT_DOUBLE_EQ(clone->get("k8"), 0.05); + + EXPECT_ANY_THROW(governor->set("t3", 0.0)); + EXPECT_ANY_THROW(governor->set("uo", -0.1)); + EXPECT_ANY_THROW(governor->set("uc", 0.1)); + EXPECT_ANY_THROW(governor->set("k4", -0.1)); + governor->set("pmax", 0.05); + EXPECT_ANY_THROW(governor->dynInitializeA(0.0, 0)); +} + TEST(GovernorModelTests, Tgov1MatchesAndesInitializationAndPerturbedEquations) { governors::GovernorTgov1 governor; @@ -224,3 +358,68 @@ TEST_F(GovernorTests, Tgov1AnalyticJacobianMatchesFiniteDifferences) ASSERT_EQ(gds->dynInitialize(), 0); EXPECT_EQ(runJacobianCheck(gds, cDaeSolverMode), 0); } + +TEST_F(GovernorTests, IeeeG1AnalyticJacobianMatchesFiniteDifferences) +{ + const std::string fileName = std::string(GOVERNOR_TEST_DIRECTORY "test_gov_stability.xml"); + + GridDynSimulation::resetObjectCounters(); + gds = readSimXMLFile(fileName); + auto* generator = dynamic_cast(gds->findByUserID("gen", 2)); + ASSERT_NE(generator, nullptr); + + auto* governor = new governors::GovernorIeeeG1(); + configureIeeeG1(*governor); + governor->set("pmax", 2.0); + governor->set("pmin", 0.0); + for (const auto* coefficient : {"k2", "k4", "k6", "k8"}) { + governor->set(coefficient, 0.0); + } + generator->add(governor); + + ASSERT_EQ(gds->dynInitialize(), 0); + EXPECT_EQ(runResidualCheck(gds, cDaeSolverMode), 0); + EXPECT_EQ(runJacobianCheck(gds, cDaeSolverMode), 0); +} + +TEST_F(GovernorTests, IeeeG1CouplesMixedGeneratorModelsExactlyOnce) +{ + const std::string fileName = std::string(GOVERNOR_TEST_DIRECTORY "test_gov_stability.xml"); + + GridDynSimulation::resetObjectCounters(); + gds = readSimXMLFile(fileName); + auto* primary = dynamic_cast(gds->findByUserID("gen", 1)); + auto* secondary = dynamic_cast(gds->findByUserID("gen", 2)); + ASSERT_NE(primary, nullptr); + ASSERT_NE(secondary, nullptr); + + // The lossless fixture dispatches 0.3 pu from gen1 and 1.2 pu from gen2. + // Use different synchronous-machine classes to verify that the shared + // governor connection is independent of the generator-model type. + secondary->add(new genmodels::GenModel6()); + auto* unusedLocalGovernor = new governors::GovernorTgov1(); + configureTgov1(*unusedLocalGovernor); + unusedLocalGovernor->set("pmax", 2.0); + unusedLocalGovernor->set("pmin", 0.0); + secondary->add(unusedLocalGovernor); + auto* governor = new governors::GovernorIeeeG1("cross_compound_ieeeg1"); + configureIeeeG1(*governor); + governor->set("pmax", 2.0); + governor->set("pmin", 0.0); + for (const auto* coefficient : {"k1", "k2", "k3", "k4", "k5", "k6", "k7", "k8"}) { + governor->set(coefficient, 0.0); + } + governor->set("k1", 0.2); + governor->set("k2", 0.8); + primary->add(governor); + secondary->setMechanicalPowerSource(governor, governors::GovernorIeeeG1::lpOutput); + + ASSERT_EQ(gds->dynInitialize(), 0); + EXPECT_EQ(primary->getMechanicalPowerSource(), governor); + EXPECT_EQ(primary->getMechanicalPowerOutput(), governors::GovernorIeeeG1::hpOutput); + EXPECT_EQ(secondary->getMechanicalPowerSource(), governor); + EXPECT_EQ(secondary->getMechanicalPowerOutput(), governors::GovernorIeeeG1::lpOutput); + EXPECT_EQ(secondary->find("governor"), unusedLocalGovernor); + EXPECT_EQ(runResidualCheck(gds, cDaeSolverMode), 0); + EXPECT_EQ(runJacobianCheck(gds, cDaeSolverMode), 0); +} diff --git a/test/test_files/andes_tests/ieee14_ieeeg1.dyr b/test/test_files/andes_tests/ieee14_ieeeg1.dyr new file mode 100644 index 00000000..f6064e80 --- /dev/null +++ b/test/test_files/andes_tests/ieee14_ieeeg1.dyr @@ -0,0 +1,2 @@ + 2 'IEEEG1' 1 0 0 18.0 0.12 0.03 0.22 0.42 -0.31 1.25 0.15 + 0.14 0.20 0.0 0.25 0.30 0.0 0.36 0.10 0.0 0.47 0.40 0.0 / diff --git a/test/test_files/andes_tests/ieee14_ieeeg1_cross.dyr b/test/test_files/andes_tests/ieee14_ieeeg1_cross.dyr new file mode 100644 index 00000000..77e9259e --- /dev/null +++ b/test/test_files/andes_tests/ieee14_ieeeg1_cross.dyr @@ -0,0 +1,2 @@ + 1 'IEEEG1' 1 2 1 20.0 0.20 0.05 0.10 0.30 -0.25 2.0 0.0 + 0.40 0.20 0.80 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 /