diff --git a/docs/developer-guide/andes-compatibility.md b/docs/developer-guide/andes-compatibility.md index 4a1d7691..2e754b99 100644 --- a/docs/developer-guide/andes-compatibility.md +++ b/docs/developer-guide/andes-compatibility.md @@ -64,31 +64,94 @@ GridDyn's current DYR reader recognizes only `GENROU`, `ESDC1A`, `EXDC2`, `TGOV1`, and `SEXS`. Recognition alone is not compatibility: the existing adapters still have attachment, parameter, equation, and validation gaps. -| PSS/e DYR record(s) | ANDES destination | Existing GridDyn candidate | Current compatibility and required action | -| ------------------------------------------------------------------------- | ----------------------- | ----------------------------------------------------------------------------- | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | -| `GENCLS` | `GENCLS` | `GenModelClassical` | **Partial.** Add a DYR adapter, verify machine/system-base conversion, initialization, and rotor trajectory. | -| `GENSAL` | `GENROU` | No verified exact equivalent | **Planned.** ANDES explicitly marks this as an approximation and fills GENROU q-axis values from d-axis data. Decide whether compatibility mode should reproduce that conversion before considering a distinct salient-pole model. | -| `GENROU` | `GENROU` | Registered `GenModelGENROU`; DYR uses the dedicated model | **Partial.** The adapter and dedicated model match all eight initialized machine states, field voltage, and mechanical power for five IEEE 14-bus machines, including quadratic saturation. Exact alphanumeric machine-ID lookup, arbitrary machine-base cases, native ANDES import, and disturbed-trajectory validation remain. | -| `IEEEVC` | `IEEEVC` | None | **No direct analogue.** Add a voltage-compensator input to the exciter path, including `rc`/`xc` terminal-current compensation. | -| `SEXS` | `SEXS` | `ExciterSEXS` | **Partial.** The adapter exists; verify lead-lag convention, limits, initialization, and trajectories. | -| `ESDC1A` | `ESDC1A` | `ExciterDC1A` / `ExciterIEEEtype1` | **Partial.** The current adapter ignores `TR`, switch behavior, and `E1`/`SE1`/`E2`/`SE2`; its fallback to type 1 when `TB` is zero must be compared with ANDES. | -| `EXDC2` | `EXDC2` | `ExciterDC2A` | **Partial.** The adapter exists but ignores `TR`, switch behavior, and saturation points. Confirm whether `DC2A` equations exactly represent ANDES `EXDC2`. | -| `IEEET1`, `IEEET3`, `ESDC2A` | Same named ANDES models | `ExciterIEEEtype1`, `ExciterIEEEtype2`, and `ExciterDC2A` are only candidates | **Planned.** Perform equation and limiter audits before selecting an analogue; similar names are insufficient. | -| `ESST3A` | `ESST3A` | `ExciterESST3A` | **Implemented.** Exact PSS/e field mapping, ANDES-equation initialization, residual, and Jacobian coverage with GENROU; reduced-order synchronous generators use the documented controller-signal approximations. | -| `EXST1` | `EXST1` | `ExciterEXST1` | **Implemented.** PSS/e field order, captured initialization and perturbed-equation values, corrected regulator-output limiter/root tests, analytic-Jacobian finite-difference checks, and GENROU attachment are covered. The documented limiter-selector divergence from frozen ANDES is intentional. Zero `TR`, `TB`, or `TA` records and a captured disturbed trajectory remain unsupported. | -| `ESST1A`, `ESAC1A`, `AC8B`, `EXAC1`, `EXAC2`, `EXAC4`, `IEEEX1`, `ESST4B` | Same named ANDES models | None exact | **No direct analogue.** Implement model-specific exciter blocks and DYR schemas, then add initialization and trajectory tests. | -| `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` | `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. | -| `ST2CUT` | `ST2CUT` | `StabilizerST2CUT` | **Partial.** Exact frozen-ANDES local modes 0/1/3/4/5, dual transducers, washout/lag, three lead-lag stages, output limits, voltage gating, DYR mapping, generator/exciter coupling, and dynamic load-step coverage are implemented. Remote `BUSR` inputs and ANDES modes 2/6 require cross-bus/frequency-derivative measurement routing and are rejected rather than approximated. | -| `IEEEST` | `IEEEST` | `StabilizerIEEEST` | **Partial.** Exact frozen-ANDES local modes 0/1/3/4/5, filter/lead-lag/washout zero-bypass semantics, output limits, voltage gate, DYR mapping, and GENROU/exciter coupling are implemented. Modes 2/6 and nonzero `BUSR` require unavailable cross-bus/frequency-derivative measurement routing and are rejected rather than approximated. A captured ANDES trajectory remains open. | -| `REGCA1`, `REECA1`, `REECB1`, `REPCA1` | Same named ANDES models | `GenModelInverter` is not equivalent | **No direct analogue.** Add the coordinated renewable generator, electrical-control, and plant-control chain rather than flattening these records into the generic inverter. | -| `WTDTA1`, `WTARA1`, `WTPTA1`, `WTTQA1` | Same named ANDES models | None | **No direct analogue.** Add drive-train, aerodynamic, pitch, and torque-control submodels with their shared interfaces. | -| `Toggle`, `Fault` | ANDES event models | GridDyn event/action and fault mechanisms | **Partial conceptually.** Define DYR record schemas and translate target resolution, timing, status changes, and fault clearing semantics; no adapter exists. | +| PSS/e DYR record(s) | ANDES destination | Existing GridDyn candidate | Current compatibility and required action | +| ---------------------------------------------- | ----------------------- | ----------------------------------------------------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ | +| `GENCLS` | `GENCLS` | `GenModelClassical` | **Partial.** Add a DYR adapter, verify machine/system-base conversion, initialization, and rotor trajectory. | +| `GENSAL` | `GENROU` | No verified exact equivalent | **Planned.** ANDES explicitly marks this as an approximation and fills GENROU q-axis values from d-axis data. Decide whether compatibility mode should reproduce that conversion before considering a distinct salient-pole model. | +| `GENROU` | `GENROU` | Registered `GenModelGENROU`; DYR uses the dedicated model | **Partial.** The adapter and dedicated model match all eight initialized machine states, field voltage, and mechanical power for five IEEE 14-bus machines, including quadratic saturation. Exact alphanumeric machine-ID lookup, arbitrary machine-base cases, native ANDES import, and disturbed-trajectory validation remain. | +| `IEEEVC` | `IEEEVC` | None | **No direct analogue.** Add a voltage-compensator input to the exciter path, including `rc`/`xc` terminal-current compensation. | +| `SEXS` | `SEXS` | `ExciterSEXS` | **Partial.** The adapter exists; verify lead-lag convention, limits, initialization, and trajectories. | +| `ESDC1A` | `ESDC1A` | `ExciterDC1A` / `ExciterIEEEtype1` | **Partial.** The current adapter ignores `TR`, switch behavior, and `E1`/`SE1`/`E2`/`SE2`; its fallback to type 1 when `TB` is zero must be compared with ANDES. | +| `EXDC2` | `EXDC2` | `ExciterDC2A` | **Partial.** The adapter exists but ignores `TR`, switch behavior, and saturation points. Confirm whether `DC2A` equations exactly represent ANDES `EXDC2`. | +| `IEEET1`, `IEEET3`, `ESDC2A` | Same named ANDES models | `ExciterIEEEtype1`, `ExciterIEEEtype2`, and `ExciterDC2A` are only candidates | **Planned.** Perform equation and limiter audits before selecting an analogue; similar names are insufficient. | +| `ESST3A` | `ESST3A` | `ExciterESST3A` | **Implemented.** Exact PSS/e field mapping, ANDES-equation initialization, residual, and Jacobian coverage with GENROU; reduced-order synchronous generators use the documented controller-signal approximations. | +| `EXST1` | `EXST1` | `ExciterEXST1` | **Implemented.** PSS/e field order, captured initialization and perturbed-equation values, corrected regulator-output limiter/root tests, analytic-Jacobian finite-difference checks, and GENROU attachment are covered. The documented limiter-selector divergence from frozen ANDES is intentional. Zero `TR`, `TB`, or `TA` records and a captured disturbed trajectory remain unsupported. | +| `EXAC1`, `EXAC2`, `EXAC4` | Same named ANDES models | `ExciterEXAC1`, `ExciterEXAC2`, `ExciterEXAC4` | **Implemented.** Exact PSS/e DYR field mapping, initialization, residual/Jacobian, limiter/root, GENROU attachment, and IEEE-14 load-step coverage are present. EXAC1/EXAC2 intentionally use the sensed-voltage transducer output instead of the disconnected frozen-ANDES path documented below. Positive `TR`, `TB`, `TA`, `TE` (where applicable), and `TF` are required; zero-time algebraic bypass records remain unsupported. | +| `ESST1A`, `ESAC1A`, `AC8B`, `IEEEX1`, `ESST4B` | Same named ANDES models | None exact | **No direct analogue.** Implement model-specific exciter blocks and DYR schemas, then add initialization and trajectory tests. | +| `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` | `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. | +| `ST2CUT` | `ST2CUT` | `StabilizerST2CUT` | **Partial.** Exact frozen-ANDES local modes 0/1/3/4/5, dual transducers, washout/lag, three lead-lag stages, output limits, voltage gating, DYR mapping, generator/exciter coupling, and dynamic load-step coverage are implemented. Remote `BUSR` inputs and ANDES modes 2/6 require cross-bus/frequency-derivative measurement routing and are rejected rather than approximated. | +| `IEEEST` | `IEEEST` | `StabilizerIEEEST` | **Partial.** Exact frozen-ANDES local modes 0/1/3/4/5, filter/lead-lag/washout zero-bypass semantics, output limits, voltage gate, DYR mapping, and GENROU/exciter coupling are implemented. Modes 2/6 and nonzero `BUSR` require unavailable cross-bus/frequency-derivative measurement routing and are rejected rather than approximated. A captured ANDES trajectory remains open. | +| `REGCA1`, `REECA1`, `REECB1`, `REPCA1` | Same named ANDES models | `GenModelInverter` is not equivalent | **No direct analogue.** Add the coordinated renewable generator, electrical-control, and plant-control chain rather than flattening these records into the generic inverter. | +| `WTDTA1`, `WTARA1`, `WTPTA1`, `WTTQA1` | Same named ANDES models | None | **No direct analogue.** Add drive-train, aerodynamic, pitch, and torque-control submodels with their shared interfaces. | +| `Toggle`, `Fault` | ANDES event models | GridDyn event/action and fault mechanisms | **Partial conceptually.** Define DYR record schemas and translate target resolution, timing, status changes, and fault clearing semantics; no adapter exists. | + +### EXAC1 and EXAC2 voltage-transducer compatibility decision + +The frozen ANDES v2.0.0 `EXAC1Model` constructs the voltage-transducer lag +`LG = Lag(v, TR)`, but its regulator-input equation uses raw terminal voltage +`v` rather than `LG_y`. Consequently `TR` creates an unconsumed state and has +no effect on the frozen ANDES EXAC1 response. `EXAC2Model` derives from the +same implementation and inherits this behavior. This appears inconsistent +with the PSS/E/IEEE AC1 block diagram, where the sensed terminal voltage is +the transducer output. + +GridDyn will implement the physically intended connection for both models: + +@f[ +T_R\dot V_m=V_t-V_m,\qquad V_i=V_{ref}-V_m-V_F. +@f] + +`ExciterEXAC1` and `ExciterEXAC2` Doxygen comments and tests name this +intentional difference. When `TR` is nonzero, tests derived solely from frozen +ANDES values must be adjusted or separately labeled; PSS/E or a second +independent implementation is required to validate the corrected trajectory. +This is a documented compatibility correction, not a claim that ANDES +behavior is exact. + +#### Copy-ready ANDES issue draft: EXAC1/EXAC2 `TR` transducer is disconnected + +```text +Title: EXAC1 and EXAC2 do not apply the TR voltage-transducer output + +In ANDES v2.0.0, `andes/models/exciter/exac1.py` creates the EXAC1 voltage +transducer as + + self.LG = Lag(self.v, T=self.TR, K=1, info='Voltage transducer') + +but the regulator summing input is defined using raw terminal voltage: + + self.vi.e_str = 'ue * (-v + vref - WF_y - vi)' + +`LG_y` is not referenced by the EXAC1 regulator path. Therefore changing TR +adds an unconsumed lag state but does not change the EXAC1 response. EXAC2 +inherits EXAC1Model, so it has the same issue. + +Expected AC1/PSS/E block-diagram connection: + + TR * d(LG_y)/dt = v - LG_y + vi = vref - LG_y - WF_y + +Suggested minimal correction: + + self.vi.e_str = 'ue * (-LG_y + vref - WF_y - vi)' + self.vi.v_str = '-LG_y + vref' + +Suggested regression: run two otherwise identical EXAC1 cases with TR = 0 +and TR > 0 after a terminal-voltage or reference-voltage step. Before the +fix, their regulator trajectories are identical apart from the unused LG +state; after the fix, the nonzero-TR case should show the expected sensing +delay. Repeat the same regression for EXAC2. + +This report concerns the regulator sensing path only. It does not propose a +change to EXAC1/EXAC2 saturation, FEX, field-current feedback, or limiter +equations. +``` ### DYR reader infrastructure plan diff --git a/src/fileInput/gridDynReadDYR.cpp b/src/fileInput/gridDynReadDYR.cpp index ba903da0..9ab4c73c 100644 --- a/src/fileInput/gridDynReadDYR.cpp +++ b/src/fileInput/gridDynReadDYR.cpp @@ -34,6 +34,9 @@ namespace { void loadESDC1A(CoreObject* parentObject, stringVec& tokens); void loadESST3A(CoreObject* parentObject, stringVec& tokens); void loadEXST1(CoreObject* parentObject, stringVec& tokens); + void loadEXAC1(CoreObject* parentObject, stringVec& tokens); + void loadEXAC2(CoreObject* parentObject, stringVec& tokens); + void loadEXAC4(CoreObject* parentObject, stringVec& tokens); void loadTGOV1(CoreObject* parentObject, stringVec& tokens); void loadIEEEG1(CoreObject* parentObject, stringVec& tokens); void loadIEEEST(CoreObject* parentObject, stringVec& tokens); @@ -86,6 +89,12 @@ void loadDyr(CoreObject* parentObject, loadESST3A(parentObject, lineTokens); } else if (type == "'EXST1'") { loadEXST1(parentObject, lineTokens); + } else if (type == "'EXAC1'") { + loadEXAC1(parentObject, lineTokens); + } else if (type == "'EXAC2'") { + loadEXAC2(parentObject, lineTokens); + } else if (type == "'EXAC4'") { + loadEXAC4(parentObject, lineTokens); } else if (type == "'EXDC2'") { loadEXDC2(parentObject, lineTokens); } else if (type == "'TGOV1'") { @@ -235,6 +244,96 @@ namespace { gen->add(exciterModel); } + void loadEXAC1(CoreObject* parentObject, stringVec& tokens) + { + const int busId = std::stoi(tokens[0]); + const auto* bus = static_cast(parentObject->findByUserID("bus", busId)); + auto* gen = bus->getGen(std::stoi(tokens[2]) - 1); + const auto params = gmlc::utilities::str2vector(tokens, kNullVal); + auto* exciter = + static_cast(CoreObjectFactory::instance()->createObject("exciter", "exac1")); + // Exact ANDES psse-dyr.yaml order: TR, TB, TC, KA, TA, VRMAX, + // VRMIN, TE, KF, TF, KC, KD, KE, E1, SE1, E2, SE2. + exciter->set("tr", params[3]); + exciter->set("tb", params[4]); + exciter->set("tc", params[5]); + exciter->set("ka", params[6]); + exciter->set("ta", params[7]); + exciter->set("vrmax", params[8]); + exciter->set("vrmin", params[9]); + exciter->set("te", params[10]); + exciter->set("kf", params[11]); + exciter->set("tf", params[12]); + exciter->set("kc", params[13]); + exciter->set("kd", params[14]); + exciter->set("ke", params[15]); + exciter->set("e1", params[16]); + exciter->set("se1", params[17]); + exciter->set("e2", params[18]); + exciter->set("se2", params[19]); + gen->add(exciter); + } + + void loadEXAC2(CoreObject* parentObject, stringVec& tokens) + { + const int busId = std::stoi(tokens[0]); + const auto* bus = static_cast(parentObject->findByUserID("bus", busId)); + auto* gen = bus->getGen(std::stoi(tokens[2]) - 1); + const auto params = gmlc::utilities::str2vector(tokens, kNullVal); + auto* exciter = + static_cast(CoreObjectFactory::instance()->createObject("exciter", "exac2")); + // Exact ANDES psse-dyr.yaml order: TR, TB, TC, KA, TA, VAMAX, + // VAMIN, KB, VRMAX, VRMIN, TE, KL, KH, KF, TF, KC, KD, KE, VLR, + // E1, SE1, E2, SE2. + exciter->set("tr", params[3]); + exciter->set("tb", params[4]); + exciter->set("tc", params[5]); + exciter->set("ka", params[6]); + exciter->set("ta", params[7]); + exciter->set("vamax", params[8]); + exciter->set("vamin", params[9]); + exciter->set("kb", params[10]); + exciter->set("vrmax", params[11]); + exciter->set("vrmin", params[12]); + exciter->set("te", params[13]); + exciter->set("kl", params[14]); + exciter->set("kh", params[15]); + exciter->set("kf", params[16]); + exciter->set("tf", params[17]); + exciter->set("kc", params[18]); + exciter->set("kd", params[19]); + exciter->set("ke", params[20]); + exciter->set("vlr", params[21]); + exciter->set("e1", params[22]); + exciter->set("se1", params[23]); + exciter->set("e2", params[24]); + exciter->set("se2", params[25]); + gen->add(exciter); + } + + void loadEXAC4(CoreObject* parentObject, stringVec& tokens) + { + const int busId = std::stoi(tokens[0]); + const auto* bus = static_cast(parentObject->findByUserID("bus", busId)); + auto* gen = bus->getGen(std::stoi(tokens[2]) - 1); + const auto params = gmlc::utilities::str2vector(tokens, kNullVal); + auto* exciter = + static_cast(CoreObjectFactory::instance()->createObject("exciter", "exac4")); + // Exact ANDES psse-dyr.yaml order: TR, VIMAX, VIMIN, TC, TB, KA, + // TA, VRMAX, VRMIN, KC. + exciter->set("tr", params[3]); + exciter->set("vimax", params[4]); + exciter->set("vimin", params[5]); + exciter->set("tc", params[6]); + exciter->set("tb", params[7]); + exciter->set("ka", params[8]); + exciter->set("ta", params[9]); + exciter->set("vrmax", params[10]); + exciter->set("vrmin", params[11]); + exciter->set("kc", params[12]); + gen->add(exciter); + } + void loadEXDC2(CoreObject* parentObject, stringVec& tokens) { const int busId = std::stoi(tokens[0]); diff --git a/src/griddyn/CMakeLists.txt b/src/griddyn/CMakeLists.txt index ab08cb71..3fd7acc6 100644 --- a/src/griddyn/CMakeLists.txt +++ b/src/griddyn/CMakeLists.txt @@ -86,6 +86,9 @@ set(exciter_headers exciters/ExciterDC1A.h exciters/ExciterDC2A.h exciters/ExciterESST3A.h + exciters/ExciterEXAC1.h + exciters/ExciterEXAC2.h + exciters/ExciterEXAC4.h exciters/ExciterEXST1.h exciters/ExciterIEEEtype1.h exciters/ExciterIEEEtype2.h @@ -98,6 +101,9 @@ set(exciter_sources exciters/ExciterDC1A.cpp exciters/ExciterDC2A.cpp exciters/ExciterESST3A.cpp + exciters/ExciterEXAC1.cpp + exciters/ExciterEXAC2.cpp + exciters/ExciterEXAC4.cpp exciters/ExciterEXST1.cpp exciters/ExciterIEEEtype1.cpp exciters/ExciterIEEEtype2.cpp diff --git a/src/griddyn/exciters/Exciter.cpp b/src/griddyn/exciters/Exciter.cpp index 2566e6e5..25d02d3d 100644 --- a/src/griddyn/exciters/Exciter.cpp +++ b/src/griddyn/exciters/Exciter.cpp @@ -13,6 +13,9 @@ #include "ExciterDC1A.h" #include "ExciterDC2A.h" #include "ExciterESST3A.h" +#include "ExciterEXAC1.h" +#include "ExciterEXAC2.h" +#include "ExciterEXAC4.h" #include "ExciterEXST1.h" #include "ExciterIEEEtype1.h" #include "ExciterIEEEtype2.h" @@ -36,6 +39,9 @@ namespace exciters { static ChildTypeFactory gfeDc1a("exciter", "dc1a"); // NOLINT static ChildTypeFactory gfeDc2a("exciter", "dc2a"); // NOLINT static ChildTypeFactory gfeEsst3a("exciter", "esst3a"); // NOLINT + static ChildTypeFactory gfeExac1("exciter", "exac1"); // NOLINT + static ChildTypeFactory gfeExac2("exciter", "exac2"); // NOLINT + static ChildTypeFactory gfeExac4("exciter", "exac4"); // NOLINT static ChildTypeFactory gfeExst1("exciter", "exst1"); // NOLINT static ChildTypeFactory gfeType1("exciter", "type1"); // NOLINT static TypeFactory gfeDefault( // NOLINT diff --git a/src/griddyn/exciters/ExciterEXAC1.cpp b/src/griddyn/exciters/ExciterEXAC1.cpp new file mode 100644 index 00000000..3fc5bd4a --- /dev/null +++ b/src/griddyn/exciters/ExciterEXAC1.cpp @@ -0,0 +1,610 @@ +/* + * 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 "ExciterEXAC1.h" + +#include "core/CoreExceptions.h" +#include "core/CoreObjectTemplates.hpp" +#include "utilities/MatrixData.hpp" +#include +#include +#include +#include + +namespace griddyn::exciters { +// The published transfer-function equations are intentionally kept in their +// conventional form. Parentheses at every multiplication obscure them. +// NOLINTBEGIN(readability-math-missing-parentheses) +namespace { + constexpr index_t voltageMeasurementState = 0; + constexpr index_t leadLagState = 1; + constexpr index_t regulatorState = 2; + constexpr index_t exciterState = 3; + constexpr index_t washoutState = 4; + constexpr double limitTolerance = 1e-7; + // These decimal values are the specified PSS/E FEX curve coefficients, + // not approximations of unrelated mathematical constants. + constexpr double lowCurrentSlope = 0.577; // NOLINT(modernize-use-std-numbers) + constexpr double highCurrentSlope = 1.732; // NOLINT(modernize-use-std-numbers) + + struct RectifierEvaluation { + double mFactor; + double mDerivative; + }; + + RectifierEvaluation rectifier(double normalizedCurrent) + { + if (normalizedCurrent <= 0.0) { + return {.mFactor = 1.0, .mDerivative = 0.0}; + } + if (normalizedCurrent <= 0.433) { + return {.mFactor = 1.0 - lowCurrentSlope * normalizedCurrent, + .mDerivative = -lowCurrentSlope}; + } + if (normalizedCurrent <= 0.75) { + const double factor = + std::sqrt(std::max(0.0, 0.75 - normalizedCurrent * normalizedCurrent)); + return {.mFactor = factor, + .mDerivative = (factor > 0.0) ? -normalizedCurrent / factor : 0.0}; + } + if (normalizedCurrent <= 1.0) { + return {.mFactor = highCurrentSlope * (1.0 - normalizedCurrent), + .mDerivative = -highCurrentSlope}; + } + return {.mFactor = 0.0, .mDerivative = 0.0}; + } +} // namespace + +ExciterEXAC1::ExciterEXAC1(const std::string& objName): Exciter(objName) +{ + m_inputSize = exciterInputCount; + Ka = 80.0; + Ta = 0.04; + Vrmax = 8.0; + Vrmin = 0.0; +} + +CoreObject* ExciterEXAC1::clone(CoreObject* obj) const +{ + auto* clone = cloneBase(this, obj); + if (clone == nullptr) { + return obj; + } + clone->Tr = Tr; + clone->Tb = Tb; + clone->Tc = Tc; + clone->Te = Te; + clone->Kf = Kf; + clone->Tf = Tf; + clone->Kc = Kc; + clone->Kd = Kd; + clone->Ke = Ke; + clone->E1 = E1; + clone->Se1 = Se1; + clone->E2 = E2; + clone->Se2 = Se2; + clone->saturation = saturation; + return clone; +} + +void ExciterEXAC1::dynObjectInitializeA(CoreTime /*time0*/, std::uint32_t /*flags*/) +{ + if (!std::isfinite(Tr) || !std::isfinite(Tb) || !std::isfinite(Tc) || !std::isfinite(Ka) || + !std::isfinite(Ta) || !std::isfinite(Vrmax) || !std::isfinite(Vrmin) || + !std::isfinite(Te) || !std::isfinite(Kf) || !std::isfinite(Tf) || !std::isfinite(Kc) || + !std::isfinite(Kd) || !std::isfinite(Ke) || !std::isfinite(E1) || !std::isfinite(Se1) || + !std::isfinite(E2) || !std::isfinite(Se2) || (Tr <= 0.0) || (Tb <= 0.0) || (Ta <= 0.0) || + (Te <= 0.0) || (Tf <= 0.0) || (Ka <= 0.0) || (Vrmax < Vrmin) || + (regulatorUpperLimit() < regulatorLowerLimit())) { + throw InvalidParameterValue("EXAC1 gains, time constants, or limits"); + } + // ANDES disables ExcQuadSat when SE2 is zero. Saturation::QUADRATIC + // intentionally does not infer that convention, because its generic + // two-point fit permits no zero denominator. + if (Se2 == 0.0) { + saturation.setType(utilities::Saturation::SaturationType::NONE); + } else { + if ((E1 <= 0.0) || (E2 <= 0.0) || (Se1 <= 0.0) || (Se2 <= 0.0)) { + throw InvalidParameterValue("EXAC1 enabled saturation points"); + } + saturation.setType(utilities::Saturation::SaturationType::QUADRATIC); + saturation.setParam(E1, E1 * Se1, E2, E2 * Se2); + } + offsets.local().local.algSize = 1; + offsets.local().local.diffSize = 5; + offsets.local().local.algRoots = 1; + offsets.local().local.jacSize = 34; +} + +void ExciterEXAC1::dynObjectInitializeB(const IOdata& inputs, + const IOdata& desiredOutput, + IOdata& fieldSet) +{ + const double voltage = inputs[exciterVoltageInLocation]; + const double fieldCurrent = inputs[exciterXadIfdInLocation]; + if (!std::isfinite(voltage) || + ((!std::isfinite(fieldCurrent) || (std::abs(fieldCurrent) > 1e20)) && + ((Kc != 0.0) || (Kd != 0.0)))) { + throw InvalidParameterValue("EXAC1 requires compatible synchronous-machine signals"); + } + if (desiredOutput.empty() || !std::isfinite(desiredOutput[0])) { + throw InvalidParameterValue("EXAC1 initial field voltage"); + } + const double fieldVoltage = desiredOutput[0]; + double exciterVoltage = std::max(0.01, std::abs(fieldVoltage)); + for (int count = 0; count < 20; ++count) { + const double current = + ((Kc != 0.0) && (exciterVoltage != 0.0)) ? Kc * fieldCurrent / exciterVoltage : 0.0; + const auto fex = rectifier(current); + const double mismatch = exciterVoltage * fex.mFactor - fieldVoltage; + const double slope = fex.mFactor - fex.mDerivative * current; + if (std::abs(slope) < 1e-12) { + break; + } + exciterVoltage = std::max(1e-8, exciterVoltage - mismatch / slope); + if (std::abs(mismatch) < 1e-12) { + break; + } + } + const double finalCurrent = + ((Kc != 0.0) && (exciterVoltage != 0.0)) ? Kc * fieldCurrent / exciterVoltage : 0.0; + if (std::abs(exciterVoltage * rectifier(finalCurrent).mFactor - fieldVoltage) > 1e-7) { + throw InvalidParameterValue( + "EXAC1 initial field voltage is incompatible with rectifier loading"); + } + double* state = m_state.data() + 1; + state[voltageMeasurementState] = voltage; + state[exciterState] = exciterVoltage; + state[washoutState] = vfe(inputs, state); + state[regulatorState] = initialRegulatorState(state[washoutState]); + if ((state[regulatorState] < regulatorLowerLimit() - limitTolerance) || + (state[regulatorState] > regulatorUpperLimit() + limitTolerance)) { + throw InvalidParameterValue("EXAC1 initial regulator output outside limits"); + } + state[leadLagState] = state[regulatorState] / Ka; + m_state[0] = fieldVoltage; + const double setpointInput = inputs[exciterVsetInLocation] - 1.0; + vBias = voltage + referenceOffset(state[washoutState]) - Vref - setpointInput - + inputs[exciterVssInLocation]; + fieldSet[exciterVsetInLocation] = Vref; + std::fill(m_dstate_dt.begin(), m_dstate_dt.end(), 0.0); + updateLimitFlags(state); +} + +double ExciterEXAC1::referenceInput(const IOdata& inputs) const +{ + return Vref + vBias + inputs[exciterVsetInLocation] - 1.0 + inputs[exciterVssInLocation]; +} + +double ExciterEXAC1::vfe(const IOdata& inputs, const double state[]) const +{ + const double fieldFeedback = Ke * state[exciterState] + saturation(state[exciterState]); + return (Kd == 0.0) ? fieldFeedback : fieldFeedback + Kd * inputs[exciterXadIfdInLocation]; +} + +double ExciterEXAC1::rectifierFactor(const IOdata& inputs, double exciterVoltage) const +{ + if (Kc == 0.0) { + return 1.0; + } + return rectifier((exciterVoltage != 0.0) ? + Kc * inputs[exciterXadIfdInLocation] / exciterVoltage : + 0.0) + .mFactor; +} + +double ExciterEXAC1::fieldVoltage(const IOdata& inputs, const double state[]) const +{ + return state[exciterState] * rectifierFactor(inputs, state[exciterState]); +} + +double ExciterEXAC1::regulatorTarget(const IOdata& /*inputs*/, const double state[]) const +{ + return state[regulatorState]; +} + +double ExciterEXAC1::regulatorUpperLimit() const +{ + return Vrmax; +} +double ExciterEXAC1::regulatorLowerLimit() const +{ + return Vrmin; +} +double ExciterEXAC1::initialRegulatorState(double vfeValue) const +{ + return vfeValue; +} +double ExciterEXAC1::referenceOffset(double vfeValue) const +{ + return vfeValue / Ka; +} +void ExciterEXAC1::regulatorTargetDerivatives(const IOdata& /*inputs*/, + const double state[], + double& regulatorDerivative, + double& exciterDerivative, + double& fieldCurrentDerivative) const +{ + static_cast(state); + regulatorDerivative = 1.0; + exciterDerivative = 0.0; + fieldCurrentDerivative = 0.0; +} + +int ExciterEXAC1::regulatorLimitStatus(const double state[]) const +{ + if (state[regulatorState] >= regulatorUpperLimit()) { + return 1; + } + return (state[regulatorState] <= regulatorLowerLimit()) ? -1 : 0; +} + +bool ExciterEXAC1::updateLimitFlags(const double state[]) +{ + const int status = regulatorLimitStatus(state); + const bool changed = (opFlags[REGULATOR_LIMITED] != (status != 0)) || + (opFlags[REGULATOR_LIMIT_HIGH] != (status > 0)); + opFlags.set(REGULATOR_LIMITED, status != 0); + opFlags.set(REGULATOR_LIMIT_HIGH, status > 0); + return changed; +} + +void ExciterEXAC1::residual(const IOdata& inputs, + const StateData& stateData, + double resid[], + const SolverMode& sMode) +{ + const auto locations = offsets.getLocations(stateData, resid, sMode, this); + if (hasAlgebraic(sMode)) { + locations.destLoc[0] = + fieldVoltage(inputs, locations.diffStateLoc) - locations.algStateLoc[0]; + } + if (hasDifferential(sMode)) { + derivative(inputs, stateData, resid, sMode); + for (index_t index = 0; index < locations.diffSize; ++index) { + locations.destDiffLoc[index] -= locations.dstateLoc[index]; + } + } +} + +void ExciterEXAC1::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* derivativeValues = locations.destDiffLoc; + const double fieldFeedback = vfe(inputs, state); + const double input = referenceInput(inputs) - state[voltageMeasurementState] - + Kf * (fieldFeedback - state[washoutState]) / Tf; + const double leadOutput = state[leadLagState] + Tc * (input - state[leadLagState]) / Tb; + const double regulatorDerivative = (Ka * leadOutput - state[regulatorState]) / Ta; + const int status = regulatorLimitStatus(state); + derivativeValues[voltageMeasurementState] = + (inputs[exciterVoltageInLocation] - state[voltageMeasurementState]) / Tr; + derivativeValues[leadLagState] = (input - state[leadLagState]) / Tb; + derivativeValues[regulatorState] = ((status > 0) && (regulatorDerivative > 0.0)) || + ((status < 0) && (regulatorDerivative < 0.0)) ? + 0.0 : + regulatorDerivative; + derivativeValues[exciterState] = (regulatorTarget(inputs, state) - fieldFeedback) / Te; + derivativeValues[washoutState] = (fieldFeedback - state[washoutState]) / Tf; +} + +void ExciterEXAC1::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; + const double exciterVoltage = state[exciterState]; + const double normalizedCurrent = + (exciterVoltage != 0.0) ? Kc * inputs[exciterXadIfdInLocation] / exciterVoltage : 0.0; + const auto fex = rectifier(normalizedCurrent); + if (hasAlgebraic(sMode)) { + matrixData.assign(algOffset, algOffset, -1.0); + matrixData.assign(algOffset, + diffOffset + exciterState, + fex.mFactor - fex.mDerivative * normalizedCurrent); + matrixData.assignCheckCol(algOffset, + inputLocs[exciterXadIfdInLocation], + fex.mDerivative * Kc); + } + if (!hasDifferential(sMode)) { + return; + } + const double saturationSlope = Ke + saturation.deriv(exciterVoltage); + const double feedbackGain = Kf / Tf; + const double leadRatio = Tc / Tb; + const double regulatorDerivative = + (Ka * + (state[leadLagState] + + leadRatio * + (referenceInput(inputs) - state[voltageMeasurementState] - + feedbackGain * (vfe(inputs, state) - state[washoutState]) - + state[leadLagState])) - + state[regulatorState]) / + Ta; + const int status = regulatorLimitStatus(state); + const bool frozen = ((status > 0) && (regulatorDerivative > 0.0)) || + ((status < 0) && (regulatorDerivative < 0.0)); + + matrixData.assign(diffOffset + voltageMeasurementState, + diffOffset + voltageMeasurementState, + -1.0 / Tr - stateData.cj); + matrixData.assignCheckCol(diffOffset + voltageMeasurementState, + inputLocs[exciterVoltageInLocation], + 1.0 / Tr); + matrixData.assign(diffOffset + leadLagState, diffOffset + voltageMeasurementState, -1.0 / Tb); + matrixData.assign(diffOffset + leadLagState, + diffOffset + leadLagState, + -1.0 / Tb - stateData.cj); + matrixData.assign(diffOffset + leadLagState, + diffOffset + exciterState, + -feedbackGain * saturationSlope / Tb); + matrixData.assign(diffOffset + leadLagState, diffOffset + washoutState, feedbackGain / Tb); + matrixData.assignCheckCol(diffOffset + leadLagState, + inputLocs[exciterVsetInLocation], + 1.0 / Tb); + matrixData.assignCheckCol(diffOffset + leadLagState, inputLocs[exciterVssInLocation], 1.0 / Tb); + matrixData.assignCheckCol(diffOffset + leadLagState, + inputLocs[exciterXadIfdInLocation], + -feedbackGain * Kd / Tb); + + if (frozen) { + matrixData.assign(diffOffset + regulatorState, diffOffset + regulatorState, -stateData.cj); + } else { + matrixData.assign(diffOffset + regulatorState, + diffOffset + voltageMeasurementState, + -Ka * leadRatio / Ta); + matrixData.assign(diffOffset + regulatorState, + diffOffset + leadLagState, + Ka * (1.0 - leadRatio) / Ta); + matrixData.assign(diffOffset + regulatorState, + diffOffset + regulatorState, + -1.0 / Ta - stateData.cj); + matrixData.assign(diffOffset + regulatorState, + diffOffset + exciterState, + -Ka * leadRatio * feedbackGain * saturationSlope / Ta); + matrixData.assign(diffOffset + regulatorState, + diffOffset + washoutState, + Ka * leadRatio * feedbackGain / Ta); + matrixData.assignCheckCol(diffOffset + regulatorState, + inputLocs[exciterVsetInLocation], + Ka * leadRatio / Ta); + matrixData.assignCheckCol(diffOffset + regulatorState, + inputLocs[exciterVssInLocation], + Ka * leadRatio / Ta); + matrixData.assignCheckCol(diffOffset + regulatorState, + inputLocs[exciterXadIfdInLocation], + -Ka * leadRatio * feedbackGain * Kd / Ta); + } + + double regulatorGain = 0.0; + double targetExciterGain = 0.0; + double targetCurrentGain = 0.0; + regulatorTargetDerivatives(inputs, state, regulatorGain, targetExciterGain, targetCurrentGain); + matrixData.assign(diffOffset + exciterState, diffOffset + regulatorState, regulatorGain / Te); + matrixData.assign(diffOffset + exciterState, + diffOffset + exciterState, + (targetExciterGain - saturationSlope) / Te - stateData.cj); + matrixData.assignCheckCol(diffOffset + exciterState, + inputLocs[exciterXadIfdInLocation], + (targetCurrentGain - Kd) / Te); + matrixData.assign(diffOffset + washoutState, diffOffset + exciterState, saturationSlope / Tf); + matrixData.assign(diffOffset + washoutState, + diffOffset + washoutState, + -1.0 / Tf - stateData.cj); + matrixData.assignCheckCol(diffOffset + washoutState, + inputLocs[exciterXadIfdInLocation], + Kd / Tf); +} + +void ExciterEXAC1::timestep(CoreTime time, const IOdata& inputs, const SolverMode& /*sMode*/) +{ + derivative(inputs, emptyStateData, m_dstate_dt.data(), cLocalSolverMode); + const double timeStep = time - prevTime; + double* state = m_state.data() + 1; + const double* derivatives = m_dstate_dt.data() + 1; + for (index_t index = 0; index < 5; ++index) { + state[index] += timeStep * derivatives[index]; + } + m_state[0] = fieldVoltage(inputs, state); + updateLimitFlags(state); + prevTime = time; +} + +void ExciterEXAC1::rootTest(const IOdata& /*inputs*/, + const StateData& stateData, + double roots[], + const SolverMode& sMode) +{ + const auto locations = offsets.getLocations(stateData, sMode, this); + const index_t rootOffset = offsets.getRootOffset(sMode); + const double regulator = locations.diffStateLoc[regulatorState]; + roots[rootOffset] = + std::min(regulatorUpperLimit() - regulator, regulator - regulatorLowerLimit()); +} + +void ExciterEXAC1::rootTrigger(CoreTime /*time*/, + const IOdata& /*inputs*/, + const std::vector& rootMask, + const SolverMode& sMode) +{ + const index_t rootOffset = offsets.getRootOffset(sMode); + if ((rootMask[rootOffset] != 0) && updateLimitFlags(m_state.data() + 1)) { + alert(this, JAC_COUNT_CHANGE); + } +} + +ChangeCode ExciterEXAC1::rootCheck(const IOdata& /*inputs*/, + const StateData& /*stateData*/, + const SolverMode& /*sMode*/, + CheckLevel /*level*/) +{ + if (updateLimitFlags(m_state.data() + 1)) { + alert(this, JAC_COUNT_CHANGE); + return ChangeCode::JACOBIAN_CHANGE; + } + return ChangeCode::NO_CHANGE; +} + +void ExciterEXAC1::set(std::string_view param, std::string_view val) +{ + Exciter::set(param, val); +} + +void ExciterEXAC1::set(std::string_view param, double val, units::unit unitType) +{ + const auto positive = [val](const char* label) { + if (!std::isfinite(val) || (val <= 0.0)) { + throw InvalidParameterValue(std::string("EXAC1 ") + label + + " must be positive and finite"); + } + }; + const auto finite = [val](const char* label) { + if (!std::isfinite(val)) { + throw InvalidParameterValue(std::string("EXAC1 ") + label + " must be finite"); + } + }; + if (param == "tr") { + positive("TR"); + Tr = val; + } else if (param == "tb") { + positive("TB"); + Tb = val; + } else if (param == "tc") { + finite("TC"); + Tc = val; + } else if (param == "te") { + positive("TE"); + Te = val; + } else if (param == "tf") { + positive("TF"); + Tf = val; + } else if (param == "kf") { + finite("KF"); + Kf = val; + } else if (param == "kc") { + finite("KC"); + Kc = val; + } else if (param == "kd") { + finite("KD"); + Kd = val; + } else if (param == "ke") { + finite("KE"); + Ke = val; + } else if (param == "e1") { + finite("E1"); + E1 = val; + } else if (param == "se1") { + finite("SE1"); + Se1 = val; + } else if (param == "e2") { + finite("E2"); + E2 = val; + } else if (param == "se2") { + finite("SE2"); + Se2 = val; + } else { + Exciter::set(param, val, unitType); + } +} + +double ExciterEXAC1::get(std::string_view param, units::unit unitType) const +{ + if (param == "ka") { + return Ka; + } + if (param == "ta") { + return Ta; + } + if ((param == "vrmax") || (param == "urmax")) { + return Vrmax; + } + if ((param == "vrmin") || (param == "urmin")) { + return Vrmin; + } + if (param == "tr") { + return Tr; + } + if (param == "tb") { + return Tb; + } + if (param == "tc") { + return Tc; + } + if (param == "te") { + return Te; + } + if (param == "tf") { + return Tf; + } + if (param == "kf") { + return Kf; + } + if (param == "kc") { + return Kc; + } + if (param == "kd") { + return Kd; + } + if (param == "ke") { + return Ke; + } + if (param == "e1") { + return E1; + } + if (param == "se1") { + return Se1; + } + if (param == "e2") { + return E2; + } + if (param == "se2") { + return Se2; + } + return Exciter::get(param, unitType); +} + +stringVec ExciterEXAC1::localStateNames() const +{ + return {"efd", "vmeas", "ll", "va", "ve", "wf"}; +} + +index_t ExciterEXAC1::findIndex(std::string_view field, const SolverMode& sMode) const +{ + if ((field == "efd") || (field == "field")) { + return offsets.getAlgOffset(sMode); + } + const index_t offset = offsets.getDiffOffset(sMode); + if (field == "vmeas") { + return offset + voltageMeasurementState; + } + if ((field == "ll") || (field == "leadlag")) { + return offset + leadLagState; + } + if ((field == "va") || (field == "regulator")) { + return offset + regulatorState; + } + if ((field == "ve") || (field == "exciter")) { + return offset + exciterState; + } + if ((field == "wf") || (field == "washout")) { + return offset + washoutState; + } + return kInvalidLocation; +} +// NOLINTEND(readability-math-missing-parentheses) +} // namespace griddyn::exciters diff --git a/src/griddyn/exciters/ExciterEXAC1.h b/src/griddyn/exciters/ExciterEXAC1.h new file mode 100644 index 00000000..d8471483 --- /dev/null +++ b/src/griddyn/exciters/ExciterEXAC1.h @@ -0,0 +1,121 @@ +/* + * 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 "../Exciter.h" +#include "utilities/Saturation.h" +#include +#include + +namespace griddyn::exciters { +/** PSS/E EXAC1 alternating-current excitation system. + * + * This class implements the EXAC1 record in the PSS/E Model Library and the + * frozen ANDES 2.0.0 parameter contract. With the measured voltage + * @f$v_m@f$, lead-lag state @f$x_L@f$, limited regulator state @f$v_A@f@, + * exciter voltage @f$v_E@f@, and washout state @f$x_F@f@, the equations are + * @f[ + * T_R\dot v_m=V_t-v_m,\quad T_B\dot x_L=v_i-x_L,\quad + * y_L=x_L+\frac{T_C}{T_B}(v_i-x_L), + * @f] + * @f[ + * T_A\dot v_A=K_Ay_L-v_A,\quad + * T_E\dot v_E=v_R-v_{FE},\quad T_F\dot x_F=v_{FE}-x_F, + * @f] + * where @f$v_i=V_{ref}-v_m-K_F(v_{FE}-x_F)/T_F@f$, + * @f$v_{FE}=K_Ev_E+S_E(v_E)+K_DX_{ad}I_{fd}@f$, and + * @f$E_{fd}=v_EF_{EX}(K_CX_{ad}I_{fd}/v_E)@f$. The regulator state has an + * anti-windup limiter @f$[V_{RMIN},V_{RMAX}]@f$. + * + * Frozen ANDES constructs the @f$T_R@f$ transducer but subtracts raw terminal + * voltage from the regulator input. GridDyn intentionally uses @f$v_m@f$ as + * required by the published AC1 block diagram; see + * docs/developer-guide/andes-compatibility.md for the compatibility note and + * an issue draft. References: IEEE Std 421.5-2016 and the PSS/E Model + * Library EXAC1 description. + */ +class ExciterEXAC1: public Exciter { + public: + enum EXACFlags { + REGULATOR_LIMITED = OBJECT_FLAG5, + REGULATOR_LIMIT_HIGH = OBJECT_FLAG6, + }; + + protected: + model_parameter Tr = 0.01; + model_parameter Tb = 1.0; + model_parameter Tc = 1.0; + model_parameter Te = 0.8; + model_parameter Kf = 0.1; + model_parameter Tf = 1.0; + model_parameter Kc = 0.1; + model_parameter Kd = 0.0; + model_parameter Ke = 1.0; + model_parameter E1 = 0.0; + model_parameter Se1 = 0.0; + model_parameter E2 = 1.0; + model_parameter Se2 = 0.0; + utilities::Saturation saturation{utilities::Saturation::SaturationType::QUADRATIC}; + + public: + explicit ExciterEXAC1(const std::string& objName = "exciterEXAC1_#"); + 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; + 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; + stringVec localStateNames() const override; + index_t findIndex(std::string_view field, const SolverMode& sMode) 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 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; + + protected: + virtual double regulatorTarget(const IOdata& inputs, const double state[]) const; + virtual double regulatorUpperLimit() const; + virtual double regulatorLowerLimit() const; + virtual double initialRegulatorState(double vfe) const; + virtual double referenceOffset(double vfe) const; + virtual void regulatorTargetDerivatives(const IOdata& inputs, + const double state[], + double& regulatorDerivative, + double& exciterDerivative, + double& fieldCurrentDerivative) const; + double referenceInput(const IOdata& inputs) const; + double vfe(const IOdata& inputs, const double state[]) const; + double rectifierFactor(const IOdata& inputs, double exciterVoltage) const; + double fieldVoltage(const IOdata& inputs, const double state[]) const; + int regulatorLimitStatus(const double state[]) const; + bool updateLimitFlags(const double state[]); +}; +} // namespace griddyn::exciters diff --git a/src/griddyn/exciters/ExciterEXAC2.cpp b/src/griddyn/exciters/ExciterEXAC2.cpp new file mode 100644 index 00000000..d6c4c1c7 --- /dev/null +++ b/src/griddyn/exciters/ExciterEXAC2.cpp @@ -0,0 +1,208 @@ +/* + * 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 "ExciterEXAC2.h" + +#include "core/CoreExceptions.h" +#include "core/CoreObjectTemplates.hpp" +#include +#include +#include + +namespace griddyn::exciters { +namespace { + // These decimal values are the specified PSS/E FEX curve coefficients, + // not approximations of unrelated mathematical constants. + constexpr double lowCurrentSlope = 0.577; // NOLINT(modernize-use-std-numbers) + constexpr double highCurrentSlope = 1.732; // NOLINT(modernize-use-std-numbers) +} // namespace + +ExciterEXAC2::ExciterEXAC2(const std::string& objName): ExciterEXAC1(objName) {} + +CoreObject* ExciterEXAC2::clone(CoreObject* obj) const +{ + auto* clone = cloneBase(this, obj); + if (clone == nullptr) { + return obj; + } + ExciterEXAC1::clone(clone); + clone->Vamax = Vamax; + clone->Vamin = Vamin; + clone->Vlr = Vlr; + clone->Kl = Kl; + clone->Kh = Kh; + clone->Kb = Kb; + clone->vlr0 = vlr0; + return clone; +} + +void ExciterEXAC2::dynObjectInitializeA(CoreTime time0, std::uint32_t flags) +{ + if (!std::isfinite(Vamax) || !std::isfinite(Vamin) || !std::isfinite(Vlr) || + !std::isfinite(Kl) || !std::isfinite(Kh) || !std::isfinite(Kb) || (Kb <= 0.0) || + (Kl <= 0.0) || (Vamax < Vamin)) { + throw InvalidParameterValue("EXAC2 gains or limits"); + } + ExciterEXAC1::dynObjectInitializeA(time0, flags); +} + +void ExciterEXAC2::dynObjectInitializeB(const IOdata& inputs, + const IOdata& desiredOutput, + IOdata& fieldSet) +{ + // VLR0 is a PostInitService in ANDES: it is calculated once from the + // initialized field-feedback voltage and remains constant thereafter. + const double fieldVoltage = desiredOutput.empty() ? 0.0 : desiredOutput[0]; + double exciterVoltage = std::max(1e-8, std::abs(fieldVoltage)); + const double fieldCurrent = inputs[exciterXadIfdInLocation]; + for (int count = 0; count < 20; ++count) { + const double normalized = (Kc != 0.0) ? Kc * fieldCurrent / exciterVoltage : 0.0; + double factor = 1.0; + double slope = 0.0; + if ((normalized > 0.0) && (normalized <= 0.433)) { + factor = 1.0 - (lowCurrentSlope * normalized); + slope = -lowCurrentSlope; + } else if ((normalized > 0.433) && (normalized <= 0.75)) { + factor = std::sqrt(std::max(0.0, 0.75 - (normalized * normalized))); + slope = (factor > 0.0) ? -normalized / factor : 0.0; + } else if ((normalized > 0.75) && (normalized <= 1.0)) { + factor = highCurrentSlope * (1.0 - normalized); + slope = -highCurrentSlope; + } else if (normalized > 1.0) { + factor = 0.0; + } + const double error = (exciterVoltage * factor) - fieldVoltage; + const double derivative = factor - (slope * normalized); + if (std::abs(derivative) < 1e-12) { + break; + } + exciterVoltage = std::max(1e-8, exciterVoltage - (error / derivative)); + if (std::abs(error) < 1e-12) { + break; + } + } + const double feedback = (Ke * exciterVoltage) + saturation(exciterVoltage) + + ((Kd == 0.0) ? 0.0 : Kd * fieldCurrent); + vlr0 = std::max(Vlr, feedback + (feedback / (Kl * Kb))); + ExciterEXAC1::dynObjectInitializeB(inputs, desiredOutput, fieldSet); +} + +void ExciterEXAC2::set(std::string_view param, std::string_view val) +{ + ExciterEXAC1::set(param, val); +} + +void ExciterEXAC2::set(std::string_view param, double val, units::unit unitType) +{ + const auto finite = [val](const char* label) { + if (!std::isfinite(val)) { + throw InvalidParameterValue(std::string("EXAC2 ") + label + " must be finite"); + } + }; + if (param == "vamax") { + finite("VAMAX"); + Vamax = val; + } else if (param == "vamin") { + finite("VAMIN"); + Vamin = val; + } else if (param == "vlr") { + finite("VLR"); + Vlr = val; + } else if (param == "kl") { + if (!std::isfinite(val) || (val <= 0.0)) { + throw InvalidParameterValue("EXAC2 KL must be positive and finite"); + } + Kl = val; + } else if (param == "kh") { + finite("KH"); + Kh = val; + } else if (param == "kb") { + if (!std::isfinite(val) || (val <= 0.0)) { + throw InvalidParameterValue("EXAC2 KB must be positive and finite"); + } + Kb = val; + } else { + ExciterEXAC1::set(param, val, unitType); + } +} + +double ExciterEXAC2::get(std::string_view param, units::unit unitType) const +{ + if (param == "vamax") { + return Vamax; + } + if (param == "vamin") { + return Vamin; + } + if (param == "vlr") { + return Vlr; + } + if (param == "kl") { + return Kl; + } + if (param == "kh") { + return Kh; + } + if (param == "kb") { + return Kb; + } + return ExciterEXAC1::get(param, unitType); +} + +double ExciterEXAC2::regulatorTarget(const IOdata& inputs, const double state[]) const +{ + const double feedback = vfe(inputs, state); + const double highGate = state[2] - (Kh * feedback); + const double lowGate = Kl * (vlr0 - feedback); + return std::clamp(Kb * std::min(highGate, lowGate), + static_cast(Vrmin), + static_cast(Vrmax)); +} + +double ExciterEXAC2::regulatorUpperLimit() const +{ + return Vamax; +} +double ExciterEXAC2::regulatorLowerLimit() const +{ + return Vamin; +} +double ExciterEXAC2::initialRegulatorState(double vfeValue) const +{ + return (vfeValue * Kl) + (vfeValue / Kb); +} +double ExciterEXAC2::referenceOffset(double vfeValue) const +{ + return ((vfeValue * Kl) + (vfeValue / Kb)) / Ka; +} + +void ExciterEXAC2::regulatorTargetDerivatives(const IOdata& inputs, + const double state[], + double& regulatorDerivative, + double& exciterDerivative, + double& fieldCurrentDerivative) const +{ + const double feedbackSlope = Ke + saturation.deriv(state[3]); + const double feedback = vfe(inputs, state); + const double highGate = state[2] - (Kh * feedback); + const double lowGate = Kl * (vlr0 - feedback); + const double unbounded = Kb * std::min(highGate, lowGate); + regulatorDerivative = 0.0; + exciterDerivative = 0.0; + fieldCurrentDerivative = 0.0; + if ((unbounded <= Vrmin) || (unbounded >= Vrmax)) { + return; + } + if (highGate <= lowGate) { + regulatorDerivative = Kb; + exciterDerivative = -Kb * Kh * feedbackSlope; + fieldCurrentDerivative = -Kb * Kh * Kd; + } else { + exciterDerivative = -Kb * Kl * feedbackSlope; + fieldCurrentDerivative = -Kb * Kl * Kd; + } +} +} // namespace griddyn::exciters diff --git a/src/griddyn/exciters/ExciterEXAC2.h b/src/griddyn/exciters/ExciterEXAC2.h new file mode 100644 index 00000000..27d8101d --- /dev/null +++ b/src/griddyn/exciters/ExciterEXAC2.h @@ -0,0 +1,59 @@ +/* + * 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 "ExciterEXAC1.h" +#include + +namespace griddyn::exciters { +/** PSS/E EXAC2 excitation system. + * + * EXAC2 shares the EXAC1 transducer, lead-lag, field integrator, washout, + * saturation, and rectifier equations. It replaces the regulator path by + * @f[ + * v_H=v_A-K_Hv_{FE},\quad v_L=K_L(V_{LR0}-v_{FE}),\quad + * v_R=\operatorname{limit}\{K_B\min(v_H,v_L)\}. + * @f] + * The @f$K_A@f$ lag is anti-windup limited by @f$[V_{AMIN},V_{AMAX}]@f$ and + * @f$V_{LR0}=\max(V_{LR},v_{FE0}+v_{FE0}/(K_LK_B))@f$ is retained after + * initialization. The same intentional sensed-voltage correction as EXAC1 + * applies; see ExciterEXAC1 and the ANDES compatibility guide. + */ +class ExciterEXAC2 final: public ExciterEXAC1 { + protected: + model_parameter Vamax = 8.0; + model_parameter Vamin = 0.0; + model_parameter Vlr = 0.0; + model_parameter Kl = 1.0; + model_parameter Kh = 1.0; + model_parameter Kb = 1.0; + double vlr0 = 0.0; + + public: + explicit ExciterEXAC2(const std::string& objName = "exciterEXAC2_#"); + 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; + 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; + + protected: + double regulatorTarget(const IOdata& inputs, const double state[]) const override; + double regulatorUpperLimit() const override; + double regulatorLowerLimit() const override; + double initialRegulatorState(double vfe) const override; + double referenceOffset(double vfe) const override; + void regulatorTargetDerivatives(const IOdata& inputs, + const double state[], + double& regulatorDerivative, + double& exciterDerivative, + double& fieldCurrentDerivative) const override; +}; +} // namespace griddyn::exciters diff --git a/src/griddyn/exciters/ExciterEXAC4.cpp b/src/griddyn/exciters/ExciterEXAC4.cpp new file mode 100644 index 00000000..c3f46eae --- /dev/null +++ b/src/griddyn/exciters/ExciterEXAC4.cpp @@ -0,0 +1,365 @@ +/* + * 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 "ExciterEXAC4.h" + +#include "core/CoreExceptions.h" +#include "core/CoreObjectTemplates.hpp" +#include "utilities/MatrixData.hpp" +#include +#include +#include +#include + +namespace griddyn::exciters { +// NOLINTBEGIN(readability-math-missing-parentheses) +namespace { + constexpr index_t voltageMeasurementState = 0; + constexpr index_t leadLagState = 1; + constexpr index_t regulatorState = 2; +} // namespace + +ExciterEXAC4::ExciterEXAC4(const std::string& objName): Exciter(objName) +{ + m_inputSize = exciterInputCount; + Ka = 80.0; + Ta = 0.04; + Vrmax = 8.0; + Vrmin = 0.0; +} + +CoreObject* ExciterEXAC4::clone(CoreObject* obj) const +{ + auto* clone = cloneBase(this, obj); + if (clone == nullptr) { + return obj; + } + clone->Tr = Tr; + clone->Vimax = Vimax; + clone->Vimin = Vimin; + clone->Tc = Tc; + clone->Tb = Tb; + clone->Kc = Kc; + clone->vref0 = vref0; + return clone; +} + +void ExciterEXAC4::dynObjectInitializeA(CoreTime /*time0*/, std::uint32_t /*flags*/) +{ + if (!std::isfinite(Tr) || !std::isfinite(Vimax) || !std::isfinite(Vimin) || + !std::isfinite(Tc) || !std::isfinite(Tb) || !std::isfinite(Ka) || !std::isfinite(Ta) || + !std::isfinite(Vrmax) || !std::isfinite(Vrmin) || !std::isfinite(Kc) || (Tr <= 0.0) || + (Tb <= 0.0) || (Ta <= 0.0) || (Ka <= 0.0) || (Vimax < Vimin) || (Vrmax < Vrmin)) { + throw InvalidParameterValue("EXAC4 gains, time constants, or limits"); + } + offsets.local().local.algSize = 1; + offsets.local().local.diffSize = 3; + offsets.local().local.algRoots = 2; + offsets.local().local.jacSize = 19; +} + +void ExciterEXAC4::dynObjectInitializeB(const IOdata& inputs, + const IOdata& desiredOutput, + IOdata& fieldSet) +{ + if (!std::isfinite(inputs[exciterVoltageInLocation]) || + ((Kc != 0.0) && !std::isfinite(inputs[exciterXadIfdInLocation])) || desiredOutput.empty() || + !std::isfinite(desiredOutput[0])) { + throw InvalidParameterValue("EXAC4 initial generator signals or field voltage"); + } + const double fieldVoltage = desiredOutput[0]; + const double fieldCurrent = inputs[exciterXadIfdInLocation]; + const double regulator = fieldVoltage; + const double lower = Vrmin - Kc * fieldCurrent; + const double upper = Vrmax - Kc * fieldCurrent; + if ((regulator < lower - 1e-7) || (regulator > upper + 1e-7)) { + throw InvalidParameterValue("EXAC4 initial field voltage outside limits"); + } + double* state = m_state.data() + 1; + state[voltageMeasurementState] = inputs[exciterVoltageInLocation]; + state[regulatorState] = regulator; + state[leadLagState] = regulator / Ka; + m_state[0] = fieldVoltage; + vref0 = inputs[exciterVoltageInLocation] + fieldVoltage / Ka; + fieldSet[exciterVsetInLocation] = Vref; + std::fill(m_dstate_dt.begin(), m_dstate_dt.end(), 0.0); + updateLimitFlags(inputs, state); +} + +double ExciterEXAC4::unlimitedInput(const double state[]) const +{ + return vref0 - state[voltageMeasurementState]; +} +double ExciterEXAC4::limitedInput(const double state[]) const +{ + return std::clamp(unlimitedInput(state), + static_cast(Vimin), + static_cast(Vimax)); +} +int ExciterEXAC4::inputLimitStatus(const double state[]) const +{ + const double value = unlimitedInput(state); + if (value >= Vimax) { + return 1; + } + return (value <= Vimin) ? -1 : 0; +} +int ExciterEXAC4::outputLimitStatus(const IOdata& inputs, const double state[]) const +{ + const double value = state[regulatorState]; + const double upper = Vrmax - Kc * inputs[exciterXadIfdInLocation]; + const double lower = Vrmin - Kc * inputs[exciterXadIfdInLocation]; + if (value >= upper) { + return 1; + } + return (value <= lower) ? -1 : 0; +} +double ExciterEXAC4::fieldVoltage(const IOdata& inputs, const double state[]) const +{ + const int status = outputLimitStatus(inputs, state); + if (status > 0) { + return Vrmax - Kc * inputs[exciterXadIfdInLocation]; + } + if (status < 0) { + return Vrmin - Kc * inputs[exciterXadIfdInLocation]; + } + return state[regulatorState]; +} +bool ExciterEXAC4::updateLimitFlags(const IOdata& inputs, const double state[]) +{ + const int input = inputLimitStatus(state); + const int output = outputLimitStatus(inputs, state); + const bool changed = (opFlags[INPUT_LIMITED] != (input != 0)) || + (opFlags[INPUT_LIMIT_HIGH] != (input > 0)) || (opFlags[OUTPUT_LIMITED] != (output != 0)) || + (opFlags[OUTPUT_LIMIT_HIGH] != (output > 0)); + opFlags.set(INPUT_LIMITED, input != 0); + opFlags.set(INPUT_LIMIT_HIGH, input > 0); + opFlags.set(OUTPUT_LIMITED, output != 0); + opFlags.set(OUTPUT_LIMIT_HIGH, output > 0); + return changed; +} + +void ExciterEXAC4::residual(const IOdata& inputs, + const StateData& stateData, + double resid[], + const SolverMode& sMode) +{ + const auto locations = offsets.getLocations(stateData, resid, sMode, this); + if (hasAlgebraic(sMode)) { + locations.destLoc[0] = + fieldVoltage(inputs, locations.diffStateLoc) - locations.algStateLoc[0]; + } + if (hasDifferential(sMode)) { + derivative(inputs, stateData, resid, sMode); + for (index_t ii = 0; ii < locations.diffSize; ++ii) { + locations.destDiffLoc[ii] -= locations.dstateLoc[ii]; + } + } +} + +void ExciterEXAC4::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* derivativeValues = locations.destDiffLoc; + const double input = limitedInput(state); + const double lead = state[leadLagState] + Tc * (input - state[leadLagState]) / Tb; + derivativeValues[voltageMeasurementState] = + (inputs[exciterVoltageInLocation] - state[voltageMeasurementState]) / Tr; + derivativeValues[leadLagState] = (input - state[leadLagState]) / Tb; + derivativeValues[regulatorState] = (Ka * lead - state[regulatorState]) / Ta; +} + +void ExciterEXAC4::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 alg = locations.algOffset; + const index_t diff = locations.diffOffset; + const double* state = locations.diffStateLoc; + if (hasAlgebraic(sMode)) { + matrixData.assign(alg, alg, -1.0); + if (outputLimitStatus(inputs, state) == 0) { + matrixData.assign(alg, diff + regulatorState, 1.0); + } else { + matrixData.assignCheckCol(alg, inputLocs[exciterXadIfdInLocation], -Kc); + } + } + if (!hasDifferential(sMode)) { + return; + } + matrixData.assign(diff + voltageMeasurementState, + diff + voltageMeasurementState, + -1.0 / Tr - stateData.cj); + matrixData.assignCheckCol(diff + voltageMeasurementState, + inputLocs[exciterVoltageInLocation], + 1.0 / Tr); + const bool limited = inputLimitStatus(state) != 0; + const double ratio = Tc / Tb; + matrixData.assign(diff + leadLagState, diff + leadLagState, -1.0 / Tb - stateData.cj); + if (!limited) { + matrixData.assign(diff + leadLagState, diff + voltageMeasurementState, -1.0 / Tb); + } + matrixData.assign(diff + regulatorState, diff + leadLagState, Ka * (1.0 - ratio) / Ta); + matrixData.assign(diff + regulatorState, diff + regulatorState, -1.0 / Ta - stateData.cj); + if (!limited) { + matrixData.assign(diff + regulatorState, diff + voltageMeasurementState, -Ka * ratio / Ta); + } +} + +void ExciterEXAC4::timestep(CoreTime time, const IOdata& inputs, const SolverMode& /*sMode*/) +{ + derivative(inputs, emptyStateData, m_dstate_dt.data(), cLocalSolverMode); + const double timeStep = time - prevTime; + double* state = m_state.data() + 1; + const double* derivativeValues = m_dstate_dt.data() + 1; + for (index_t ii = 0; ii < 3; ++ii) { + state[ii] += timeStep * derivativeValues[ii]; + } + m_state[0] = fieldVoltage(inputs, state); + updateLimitFlags(inputs, state); + prevTime = time; +} + +void ExciterEXAC4::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 root = offsets.getRootOffset(sMode); + const double input = unlimitedInput(state); + roots[root] = std::min(Vimax - input, input - Vimin); + const double output = state[regulatorState]; + roots[root + 1] = std::min(Vrmax - Kc * inputs[exciterXadIfdInLocation] - output, + output - (Vrmin - Kc * inputs[exciterXadIfdInLocation])); +} +void ExciterEXAC4::rootTrigger(CoreTime /*time*/, + const IOdata& inputs, + const std::vector& rootMask, + const SolverMode& sMode) +{ + const index_t root = offsets.getRootOffset(sMode); + if (((rootMask[root] != 0) || (rootMask[root + 1] != 0)) && + updateLimitFlags(inputs, m_state.data() + 1)) { + alert(this, JAC_COUNT_CHANGE); + } +} +ChangeCode ExciterEXAC4::rootCheck(const IOdata& inputs, + const StateData& /*stateData*/, + const SolverMode& /*sMode*/, + CheckLevel /*level*/) +{ + if (updateLimitFlags(inputs, m_state.data() + 1)) { + alert(this, JAC_COUNT_CHANGE); + return ChangeCode::JACOBIAN_CHANGE; + } + return ChangeCode::NO_CHANGE; +} +void ExciterEXAC4::set(std::string_view param, std::string_view val) +{ + Exciter::set(param, val); +} +void ExciterEXAC4::set(std::string_view param, double val, units::unit unitType) +{ + const auto finite = [val](const char* label) { + if (!std::isfinite(val)) { + throw InvalidParameterValue(std::string("EXAC4 ") + label + " must be finite"); + } + }; + if (param == "tr") { + if (!std::isfinite(val) || (val <= 0.0)) { + throw InvalidParameterValue("EXAC4 TR must be positive and finite"); + } + Tr = val; + } else if (param == "tb") { + if (!std::isfinite(val) || (val <= 0.0)) { + throw InvalidParameterValue("EXAC4 TB must be positive and finite"); + } + Tb = val; + } else if (param == "tc") { + finite("TC"); + Tc = val; + } else if (param == "vimax") { + finite("VIMAX"); + Vimax = val; + } else if (param == "vimin") { + finite("VIMIN"); + Vimin = val; + } else if (param == "kc") { + finite("KC"); + Kc = val; + } else { + Exciter::set(param, val, unitType); + } +} +double ExciterEXAC4::get(std::string_view param, units::unit unitType) const +{ + if (param == "ka") { + return Ka; + } + if (param == "ta") { + return Ta; + } + if ((param == "vrmax") || (param == "urmax")) { + return Vrmax; + } + if ((param == "vrmin") || (param == "urmin")) { + return Vrmin; + } + if (param == "tr") { + return Tr; + } + if (param == "tb") { + return Tb; + } + if (param == "tc") { + return Tc; + } + if (param == "vimax") { + return Vimax; + } + if (param == "vimin") { + return Vimin; + } + if (param == "kc") { + return Kc; + } + return Exciter::get(param, unitType); +} +stringVec ExciterEXAC4::localStateNames() const +{ + return {"efd", "vmeas", "ll", "vr"}; +} +index_t ExciterEXAC4::findIndex(std::string_view field, const SolverMode& sMode) const +{ + if ((field == "efd") || (field == "field")) { + return offsets.getAlgOffset(sMode); + } + const index_t offset = offsets.getDiffOffset(sMode); + if (field == "vmeas") { + return offset + voltageMeasurementState; + } + if ((field == "ll") || (field == "leadlag")) { + return offset + leadLagState; + } + if ((field == "vr") || (field == "regulator")) { + return offset + regulatorState; + } + return kInvalidLocation; +} +// NOLINTEND(readability-math-missing-parentheses) +} // namespace griddyn::exciters diff --git a/src/griddyn/exciters/ExciterEXAC4.h b/src/griddyn/exciters/ExciterEXAC4.h new file mode 100644 index 00000000..a1049293 --- /dev/null +++ b/src/griddyn/exciters/ExciterEXAC4.h @@ -0,0 +1,99 @@ +/* + * 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 "../Exciter.h" +#include +#include + +namespace griddyn::exciters { +/** PSS/E EXAC4 fast alternating-current excitation system. + * + * The EXAC4 record uses a sensed terminal-voltage lag, a limited input and a + * lead-lag compensator followed by a regulator: + * @f[ + * T_R\dot v_m=V_t-v_m,\quad v_i=V_{ref0}-v_m,\quad + * u=\operatorname{limit}_{[V_{IMIN},V_{IMAX}]}(v_i), + * @f] + * @f[ + * T_B\dot x_L=u-x_L,\quad y_L=x_L+\frac{T_C}{T_B}(u-x_L),\quad + * T_A\dot v_R=K_Ay_L-v_R, + * @f] + * @f[ + * E_{fd}=\operatorname{limit}_{[V_{RMIN}-K_CX_{ad}I_{fd},\ + * V_{RMAX}-K_CX_{ad}I_{fd}]}(v_R). + * @f] + * @f$V_{ref0}=V_{t0}+E_{fd0}/K_A@f$ is retained after initialization, as in + * frozen ANDES 2.0.0. References: IEEE Std 421.5-2016 and the PSS/E Model + * Library EXAC4 description. + */ +class ExciterEXAC4 final: public Exciter { + public: + enum EXAC4Flags { + INPUT_LIMITED = OBJECT_FLAG5, + INPUT_LIMIT_HIGH = OBJECT_FLAG6, + OUTPUT_LIMITED = OBJECT_FLAG7, + OUTPUT_LIMIT_HIGH = OBJECT_FLAG8 + }; + + private: + model_parameter Tr = 0.01; + model_parameter Vimax = 5.0; + model_parameter Vimin = -0.1; + model_parameter Tc = 1.0; + model_parameter Tb = 1.0; + model_parameter Kc = 0.0; + double vref0 = 1.0; + + public: + explicit ExciterEXAC4(const std::string& objName = "exciterEXAC4_#"); + 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; + 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; + stringVec localStateNames() const override; + index_t findIndex(std::string_view field, const SolverMode& sMode) 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 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; + + private: + double unlimitedInput(const double state[]) const; + double limitedInput(const double state[]) const; + int inputLimitStatus(const double state[]) const; + int outputLimitStatus(const IOdata& inputs, const double state[]) const; + double fieldVoltage(const IOdata& inputs, const double state[]) const; + bool updateLimitFlags(const IOdata& inputs, const double state[]); +}; +} // namespace griddyn::exciters diff --git a/test/andesTests/testAndesDyrReader.cpp b/test/andesTests/testAndesDyrReader.cpp index 2b9a294a..3daf7422 100644 --- a/test/andesTests/testAndesDyrReader.cpp +++ b/test/andesTests/testAndesDyrReader.cpp @@ -12,6 +12,9 @@ #include "griddyn/GridSubModel.h" #include "griddyn/events/Event.h" #include "griddyn/exciters/ExciterESST3A.h" +#include "griddyn/exciters/ExciterEXAC1.h" +#include "griddyn/exciters/ExciterEXAC2.h" +#include "griddyn/exciters/ExciterEXAC4.h" #include "griddyn/exciters/ExciterEXST1.h" #include "griddyn/generators/DynamicGenerator.h" #include "griddyn/genmodels/GenModelGENROU.h" @@ -364,6 +367,35 @@ TEST(AndesDyrReaderTests, MapsExst1ParametersAndCouplesToGenrou) EXPECT_EQ(runJacobianCheck(simulation, griddyn::cDaeSolverMode, false), 0); } +TEST(AndesDyrReaderTests, MapsExacParameterRecordsAndCouplesToGenrou) +{ + const std::array, 3> records{{ + {"ieee14_exac1.dyr", "exac1"}, + {"ieee14_exac2.dyr", "exac2"}, + {"ieee14_exac4.dyr", "exac4"}, + }}; + for (const auto& [record, model] : records) { + 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(record)); + auto* bus = dynamic_cast(simulation->findByUserID("bus", 2)); + ASSERT_NE(bus, nullptr); + auto* generator = bus->getGen(0); + ASSERT_NE(generator, nullptr); + auto* exciter = dynamic_cast(generator->find("exciter")); + ASSERT_NE(exciter, nullptr); + EXPECT_DOUBLE_EQ(exciter->get("tr"), 0.031) << model; + EXPECT_DOUBLE_EQ(exciter->get("ka"), 41.0) << model; + EXPECT_DOUBLE_EQ(exciter->get("ta"), 0.034) << model; + EXPECT_DOUBLE_EQ(exciter->get("vrmax"), 7.2) << model; + EXPECT_DOUBLE_EQ(exciter->get("vrmin"), -4.3) << model; + ASSERT_EQ(simulation->dynInitialize(), 0) << model; + EXPECT_EQ(runResidualCheck(simulation, griddyn::cDaeSolverMode, false), 0) << model; + EXPECT_EQ(runJacobianCheck(simulation, griddyn::cDaeSolverMode, false), 0) << model; + } +} + TEST(AndesDyrReaderTests, MapsSt2cutParametersAndCouplesToExciters) { auto simulation = std::make_unique(); @@ -461,6 +493,14 @@ TEST(AndesDynamicTests, Exst1RespondsToLoadStep) EXPECT_FALSE(finalState.empty()); } +TEST(AndesDynamicTests, ExacExcitersRespondToLoadStep) +{ + for (const auto record : {"ieee14_exac1.dyr", "ieee14_exac2.dyr", "ieee14_exac4.dyr"}) { + const auto finalState = runLoadStepCase({record}); + EXPECT_FALSE(finalState.empty()) << record; + } +} + TEST(AndesDynamicTests, St2cutRespondsToLoadStep) { const auto finalState = diff --git a/test/componentTests/testExciters.cpp b/test/componentTests/testExciters.cpp index f05b7087..27ac5643 100644 --- a/test/componentTests/testExciters.cpp +++ b/test/componentTests/testExciters.cpp @@ -11,6 +11,8 @@ #include "gmlc/utilities/vectorOps.hpp" #include "griddyn/Generator.h" #include "griddyn/exciters/ExciterESST3A.h" +#include "griddyn/exciters/ExciterEXAC1.h" +#include "griddyn/exciters/ExciterEXAC4.h" #include "griddyn/exciters/ExciterEXST1.h" #include "solvers/SolverMode.hpp" #include @@ -86,7 +88,8 @@ void verifyStabilityCase(ExciterTests& fixture, for (const auto& exciterName : exciterList) { if (exciterName.starts_with("fmi") || (exciterName == "esst3a") || - (exciterName == "exst1")) { + (exciterName == "exst1") || (exciterName == "exac1") || (exciterName == "exac2") || + (exciterName == "exac4")) { continue; } if (std::find(skippedExcters.begin(), skippedExcters.end(), exciterName) != @@ -264,6 +267,103 @@ TEST(ExciterModelTests, Exst1MatchesAndesInitializationAndPerturbedEquations) EXPECT_NEAR(residual[0], 0.1, 1e-12); } +TEST(ExciterModelTests, Exac1InitializesCorrectedTransducerAndAntiWindup) +{ + exciters::ExciterEXAC1 exciter; + exciter.set("tr", 0.1); + exciter.set("tb", 0.2); + exciter.set("tc", 0.05); + exciter.set("ka", 10.0); + exciter.set("ta", 0.1); + exciter.set("te", 0.5); + exciter.set("kf", 0.1); + exciter.set("tf", 1.0); + exciter.set("kc", 0.0); + exciter.set("kd", 0.0); + exciter.set("ke", 1.0); + exciter.set("vrmax", 0.5); + exciter.set("vrmin", -0.5); + exciter.dynInitializeA(0.0, 0); + exciter.setRootOffset(0, cLocalSolverMode); + + IOdata inputs(exciterInputCount, 0.0); + inputs[exciterVoltageInLocation] = 1.0; + inputs[exciterVsetInLocation] = 1.0; + inputs[exciterXadIfdInLocation] = 0.0; + inputs[exciterVssInLocation] = 0.0; + IOdata fieldSet(4, 0.0); + // Use an in-limit operating point for initialization, then perturb all + // blocks independently. The sensed-voltage state is intentionally used. + exciter.dynInitializeB(inputs, {0.4}, fieldSet); + const auto& initialized = exciter.getStates(); + ASSERT_EQ(initialized.size(), 6U); + EXPECT_NEAR(initialized[1], 1.0, 1e-12); + EXPECT_NEAR(initialized[4], 0.4, 1e-12); + EXPECT_NEAR(initialized[5], 0.4, 1e-12); + + std::vector state{0.9, 0.98, 0.1, 0.3, 1.1, 1.0}; + std::vector stateDerivative(state.size(), 0.0); + exciter.setState(0.0, state.data(), stateDerivative.data(), cLocalSolverMode); + inputs[exciterVoltageInLocation] = 1.05; + inputs[exciterVsetInLocation] = 1.02; + std::vector derivative(state.size(), 0.0); + exciter.derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode); + EXPECT_NEAR(derivative[1], 0.7, 1e-12); + EXPECT_NEAR(derivative[2], -0.15, 1e-12); + EXPECT_NEAR(derivative[3], 6.25, 1e-12); + EXPECT_NEAR(derivative[4], -1.6, 1e-12); + EXPECT_NEAR(derivative[5], 0.1, 1e-12); + std::vector residual(state.size(), 0.0); + exciter.residual(inputs, emptyStateData, residual.data(), cLocalSolverMode); + EXPECT_NEAR(residual[0], 0.2, 1e-12); + + state[3] = 0.6; + exciter.setState(0.0, state.data(), stateDerivative.data(), cLocalSolverMode); + exciter.derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode); + EXPECT_DOUBLE_EQ(derivative[3], 0.0); + std::vector roots(1, 0.0); + exciter.rootTest(inputs, emptyStateData, roots.data(), cLocalSolverMode); + EXPECT_LT(roots[0], 0.0); +} + +TEST(ExciterModelTests, Exac4HardLimitsAndFactoryClone) +{ + auto factory = CoreObjectFactory::instance(); + std::unique_ptr object(factory->createObject("exciter", "exac4")); + auto* exciter = dynamic_cast(object.get()); + ASSERT_NE(exciter, nullptr); + exciter->set("tr", 0.02); + exciter->set("tb", 0.2); + exciter->set("tc", 0.1); + exciter->set("ka", 20.0); + exciter->set("ta", 0.05); + exciter->set("vimax", 0.1); + exciter->set("vimin", -0.1); + exciter->set("vrmax", 0.5); + exciter->set("vrmin", -0.5); + exciter->dynInitializeA(0.0, 0); + exciter->setRootOffset(0, cLocalSolverMode); + IOdata inputs(exciterInputCount, 0.0); + inputs[exciterVoltageInLocation] = 1.0; + inputs[exciterVsetInLocation] = 1.0; + IOdata fieldSet(4, 0.0); + exciter->dynInitializeB(inputs, {0.25}, fieldSet); + std::vector state{0.25, 0.8, 0.0, 0.75}; + std::vector dstate(state.size(), 0.0); + exciter->setState(0.0, state.data(), dstate.data(), cLocalSolverMode); + std::vector residual(state.size(), 0.0); + exciter->residual(inputs, emptyStateData, residual.data(), cLocalSolverMode); + EXPECT_NEAR(residual[0], 0.25, 1e-12); + std::vector roots(2, 0.0); + exciter->rootTest(inputs, emptyStateData, roots.data(), cLocalSolverMode); + EXPECT_LT(roots[0], 0.0); + EXPECT_LT(roots[1], 0.0); + std::unique_ptr cloned(exciter->clone()); + auto* clonedExciter = dynamic_cast(cloned.get()); + ASSERT_NE(clonedExciter, nullptr); + EXPECT_DOUBLE_EQ(clonedExciter->get("vimax"), 0.1); +} + TEST(ExciterModelTests, Exst1HardLimitersAndRootsUseRegulatorOutput) { exciters::ExciterEXST1 exciter; @@ -494,7 +594,8 @@ TEST_F(ExciterTests, ExciterTest2AlgDiffTests) // exclist.insert(exclist.begin(), "none"); for (auto& excname : exclist) { - if (excname.starts_with("fmi") || (excname == "esst3a") || (excname == "exst1")) { + if (excname.starts_with("fmi") || (excname == "esst3a") || (excname == "exst1") || + (excname == "exac1") || (excname == "exac2") || (excname == "exac4")) { continue; } gds = readSimXMLFile(fileName); @@ -541,7 +642,8 @@ TEST_F(ExciterTests, ExciterAlgDiffJacobianTests) // exclist.insert(exclist.begin(), "none"); for (auto& excname : exclist) { - if (excname.starts_with("fmi") || (excname == "esst3a") || (excname == "exst1")) { + if (excname.starts_with("fmi") || (excname == "esst3a") || (excname == "exst1") || + (excname == "exac1") || (excname == "exac2") || (excname == "exac4")) { continue; } gds = readSimXMLFile(fileName); diff --git a/test/test_files/andes_tests/ieee14_exac1.dyr b/test/test_files/andes_tests/ieee14_exac1.dyr new file mode 100644 index 00000000..9bf5ef2f --- /dev/null +++ b/test/test_files/andes_tests/ieee14_exac1.dyr @@ -0,0 +1,4 @@ + 2 'EXAC1' 1 0.031 0.23 0.12 41.0 + 0.034 7.2 -4.3 0.8 0.08 + 1.3 0.0 0.0 1.0 0.0 + 0.0 1.0 0.0 / diff --git a/test/test_files/andes_tests/ieee14_exac2.dyr b/test/test_files/andes_tests/ieee14_exac2.dyr new file mode 100644 index 00000000..2e639b0c --- /dev/null +++ b/test/test_files/andes_tests/ieee14_exac2.dyr @@ -0,0 +1,5 @@ + 2 'EXAC2' 1 0.031 0.23 0.12 41.0 + 0.034 7.2 -4.3 1.0 7.2 + -4.3 0.8 1.0 0.0 0.08 + 1.3 0.0 0.0 1.0 0.0 + 0.0 0.0 1.0 0.0 / diff --git a/test/test_files/andes_tests/ieee14_exac4.dyr b/test/test_files/andes_tests/ieee14_exac4.dyr new file mode 100644 index 00000000..661974e4 --- /dev/null +++ b/test/test_files/andes_tests/ieee14_exac4.dyr @@ -0,0 +1,3 @@ + 2 'EXAC4' 1 0.031 0.41 -0.37 0.12 + 0.23 41.0 0.034 7.2 -4.3 + 0.0 /