From 6fa4f0f07d20ee144dd5af34f719565012146bd4 Mon Sep 17 00:00:00 2001 From: Philip Top Date: Sun, 6 Sep 2026 20:55:31 -0700 Subject: [PATCH 1/3] add AC8B, ESST2A, IEEET3, and AC7B exciters --- src/fileInput/gridDynReadDYR.cpp | 97 +++ src/griddyn/CMakeLists.txt | 9 + src/griddyn/exciters/Exciter.cpp | 8 + src/griddyn/exciters/ExciterAC7B.cpp | 746 ++++++++++++++++++ src/griddyn/exciters/ExciterAC7B.h | 155 ++++ src/griddyn/exciters/ExciterAC8B.cpp | 680 ++++++++++++++++ src/griddyn/exciters/ExciterAC8B.h | 147 ++++ src/griddyn/exciters/ExciterESST2A.cpp | 520 ++++++++++++ src/griddyn/exciters/ExciterESST2A.h | 126 +++ src/griddyn/exciters/ExciterIEEET3.cpp | 482 +++++++++++ src/griddyn/exciters/ExciterIEEET3.h | 119 +++ src/griddyn/exciters/ExciterSignalHelper.h | 203 +++++ src/utilities/Saturation.cpp | 8 +- test/andesTests/testAndesDyrReader.cpp | 114 +++ test/componentTests/testExciters.cpp | 401 +++++++++- test/libraryTests/testSaturation.cpp | 13 + test/test_files/andes_tests/ieee14_ac7b.dyr | 1 + test/test_files/andes_tests/ieee14_ac8b.dyr | 1 + test/test_files/andes_tests/ieee14_esst2a.dyr | 1 + test/test_files/andes_tests/ieee14_ieeet3.dyr | 1 + 20 files changed, 3826 insertions(+), 6 deletions(-) create mode 100644 src/griddyn/exciters/ExciterAC7B.cpp create mode 100644 src/griddyn/exciters/ExciterAC7B.h create mode 100644 src/griddyn/exciters/ExciterAC8B.cpp create mode 100644 src/griddyn/exciters/ExciterAC8B.h create mode 100644 src/griddyn/exciters/ExciterESST2A.cpp create mode 100644 src/griddyn/exciters/ExciterESST2A.h create mode 100644 src/griddyn/exciters/ExciterIEEET3.cpp create mode 100644 src/griddyn/exciters/ExciterIEEET3.h create mode 100644 src/griddyn/exciters/ExciterSignalHelper.h create mode 100644 test/test_files/andes_tests/ieee14_ac7b.dyr create mode 100644 test/test_files/andes_tests/ieee14_ac8b.dyr create mode 100644 test/test_files/andes_tests/ieee14_esst2a.dyr create mode 100644 test/test_files/andes_tests/ieee14_ieeet3.dyr diff --git a/src/fileInput/gridDynReadDYR.cpp b/src/fileInput/gridDynReadDYR.cpp index c79c2fdf..466f158e 100644 --- a/src/fileInput/gridDynReadDYR.cpp +++ b/src/fileInput/gridDynReadDYR.cpp @@ -60,8 +60,12 @@ namespace { void loadESDC1A(CoreObject* parentObject, stringVec& tokens); void loadESDC2A(CoreObject* parentObject, stringVec& tokens); void loadIEEET1(CoreObject* parentObject, stringVec& tokens); + void loadIEEET3(CoreObject* parentObject, stringVec& tokens); void loadIEEEX1(CoreObject* parentObject, stringVec& tokens); + void loadAC7B(CoreObject* parentObject, stringVec& tokens); + void loadAC8B(CoreObject* parentObject, stringVec& tokens); void loadESST1A(CoreObject* parentObject, stringVec& tokens); + void loadESST2A(CoreObject* parentObject, stringVec& tokens); void loadESST3A(CoreObject* parentObject, stringVec& tokens); void loadESST4B(CoreObject* parentObject, stringVec& tokens); void loadEXPIC1(CoreObject* parentObject, stringVec& tokens); @@ -147,10 +151,18 @@ void loadDyr(CoreObject* parentObject, loadESDC2A(parentObject, lineTokens); } else if (type == "'IEEET1'") { loadIEEET1(parentObject, lineTokens); + } else if (type == "'IEEET3'") { + loadIEEET3(parentObject, lineTokens); } else if (type == "'IEEEX1'") { loadIEEEX1(parentObject, lineTokens); + } else if (type == "'AC7B'") { + loadAC7B(parentObject, lineTokens); + } else if (type == "'AC8B'") { + loadAC8B(parentObject, lineTokens); } else if (type == "'ESST1A'") { loadESST1A(parentObject, lineTokens); + } else if (type == "'ESST2A'") { + loadESST2A(parentObject, lineTokens); } else if (type == "'ESST3A'") { loadESST3A(parentObject, lineTokens); } else if (type == "'ESST4B'") { @@ -466,6 +478,24 @@ namespace { gen->add(exciterModel); } + void loadIEEET3(CoreObject* parentObject, stringVec& tokens) + { + if (tokens.size() != 15U) { + throw InvalidParameterValue("IEEET3 DYR record must contain 15 fields"); + } + auto* gen = requireDyrGenerator(parentObject, tokens, "IEEET3"); + const auto params = gmlc::utilities::str2vector(tokens, kNullVal); + auto* model = + static_cast(CoreObjectFactory::instance()->createObject("exciter", "ieeet3")); + // PSS/E order verified against the ANDES psse-dyr.yaml IEEET3 schema. + static constexpr std::array names{ + "tr", "ka", "ta", "vrmax", "vrmin", "vbmax", "ke", "te", "kf", "tf", "kp", "ki"}; + for (std::size_t index = 0; index < names.size(); ++index) { + model->set(names[index], params[index + 3]); + } + gen->add(model); + } + void loadIEEEX1(CoreObject* parentObject, stringVec& tokens) { if (tokens.size() != 19U) { @@ -498,6 +528,54 @@ namespace { gen->add(exciterModel); } + void loadAC7B(CoreObject* parentObject, stringVec& tokens) + { + if (tokens.size() != 30U) { + throw InvalidParameterValue("AC7B DYR record must contain 30 fields"); + } + auto* gen = requireDyrGenerator(parentObject, tokens, "AC7B"); + const auto params = gmlc::utilities::str2vector(tokens, kNullVal); + auto* model = + static_cast(CoreObjectFactory::instance()->createObject("exciter", "ac7b")); + // AC7B uses the 27-parameter order of the IEEE/OpenIPSL block: + // TR, KPR, KIR, KDR, TDR, VRMAX, VRMIN, KPA, KIA, VAMAX, VAMIN, + // KP, KL, TE, KC, KD, KE, KF1, KF2, KF3, TF3, VEMIN, VFEMAX, + // E1, SE1, E2, SE2. + static constexpr std::array names{"tr", "kpr", "kir", "kdr", + "tdr", "vrmax", "vrmin", "kpa", + "kia", "vamax", "vamin", "kp", + "kl", "te", "kc", "kd", + "ke", "kf1", "kf2", "kf3", + "tf3", "vemin", "vfemax", "e1", + "se1", "e2", "se2"}; + for (std::size_t index = 0; index < names.size(); ++index) { + model->set(names[index], params[index + 3]); + } + gen->add(model); + } + + void loadAC8B(CoreObject* parentObject, stringVec& tokens) + { + if (tokens.size() != 24U) { + throw InvalidParameterValue("AC8B DYR record must contain 24 fields"); + } + auto* gen = requireDyrGenerator(parentObject, tokens, "AC8B"); + const auto params = gmlc::utilities::str2vector(tokens, kNullVal); + auto* model = + static_cast(CoreObjectFactory::instance()->createObject("exciter", "ac8b")); + // Exact ANDES psse-dyr.yaml order. + static constexpr std::array names{"tr", "kpr", "kir", "kdr", + "tdr", "vpmax", "vpmin", "vrmax", + "vrmin", "vfemax", "vemin", "ta", + "ka", "te", "kc", "kd", + "ke", "e1", "se1", "e2", + "se2"}; + for (std::size_t index = 0; index < names.size(); ++index) { + model->set(names[index], params[index + 3]); + } + gen->add(model); + } + void loadESST1A(CoreObject* parentObject, stringVec& tokens) { if (tokens.size() != 23U) { @@ -533,6 +611,25 @@ namespace { gen->add(model); } + void loadESST2A(CoreObject* parentObject, stringVec& tokens) + { + if (tokens.size() != 16U) { + throw InvalidParameterValue("ESST2A DYR record must contain 16 fields"); + } + auto* gen = requireDyrGenerator(parentObject, tokens, "ESST2A"); + const auto params = gmlc::utilities::str2vector(tokens, kNullVal); + auto* model = + static_cast(CoreObjectFactory::instance()->createObject("exciter", "esst2a")); + // GridDyn's compact 13-parameter layout for the OpenIPSL ESST2A core. + // Separate VUEL/VOEL selector fields are intentionally not accepted. + static constexpr std::array names{ + "tr", "ka", "ta", "vrmax", "vrmin", "kp", "ki", "kc", "kf", "tf", "ke", "te", "efdmax"}; + for (std::size_t index = 0; index < names.size(); ++index) { + model->set(names[index], params[index + 3]); + } + gen->add(model); + } + void loadESST3A(CoreObject* parentObject, stringVec& tokens) { auto* gen = requireDyrGenerator(parentObject, tokens, "ESST3A"); diff --git a/src/griddyn/CMakeLists.txt b/src/griddyn/CMakeLists.txt index 864ac122..bfe11b88 100644 --- a/src/griddyn/CMakeLists.txt +++ b/src/griddyn/CMakeLists.txt @@ -91,10 +91,14 @@ set(exciter_headers stabilizers/StabilizerST2CUT.h exciters/ExciterDC1A.h exciters/ExciterDC2A.h + exciters/ExciterAC7B.h + exciters/ExciterAC8B.h exciters/ExciterESST1A.h + exciters/ExciterESST2A.h exciters/ExciterESST3A.h exciters/ExciterESST4B.h exciters/ExciterEXPIC1.h + exciters/ExciterSignalHelper.h exciters/StaticExciterRectifier.h exciters/ExciterEXAC1.h exciters/ExciterESAC1A.h @@ -104,6 +108,7 @@ set(exciter_headers exciters/ExciterEXST1.h exciters/ExciterIEEEtype1.h exciters/ExciterIEEEtype2.h + exciters/ExciterIEEET3.h exciters/ExciterIEEEX1.h exciters/ExciterSEXS.h exciters/ExciterSCRX.h @@ -114,7 +119,10 @@ set(exciter_sources exciters/Exciter.cpp exciters/ExciterDC1A.cpp exciters/ExciterDC2A.cpp + exciters/ExciterAC7B.cpp + exciters/ExciterAC8B.cpp exciters/ExciterESST1A.cpp + exciters/ExciterESST2A.cpp exciters/ExciterESST3A.cpp exciters/ExciterESST4B.cpp exciters/ExciterEXPIC1.cpp @@ -127,6 +135,7 @@ set(exciter_sources exciters/ExciterEXST1.cpp exciters/ExciterIEEEtype1.cpp exciters/ExciterIEEEtype2.cpp + exciters/ExciterIEEET3.cpp exciters/ExciterIEEEX1.cpp exciters/ExciterSEXS.cpp exciters/ExciterSCRX.cpp diff --git a/src/griddyn/exciters/Exciter.cpp b/src/griddyn/exciters/Exciter.cpp index a43692f0..a7374b3b 100644 --- a/src/griddyn/exciters/Exciter.cpp +++ b/src/griddyn/exciters/Exciter.cpp @@ -10,11 +10,14 @@ #include "../GridComponentHelperClasses.h" #include "../GridPrimary.h" #include "../gridDynDefinitions.hpp" +#include "ExciterAC7B.h" +#include "ExciterAC8B.h" #include "ExciterDC1A.h" #include "ExciterDC2A.h" #include "ExciterESAC1A.h" #include "ExciterESAC6A.h" #include "ExciterESST1A.h" +#include "ExciterESST2A.h" #include "ExciterESST3A.h" #include "ExciterESST4B.h" #include "ExciterEXAC1.h" @@ -22,6 +25,7 @@ #include "ExciterEXAC4.h" #include "ExciterEXPIC1.h" #include "ExciterEXST1.h" +#include "ExciterIEEET3.h" #include "ExciterIEEEX1.h" #include "ExciterIEEEtype1.h" #include "ExciterIEEEtype2.h" @@ -49,7 +53,10 @@ namespace exciters { static ChildTypeFactory gfeEsdc1a("exciter", "esdc1a"); // NOLINT static ChildTypeFactory gfeEsdc2a("exciter", "esdc2a"); // NOLINT static ChildTypeFactory gfeExdc2("exciter", "exdc2"); // NOLINT + static ChildTypeFactory gfeAc7b("exciter", "ac7b"); // NOLINT + static ChildTypeFactory gfeAc8b("exciter", "ac8b"); // NOLINT static ChildTypeFactory gfeEsst1a("exciter", "esst1a"); // NOLINT + static ChildTypeFactory gfeEsst2a("exciter", "esst2a"); // NOLINT static ChildTypeFactory gfeEsst3a("exciter", "esst3a"); // NOLINT static ChildTypeFactory gfeEsst4b("exciter", "esst4b"); // NOLINT static ChildTypeFactory gfeExpic1("exciter", "expic1"); // NOLINT @@ -61,6 +68,7 @@ namespace exciters { static ChildTypeFactory gfeExst1("exciter", "exst1"); // NOLINT static ChildTypeFactory gfeType1("exciter", "type1"); // NOLINT ChildTypeFactory gFeIeeet1("exciter", "ieeet1"); // NOLINT + ChildTypeFactory gFeIeeet3("exciter", "ieeet3"); // NOLINT ChildTypeFactory gFeIeeex1("exciter", "ieeex1"); // NOLINT static TypeFactory gfeDefault( // NOLINT "exciter", diff --git a/src/griddyn/exciters/ExciterAC7B.cpp b/src/griddyn/exciters/ExciterAC7B.cpp new file mode 100644 index 00000000..e2e2a1c8 --- /dev/null +++ b/src/griddyn/exciters/ExciterAC7B.cpp @@ -0,0 +1,746 @@ +/* + * 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 "ExciterAC7B.h" + +#include "ExciterSignalHelper.h" +#include "StaticExciterRectifier.h" +#include "core/CoreExceptions.h" +#include "core/CoreObjectTemplates.hpp" +#include "utilities/MatrixData.hpp" +#include +#include +#include +#include +#include +#include +#include + +namespace griddyn::exciters { +// NOLINTBEGIN(readability-math-missing-parentheses) +namespace { + constexpr index_t voltageMeasurementState = 0; + constexpr index_t pidIntegralState = 1; + constexpr index_t pidDerivativeState = 2; + constexpr index_t piIntegralState = 3; + constexpr index_t feedbackState = 4; + constexpr index_t exciterVoltageState = 5; + constexpr double initializationTolerance = 1e-7; + + bool finiteMachineSignal(double value) + { + return std::isfinite(value) && (std::abs(value) < 1e20); + } +} // namespace + +ExciterAC7B::ExciterAC7B(const std::string& objName): Exciter(objName) +{ + m_inputSize = exciterInputCount; + Vrmax = 5.79; + Vrmin = -5.79; +} + +CoreObject* ExciterAC7B::clone(CoreObject* obj) const +{ + auto* result = cloneBase(this, obj); + if (result == nullptr) { + return obj; + } + result->Tr = Tr; + result->Kpr = Kpr; + result->Kir = Kir; + result->Kdr = Kdr; + result->Tdr = Tdr; + result->Kpa = Kpa; + result->Kia = Kia; + result->Vamax = Vamax; + result->Vamin = Vamin; + result->Kp = Kp; + result->Kl = Kl; + result->Te = Te; + result->Kc = Kc; + result->Kd = Kd; + result->Ke = Ke; + result->Kf1 = Kf1; + result->Kf2 = Kf2; + result->Kf3 = Kf3; + result->Tf3 = Tf3; + result->Vemax = Vemax; + result->Vemin = Vemin; + result->E1 = E1; + result->Se1 = Se1; + result->E2 = E2; + result->Se2 = Se2; + result->saturation = saturation; + return result; +} + +void ExciterAC7B::dynObjectInitializeA(CoreTime /*time0*/, std::uint32_t /*flags*/) +{ + const std::array parameters{Tr, Kpr, Kir, Kdr, Tdr, Vrmax, Vrmin, + Kpa, Kia, Vamax, Vamin, Kp, Kl, Te, + Kc, Kd, Ke, Kf1, Kf2, Kf3, Tf3, + Vemax, Vemin, E1, Se1, E2, Se2}; + const bool disabledSaturation = (Se1 == 0.0) && (Se2 == 0.0); + const bool invalidSaturation = !disabledSaturation && + ((E1 <= 0.0) || (E2 <= 0.0) || (E1 == E2) || (Se1 <= 0.0) || (Se2 <= 0.0)); + if (std::any_of(parameters.begin(), + parameters.end(), + [](double value) { return !std::isfinite(value); }) || + (Tr < 0.0) || (Tdr < 0.0) || (Te <= 0.0) || (Tf3 < 0.0) || (Vrmax < Vrmin) || + (Vamax < Vamin) || (Vemax < Vemin) || invalidSaturation) { + throw InvalidParameterValue("AC7B gains, time constants, limits, or saturation points"); + } + saturation.setType(disabledSaturation ? + utilities::Saturation::SaturationType::NONE : + utilities::Saturation::SaturationType::CUTOFF_SCALED_QUADRATIC); + if (!disabledSaturation) { + saturation.setParam(E1, Se1, E2, Se2); + } + offsets.local().local.algSize = 1; + offsets.local().local.diffSize = stateCount; + offsets.local().local.jacSize = 120; + offsets.local().local.algRoots = 3; +} + +double ExciterAC7B::solveExciterVoltage(double fieldVoltage, double fieldCurrent) const +{ + const double loading = Kc * fieldCurrent; + const double middleMagnitude = + std::sqrt(((fieldVoltage * fieldVoltage) + (loading * loading)) / 0.75); + const double middleSign = (fieldVoltage != 0.0) ? fieldVoltage : loading; + const std::array candidates{fieldVoltage, + fieldVoltage + (0.577 * loading), + std::copysign(middleMagnitude, middleSign), + (fieldVoltage / 1.732) + loading, + loading}; + double bestVoltage = candidates.front(); + double bestMismatch = std::numeric_limits::infinity(); + for (double candidate : candidates) { + const double normalizedCurrent = (std::abs(candidate) > 1e-14) ? loading / candidate : 0.0; + const double output = candidate * detail::computeRectifierFactor(normalizedCurrent).factor; + const double mismatch = std::abs(output - fieldVoltage); + if (mismatch < bestMismatch) { + bestMismatch = mismatch; + bestVoltage = candidate; + } + } + if (!std::isfinite(bestVoltage) || + (bestMismatch > 1e-9 * std::max(1.0, std::abs(fieldVoltage)))) { + throw InvalidParameterValue("AC7B initial rectifier equation has no finite solution"); + } + return bestVoltage; +} + +double ExciterAC7B::solveFieldFeedbackLimit(double fieldCurrent, double initialVoltage) const +{ + const double target = Vemax - (Kd * fieldCurrent); + double voltage = initialVoltage; + if (!std::isfinite(voltage)) { + voltage = (std::abs(Ke) > 1e-12) ? target / Ke : Vemax; + } + for (int iteration = 0; iteration < 40; ++iteration) { + const auto sat = saturation.evaluate(voltage); + const double mismatch = ((Ke + sat.value) * voltage) - target; + if (std::abs(mismatch) <= 1e-11 * std::max(1.0, std::abs(target))) { + return voltage; + } + const double slope = Ke + sat.value + (voltage * sat.derivative); + if (!std::isfinite(slope) || (std::abs(slope) < 1e-12)) { + break; + } + voltage -= mismatch / slope; + } + throw InvalidParameterValue("AC7B exciter upper-limit equation has no finite solution"); +} + +ExciterAC7B::Evaluation + ExciterAC7B::evaluate(const IOdata& inputs, double fieldVoltage, const double state[]) const +{ + using detail::addSignals; + using detail::algebraicSignal; + using detail::clampSignal; + using detail::constantSignal; + using detail::divideSignals; + using detail::inputSignal; + using detail::multiplySignals; + using detail::rectifierFactorSignal; + using detail::scaleSignal; + using detail::stateSignal; + using detail::subtractSignals; + using Signal = detail::ExciterSignal; + + Evaluation evaluation; + std::array rates{}; + const Signal terminalVoltage = inputSignal(inputs, exciterVoltageInLocation); + const Signal measuredVoltage = (Tr > 0.0) ? + stateSignal(state[voltageMeasurementState], voltageMeasurementState) : + terminalVoltage; + if (Tr > 0.0) { + rates[voltageMeasurementState] = + scaleSignal(subtractSignals(terminalVoltage, measuredVoltage), 1.0 / Tr); + } + + const Signal exciterVoltage = + stateSignal(state[exciterVoltageState], exciterVoltageState); + const auto saturationData = saturation.evaluate(exciterVoltage.value); + const Signal saturationSignal = + detail::applyFunction(exciterVoltage, saturationData.value, saturationData.derivative); + const Signal fieldCurrent = inputSignal(inputs, exciterXadIfdInLocation); + const Signal saturationGain = addSignals(constantSignal(Ke), saturationSignal); + const Signal fieldFeedback = + addSignals(multiplySignals(saturationGain, exciterVoltage), scaleSignal(fieldCurrent, Kd)); + const Signal fieldOutput = + multiplySignals(rectifierFactorSignal( + divideSignals(scaleSignal(fieldCurrent, Kc), exciterVoltage)), + exciterVoltage); + + const Signal feedbackFilter = stateSignal(state[feedbackState], feedbackState); + Signal derivativeFeedback = constantSignal(0.0); + if (Tf3 > 0.0) { + rates[feedbackState] = + scaleSignal(subtractSignals(fieldFeedback, feedbackFilter), 1.0 / Tf3); + derivativeFeedback = scaleSignal(subtractSignals(fieldFeedback, feedbackFilter), Kf3 / Tf3); + } + + Signal input = addSignals(constantSignal(Vref + vBias - 1.0), + inputSignal(inputs, exciterVsetInLocation)); + input = addSignals(input, inputSignal(inputs, exciterVssInLocation)); + input = subtractSignals(subtractSignals(input, measuredVoltage), derivativeFeedback); + + const Signal derivativeState = + stateSignal(state[pidDerivativeState], pidDerivativeState); + Signal derivativeOutput = constantSignal(0.0); + if (Tdr > 0.0) { + rates[pidDerivativeState] = scaleSignal(subtractSignals(input, derivativeState), 1.0 / Tdr); + derivativeOutput = scaleSignal(subtractSignals(input, derivativeState), Kdr / Tdr); + } + const Signal pidIntegral = stateSignal(state[pidIntegralState], pidIntegralState); + const Signal pidDrive = + addSignals(addSignals(scaleSignal(input, Kpr), pidIntegral), derivativeOutput); + const Signal pidOutput = clampSignal(pidDrive, Vrmin, Vrmax); + const Signal pidIntegratorDrive = scaleSignal(input, Kir); + if (!opFlags[PID_LIMITED]) { + rates[pidIntegralState] = pidIntegratorDrive; + } + + const Signal regulatorFeedback = + addSignals(scaleSignal(algebraicSignal(fieldVoltage), Kf1), + scaleSignal(fieldFeedback, Kf2)); + const Signal piError = subtractSignals(pidOutput, regulatorFeedback); + const Signal piIntegral = stateSignal(state[piIntegralState], piIntegralState); + const Signal piDrive = addSignals(scaleSignal(piError, Kpa), piIntegral); + const Signal piOutput = clampSignal(piDrive, Vamin, Vamax); + const Signal piIntegratorDrive = scaleSignal(piError, Kia); + if (!opFlags[PI_LIMITED]) { + rates[piIntegralState] = piIntegratorDrive; + } + + const Signal lowerControlLimit = scaleSignal(fieldFeedback, -Kl); + Signal exciterControl = multiplySignals(scaleSignal(terminalVoltage, Kp), piOutput); + if (exciterControl.value < lowerControlLimit.value) { + exciterControl = lowerControlLimit; + } + // VE limits constrain the integrator state. The input limiter only + // imposes I_C >= -K_L V_FE; it does not clamp I_C to the VE bounds. + const Signal exciterDrive = subtractSignals(exciterControl, fieldFeedback); + if (!opFlags[EXCITER_LIMITED]) { + rates[exciterVoltageState] = scaleSignal(exciterDrive, 1.0 / Te); + } + + const Signal algebraicResidual = + subtractSignals(fieldOutput, algebraicSignal(fieldVoltage)); + evaluation.fieldOutput = fieldOutput.value; + evaluation.pidDrive = pidDrive.value; + evaluation.pidIntegratorDrive = pidIntegratorDrive.value; + evaluation.piDrive = piDrive.value; + evaluation.piIntegratorDrive = piIntegratorDrive.value; + evaluation.fieldFeedback = fieldFeedback.value; + evaluation.exciterDrive = exciterDrive.value; + evaluation.algebraicDerivative = algebraicResidual.algebraic; + evaluation.fieldStateDerivatives = algebraicResidual.state; + evaluation.fieldInputDerivatives = algebraicResidual.input; + for (index_t row = 0; row < stateCount; ++row) { + evaluation.rates[row] = rates[row].value; + evaluation.rateAlgebraicDerivatives[row] = rates[row].algebraic; + evaluation.rateStateDerivatives[row] = rates[row].state; + evaluation.rateInputDerivatives[row] = rates[row].input; + } + return evaluation; +} + +void ExciterAC7B::dynObjectInitializeB(const IOdata& inputs, + const IOdata& desiredOutput, + IOdata& fieldSet) +{ + if (inputs.size() < exciterInputCount || desiredOutput.empty() || + !std::isfinite(inputs[exciterVoltageInLocation]) || + !std::isfinite(inputs[exciterVsetInLocation]) || + !std::isfinite(inputs[exciterVssInLocation]) || !std::isfinite(desiredOutput[0]) || + (((Kc != 0.0) || (Kd != 0.0)) && !finiteMachineSignal(inputs[exciterXadIfdInLocation]))) { + throw InvalidParameterValue("AC7B initial signals"); + } + const double fieldVoltage = desiredOutput[0]; + const double fieldCurrent = inputs[exciterXadIfdInLocation]; + const double exciterVoltage = solveExciterVoltage(fieldVoltage, fieldCurrent); + const double fieldFeedback = + (Ke + saturation(exciterVoltage)) * exciterVoltage + (Kd * fieldCurrent); + const double terminalVoltage = std::max(1e-8, inputs[exciterVoltageInLocation]); + const double piOutput = fieldFeedback / (Kp * terminalVoltage); + const double pidOutput = (Kf1 * fieldVoltage) + (Kf2 * fieldFeedback); + const double rectifierMismatch = + (exciterVoltage * + detail::computeRectifierFactor( + (std::abs(exciterVoltage) > 1e-14) ? Kc * fieldCurrent / exciterVoltage : 0.0) + .factor) - + fieldVoltage; + if ((pidOutput < Vrmin - initializationTolerance) || + (pidOutput > Vrmax + initializationTolerance) || + (piOutput < Vamin - initializationTolerance) || + (piOutput > Vamax + initializationTolerance) || + (exciterVoltage < Vemin - initializationTolerance) || + (fieldFeedback > Vemax + initializationTolerance) || + (std::abs(rectifierMismatch) > initializationTolerance)) { + throw InvalidParameterValue( + "AC7B initial state outside limits: PID=" + std::to_string(pidOutput) + + ", PI=" + std::to_string(piOutput) + ", VE=" + std::to_string(exciterVoltage) + + ", VFE=" + std::to_string(fieldFeedback)); + } + + m_state[0] = fieldVoltage; + double* state = m_state.data() + 1; + state[voltageMeasurementState] = inputs[exciterVoltageInLocation]; + state[pidDerivativeState] = 0.0; + state[pidIntegralState] = pidOutput; + state[piIntegralState] = piOutput; + state[feedbackState] = fieldFeedback; + state[exciterVoltageState] = exciterVoltage; + vBias = state[voltageMeasurementState] - Vref - (inputs[exciterVsetInLocation] - 1.0) - + inputs[exciterVssInLocation]; + fieldSet.resize(std::max(fieldSet.size(), 2U)); + fieldSet[exciterVsetInLocation] = Vref; + std::fill(m_dstate_dt.begin(), m_dstate_dt.end(), 0.0); + updateLimitFlags(inputs, state, true); +} + +void ExciterAC7B::residual(const IOdata& inputs, + const StateData& stateData, + double resid[], + const SolverMode& sMode) +{ + const auto loc = offsets.getLocations(stateData, resid, sMode, this); + const auto evaluation = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc); + if (hasAlgebraic(sMode)) { + loc.destLoc[0] = evaluation.fieldOutput - loc.algStateLoc[0]; + } + if (hasDifferential(sMode)) { + for (index_t index = 0; index < loc.diffSize; ++index) { + loc.destDiffLoc[index] = evaluation.rates[index] - loc.dstateLoc[index]; + } + } +} + +void ExciterAC7B::derivative(const IOdata& inputs, + const StateData& stateData, + double deriv[], + const SolverMode& sMode) +{ + if (!hasDifferential(sMode)) { + return; + } + const auto loc = offsets.getLocations(stateData, deriv, sMode, this); + const auto evaluation = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc); + std::copy_n(evaluation.rates.begin(), loc.diffSize, loc.destDiffLoc); +} + +void ExciterAC7B::algebraicUpdate(const IOdata& inputs, + const StateData& stateData, + double update[], + const SolverMode& sMode, + double /*alpha*/) +{ + if (hasAlgebraic(sMode)) { + const auto loc = offsets.getLocations(stateData, update, sMode, this); + loc.destLoc[0] = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc).fieldOutput; + } +} + +void ExciterAC7B::jacobianElements(const IOdata& inputs, + const StateData& stateData, + MatrixData& matrixData, + const IOlocs& inputLocs, + const SolverMode& sMode) +{ + const auto loc = offsets.getLocations(stateData, sMode, this); + const auto evaluation = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc); + if (hasAlgebraic(sMode)) { + matrixData.assign(loc.algOffset, loc.algOffset, evaluation.algebraicDerivative); + for (index_t column = 0; column < stateCount; ++column) { + matrixData.assign(loc.algOffset, + loc.diffOffset + column, + evaluation.fieldStateDerivatives[column]); + } + for (index_t input = 0; input < exciterInputCount; ++input) { + matrixData.assignCheckCol(loc.algOffset, + inputLocs[input], + evaluation.fieldInputDerivatives[input]); + } + } + if (!hasDifferential(sMode)) { + return; + } + for (index_t row = 0; row < stateCount; ++row) { + matrixData.assign(loc.diffOffset + row, + loc.algOffset, + evaluation.rateAlgebraicDerivatives[row]); + for (index_t column = 0; column < stateCount; ++column) { + double value = evaluation.rateStateDerivatives[row][column]; + if (row == column) { + value -= stateData.cj; + } + matrixData.assign(loc.diffOffset + row, loc.diffOffset + column, value); + } + for (index_t input = 0; input < exciterInputCount; ++input) { + matrixData.assignCheckCol(loc.diffOffset + row, + inputLocs[input], + evaluation.rateInputDerivatives[row][input]); + } + } +} + +void ExciterAC7B::timestep(CoreTime time, const IOdata& inputs, const SolverMode& /*sMode*/) +{ + derivative(inputs, emptyStateData, m_dstate_dt.data(), cLocalSolverMode); + const double timeStep = time - prevTime; + for (index_t index = 0; index < stateCount; ++index) { + m_state[index + 1] += timeStep * m_dstate_dt[index + 1]; + } + updateLimitFlags(inputs, m_state.data() + 1, true); + m_state[0] = evaluate(inputs, m_state[0], m_state.data() + 1).fieldOutput; + prevTime = time; +} + +bool ExciterAC7B::updateLimitFlags(const IOdata& inputs, double state[], bool projectStates) +{ + constexpr double tolerance = 1e-7; + const auto evaluation = evaluate(inputs, m_state[0], state); + bool changed = false; + const auto setLimited = [this, + &changed](int limitedFlag, int highFlag, bool limited, bool high) { + if (opFlags[limitedFlag] != limited) { + opFlags.set(limitedFlag, limited); + changed = true; + } + if (limited && (opFlags[highFlag] != high)) { + opFlags.set(highFlag, high); + changed = true; + } + }; + + if (opFlags[PID_LIMITED]) { + const bool release = opFlags[PID_LIMIT_HIGH] ? + ((evaluation.pidDrive < Vrmax - tolerance) || (evaluation.pidIntegratorDrive < 0.0)) : + ((evaluation.pidDrive > Vrmin + tolerance) || (evaluation.pidIntegratorDrive > 0.0)); + if (release) { + setLimited(PID_LIMITED, PID_LIMIT_HIGH, false, false); + } + } else if ((evaluation.pidDrive >= Vrmax - tolerance) && + (evaluation.pidIntegratorDrive > 0.0)) { + setLimited(PID_LIMITED, PID_LIMIT_HIGH, true, true); + } else if ((evaluation.pidDrive <= Vrmin + tolerance) && + (evaluation.pidIntegratorDrive < 0.0)) { + setLimited(PID_LIMITED, PID_LIMIT_HIGH, true, false); + } + + if (opFlags[PI_LIMITED]) { + const bool release = opFlags[PI_LIMIT_HIGH] ? + ((evaluation.piDrive < Vamax - tolerance) || (evaluation.piIntegratorDrive < 0.0)) : + ((evaluation.piDrive > Vamin + tolerance) || (evaluation.piIntegratorDrive > 0.0)); + if (release) { + setLimited(PI_LIMITED, PI_LIMIT_HIGH, false, false); + } + } else if ((evaluation.piDrive >= Vamax - tolerance) && (evaluation.piIntegratorDrive > 0.0)) { + setLimited(PI_LIMITED, PI_LIMIT_HIGH, true, true); + } else if ((evaluation.piDrive <= Vamin + tolerance) && (evaluation.piIntegratorDrive < 0.0)) { + setLimited(PI_LIMITED, PI_LIMIT_HIGH, true, false); + } + + const double fieldCurrent = inputs[exciterXadIfdInLocation]; + const bool exciterAbove = evaluation.fieldFeedback > Vemax + tolerance; + const bool exciterBelow = state[exciterVoltageState] < Vemin - tolerance; + if (projectStates) { + state[exciterVoltageState] = + std::max(state[exciterVoltageState], static_cast(Vemin)); + const double voltage = state[exciterVoltageState]; + const double feedback = (Ke + saturation(voltage)) * voltage + (Kd * fieldCurrent); + if (feedback > Vemax) { + state[exciterVoltageState] = solveFieldFeedbackLimit(fieldCurrent, voltage); + } + } + if (opFlags[EXCITER_LIMITED]) { + const bool release = opFlags[EXCITER_LIMIT_HIGH] ? + ((evaluation.fieldFeedback < Vemax - tolerance) || (evaluation.exciterDrive < 0.0)) : + ((state[exciterVoltageState] > Vemin + tolerance) || (evaluation.exciterDrive > 0.0)); + if (release) { + setLimited(EXCITER_LIMITED, EXCITER_LIMIT_HIGH, false, false); + } + } else if (exciterAbove || + ((evaluation.fieldFeedback >= Vemax - tolerance) && + (evaluation.exciterDrive > 0.0))) { + setLimited(EXCITER_LIMITED, EXCITER_LIMIT_HIGH, true, true); + } else if (exciterBelow || + ((state[exciterVoltageState] <= Vemin + tolerance) && + (evaluation.exciterDrive < 0.0))) { + setLimited(EXCITER_LIMITED, EXCITER_LIMIT_HIGH, true, false); + } + return changed; +} + +void ExciterAC7B::rootTest(const IOdata& inputs, + const StateData& stateData, + double roots[], + const SolverMode& sMode) +{ + const auto loc = offsets.getLocations(stateData, sMode, this); + const auto evaluation = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc); + const index_t root = offsets.getRootOffset(sMode); + if (opFlags[PID_LIMITED]) { + roots[root] = opFlags[PID_LIMIT_HIGH] ? + std::max(Vrmax - evaluation.pidDrive, -evaluation.pidIntegratorDrive) : + std::max(evaluation.pidDrive - Vrmin, evaluation.pidIntegratorDrive); + } else { + roots[root] = std::min(Vrmax - evaluation.pidDrive, evaluation.pidDrive - Vrmin); + opFlags.set(PID_LIMIT_HIGH, evaluation.pidDrive >= Vrmax); + } + if (opFlags[PI_LIMITED]) { + roots[root + 1] = opFlags[PI_LIMIT_HIGH] ? + std::max(Vamax - evaluation.piDrive, -evaluation.piIntegratorDrive) : + std::max(evaluation.piDrive - Vamin, evaluation.piIntegratorDrive); + } else { + roots[root + 1] = std::min(Vamax - evaluation.piDrive, evaluation.piDrive - Vamin); + opFlags.set(PI_LIMIT_HIGH, evaluation.piDrive >= Vamax); + } + if (opFlags[EXCITER_LIMITED]) { + roots[root + 2] = opFlags[EXCITER_LIMIT_HIGH] ? + std::max(Vemax - evaluation.fieldFeedback, -evaluation.exciterDrive) : + std::max(loc.diffStateLoc[exciterVoltageState] - Vemin, evaluation.exciterDrive); + } else { + roots[root + 2] = std::min(Vemax - evaluation.fieldFeedback, + loc.diffStateLoc[exciterVoltageState] - Vemin); + opFlags.set(EXCITER_LIMIT_HIGH, evaluation.fieldFeedback >= Vemax); + } +} + +void ExciterAC7B::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) || (rootMask[root + 2] != 0)) { + updateLimitFlags(inputs, m_state.data() + 1, true); + alert(this, JAC_COUNT_CHANGE); + } +} + +ChangeCode ExciterAC7B::rootCheck(const IOdata& inputs, + const StateData& /*stateData*/, + const SolverMode& /*sMode*/, + CheckLevel /*level*/) +{ + if (updateLimitFlags(inputs, m_state.data() + 1, true)) { + alert(this, JAC_COUNT_CHANGE); + return ChangeCode::JACOBIAN_CHANGE; + } + return ChangeCode::NO_CHANGE; +} + +stringVec ExciterAC7B::localStateNames() const +{ + return {"efd", "vmeas", "pidint", "pidder", "piint", "xf3", "ve"}; +} + +index_t ExciterAC7B::findIndex(std::string_view field, const SolverMode& sMode) const +{ + if ((field == "efd") || (field == "field")) { + return offsets.getAlgOffset(sMode); + } + const auto offset = offsets.getDiffOffset(sMode); + if ((field == "vmeas") || (field == "vc")) { + return offset + voltageMeasurementState; + } + if (field == "pidint") { + return offset + pidIntegralState; + } + if (field == "pidder") { + return offset + pidDerivativeState; + } + if (field == "piint") { + return offset + piIntegralState; + } + if (field == "xf3") { + return offset + feedbackState; + } + if ((field == "ve") || (field == "exciter")) { + return offset + exciterVoltageState; + } + return kInvalidLocation; +} + +void ExciterAC7B::set(std::string_view param, std::string_view val) +{ + Exciter::set(param, val); +} + +void ExciterAC7B::set(std::string_view param, double val, units::unit unitType) +{ + if (param == "tr") { + Tr = val; + } else if (param == "kpr") { + Kpr = val; + } else if (param == "kir") { + Kir = val; + } else if ((param == "kdr") || (param == "kd_pid")) { + Kdr = val; + } else if ((param == "tdr") || (param == "td")) { + Tdr = val; + } else if (param == "kpa") { + Kpa = val; + } else if (param == "kia") { + Kia = val; + } else if ((param == "vamax") || (param == "va_max")) { + Vamax = val; + } else if ((param == "vamin") || (param == "va_min")) { + Vamin = val; + } else if (param == "kp") { + Kp = val; + } else if (param == "kl") { + Kl = val; + } else if (param == "te") { + Te = val; + } else if (param == "kc") { + Kc = val; + } else if (param == "kd") { + Kd = val; + } else if (param == "ke") { + Ke = val; + } else if (param == "kf1") { + Kf1 = val; + } else if (param == "kf2") { + Kf2 = val; + } else if (param == "kf3") { + Kf3 = val; + } else if (param == "tf3") { + Tf3 = val; + } else if ((param == "vfemax") || (param == "vemax")) { + Vemax = val; + } else if (param == "vemin") { + Vemin = val; + } else if (param == "e1") { + E1 = val; + } else if (param == "se1") { + Se1 = val; + } else if (param == "e2") { + E2 = val; + } else if (param == "se2") { + Se2 = val; + } else { + Exciter::set(param, val, unitType); + } +} + +double ExciterAC7B::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 == "kpr") { + return Kpr; + } + if (param == "kir") { + return Kir; + } + if ((param == "kdr") || (param == "kd_pid")) { + return Kdr; + } + if ((param == "tdr") || (param == "td")) { + return Tdr; + } + if (param == "kpa") { + return Kpa; + } + if (param == "kia") { + return Kia; + } + if ((param == "vamax") || (param == "va_max")) { + return Vamax; + } + if ((param == "vamin") || (param == "va_min")) { + return Vamin; + } + if (param == "kp") { + return Kp; + } + if (param == "kl") { + return Kl; + } + if (param == "te") { + return Te; + } + if (param == "kc") { + return Kc; + } + if (param == "kd") { + return Kd; + } + if (param == "ke") { + return Ke; + } + if (param == "kf1") { + return Kf1; + } + if (param == "kf2") { + return Kf2; + } + if (param == "kf3") { + return Kf3; + } + if (param == "tf3") { + return Tf3; + } + if ((param == "vfemax") || (param == "vemax")) { + return Vemax; + } + if (param == "vemin") { + return Vemin; + } + 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); +} +// NOLINTEND(readability-math-missing-parentheses) +} // namespace griddyn::exciters diff --git a/src/griddyn/exciters/ExciterAC7B.h b/src/griddyn/exciters/ExciterAC7B.h new file mode 100644 index 00000000..74630e99 --- /dev/null +++ b/src/griddyn/exciters/ExciterAC7B.h @@ -0,0 +1,155 @@ +/* + * 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 +#include +#include +#include + +namespace griddyn::exciters { +/** IEEE/PSS/E AC7B alternator-rectifier excitation system. + * + * AC7B uses the AC8B-style PID voltage regulator and rotating exciter, with an + * additional inner PI loop and excitation-system stabilizer feedback. The + * implemented equations are + * \f[ + * V_i=V_{ref}+V_{bias}+V_{set}-1+V_{SS}-V_m- + * K_{F3}(V_{FE}-x_{F3})/T_{F3}, + * \f] + * \f[ + * T_R\dot V_m=V_T-V_m,\quad T_{DR}\dot x_D=V_i-x_D,\\ + * V_P^*=K_{PR}V_i+x_I+K_{DR}(V_i-x_D)/T_{DR},\\ + * V_P=\operatorname{lim}(V_P^*,V_{RMIN},V_{RMAX}),\quad \dot x_I=K_{IR}V_i,\\ + * e_A=V_P-K_{F1}E_{fd}-K_{F2}V_{FE},\\ + * V_A^*=K_{PA}e_A+x_A,\quad + * V_A=\operatorname{lim}(V_A^*,V_{AMIN},V_{AMAX}),\quad \dot x_A=K_{IA}e_A, + * \f] + * \f[ + * I_C=\max(K_PV_TV_A,-K_LV_{FE}),\quad + * T_E\dot V_E=I_C-\{(K_E+S_E(V_E))V_E+K_DI_f\},\quad + * E_{fd}=V_EF_{EX}(K_CI_f/V_E). + * \f] + * where \f$T_{F3}\dot x_{F3}=V_{FE}-x_{F3}\f$ and + * \f$S_E(V_E)=B(V_E-A)^2/V_E\f$ above its fitted cutoff (zero below). + * The PID and PI integrators are held only when their complete pre-limit + * outputs are saturated outward. The \f$V_E\f$ state is held at + * \f$V_{EMIN}\f$ or the implicit \f$V_{FE}=V_{FEMAX}\f$ boundary; the + * control input \f$I_C\f$ is not clipped to those state bounds. + * + * @par Equation source + * OpenIPSL `Electrical.Controls.PSSE.ES.AC7B` and its shared rotating + * exciter, rectifier, and no-windup PID/PI blocks, commit 8155c73f. + */ +class ExciterAC7B final: public Exciter { + public: + enum AC7BFlags { + PID_LIMITED = OBJECT_FLAG5, + PID_LIMIT_HIGH = OBJECT_FLAG6, + PI_LIMITED = OBJECT_FLAG7, + PI_LIMIT_HIGH = OBJECT_FLAG8, + EXCITER_LIMITED = OBJECT_FLAG9, + EXCITER_LIMIT_HIGH = OBJECT_FLAG10, + }; + + explicit ExciterAC7B(const std::string& objName = "exciterAC7B_#"); + 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 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; + 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; + + private: + static constexpr index_t stateCount = 6; + struct Evaluation { + double fieldOutput = 0.0; + double algebraicDerivative = 0.0; + std::array fieldStateDerivatives{}; + std::array fieldInputDerivatives{}; + std::array rates{}; + std::array rateAlgebraicDerivatives{}; + std::array, stateCount> rateStateDerivatives{}; + std::array, stateCount> rateInputDerivatives{}; + double pidDrive = 0.0; + double pidIntegratorDrive = 0.0; + double piDrive = 0.0; + double piIntegratorDrive = 0.0; + double fieldFeedback = 0.0; + double exciterDrive = 0.0; + }; + Evaluation evaluate(const IOdata& inputs, double fieldVoltage, const double state[]) const; + double solveExciterVoltage(double fieldVoltage, double fieldCurrent) const; + double solveFieldFeedbackLimit(double fieldCurrent, double initialVoltage) const; + bool updateLimitFlags(const IOdata& inputs, double state[], bool projectStates); + + model_parameter Tr = 0.0; + model_parameter Kpr = 4.24; + model_parameter Kir = 4.24; + model_parameter Kdr = 0.0; + model_parameter Tdr = 0.0; + model_parameter Kpa = 65.36; + model_parameter Kia = 59.69; + model_parameter Vamax = 1.0; + model_parameter Vamin = -0.95; + model_parameter Kp = 4.96; + model_parameter Kl = 10.0; + model_parameter Te = 1.1; + model_parameter Kc = 0.18; + model_parameter Kd = 0.02; + model_parameter Ke = 1.0; + model_parameter Kf1 = 0.212; + model_parameter Kf2 = 0.0; + model_parameter Kf3 = 0.0; + model_parameter Tf3 = 0.0; + model_parameter Vemax = 6.9; + model_parameter Vemin = -99.0; + model_parameter E1 = 6.3; + model_parameter Se1 = 0.44; + model_parameter E2 = 4.725; + model_parameter Se2 = 0.075; + utilities::Saturation saturation{ + utilities::Saturation::SaturationType::CUTOFF_SCALED_QUADRATIC}; +}; +} // namespace griddyn::exciters diff --git a/src/griddyn/exciters/ExciterAC8B.cpp b/src/griddyn/exciters/ExciterAC8B.cpp new file mode 100644 index 00000000..33e169c2 --- /dev/null +++ b/src/griddyn/exciters/ExciterAC8B.cpp @@ -0,0 +1,680 @@ +/* + * 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 "ExciterAC8B.h" + +#include "ExciterSignalHelper.h" +#include "StaticExciterRectifier.h" +#include "core/CoreExceptions.h" +#include "core/CoreObjectTemplates.hpp" +#include "utilities/MatrixData.hpp" +#include +#include +#include +#include +#include +#include +#include + +namespace griddyn::exciters { +// NOLINTBEGIN(readability-math-missing-parentheses) +namespace { + constexpr index_t voltageMeasurementState = 0; + constexpr index_t pidIntegralState = 1; + constexpr index_t pidDerivativeState = 2; + constexpr index_t regulatorState = 3; + constexpr index_t exciterVoltageState = 4; + constexpr double initializationTolerance = 1e-7; + + bool finiteMachineSignal(double value) + { + return std::isfinite(value) && (std::abs(value) < 1e20); + } +} // namespace + +ExciterAC8B::ExciterAC8B(const std::string& objName): Exciter(objName) +{ + m_inputSize = exciterInputCount; + Ka = 40.0; + Ta = 0.04; + Vrmax = 7.3; + Vrmin = 1.0; +} + +CoreObject* ExciterAC8B::clone(CoreObject* obj) const +{ + auto* result = cloneBase(this, obj); + if (result == nullptr) { + return obj; + } + result->Tr = Tr; + result->Kpr = Kpr; + result->Kir = Kir; + result->Kdr = Kdr; + result->Tdr = Tdr; + result->Vpmax = Vpmax; + result->Vpmin = Vpmin; + result->Vemax = Vemax; + result->Vemin = Vemin; + result->Te = Te; + result->Kc = Kc; + result->Kd = Kd; + result->Ke = Ke; + result->E1 = E1; + result->Se1 = Se1; + result->E2 = E2; + result->Se2 = Se2; + result->saturation = saturation; + return result; +} + +void ExciterAC8B::dynObjectInitializeA(CoreTime /*time0*/, std::uint32_t /*flags*/) +{ + const std::array parameters{Tr, Kpr, Kir, Kdr, Tdr, Vpmax, Vpmin, + Vrmax, Vrmin, Vemax, Vemin, Ta, Ka, Te, + Kc, Kd, Ke, E1, Se1, E2}; + const bool disabledSaturation = (Se1 == 0.0) && (Se2 == 0.0); + const bool invalidSaturation = !disabledSaturation && + ((E1 <= 0.0) || (E2 <= 0.0) || (E1 == E2) || (Se1 <= 0.0) || (Se2 <= 0.0)); + if (std::any_of(parameters.begin(), + parameters.end(), + [](double value) { return !std::isfinite(value); }) || + !std::isfinite(Se2) || (Tr < 0.0) || (Tdr < 0.0) || (Ta <= 0.0) || (Te <= 0.0) || + (Ka <= 0.0) || (Vpmax < Vpmin) || (Vrmax < Vrmin) || (Vemax < Vemin) || invalidSaturation) { + throw InvalidParameterValue("AC8B gains, time constants, limits, or saturation points"); + } + saturation.setType(disabledSaturation ? + utilities::Saturation::SaturationType::NONE : + utilities::Saturation::SaturationType::CUTOFF_SCALED_QUADRATIC); + if (!disabledSaturation) { + saturation.setParam(E1, Se1, E2, Se2); + } + offsets.local().local.algSize = 1; + offsets.local().local.diffSize = stateCount; + offsets.local().local.jacSize = 96; + offsets.local().local.algRoots = 3; +} + +double ExciterAC8B::solveExciterVoltage(double fieldVoltage, double fieldCurrent) const +{ + const double loading = Kc * fieldCurrent; + const double middleMagnitude = + std::sqrt(((fieldVoltage * fieldVoltage) + (loading * loading)) / 0.75); + const double middleSign = (fieldVoltage != 0.0) ? fieldVoltage : loading; + const std::array candidates{fieldVoltage, + fieldVoltage + (0.577 * loading), + std::copysign(middleMagnitude, middleSign), + (fieldVoltage / 1.732) + loading, + loading}; + double bestVoltage = candidates.front(); + double bestMismatch = std::numeric_limits::infinity(); + for (double candidate : candidates) { + const double normalizedCurrent = (std::abs(candidate) > 1e-14) ? loading / candidate : 0.0; + const double output = candidate * detail::computeRectifierFactor(normalizedCurrent).factor; + const double mismatch = std::abs(output - fieldVoltage); + if (mismatch < bestMismatch) { + bestMismatch = mismatch; + bestVoltage = candidate; + } + } + if (!std::isfinite(bestVoltage) || + (bestMismatch > 1e-9 * std::max(1.0, std::abs(fieldVoltage)))) { + throw InvalidParameterValue("AC8B initial rectifier equation has no finite solution"); + } + return bestVoltage; +} + +double ExciterAC8B::solveFieldFeedbackLimit(double fieldCurrent, double initialVoltage) const +{ + const double target = Vemax - (Kd * fieldCurrent); + double voltage = initialVoltage; + if (!std::isfinite(voltage)) { + voltage = (std::abs(Ke) > 1e-12) ? target / Ke : Vemax; + } + for (int iteration = 0; iteration < 40; ++iteration) { + const auto sat = saturation.evaluate(voltage); + const double mismatch = ((Ke + sat.value) * voltage) - target; + if (std::abs(mismatch) <= 1e-11 * std::max(1.0, std::abs(target))) { + return voltage; + } + const double slope = Ke + sat.value + (voltage * sat.derivative); + if (!std::isfinite(slope) || (std::abs(slope) < 1e-12)) { + break; + } + voltage -= mismatch / slope; + } + throw InvalidParameterValue("AC8B exciter upper-limit equation has no finite solution"); +} + +ExciterAC8B::Evaluation + ExciterAC8B::evaluate(const IOdata& inputs, double fieldVoltage, const double state[]) const +{ + using detail::addSignals; + using detail::algebraicSignal; + using detail::clampSignal; + using detail::constantSignal; + using detail::divideSignals; + using detail::inputSignal; + using detail::multiplySignals; + using detail::rectifierFactorSignal; + using detail::scaleSignal; + using detail::stateSignal; + using detail::subtractSignals; + using Signal = detail::ExciterSignal; + + Evaluation evaluation; + std::array rates{}; + const Signal measuredVoltage = (Tr > 0.0) ? + stateSignal(state[voltageMeasurementState], voltageMeasurementState) : + inputSignal(inputs, exciterVoltageInLocation); + if (Tr > 0.0) { + rates[voltageMeasurementState] = + scaleSignal(subtractSignals(inputSignal(inputs, exciterVoltageInLocation), + measuredVoltage), + 1.0 / Tr); + } + + Signal input = addSignals(constantSignal(Vref + vBias - 1.0), + inputSignal(inputs, exciterVsetInLocation)); + input = addSignals(input, inputSignal(inputs, exciterVssInLocation)); + input = subtractSignals(input, measuredVoltage); + + const Signal derivativeState = + stateSignal(state[pidDerivativeState], pidDerivativeState); + Signal derivativeOutput = constantSignal(0.0); + if (Tdr > 0.0) { + rates[pidDerivativeState] = scaleSignal(subtractSignals(input, derivativeState), 1.0 / Tdr); + derivativeOutput = scaleSignal(subtractSignals(input, derivativeState), Kdr / Tdr); + } + + const Signal pidIntegral = stateSignal(state[pidIntegralState], pidIntegralState); + const Signal pidDrive = + addSignals(addSignals(scaleSignal(input, Kpr), pidIntegral), derivativeOutput); + const Signal pidOutput = clampSignal(pidDrive, Vpmin, Vpmax); + const Signal integratorDrive = scaleSignal(input, Kir); + if (!opFlags[PID_LIMITED]) { + rates[pidIntegralState] = integratorDrive; + } + + const Signal regulator = stateSignal(state[regulatorState], regulatorState); + const Signal regulatorDrive = subtractSignals(scaleSignal(pidOutput, Ka), regulator); + if (!opFlags[REGULATOR_LIMITED]) { + rates[regulatorState] = scaleSignal(regulatorDrive, 1.0 / Ta); + } + + const Signal exciterVoltage = + stateSignal(state[exciterVoltageState], exciterVoltageState); + const auto saturationData = saturation.evaluate(exciterVoltage.value); + const Signal saturationSignal = + detail::applyFunction(exciterVoltage, saturationData.value, saturationData.derivative); + const Signal fieldCurrent = inputSignal(inputs, exciterXadIfdInLocation); + const Signal saturationGain = addSignals(constantSignal(Ke), saturationSignal); + const Signal fieldFeedback = + addSignals(multiplySignals(saturationGain, exciterVoltage), scaleSignal(fieldCurrent, Kd)); + // VEMIN and VFEMAX bound the rotating-exciter state, not its control input. + // Therefore T_E dV_E/dt = V_R - V_FE uses the unclipped regulator output. + const Signal exciterDrive = subtractSignals(regulator, fieldFeedback); + if (!opFlags[EXCITER_LIMITED]) { + rates[exciterVoltageState] = scaleSignal(exciterDrive, 1.0 / Te); + } + + const Signal normalizedCurrent = divideSignals(scaleSignal(fieldCurrent, Kc), exciterVoltage); + const Signal fieldOutput = + multiplySignals(rectifierFactorSignal(normalizedCurrent), exciterVoltage); + const Signal algebraicResidual = + subtractSignals(fieldOutput, algebraicSignal(fieldVoltage)); + + evaluation.fieldOutput = fieldOutput.value; + evaluation.pidDrive = pidDrive.value; + evaluation.pidIntegratorDrive = integratorDrive.value; + evaluation.regulatorDrive = regulatorDrive.value; + evaluation.fieldFeedback = fieldFeedback.value; + evaluation.exciterDrive = exciterDrive.value; + evaluation.algebraicDerivative = algebraicResidual.algebraic; + evaluation.fieldStateDerivatives = algebraicResidual.state; + evaluation.fieldInputDerivatives = algebraicResidual.input; + for (index_t row = 0; row < stateCount; ++row) { + evaluation.rates[row] = rates[row].value; + evaluation.rateAlgebraicDerivatives[row] = rates[row].algebraic; + evaluation.rateStateDerivatives[row] = rates[row].state; + evaluation.rateInputDerivatives[row] = rates[row].input; + } + return evaluation; +} + +void ExciterAC8B::dynObjectInitializeB(const IOdata& inputs, + const IOdata& desiredOutput, + IOdata& fieldSet) +{ + if (inputs.size() < exciterInputCount || desiredOutput.empty() || + !std::isfinite(inputs[exciterVoltageInLocation]) || + !std::isfinite(inputs[exciterVsetInLocation]) || + !std::isfinite(inputs[exciterVssInLocation]) || !std::isfinite(desiredOutput[0]) || + (((Kc != 0.0) || (Kd != 0.0)) && !finiteMachineSignal(inputs[exciterXadIfdInLocation]))) { + throw InvalidParameterValue("AC8B initial signals"); + } + const double fieldVoltage = desiredOutput[0]; + const double fieldCurrent = inputs[exciterXadIfdInLocation]; + const double exciterVoltage = solveExciterVoltage(fieldVoltage, fieldCurrent); + const double fieldFeedback = + (Ke + saturation(exciterVoltage)) * exciterVoltage + (Kd * fieldCurrent); + const double pidOutput = fieldFeedback / Ka; + const double rectifierMismatch = + (exciterVoltage * + detail::computeRectifierFactor( + (std::abs(exciterVoltage) > 1e-14) ? Kc * fieldCurrent / exciterVoltage : 0.0) + .factor) - + fieldVoltage; + if ((pidOutput < Vpmin - initializationTolerance) || + (pidOutput > Vpmax + initializationTolerance) || + (fieldFeedback < Vrmin - initializationTolerance) || + (fieldFeedback > Vrmax + initializationTolerance) || + (exciterVoltage < Vemin - initializationTolerance) || + (fieldFeedback > Vemax + initializationTolerance) || + (std::abs(rectifierMismatch) > initializationTolerance)) { + throw InvalidParameterValue("AC8B initial regulator output outside limits"); + } + + m_state[0] = fieldVoltage; + double* state = m_state.data() + 1; + state[voltageMeasurementState] = inputs[exciterVoltageInLocation]; + state[pidDerivativeState] = 0.0; + state[pidIntegralState] = pidOutput; + state[regulatorState] = fieldFeedback; + state[exciterVoltageState] = exciterVoltage; + vBias = state[voltageMeasurementState] - Vref - (inputs[exciterVsetInLocation] - 1.0) - + inputs[exciterVssInLocation]; + fieldSet.resize(std::max(fieldSet.size(), 2U)); + fieldSet[exciterVsetInLocation] = Vref; + std::fill(m_dstate_dt.begin(), m_dstate_dt.end(), 0.0); + updateLimitFlags(inputs, state, true); +} + +void ExciterAC8B::residual(const IOdata& inputs, + const StateData& stateData, + double resid[], + const SolverMode& sMode) +{ + const auto loc = offsets.getLocations(stateData, resid, sMode, this); + const auto evaluation = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc); + if (hasAlgebraic(sMode)) { + loc.destLoc[0] = evaluation.fieldOutput - loc.algStateLoc[0]; + } + if (hasDifferential(sMode)) { + for (index_t index = 0; index < loc.diffSize; ++index) { + loc.destDiffLoc[index] = evaluation.rates[index] - loc.dstateLoc[index]; + } + } +} + +void ExciterAC8B::derivative(const IOdata& inputs, + const StateData& stateData, + double deriv[], + const SolverMode& sMode) +{ + if (!hasDifferential(sMode)) { + return; + } + const auto loc = offsets.getLocations(stateData, deriv, sMode, this); + const auto evaluation = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc); + std::copy_n(evaluation.rates.begin(), loc.diffSize, loc.destDiffLoc); +} + +void ExciterAC8B::algebraicUpdate(const IOdata& inputs, + const StateData& stateData, + double update[], + const SolverMode& sMode, + double /*alpha*/) +{ + if (hasAlgebraic(sMode)) { + const auto loc = offsets.getLocations(stateData, update, sMode, this); + loc.destLoc[0] = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc).fieldOutput; + } +} + +void ExciterAC8B::jacobianElements(const IOdata& inputs, + const StateData& stateData, + MatrixData& matrixData, + const IOlocs& inputLocs, + const SolverMode& sMode) +{ + const auto loc = offsets.getLocations(stateData, sMode, this); + const auto evaluation = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc); + if (hasAlgebraic(sMode)) { + matrixData.assign(loc.algOffset, loc.algOffset, evaluation.algebraicDerivative); + for (index_t column = 0; column < stateCount; ++column) { + matrixData.assign(loc.algOffset, + loc.diffOffset + column, + evaluation.fieldStateDerivatives[column]); + } + for (index_t input = 0; input < exciterInputCount; ++input) { + matrixData.assignCheckCol(loc.algOffset, + inputLocs[input], + evaluation.fieldInputDerivatives[input]); + } + } + if (!hasDifferential(sMode)) { + return; + } + for (index_t row = 0; row < stateCount; ++row) { + matrixData.assign(loc.diffOffset + row, + loc.algOffset, + evaluation.rateAlgebraicDerivatives[row]); + for (index_t column = 0; column < stateCount; ++column) { + double value = evaluation.rateStateDerivatives[row][column]; + if (row == column) { + value -= stateData.cj; + } + matrixData.assign(loc.diffOffset + row, loc.diffOffset + column, value); + } + for (index_t input = 0; input < exciterInputCount; ++input) { + matrixData.assignCheckCol(loc.diffOffset + row, + inputLocs[input], + evaluation.rateInputDerivatives[row][input]); + } + } +} + +void ExciterAC8B::timestep(CoreTime time, const IOdata& inputs, const SolverMode& /*sMode*/) +{ + derivative(inputs, emptyStateData, m_dstate_dt.data(), cLocalSolverMode); + const double timeStep = time - prevTime; + for (index_t index = 0; index < stateCount; ++index) { + m_state[index + 1] += timeStep * m_dstate_dt[index + 1]; + } + updateLimitFlags(inputs, m_state.data() + 1, true); + m_state[0] = evaluate(inputs, m_state[0], m_state.data() + 1).fieldOutput; + prevTime = time; +} + +bool ExciterAC8B::updateLimitFlags(const IOdata& inputs, double state[], bool projectStates) +{ + constexpr double tolerance = 1e-7; + const auto evaluation = evaluate(inputs, m_state[0], state); + bool changed = false; + const auto setLimited = [this, + &changed](int limitedFlag, int highFlag, bool limited, bool high) { + if (opFlags[limitedFlag] != limited) { + opFlags.set(limitedFlag, limited); + changed = true; + } + if (limited && (opFlags[highFlag] != high)) { + opFlags.set(highFlag, high); + changed = true; + } + }; + + if (opFlags[PID_LIMITED]) { + const bool release = opFlags[PID_LIMIT_HIGH] ? + ((evaluation.pidDrive < Vpmax - tolerance) || (evaluation.pidIntegratorDrive < 0.0)) : + ((evaluation.pidDrive > Vpmin + tolerance) || (evaluation.pidIntegratorDrive > 0.0)); + if (release) { + setLimited(PID_LIMITED, PID_LIMIT_HIGH, false, false); + } + } else if ((evaluation.pidDrive >= Vpmax - tolerance) && + (evaluation.pidIntegratorDrive > 0.0)) { + setLimited(PID_LIMITED, PID_LIMIT_HIGH, true, true); + } else if ((evaluation.pidDrive <= Vpmin + tolerance) && + (evaluation.pidIntegratorDrive < 0.0)) { + setLimited(PID_LIMITED, PID_LIMIT_HIGH, true, false); + } + + const bool regulatorAbove = state[regulatorState] > Vrmax + tolerance; + const bool regulatorBelow = state[regulatorState] < Vrmin - tolerance; + if (projectStates) { + state[regulatorState] = std::clamp(state[regulatorState], + static_cast(Vrmin), + static_cast(Vrmax)); + } + if (opFlags[REGULATOR_LIMITED]) { + const bool release = opFlags[REGULATOR_LIMIT_HIGH] ? (evaluation.regulatorDrive < 0.0) : + (evaluation.regulatorDrive > 0.0); + if (release) { + setLimited(REGULATOR_LIMITED, REGULATOR_LIMIT_HIGH, false, false); + } + } else if (regulatorAbove || + ((state[regulatorState] >= Vrmax - tolerance) && + (evaluation.regulatorDrive > 0.0))) { + setLimited(REGULATOR_LIMITED, REGULATOR_LIMIT_HIGH, true, true); + } else if (regulatorBelow || + ((state[regulatorState] <= Vrmin + tolerance) && + (evaluation.regulatorDrive < 0.0))) { + setLimited(REGULATOR_LIMITED, REGULATOR_LIMIT_HIGH, true, false); + } + + const double fieldCurrent = inputs[exciterXadIfdInLocation]; + const bool exciterAbove = evaluation.fieldFeedback > Vemax + tolerance; + const bool exciterBelow = state[exciterVoltageState] < Vemin - tolerance; + if (projectStates) { + state[exciterVoltageState] = + std::max(state[exciterVoltageState], static_cast(Vemin)); + const double voltage = state[exciterVoltageState]; + const double feedback = (Ke + saturation(voltage)) * voltage + (Kd * fieldCurrent); + if (feedback > Vemax) { + state[exciterVoltageState] = solveFieldFeedbackLimit(fieldCurrent, voltage); + } + } + if (opFlags[EXCITER_LIMITED]) { + const bool release = opFlags[EXCITER_LIMIT_HIGH] ? + ((evaluation.fieldFeedback < Vemax - tolerance) || (evaluation.exciterDrive < 0.0)) : + ((state[exciterVoltageState] > Vemin + tolerance) || (evaluation.exciterDrive > 0.0)); + if (release) { + setLimited(EXCITER_LIMITED, EXCITER_LIMIT_HIGH, false, false); + } + } else if (exciterAbove || + ((evaluation.fieldFeedback >= Vemax - tolerance) && + (evaluation.exciterDrive > 0.0))) { + setLimited(EXCITER_LIMITED, EXCITER_LIMIT_HIGH, true, true); + } else if (exciterBelow || + ((state[exciterVoltageState] <= Vemin + tolerance) && + (evaluation.exciterDrive < 0.0))) { + setLimited(EXCITER_LIMITED, EXCITER_LIMIT_HIGH, true, false); + } + return changed; +} + +void ExciterAC8B::rootTest(const IOdata& inputs, + const StateData& stateData, + double roots[], + const SolverMode& sMode) +{ + const auto loc = offsets.getLocations(stateData, sMode, this); + const auto evaluation = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc); + const index_t root = offsets.getRootOffset(sMode); + if (opFlags[PID_LIMITED]) { + roots[root] = opFlags[PID_LIMIT_HIGH] ? + std::max(Vpmax - evaluation.pidDrive, -evaluation.pidIntegratorDrive) : + std::max(evaluation.pidDrive - Vpmin, evaluation.pidIntegratorDrive); + } else { + roots[root] = std::min(Vpmax - evaluation.pidDrive, evaluation.pidDrive - Vpmin); + opFlags.set(PID_LIMIT_HIGH, evaluation.pidDrive >= Vpmax); + } + if (opFlags[REGULATOR_LIMITED]) { + roots[root + 1] = + opFlags[REGULATOR_LIMIT_HIGH] ? -evaluation.regulatorDrive : evaluation.regulatorDrive; + } else { + roots[root + 1] = std::min(Vrmax - loc.diffStateLoc[regulatorState], + loc.diffStateLoc[regulatorState] - Vrmin); + opFlags.set(REGULATOR_LIMIT_HIGH, loc.diffStateLoc[regulatorState] >= Vrmax); + } + if (opFlags[EXCITER_LIMITED]) { + roots[root + 2] = opFlags[EXCITER_LIMIT_HIGH] ? + std::max(Vemax - evaluation.fieldFeedback, -evaluation.exciterDrive) : + std::max(loc.diffStateLoc[exciterVoltageState] - Vemin, evaluation.exciterDrive); + } else { + roots[root + 2] = std::min(Vemax - evaluation.fieldFeedback, + loc.diffStateLoc[exciterVoltageState] - Vemin); + opFlags.set(EXCITER_LIMIT_HIGH, evaluation.fieldFeedback >= Vemax); + } +} + +void ExciterAC8B::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) || (rootMask[root + 2] != 0)) { + updateLimitFlags(inputs, m_state.data() + 1, true); + alert(this, JAC_COUNT_CHANGE); + } +} + +ChangeCode ExciterAC8B::rootCheck(const IOdata& inputs, + const StateData& /*stateData*/, + const SolverMode& /*sMode*/, + CheckLevel /*level*/) +{ + if (updateLimitFlags(inputs, m_state.data() + 1, true)) { + alert(this, JAC_COUNT_CHANGE); + return ChangeCode::JACOBIAN_CHANGE; + } + return ChangeCode::NO_CHANGE; +} + +stringVec ExciterAC8B::localStateNames() const +{ + return {"efd", "vmeas", "pidint", "pidder", "vr", "ve"}; +} + +index_t ExciterAC8B::findIndex(std::string_view field, const SolverMode& sMode) const +{ + if ((field == "efd") || (field == "field")) { + return offsets.getAlgOffset(sMode); + } + const auto offset = offsets.getDiffOffset(sMode); + if ((field == "vmeas") || (field == "vc")) { + return offset + voltageMeasurementState; + } + if (field == "pidint") { + return offset + pidIntegralState; + } + if (field == "pidder") { + return offset + pidDerivativeState; + } + if ((field == "vr") || (field == "regulator")) { + return offset + regulatorState; + } + if ((field == "ve") || (field == "exciter")) { + return offset + exciterVoltageState; + } + return kInvalidLocation; +} + +void ExciterAC8B::set(std::string_view param, std::string_view val) +{ + Exciter::set(param, val); +} + +void ExciterAC8B::set(std::string_view param, double val, units::unit unitType) +{ + if (param == "tr") { + Tr = val; + } else if ((param == "kpr") || (param == "kp_pid") || (param == "kp")) { + Kpr = val; + } else if ((param == "kir") || (param == "ki_pid") || (param == "ki")) { + Kir = val; + } else if ((param == "kdr") || (param == "kd_pid")) { + Kdr = val; + } else if ((param == "tdr") || (param == "td")) { + Tdr = val; + } else if ((param == "vpmax") || (param == "vpidmax")) { + Vpmax = val; + } else if ((param == "vpmin") || (param == "vpidmin")) { + Vpmin = val; + } else if ((param == "vfemax") || (param == "vemax")) { + Vemax = val; + } else if (param == "vemin") { + Vemin = val; + } else if (param == "te") { + Te = val; + } else if (param == "kc") { + Kc = val; + } else if (param == "kd") { + Kd = val; + } else if (param == "ke") { + Ke = val; + } else if (param == "e1") { + E1 = val; + } else if (param == "se1") { + Se1 = val; + } else if (param == "e2") { + E2 = val; + } else if (param == "se2") { + Se2 = val; + } else { + Exciter::set(param, val, unitType); + } +} + +double ExciterAC8B::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 == "kpr") || (param == "kp_pid") || (param == "kp")) { + return Kpr; + } + if ((param == "kir") || (param == "ki_pid") || (param == "ki")) { + return Kir; + } + if ((param == "kdr") || (param == "kd_pid")) { + return Kdr; + } + if ((param == "tdr") || (param == "td")) { + return Tdr; + } + if ((param == "vpmax") || (param == "vpidmax")) { + return Vpmax; + } + if ((param == "vpmin") || (param == "vpidmin")) { + return Vpmin; + } + if ((param == "vfemax") || (param == "vemax")) { + return Vemax; + } + if (param == "vemin") { + return Vemin; + } + if (param == "te") { + return Te; + } + 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); +} +// NOLINTEND(readability-math-missing-parentheses) +} // namespace griddyn::exciters diff --git a/src/griddyn/exciters/ExciterAC8B.h b/src/griddyn/exciters/ExciterAC8B.h new file mode 100644 index 00000000..a113db1d --- /dev/null +++ b/src/griddyn/exciters/ExciterAC8B.h @@ -0,0 +1,147 @@ +/* + * 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 +#include +#include +#include + +namespace griddyn::exciters { +/** IEEE/PSS/E AC8B alternator-rectifier excitation system. + * + * AC8B is implemented as a terminal-voltage transducer, non-windup PID voltage + * regulator, non-windup regulator lag, and rotating AC exciter with + * demagnetizing feedback and rectifier commutation-voltage drop. With + * measured voltage \f$V_m\f$, PID integral \f$x_I\f$, derivative lag + * \f$x_D\f$, regulator output \f$V_R\f$, exciter internal voltage \f$V_E\f$, + * and \f$I_f=X_{ad}I_{fd}\f$, + * \f[ + * \begin{aligned} + * T_R\dot V_m &= V_T-V_m,\\ + * V_i &= V_{ref}+V_{bias}+V_{set}-1+V_{SS}-V_m,\\ + * T_D\dot x_D &= V_i-x_D,\\ + * V_P^* &= K_{PR}V_i+x_I+K_{DR}(V_i-x_D)/T_D,\\ + * V_P &= \operatorname{lim}(V_P^*,V_{PMIN},V_{PMAX}),\\ + * \dot x_I &= K_{IR}V_i\quad\hbox{unless }V_P\hbox{ is saturated outward},\\ + * T_A\dot V_R &= K_AV_P-V_R,\\ + * I_N &= K_CI_f/V_E,\qquad E_{fd}=V_EF_{EX}(I_N),\\ + * S_E(V_E)&=\begin{cases}0,&V_E& 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; + 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; + + private: + static constexpr index_t stateCount = 5; + struct Evaluation { + double fieldOutput = 0.0; + double algebraicDerivative = 0.0; + std::array fieldStateDerivatives{}; + std::array fieldInputDerivatives{}; + std::array rates{}; + std::array rateAlgebraicDerivatives{}; + std::array, stateCount> rateStateDerivatives{}; + std::array, stateCount> rateInputDerivatives{}; + double pidDrive = 0.0; + double pidIntegratorDrive = 0.0; + double regulatorDrive = 0.0; + double fieldFeedback = 0.0; + double exciterDrive = 0.0; + }; + Evaluation evaluate(const IOdata& inputs, double fieldVoltage, const double state[]) const; + double solveExciterVoltage(double fieldVoltage, double fieldCurrent) const; + double solveFieldFeedbackLimit(double fieldCurrent, double initialVoltage) const; + bool updateLimitFlags(const IOdata& inputs, double state[], bool projectStates); + + model_parameter Tr = 0.01; + model_parameter Kpr = 10.0; + model_parameter Kir = 10.0; + model_parameter Kdr = 10.0; + model_parameter Tdr = 0.2; + model_parameter Vpmax = 999.0; + model_parameter Vpmin = -999.0; + model_parameter Vemax = 999.0; + model_parameter Vemin = -999.0; + model_parameter Te = 0.8; + 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::CUTOFF_SCALED_QUADRATIC}; +}; +} // namespace griddyn::exciters diff --git a/src/griddyn/exciters/ExciterESST2A.cpp b/src/griddyn/exciters/ExciterESST2A.cpp new file mode 100644 index 00000000..a189999b --- /dev/null +++ b/src/griddyn/exciters/ExciterESST2A.cpp @@ -0,0 +1,520 @@ +/* + * 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 "ExciterESST2A.h" + +#include "ExciterSignalHelper.h" +#include "StaticExciterRectifier.h" +#include "core/CoreExceptions.h" +#include "core/CoreObjectTemplates.hpp" +#include "utilities/MatrixData.hpp" +#include +#include +#include +#include +#include +#include +#include + +namespace griddyn::exciters { +// NOLINTBEGIN(readability-math-missing-parentheses) +namespace { + constexpr index_t voltageMeasurementState = 0; + constexpr index_t regulatorState = 1; + constexpr index_t feedbackState = 2; + constexpr index_t fieldState = 3; + constexpr double initializationTolerance = 1e-7; + + bool finiteMachineSignal(double value) + { + return std::isfinite(value) && (std::abs(value) < 1e20); + } +} // namespace + +ExciterESST2A::ExciterESST2A(const std::string& objName): Exciter(objName) +{ + m_inputSize = exciterInputCount; + Ka = 240.0; + Ta = 0.01; + Vrmax = 4.5; + Vrmin = -4.5; +} + +CoreObject* ExciterESST2A::clone(CoreObject* obj) const +{ + auto* result = cloneBase(this, obj); + if (result == nullptr) { + return obj; + } + result->Tr = Tr; + result->Kp = Kp; + result->Ki = Ki; + result->Kc = Kc; + result->Kf = Kf; + result->Tf = Tf; + result->Ke = Ke; + result->Te = Te; + result->Efdmax = Efdmax; + return result; +} + +void ExciterESST2A::dynObjectInitializeA(CoreTime /*time0*/, std::uint32_t /*flags*/) +{ + const std::array parameters{ + Tr, Vrmax, Vrmin, Ka, Ta, Kp, Ki, Kc, Kf, Tf, Ke, Te, Efdmax, Vref}; + if (std::any_of(parameters.begin(), + parameters.end(), + [](double value) { return !std::isfinite(value); }) || + (Tr < 0.0) || (Ka <= 0.0) || (Ta <= 0.0) || (Vrmax < Vrmin) || (Tf <= 0.0) || (Te <= 0.0) || + (Efdmax < 0.0)) { + throw InvalidParameterValue("ESST2A gains, time constants, or limits"); + } + offsets.local().local.algSize = 1; + offsets.local().local.diffSize = stateCount; + offsets.local().local.jacSize = 80; + offsets.local().local.algRoots = 2; +} + +double ExciterESST2A::rectifierVoltage(const IOdata& inputs) const +{ + if ((Kp <= 0.0) && (Ki <= 0.0)) { + return 1.0; + } + return detail::computeRectifierData( + inputs, Kp, Ki, Kc, 0.0, 0.0, std::numeric_limits::infinity()) + .voltage; +} + +ExciterESST2A::Evaluation + ExciterESST2A::evaluate(const IOdata& inputs, double fieldVoltage, const double state[]) const +{ + using detail::addSignals; + using detail::algebraicSignal; + using detail::clampSignal; + using detail::constantSignal; + using detail::inputSignal; + using detail::scaleSignal; + using detail::stateSignal; + using detail::subtractSignals; + using Signal = detail::ExciterSignal; + + Evaluation evaluation; + std::array rates{}; + const Signal measuredVoltage = (Tr > 0.0) ? + stateSignal(state[voltageMeasurementState], voltageMeasurementState) : + inputSignal(inputs, exciterVoltageInLocation); + if (Tr > 0.0) { + rates[voltageMeasurementState] = + scaleSignal(subtractSignals(inputSignal(inputs, exciterVoltageInLocation), + measuredVoltage), + 1.0 / Tr); + } + + Signal rectifier = constantSignal(1.0); + if ((Kp > 0.0) || (Ki > 0.0)) { + const auto data = detail::computeRectifierData( + inputs, Kp, Ki, Kc, 0.0, 0.0, std::numeric_limits::infinity()); + rectifier.value = data.voltage; + const std::array signalIndices{exciterIdInLocation, + exciterIqInLocation, + exciterVdInLocation, + exciterVqInLocation, + exciterXadIfdInLocation}; + for (index_t index = 0; index < 5; ++index) { + rectifier.input[signalIndices[index]] = data.derivatives[index]; + } + } + + const Signal fieldStateSignal = stateSignal(state[fieldState], fieldState); + const Signal feedbackStateSignal = stateSignal(state[feedbackState], feedbackState); + const Signal feedback = + scaleSignal(subtractSignals(fieldStateSignal, feedbackStateSignal), Kf / Tf); + rates[feedbackState] = + scaleSignal(subtractSignals(fieldStateSignal, feedbackStateSignal), 1.0 / Tf); + + Signal input = addSignals(constantSignal(Vref + vBias - 1.0), + inputSignal(inputs, exciterVsetInLocation)); + input = addSignals(input, inputSignal(inputs, exciterVssInLocation)); + input = subtractSignals(subtractSignals(input, measuredVoltage), feedback); + const Signal regulator = stateSignal(state[regulatorState], regulatorState); + const Signal regulatorDrive = subtractSignals(scaleSignal(input, Ka), regulator); + if (!opFlags[REGULATOR_LIMITED]) { + rates[regulatorState] = scaleSignal(regulatorDrive, 1.0 / Ta); + } + + const Signal fieldDrive = subtractSignals(detail::multiplySignals(rectifier, regulator), + scaleSignal(fieldStateSignal, Ke)); + if (!opFlags[FIELD_LIMITED]) { + rates[fieldState] = scaleSignal(fieldDrive, 1.0 / Te); + } + const Signal algebraicResidual = + subtractSignals(fieldStateSignal, algebraicSignal(fieldVoltage)); + + evaluation.fieldOutput = state[fieldState]; + evaluation.regulatorDrive = regulatorDrive.value; + evaluation.fieldDrive = fieldDrive.value; + evaluation.algebraicDerivative = algebraicResidual.algebraic; + evaluation.fieldStateDerivatives = algebraicResidual.state; + evaluation.fieldInputDerivatives = algebraicResidual.input; + for (index_t row = 0; row < stateCount; ++row) { + evaluation.rates[row] = rates[row].value; + evaluation.rateAlgebraicDerivatives[row] = rates[row].algebraic; + evaluation.rateStateDerivatives[row] = rates[row].state; + evaluation.rateInputDerivatives[row] = rates[row].input; + } + return evaluation; +} + +void ExciterESST2A::dynObjectInitializeB(const IOdata& inputs, + const IOdata& desiredOutput, + IOdata& fieldSet) +{ + if (inputs.size() < exciterInputCount || desiredOutput.empty() || + !std::isfinite(inputs[exciterVoltageInLocation]) || + !std::isfinite(inputs[exciterVsetInLocation]) || + !std::isfinite(inputs[exciterVssInLocation]) || !std::isfinite(desiredOutput[0])) { + throw InvalidParameterValue("ESST2A initial voltage signals"); + } + if ((Kp > 0.0) || (Ki > 0.0)) { + for (const index_t inputIndex : {exciterIdInLocation, + exciterIqInLocation, + exciterVdInLocation, + exciterVqInLocation, + exciterXadIfdInLocation}) { + if (!finiteMachineSignal(inputs[inputIndex])) { + throw InvalidParameterValue( + "ESST2A requires synchronous-machine dq and field signals"); + } + } + } + const double fieldVoltage = desiredOutput[0]; + const double rectifier = rectifierVoltage(inputs); + if (std::abs(rectifier) <= 1e-12) { + throw InvalidParameterValue("ESST2A initial rectifier voltage"); + } + const double regulator = Ke * fieldVoltage / rectifier; + if ((regulator < Vrmin - initializationTolerance) || + (regulator > Vrmax + initializationTolerance) || + (fieldVoltage < -initializationTolerance) || + (fieldVoltage > Efdmax + initializationTolerance)) { + throw InvalidParameterValue("ESST2A initial state outside limits"); + } + m_state[0] = fieldVoltage; + double* state = m_state.data() + 1; + state[voltageMeasurementState] = inputs[exciterVoltageInLocation]; + state[regulatorState] = regulator; + state[feedbackState] = fieldVoltage; + state[fieldState] = fieldVoltage; + vBias = (regulator / Ka) + state[voltageMeasurementState] - Vref - + (inputs[exciterVsetInLocation] - 1.0) - inputs[exciterVssInLocation]; + fieldSet.resize(std::max(fieldSet.size(), 2U)); + fieldSet[exciterVsetInLocation] = Vref; + std::fill(m_dstate_dt.begin(), m_dstate_dt.end(), 0.0); + updateLimitFlags(inputs, state, true); +} + +void ExciterESST2A::residual(const IOdata& inputs, + const StateData& stateData, + double resid[], + const SolverMode& sMode) +{ + const auto loc = offsets.getLocations(stateData, resid, sMode, this); + const auto evaluation = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc); + if (hasAlgebraic(sMode)) { + loc.destLoc[0] = evaluation.fieldOutput - loc.algStateLoc[0]; + } + if (hasDifferential(sMode)) { + for (index_t index = 0; index < loc.diffSize; ++index) { + loc.destDiffLoc[index] = evaluation.rates[index] - loc.dstateLoc[index]; + } + } +} + +void ExciterESST2A::derivative(const IOdata& inputs, + const StateData& stateData, + double deriv[], + const SolverMode& sMode) +{ + if (!hasDifferential(sMode)) { + return; + } + const auto loc = offsets.getLocations(stateData, deriv, sMode, this); + const auto evaluation = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc); + std::copy_n(evaluation.rates.begin(), loc.diffSize, loc.destDiffLoc); +} + +void ExciterESST2A::algebraicUpdate(const IOdata& /*inputs*/, + const StateData& stateData, + double update[], + const SolverMode& sMode, + double /*alpha*/) +{ + if (hasAlgebraic(sMode)) { + const auto loc = offsets.getLocations(stateData, update, sMode, this); + loc.destLoc[0] = loc.diffStateLoc[fieldState]; + } +} + +void ExciterESST2A::jacobianElements(const IOdata& inputs, + const StateData& stateData, + MatrixData& matrixData, + const IOlocs& inputLocs, + const SolverMode& sMode) +{ + const auto loc = offsets.getLocations(stateData, sMode, this); + const auto evaluation = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc); + if (hasAlgebraic(sMode)) { + matrixData.assign(loc.algOffset, loc.algOffset, evaluation.algebraicDerivative); + for (index_t column = 0; column < stateCount; ++column) { + matrixData.assign(loc.algOffset, + loc.diffOffset + column, + evaluation.fieldStateDerivatives[column]); + } + } + if (!hasDifferential(sMode)) { + return; + } + for (index_t row = 0; row < stateCount; ++row) { + matrixData.assign(loc.diffOffset + row, + loc.algOffset, + evaluation.rateAlgebraicDerivatives[row]); + for (index_t column = 0; column < stateCount; ++column) { + double value = evaluation.rateStateDerivatives[row][column]; + if (row == column) { + value -= stateData.cj; + } + matrixData.assign(loc.diffOffset + row, loc.diffOffset + column, value); + } + for (index_t input = 0; input < exciterInputCount; ++input) { + matrixData.assignCheckCol(loc.diffOffset + row, + inputLocs[input], + evaluation.rateInputDerivatives[row][input]); + } + } +} + +void ExciterESST2A::timestep(CoreTime time, const IOdata& inputs, const SolverMode& /*sMode*/) +{ + derivative(inputs, emptyStateData, m_dstate_dt.data(), cLocalSolverMode); + const double timeStep = time - prevTime; + for (index_t index = 0; index < stateCount; ++index) { + m_state[index + 1] += timeStep * m_dstate_dt[index + 1]; + } + updateLimitFlags(inputs, m_state.data() + 1, true); + m_state[0] = m_state[fieldState + 1]; + prevTime = time; +} + +bool ExciterESST2A::updateLimitFlags(const IOdata& inputs, double state[], bool projectStates) +{ + constexpr double tolerance = 1e-7; + const auto evaluation = evaluate(inputs, m_state[0], state); + bool changed = false; + const auto setLimited = [this, + &changed](int limitedFlag, int highFlag, bool limited, bool high) { + if (opFlags[limitedFlag] != limited) { + opFlags.set(limitedFlag, limited); + changed = true; + } + if (limited && (opFlags[highFlag] != high)) { + opFlags.set(highFlag, high); + changed = true; + } + }; + + const bool regulatorAbove = state[regulatorState] > Vrmax + tolerance; + const bool regulatorBelow = state[regulatorState] < Vrmin - tolerance; + if (projectStates) { + state[regulatorState] = std::clamp(state[regulatorState], + static_cast(Vrmin), + static_cast(Vrmax)); + } + if (opFlags[REGULATOR_LIMITED]) { + const bool release = opFlags[REGULATOR_LIMIT_HIGH] ? (evaluation.regulatorDrive < 0.0) : + (evaluation.regulatorDrive > 0.0); + if (release) { + setLimited(REGULATOR_LIMITED, REGULATOR_LIMIT_HIGH, false, false); + } + } else if (regulatorAbove || + ((state[regulatorState] >= Vrmax - tolerance) && + (evaluation.regulatorDrive > 0.0))) { + setLimited(REGULATOR_LIMITED, REGULATOR_LIMIT_HIGH, true, true); + } else if (regulatorBelow || + ((state[regulatorState] <= Vrmin + tolerance) && + (evaluation.regulatorDrive < 0.0))) { + setLimited(REGULATOR_LIMITED, REGULATOR_LIMIT_HIGH, true, false); + } + + const bool fieldAbove = state[fieldState] > Efdmax + tolerance; + const bool fieldBelow = state[fieldState] < -tolerance; + if (projectStates) { + state[fieldState] = std::clamp(state[fieldState], 0.0, static_cast(Efdmax)); + } + if (opFlags[FIELD_LIMITED]) { + const bool release = opFlags[FIELD_LIMIT_HIGH] ? (evaluation.fieldDrive < 0.0) : + (evaluation.fieldDrive > 0.0); + if (release) { + setLimited(FIELD_LIMITED, FIELD_LIMIT_HIGH, false, false); + } + } else if (fieldAbove || + ((state[fieldState] >= Efdmax - tolerance) && (evaluation.fieldDrive > 0.0))) { + setLimited(FIELD_LIMITED, FIELD_LIMIT_HIGH, true, true); + } else if (fieldBelow || ((state[fieldState] <= tolerance) && (evaluation.fieldDrive < 0.0))) { + setLimited(FIELD_LIMITED, FIELD_LIMIT_HIGH, true, false); + } + return changed; +} + +void ExciterESST2A::rootTest(const IOdata& inputs, + const StateData& stateData, + double roots[], + const SolverMode& sMode) +{ + const auto loc = offsets.getLocations(stateData, sMode, this); + const auto evaluation = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc); + const index_t root = offsets.getRootOffset(sMode); + if (opFlags[REGULATOR_LIMITED]) { + roots[root] = + opFlags[REGULATOR_LIMIT_HIGH] ? -evaluation.regulatorDrive : evaluation.regulatorDrive; + } else { + roots[root] = std::min(Vrmax - loc.diffStateLoc[regulatorState], + loc.diffStateLoc[regulatorState] - Vrmin); + opFlags.set(REGULATOR_LIMIT_HIGH, loc.diffStateLoc[regulatorState] >= Vrmax); + } + if (opFlags[FIELD_LIMITED]) { + roots[root + 1] = + opFlags[FIELD_LIMIT_HIGH] ? -evaluation.fieldDrive : evaluation.fieldDrive; + } else { + roots[root + 1] = + std::min(Efdmax - loc.diffStateLoc[fieldState], loc.diffStateLoc[fieldState]); + opFlags.set(FIELD_LIMIT_HIGH, loc.diffStateLoc[fieldState] >= Efdmax); + } +} + +void ExciterESST2A::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, true)) { + alert(this, JAC_COUNT_CHANGE); + } +} + +ChangeCode ExciterESST2A::rootCheck(const IOdata& inputs, + const StateData& /*stateData*/, + const SolverMode& /*sMode*/, + CheckLevel /*level*/) +{ + if (updateLimitFlags(inputs, m_state.data() + 1, true)) { + alert(this, JAC_COUNT_CHANGE); + return ChangeCode::JACOBIAN_CHANGE; + } + return ChangeCode::NO_CHANGE; +} + +stringVec ExciterESST2A::localStateNames() const +{ + return {"efd", "vmeas", "va", "wf", "efstate"}; +} + +index_t ExciterESST2A::findIndex(std::string_view field, const SolverMode& sMode) const +{ + if ((field == "efd") || (field == "field")) { + return offsets.getAlgOffset(sMode); + } + const auto offset = offsets.getDiffOffset(sMode); + if ((field == "vmeas") || (field == "vc")) { + return offset + voltageMeasurementState; + } + if ((field == "va") || (field == "regulator")) { + return offset + regulatorState; + } + if ((field == "wf") || (field == "feedback")) { + return offset + feedbackState; + } + if (field == "efstate") { + return offset + fieldState; + } + return kInvalidLocation; +} + +void ExciterESST2A::set(std::string_view param, std::string_view val) +{ + Exciter::set(param, val); +} + +void ExciterESST2A::set(std::string_view param, double val, units::unit unitType) +{ + if (param == "tr") { + Tr = val; + } else if (param == "kp") { + Kp = val; + } else if (param == "ki") { + Ki = val; + } else if (param == "kc") { + Kc = val; + } else if (param == "kf") { + Kf = val; + } else if (param == "tf") { + Tf = val; + } else if (param == "ke") { + Ke = val; + } else if (param == "te") { + Te = val; + } else if ((param == "efdmax") || (param == "efmax")) { + Efdmax = val; + } else { + Exciter::set(param, val, unitType); + } +} + +double ExciterESST2A::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 == "kp") { + return Kp; + } + if (param == "ki") { + return Ki; + } + if (param == "kc") { + return Kc; + } + if (param == "kf") { + return Kf; + } + if (param == "tf") { + return Tf; + } + if (param == "ke") { + return Ke; + } + if (param == "te") { + return Te; + } + if ((param == "efdmax") || (param == "efmax")) { + return Efdmax; + } + return Exciter::get(param, unitType); +} +// NOLINTEND(readability-math-missing-parentheses) +} // namespace griddyn::exciters diff --git a/src/griddyn/exciters/ExciterESST2A.h b/src/griddyn/exciters/ExciterESST2A.h new file mode 100644 index 00000000..e02f734c --- /dev/null +++ b/src/griddyn/exciters/ExciterESST2A.h @@ -0,0 +1,126 @@ +/* + * 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 +#include +#include +#include + +namespace griddyn::exciters { +/** IEEE/PSS/E ESST2A static excitation system. + * + * ESST2A is a compound-source controlled-rectifier exciter. The source term is + * \f$V_B=V_EF_{EX}(K_CX_{ad}I_{fd}/V_E)\f$, where + * \f$V_E=|K_PV_T+jK_II_T|\f$; if both \f$K_P\f$ and \f$K_I\f$ are nonpositive, + * the OpenIPSL/PSS/E bypass uses \f$V_B=1\f$. With measured voltage + * \f$V_m\f$, regulator output \f$V_A\f$, feedback state \f$x_F\f$, and field + * state \f$E_{fd}\f$, + * \f[ + * T_R\dot V_m=V_T-V_m,\quad V_F=K_F(E_{fd}-x_F)/T_F,\quad + * V_i=V_{ref}+V_{bias}+V_{set}-1+V_{SS}-V_m-V_F, + * \f] + * \f[ + * T_A\dot V_A=K_AV_i-V_A,\quad + * T_E\dot E_{fd}=V_BV_A-K_EE_{fd},\quad + * T_F\dot x_F=E_{fd}-x_F. + * \f] + * \f$V_A\f$ is limited by \f$[V_{RMIN},V_{RMAX}]\f$ and \f$E_{fd}\f$ by + * \f$[0,E_{FDMAX}]\f$ with event-aware non-windup behavior. + * + * This is the 13-parameter OpenIPSL core model. GridDyn's standard + * stabilizer input supplies the normal summed supplementary path. Separate + * VUEL/VOEL inputs and the high-value gate are not available in the current + * exciter interface and are therefore outside this implementation. + * + * @par Equation source + * OpenIPSL `Electrical.Controls.PSSE.ES.ESST2A` and its shared rectifier, + * derivative-feedback, limited-lag, and limited-integrator blocks, commit + * 8155c73f. + */ +class ExciterESST2A final: public Exciter { + public: + enum ESST2AFlags { + REGULATOR_LIMITED = OBJECT_FLAG5, + REGULATOR_LIMIT_HIGH = OBJECT_FLAG6, + FIELD_LIMITED = OBJECT_FLAG7, + FIELD_LIMIT_HIGH = OBJECT_FLAG8, + }; + + explicit ExciterESST2A(const std::string& objName = "exciterESST2A_#"); + 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 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; + 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; + + private: + static constexpr index_t stateCount = 4; + struct Evaluation { + double fieldOutput = 0.0; + double algebraicDerivative = 0.0; + std::array fieldStateDerivatives{}; + std::array fieldInputDerivatives{}; + std::array rates{}; + std::array rateAlgebraicDerivatives{}; + std::array, stateCount> rateStateDerivatives{}; + std::array, stateCount> rateInputDerivatives{}; + double regulatorDrive = 0.0; + double fieldDrive = 0.0; + }; + Evaluation evaluate(const IOdata& inputs, double fieldVoltage, const double state[]) const; + double rectifierVoltage(const IOdata& inputs) const; + bool updateLimitFlags(const IOdata& inputs, double state[], bool projectStates); + + model_parameter Tr = 0.01; + model_parameter Kp = 0.7; + model_parameter Ki = 1.0; + model_parameter Kc = 0.03; + model_parameter Kf = 0.05; + model_parameter Tf = 0.7; + model_parameter Ke = 1.0; + model_parameter Te = 0.5; + model_parameter Efdmax = 5.0; +}; +} // namespace griddyn::exciters diff --git a/src/griddyn/exciters/ExciterIEEET3.cpp b/src/griddyn/exciters/ExciterIEEET3.cpp new file mode 100644 index 00000000..c3550139 --- /dev/null +++ b/src/griddyn/exciters/ExciterIEEET3.cpp @@ -0,0 +1,482 @@ +/* + * 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 "ExciterIEEET3.h" + +#include "ExciterSignalHelper.h" +#include "StaticExciterRectifier.h" +#include "core/CoreExceptions.h" +#include "core/CoreObjectTemplates.hpp" +#include "utilities/MatrixData.hpp" +#include +#include +#include +#include +#include +#include +#include + +namespace griddyn::exciters { +// NOLINTBEGIN(readability-math-missing-parentheses) +namespace { + constexpr index_t voltageMeasurementState = 0; + constexpr index_t regulatorState = 1; + constexpr index_t feedbackState = 2; + constexpr index_t fieldState = 3; + constexpr double initializationTolerance = 1e-7; + + double upperRegulatorLimit(double limit) + { + return (limit == 0.0) ? 999.0 : limit; + } + + bool finiteMachineSignal(double value) + { + return std::isfinite(value) && (std::abs(value) < 1e20); + } +} // namespace + +ExciterIEEET3::ExciterIEEET3(const std::string& objName): Exciter(objName) +{ + m_inputSize = exciterInputCount; + Ka = 5.0; + Ta = 0.04; + Vrmax = 7.3; + Vrmin = -7.3; +} + +CoreObject* ExciterIEEET3::clone(CoreObject* obj) const +{ + auto* result = cloneBase(this, obj); + if (result == nullptr) { + return obj; + } + result->Tr = Tr; + result->Vbmax = Vbmax; + result->Ke = Ke; + result->Te = Te; + result->Kf = Kf; + result->Tf = Tf; + result->Kp = Kp; + result->Ki = Ki; + return result; +} + +void ExciterIEEET3::dynObjectInitializeA(CoreTime /*time0*/, std::uint32_t /*flags*/) +{ + const std::array parameters{ + Tr, Ka, Ta, Vrmax, Vrmin, Vbmax, Ke, Te, Kf, Tf, Kp, Ki, Vref}; + if (std::any_of(parameters.begin(), + parameters.end(), + [](double value) { return !std::isfinite(value); }) || + (Tr < 0.0) || (Ka <= 0.0) || (Ta <= 0.0) || (upperRegulatorLimit(Vrmax) < Vrmin) || + (Vbmax <= 0.0) || (Te <= 0.0) || (Tf <= 0.0)) { + throw InvalidParameterValue("IEEET3 gains, time constants, or limits"); + } + offsets.local().local.algSize = 1; + offsets.local().local.diffSize = stateCount; + offsets.local().local.jacSize = 80; + offsets.local().local.algRoots = 1; +} + +double ExciterIEEET3::sourceVoltage(const IOdata& inputs) const +{ + return detail::computeRectifierData( + inputs, Kp, Ki, 0.0, 0.0, 0.0, std::numeric_limits::infinity()) + .voltage; +} + +ExciterIEEET3::Evaluation + ExciterIEEET3::evaluate(const IOdata& inputs, double fieldVoltage, const double state[]) const +{ + using detail::addSignals; + using detail::algebraicSignal; + using detail::clampSignal; + using detail::constantSignal; + using detail::inputSignal; + using detail::multiplySignals; + using detail::scaleSignal; + using detail::sqrtSignal; + using detail::stateSignal; + using detail::subtractSignals; + using Signal = detail::ExciterSignal; + + Evaluation evaluation; + std::array rates{}; + const Signal measuredVoltage = (Tr > 0.0) ? + stateSignal(state[voltageMeasurementState], voltageMeasurementState) : + inputSignal(inputs, exciterVoltageInLocation); + if (Tr > 0.0) { + rates[voltageMeasurementState] = + scaleSignal(subtractSignals(inputSignal(inputs, exciterVoltageInLocation), + measuredVoltage), + 1.0 / Tr); + } + + const auto source = detail::computeRectifierData( + inputs, Kp, Ki, 0.0, 0.0, 0.0, std::numeric_limits::infinity()); + Signal sourceSignal = constantSignal(source.voltage); + const std::array signalIndices{exciterIdInLocation, + exciterIqInLocation, + exciterVdInLocation, + exciterVqInLocation, + exciterXadIfdInLocation}; + for (index_t index = 0; index < 5; ++index) { + sourceSignal.input[signalIndices[index]] = source.derivatives[index]; + } + const Signal fieldCurrent = inputSignal(inputs, exciterXadIfdInLocation); + const Signal v40Input = subtractSignals(multiplySignals(sourceSignal, sourceSignal), + multiplySignals(scaleSignal(fieldCurrent, 0.78), + scaleSignal(fieldCurrent, 0.78))); + const Signal v40 = sqrtSignal(v40Input); + const Signal fieldStateSignal = stateSignal(state[fieldState], fieldState); + const Signal feedbackStateSignal = stateSignal(state[feedbackState], feedbackState); + const Signal feedback = + scaleSignal(subtractSignals(fieldStateSignal, feedbackStateSignal), Kf / Tf); + rates[feedbackState] = + scaleSignal(subtractSignals(fieldStateSignal, feedbackStateSignal), 1.0 / Tf); + + Signal input = addSignals(constantSignal(Vref + vBias - 1.0), + inputSignal(inputs, exciterVsetInLocation)); + input = addSignals(input, inputSignal(inputs, exciterVssInLocation)); + input = subtractSignals(subtractSignals(input, measuredVoltage), feedback); + const Signal regulator = stateSignal(state[regulatorState], regulatorState); + const Signal regulatorDrive = subtractSignals(scaleSignal(input, Ka), regulator); + if (!opFlags[REGULATOR_LIMITED]) { + rates[regulatorState] = scaleSignal(regulatorDrive, 1.0 / Ta); + } + const Signal vbRaw = addSignals(regulator, v40); + const Signal vb = clampSignal(vbRaw, 0.0, Vbmax); + const Signal fieldDrive = subtractSignals(vb, scaleSignal(fieldStateSignal, Ke)); + rates[fieldState] = scaleSignal(fieldDrive, 1.0 / Te); + const Signal algebraicResidual = + subtractSignals(fieldStateSignal, algebraicSignal(fieldVoltage)); + + evaluation.fieldOutput = state[fieldState]; + evaluation.regulatorDrive = regulatorDrive.value; + evaluation.algebraicDerivative = algebraicResidual.algebraic; + evaluation.fieldStateDerivatives = algebraicResidual.state; + evaluation.fieldInputDerivatives = algebraicResidual.input; + for (index_t row = 0; row < stateCount; ++row) { + evaluation.rates[row] = rates[row].value; + evaluation.rateAlgebraicDerivatives[row] = rates[row].algebraic; + evaluation.rateStateDerivatives[row] = rates[row].state; + evaluation.rateInputDerivatives[row] = rates[row].input; + } + return evaluation; +} + +void ExciterIEEET3::dynObjectInitializeB(const IOdata& inputs, + const IOdata& desiredOutput, + IOdata& fieldSet) +{ + if (inputs.size() < exciterInputCount || desiredOutput.empty() || + !std::isfinite(inputs[exciterVoltageInLocation]) || + !std::isfinite(inputs[exciterVsetInLocation]) || + !std::isfinite(inputs[exciterVssInLocation]) || !std::isfinite(desiredOutput[0])) { + throw InvalidParameterValue("IEEET3 initial voltage signals"); + } + for (const index_t inputIndex : {exciterIdInLocation, + exciterIqInLocation, + exciterVdInLocation, + exciterVqInLocation, + exciterXadIfdInLocation}) { + if (!finiteMachineSignal(inputs[inputIndex])) { + throw InvalidParameterValue("IEEET3 requires synchronous-machine dq and field signals"); + } + } + const double fieldVoltage = desiredOutput[0]; + const double source = sourceVoltage(inputs); + const double v40Argument = + (source * source) - std::pow(0.78 * inputs[exciterXadIfdInLocation], 2.0); + const double v40 = (v40Argument > 0.0) ? std::sqrt(v40Argument) : 0.0; + const double regulator = (Ke * fieldVoltage) - v40; + if ((regulator < Vrmin - initializationTolerance) || + (regulator > upperRegulatorLimit(Vrmax) + initializationTolerance)) { + throw InvalidParameterValue("IEEET3 initial regulator output outside limits"); + } + m_state[0] = fieldVoltage; + double* state = m_state.data() + 1; + state[voltageMeasurementState] = inputs[exciterVoltageInLocation]; + state[regulatorState] = regulator; + state[feedbackState] = fieldVoltage; + state[fieldState] = fieldVoltage; + vBias = (regulator / Ka) + state[voltageMeasurementState] - Vref - + (inputs[exciterVsetInLocation] - 1.0) - inputs[exciterVssInLocation]; + fieldSet.resize(std::max(fieldSet.size(), 2U)); + fieldSet[exciterVsetInLocation] = Vref; + std::fill(m_dstate_dt.begin(), m_dstate_dt.end(), 0.0); + updateLimitFlag(inputs, state, true); +} + +void ExciterIEEET3::residual(const IOdata& inputs, + const StateData& stateData, + double resid[], + const SolverMode& sMode) +{ + const auto loc = offsets.getLocations(stateData, resid, sMode, this); + const auto evaluation = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc); + if (hasAlgebraic(sMode)) { + loc.destLoc[0] = evaluation.fieldOutput - loc.algStateLoc[0]; + } + if (hasDifferential(sMode)) { + for (index_t index = 0; index < loc.diffSize; ++index) { + loc.destDiffLoc[index] = evaluation.rates[index] - loc.dstateLoc[index]; + } + } +} + +void ExciterIEEET3::derivative(const IOdata& inputs, + const StateData& stateData, + double deriv[], + const SolverMode& sMode) +{ + if (!hasDifferential(sMode)) { + return; + } + const auto loc = offsets.getLocations(stateData, deriv, sMode, this); + const auto evaluation = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc); + std::copy_n(evaluation.rates.begin(), loc.diffSize, loc.destDiffLoc); +} + +void ExciterIEEET3::algebraicUpdate(const IOdata& /*inputs*/, + const StateData& stateData, + double update[], + const SolverMode& sMode, + double /*alpha*/) +{ + if (hasAlgebraic(sMode)) { + const auto loc = offsets.getLocations(stateData, update, sMode, this); + loc.destLoc[0] = loc.diffStateLoc[fieldState]; + } +} + +void ExciterIEEET3::jacobianElements(const IOdata& inputs, + const StateData& stateData, + MatrixData& matrixData, + const IOlocs& inputLocs, + const SolverMode& sMode) +{ + const auto loc = offsets.getLocations(stateData, sMode, this); + const auto evaluation = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc); + if (hasAlgebraic(sMode)) { + matrixData.assign(loc.algOffset, loc.algOffset, evaluation.algebraicDerivative); + for (index_t column = 0; column < stateCount; ++column) { + matrixData.assign(loc.algOffset, + loc.diffOffset + column, + evaluation.fieldStateDerivatives[column]); + } + } + if (!hasDifferential(sMode)) { + return; + } + for (index_t row = 0; row < stateCount; ++row) { + matrixData.assign(loc.diffOffset + row, + loc.algOffset, + evaluation.rateAlgebraicDerivatives[row]); + for (index_t column = 0; column < stateCount; ++column) { + double value = evaluation.rateStateDerivatives[row][column]; + if (row == column) { + value -= stateData.cj; + } + matrixData.assign(loc.diffOffset + row, loc.diffOffset + column, value); + } + for (index_t input = 0; input < exciterInputCount; ++input) { + matrixData.assignCheckCol(loc.diffOffset + row, + inputLocs[input], + evaluation.rateInputDerivatives[row][input]); + } + } +} + +void ExciterIEEET3::timestep(CoreTime time, const IOdata& inputs, const SolverMode& /*sMode*/) +{ + derivative(inputs, emptyStateData, m_dstate_dt.data(), cLocalSolverMode); + const double timeStep = time - prevTime; + for (index_t index = 0; index < stateCount; ++index) { + m_state[index + 1] += timeStep * m_dstate_dt[index + 1]; + } + updateLimitFlag(inputs, m_state.data() + 1, true); + m_state[0] = m_state[fieldState + 1]; + prevTime = time; +} + +bool ExciterIEEET3::updateLimitFlag(const IOdata& inputs, double state[], bool projectState) +{ + constexpr double tolerance = 1e-7; + const double upper = upperRegulatorLimit(Vrmax); + const auto evaluation = evaluate(inputs, m_state[0], state); + bool limited = opFlags[REGULATOR_LIMITED]; + bool high = opFlags[REGULATOR_LIMIT_HIGH]; + const bool above = state[regulatorState] > upper + tolerance; + const bool below = state[regulatorState] < Vrmin - tolerance; + if (projectState) { + state[regulatorState] = + std::clamp(state[regulatorState], static_cast(Vrmin), upper); + } + if (limited) { + const bool release = + high ? (evaluation.regulatorDrive < 0.0) : (evaluation.regulatorDrive > 0.0); + if (release) { + limited = false; + } + } else if (above || + ((state[regulatorState] >= upper - tolerance) && + (evaluation.regulatorDrive > 0.0))) { + limited = true; + high = true; + } else if (below || + ((state[regulatorState] <= Vrmin + tolerance) && + (evaluation.regulatorDrive < 0.0))) { + limited = true; + high = false; + } + const bool changed = (limited != opFlags[REGULATOR_LIMITED]) || + (limited && (high != opFlags[REGULATOR_LIMIT_HIGH])); + opFlags.set(REGULATOR_LIMITED, limited); + if (limited) { + opFlags.set(REGULATOR_LIMIT_HIGH, high); + } + return changed; +} + +void ExciterIEEET3::rootTest(const IOdata& inputs, + const StateData& stateData, + double roots[], + const SolverMode& sMode) +{ + const auto loc = offsets.getLocations(stateData, sMode, this); + const auto evaluation = evaluate(inputs, loc.algStateLoc[0], loc.diffStateLoc); + const index_t root = offsets.getRootOffset(sMode); + if (opFlags[REGULATOR_LIMITED]) { + roots[root] = + opFlags[REGULATOR_LIMIT_HIGH] ? -evaluation.regulatorDrive : evaluation.regulatorDrive; + } else { + const double upper = upperRegulatorLimit(Vrmax); + roots[root] = std::min(upper - loc.diffStateLoc[regulatorState], + loc.diffStateLoc[regulatorState] - Vrmin); + opFlags.set(REGULATOR_LIMIT_HIGH, loc.diffStateLoc[regulatorState] >= upper); + } +} + +void ExciterIEEET3::rootTrigger(CoreTime /*time*/, + const IOdata& inputs, + const std::vector& rootMask, + const SolverMode& sMode) +{ + if ((rootMask[offsets.getRootOffset(sMode)] != 0) && + updateLimitFlag(inputs, m_state.data() + 1, true)) { + alert(this, JAC_COUNT_CHANGE); + } +} + +ChangeCode ExciterIEEET3::rootCheck(const IOdata& inputs, + const StateData& /*stateData*/, + const SolverMode& /*sMode*/, + CheckLevel /*level*/) +{ + if (updateLimitFlag(inputs, m_state.data() + 1, true)) { + alert(this, JAC_COUNT_CHANGE); + return ChangeCode::JACOBIAN_CHANGE; + } + return ChangeCode::NO_CHANGE; +} + +stringVec ExciterIEEET3::localStateNames() const +{ + return {"efd", "vmeas", "vr", "wf", "efstate"}; +} + +index_t ExciterIEEET3::findIndex(std::string_view field, const SolverMode& sMode) const +{ + if ((field == "efd") || (field == "field")) { + return offsets.getAlgOffset(sMode); + } + const auto offset = offsets.getDiffOffset(sMode); + if ((field == "vmeas") || (field == "vc")) { + return offset + voltageMeasurementState; + } + if ((field == "vr") || (field == "regulator")) { + return offset + regulatorState; + } + if ((field == "wf") || (field == "feedback")) { + return offset + feedbackState; + } + if (field == "efstate") { + return offset + fieldState; + } + return kInvalidLocation; +} + +void ExciterIEEET3::set(std::string_view param, std::string_view val) +{ + Exciter::set(param, val); +} + +void ExciterIEEET3::set(std::string_view param, double val, units::unit unitType) +{ + if (param == "tr") { + Tr = val; + } else if (param == "vbmax") { + Vbmax = val; + } else if (param == "ke") { + Ke = val; + } else if (param == "te") { + Te = val; + } else if (param == "kf") { + Kf = val; + } else if (param == "tf") { + Tf = val; + } else if (param == "kp") { + Kp = val; + } else if (param == "ki") { + Ki = val; + } else { + Exciter::set(param, val, unitType); + } +} + +double ExciterIEEET3::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 == "vbmax") { + return Vbmax; + } + if (param == "ke") { + return Ke; + } + if (param == "te") { + return Te; + } + if (param == "kf") { + return Kf; + } + if (param == "tf") { + return Tf; + } + if (param == "kp") { + return Kp; + } + if (param == "ki") { + return Ki; + } + return Exciter::get(param, unitType); +} +// NOLINTEND(readability-math-missing-parentheses) +} // namespace griddyn::exciters diff --git a/src/griddyn/exciters/ExciterIEEET3.h b/src/griddyn/exciters/ExciterIEEET3.h new file mode 100644 index 00000000..9a31fc18 --- /dev/null +++ b/src/griddyn/exciters/ExciterIEEET3.h @@ -0,0 +1,119 @@ +/* + * 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 +#include +#include +#include + +namespace griddyn::exciters { +/** IEEE/PSS/E Type 3 excitation system. + * + * IEEET3 is a compound-source exciter. GridDyn uses the existing synchronous + * machine dq signal interface to form + * \f$V_E=|K_P(V_d+jV_q)+jK_I(I_d+jI_q)|\f$ and then + * \f$V_{40}=\sqrt{\max(0,V_E^2-(0.78X_{ad}I_{fd})^2)}\f$. With measured + * voltage \f$V_m\f$, regulator state \f$V_R\f$, washout state \f$x_F\f$, and + * exciter output state \f$E_{fd}\f$, + * \f[ + * T_R\dot V_m=V_T-V_m,\quad V_F=K_F(E_{fd}-x_F)/T_F,\quad + * V_i=V_{ref}+V_{bias}+V_{set}-1+V_{SS}-V_m-V_F, + * \f] + * \f[ + * T_A\dot V_R=K_AV_i-V_R,\quad + * V_B=\operatorname{lim}(V_R+V_{40},0,V_{BMAX}),\quad + * T_E\dot E_{fd}=V_B-K_EE_{fd},\quad + * T_F\dot x_F=E_{fd}-x_F. + * \f] + * The regulator uses anti-windup bounds \f$[V_{RMIN},V_{RMAX}]\f$ with + * PSS/E's convention that \f$V_{RMAX}=0\f$ means an effectively open upper + * limit. The regulator state is projected to either limit and held only for + * outward drive. + * + * @par Equation source + * ANDES `IEEET3` equations and PSS/E DYR schema, commit + * eda5163c9ee8d19945a1dd5d1771fec5da608c27. + */ +class ExciterIEEET3 final: public Exciter { + public: + enum IEEET3Flags { + REGULATOR_LIMITED = OBJECT_FLAG5, + REGULATOR_LIMIT_HIGH = OBJECT_FLAG6, + }; + + explicit ExciterIEEET3(const std::string& objName = "exciterIEEET3_#"); + 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 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; + 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; + + private: + static constexpr index_t stateCount = 4; + struct Evaluation { + double fieldOutput = 0.0; + double algebraicDerivative = 0.0; + std::array fieldStateDerivatives{}; + std::array fieldInputDerivatives{}; + std::array rates{}; + std::array rateAlgebraicDerivatives{}; + std::array, stateCount> rateStateDerivatives{}; + std::array, stateCount> rateInputDerivatives{}; + double regulatorDrive = 0.0; + }; + Evaluation evaluate(const IOdata& inputs, double fieldVoltage, const double state[]) const; + double sourceVoltage(const IOdata& inputs) const; + bool updateLimitFlag(const IOdata& inputs, double state[], bool projectState); + + model_parameter Tr = 0.02; + model_parameter Vbmax = 18.0; + model_parameter Ke = 1.0; + model_parameter Te = 1.0; + model_parameter Kf = 0.1; + model_parameter Tf = 1.0; + model_parameter Kp = 4.0; + model_parameter Ki = 0.1; +}; +} // namespace griddyn::exciters diff --git a/src/griddyn/exciters/ExciterSignalHelper.h b/src/griddyn/exciters/ExciterSignalHelper.h new file mode 100644 index 00000000..954608dd --- /dev/null +++ b/src/griddyn/exciters/ExciterSignalHelper.h @@ -0,0 +1,203 @@ +/* + * 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 "StaticExciterRectifier.h" +#include +#include + +namespace griddyn::exciters::detail { +template +struct ExciterSignal { + double value = 0.0; + double algebraic = 0.0; + std::array state{}; + std::array input{}; +}; + +template +ExciterSignal constantSignal(double value) +{ + ExciterSignal result; + result.value = value; + return result; +} + +template +ExciterSignal algebraicSignal(double value) +{ + auto result = constantSignal(value); + result.algebraic = 1.0; + return result; +} + +template +ExciterSignal stateSignal(double value, index_t index) +{ + auto result = constantSignal(value); + result.state[index] = 1.0; + return result; +} + +template +ExciterSignal inputSignal(const IOdata& inputs, index_t index) +{ + auto result = constantSignal(inputs[index]); + result.input[index] = 1.0; + return result; +} + +template +ExciterSignal addSignals(const ExciterSignal& left, + const ExciterSignal& right) +{ + ExciterSignal result; + result.value = left.value + right.value; + result.algebraic = left.algebraic + right.algebraic; + for (index_t index = 0; index < StateCount; ++index) { + result.state[index] = left.state[index] + right.state[index]; + } + for (index_t index = 0; index < exciterInputCount; ++index) { + result.input[index] = left.input[index] + right.input[index]; + } + return result; +} + +template +ExciterSignal scaleSignal(const ExciterSignal& signal, double factor) +{ + ExciterSignal result; + result.value = factor * signal.value; + result.algebraic = factor * signal.algebraic; + for (index_t index = 0; index < StateCount; ++index) { + result.state[index] = factor * signal.state[index]; + } + for (index_t index = 0; index < exciterInputCount; ++index) { + result.input[index] = factor * signal.input[index]; + } + return result; +} + +template +ExciterSignal subtractSignals(const ExciterSignal& left, + const ExciterSignal& right) +{ + return addSignals(left, scaleSignal(right, -1.0)); +} + +template +ExciterSignal multiplySignals(const ExciterSignal& left, + const ExciterSignal& right) +{ + ExciterSignal result; + result.value = left.value * right.value; + result.algebraic = (left.algebraic * right.value) + (left.value * right.algebraic); + for (index_t index = 0; index < StateCount; ++index) { + result.state[index] = (left.state[index] * right.value) + (left.value * right.state[index]); + } + for (index_t index = 0; index < exciterInputCount; ++index) { + result.input[index] = (left.input[index] * right.value) + (left.value * right.input[index]); + } + return result; +} + +template +ExciterSignal divideSignals(const ExciterSignal& numerator, + const ExciterSignal& denominator) +{ + if (std::abs(denominator.value) <= 1e-12) { + return constantSignal(0.0); + } + ExciterSignal result; + result.value = numerator.value / denominator.value; + const double denominatorSquared = denominator.value * denominator.value; + result.algebraic = + ((numerator.algebraic * denominator.value) - (numerator.value * denominator.algebraic)) / + denominatorSquared; + for (index_t index = 0; index < StateCount; ++index) { + result.state[index] = ((numerator.state[index] * denominator.value) - + (numerator.value * denominator.state[index])) / + denominatorSquared; + } + for (index_t index = 0; index < exciterInputCount; ++index) { + result.input[index] = ((numerator.input[index] * denominator.value) - + (numerator.value * denominator.input[index])) / + denominatorSquared; + } + return result; +} + +template +ExciterSignal + applyFunction(const ExciterSignal& signal, double value, double derivative) +{ + ExciterSignal result; + result.value = value; + result.algebraic = derivative * signal.algebraic; + for (index_t index = 0; index < StateCount; ++index) { + result.state[index] = derivative * signal.state[index]; + } + for (index_t index = 0; index < exciterInputCount; ++index) { + result.input[index] = derivative * signal.input[index]; + } + return result; +} + +template +ExciterSignal sqrtSignal(const ExciterSignal& signal) +{ + if (signal.value <= 0.0) { + return constantSignal(0.0); + } + const double root = std::sqrt(signal.value); + return applyFunction(signal, root, 0.5 / root); +} + +template +ExciterSignal + clampSignal(const ExciterSignal& signal, double lower, double upper) +{ + if (signal.value <= lower) { + return constantSignal(lower); + } + if (signal.value >= upper) { + return constantSignal(upper); + } + return signal; +} + +template +ExciterSignal clampSignal(const ExciterSignal& signal, + const ExciterSignal& lower, + const ExciterSignal& upper) +{ + if (signal.value <= lower.value) { + return lower; + } + if (signal.value >= upper.value) { + return upper; + } + return signal; +} + +template +ExciterSignal rectifierFactorSignal(const ExciterSignal& normalizedCurrent) +{ + const auto data = computeRectifierFactor(normalizedCurrent.value); + return applyFunction(normalizedCurrent, data.factor, data.derivative); +} + +template +bool integrationBlocked(double state, + double minimum, + double maximum, + const ExciterSignal& drive) +{ + return ((state >= maximum) && (drive.value > 0.0)) || + ((state <= minimum) && (drive.value < 0.0)); +} +} // namespace griddyn::exciters::detail diff --git a/src/utilities/Saturation.cpp b/src/utilities/Saturation.cpp index 8ba3c112..9de87be7 100644 --- a/src/utilities/Saturation.cpp +++ b/src/utilities/Saturation.cpp @@ -128,15 +128,17 @@ void Saturation::setParam(double firstInput, } const double ssv = sqrt((firstSaturation * firstInput) / (secondSaturation * secondInput)); - const double fitDenominator = firstInput - ssv; + const double fitDenominator = 1.0 - ssv; if (std::abs(fitDenominator) < 1e-12) { A = 0.0; B = 0.0; break; } - A = -(secondInput * ssv - firstInput) / fitDenominator; + A = (firstInput - secondInput * ssv) / fitDenominator; const double distance = firstInput - A; - B = (std::abs(distance) < 1e-12) ? 0.0 : firstSaturation / (distance * distance); + B = (std::abs(distance) < 1e-12) ? + 0.0 : + (firstSaturation * firstInput) / (distance * distance); } break; case SaturationType::EXPONENTIAL: A = log(firstSaturation / secondSaturation) / log(firstInput / secondInput); diff --git a/test/andesTests/testAndesDyrReader.cpp b/test/andesTests/testAndesDyrReader.cpp index 17f7e5d5..ec25680d 100644 --- a/test/andesTests/testAndesDyrReader.cpp +++ b/test/andesTests/testAndesDyrReader.cpp @@ -12,11 +12,14 @@ #include "griddyn/GridDynSimulation.h" #include "griddyn/GridSubModel.h" #include "griddyn/events/Event.h" +#include "griddyn/exciters/ExciterAC7B.h" +#include "griddyn/exciters/ExciterAC8B.h" #include "griddyn/exciters/ExciterDC1A.h" #include "griddyn/exciters/ExciterDC2A.h" #include "griddyn/exciters/ExciterESAC1A.h" #include "griddyn/exciters/ExciterESAC6A.h" #include "griddyn/exciters/ExciterESST1A.h" +#include "griddyn/exciters/ExciterESST2A.h" #include "griddyn/exciters/ExciterESST3A.h" #include "griddyn/exciters/ExciterESST4B.h" #include "griddyn/exciters/ExciterEXAC1.h" @@ -24,6 +27,7 @@ #include "griddyn/exciters/ExciterEXAC4.h" #include "griddyn/exciters/ExciterEXPIC1.h" #include "griddyn/exciters/ExciterEXST1.h" +#include "griddyn/exciters/ExciterIEEET3.h" #include "griddyn/exciters/ExciterIEEEX1.h" #include "griddyn/exciters/ExciterIEEEtype1.h" #include "griddyn/exciters/ExciterSCRX.h" @@ -420,6 +424,116 @@ TEST(AndesDyrReaderTests, MapsGensaeAndEsst1aParametersInPsseDyrOrder) EXPECT_EQ(runJacobianCheck(simulation, griddyn::cDaeSolverMode, false), 0); } +TEST(AndesDyrReaderTests, MapsNewExciterFocusedChunkInDyrOrder) +{ + { + 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_ieeet3.dyr")); + auto* bus = dynamic_cast(simulation->findByUserID("bus", 1)); + ASSERT_NE(bus, nullptr); + auto* exciter = + dynamic_cast(bus->getGen(0)->find("exciter")); + ASSERT_NE(exciter, nullptr); + const std::pair expected[]{{"tr", 0.011}, + {"ka", 5.1}, + {"ta", 0.041}, + {"vrmax", 20.1}, + {"vrmin", -20.2}, + {"vbmax", 18.3}, + {"ke", 1.04}, + {"te", 1.05}, + {"kf", 0.106}, + {"tf", 1.07}, + {"kp", 4.08}, + {"ki", 0.109}}; + for (const auto& [name, value] : expected) { + EXPECT_DOUBLE_EQ(exciter->get(name), value) << name; + } + ASSERT_EQ(simulation->dynInitialize(), 0); + EXPECT_EQ(runResidualCheck(simulation, griddyn::cDaeSolverMode, false), 0); + EXPECT_EQ(runJacobianCheck(simulation, griddyn::cDaeSolverMode, false), 0); + } + { + 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_ac8b.dyr")); + auto* bus = dynamic_cast(simulation->findByUserID("bus", 1)); + ASSERT_NE(bus, nullptr); + auto* exciter = + dynamic_cast(bus->getGen(0)->find("exciter")); + ASSERT_NE(exciter, nullptr); + const std::pair expected[]{ + {"tr", 0.011}, {"kpr", 10.1}, {"kir", 11.2}, {"kdr", 12.3}, {"tdr", 0.24}, + {"vpmax", 99.5}, {"vpmin", -99.6}, {"vrmax", 99.7}, {"vrmin", -99.8}, {"vfemax", 99.9}, + {"vemin", -9.1}, {"ta", 0.42}, {"ka", 40.3}, {"te", 0.84}, {"kc", 0.015}, + {"kd", 0.016}, {"ke", 1.07}, {"e1", 3.8}, {"se1", 0.18}, {"e2", 2.6}, + {"se2", 0.05}}; + for (const auto& [name, value] : expected) { + EXPECT_DOUBLE_EQ(exciter->get(name), value) << name; + } + ASSERT_EQ(simulation->dynInitialize(), 0); + EXPECT_EQ(runResidualCheck(simulation, griddyn::cDaeSolverMode, false), 0); + EXPECT_EQ(runJacobianCheck(simulation, griddyn::cDaeSolverMode, false), 0); + } + { + 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_ac7b.dyr")); + auto* bus = dynamic_cast(simulation->findByUserID("bus", 1)); + ASSERT_NE(bus, nullptr); + auto* exciter = + dynamic_cast(bus->getGen(0)->find("exciter")); + ASSERT_NE(exciter, nullptr); + const std::pair expected[]{ + {"tr", 0.011}, {"kpr", 4.1}, {"kir", 4.2}, {"kdr", 0.13}, {"tdr", 0.14}, + {"vrmax", 20.1}, {"vrmin", -20.2}, {"kpa", 2.1}, {"kia", 1.1}, {"vamax", 20.3}, + {"vamin", -20.4}, {"kp", 1.2}, {"kl", 0.3}, {"te", 1.3}, {"kc", 0.01}, + {"kd", 0.02}, {"ke", 1.01}, {"kf1", 0.2}, {"kf2", 0.03}, {"kf3", 0.04}, + {"tf3", 0.5}, {"vemin", -5.1}, {"vfemax", 20.5}, {"e1", 3.1}, {"se1", 0.1}, + {"e2", 2.1}, {"se2", 0.02}}; + for (const auto& [name, value] : expected) { + EXPECT_DOUBLE_EQ(exciter->get(name), value) << name; + } + ASSERT_EQ(simulation->dynInitialize(), 0); + EXPECT_EQ(runResidualCheck(simulation, griddyn::cDaeSolverMode, false), 0); + EXPECT_EQ(runJacobianCheck(simulation, griddyn::cDaeSolverMode, false), 0); + } + { + 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_esst2a.dyr")); + auto* bus = dynamic_cast(simulation->findByUserID("bus", 1)); + ASSERT_NE(bus, nullptr); + auto* exciter = + dynamic_cast(bus->getGen(0)->find("exciter")); + ASSERT_NE(exciter, nullptr); + const std::pair expected[]{{"tr", 0.011}, + {"ka", 40.2}, + {"ta", 0.053}, + {"vrmax", 99.4}, + {"vrmin", -99.5}, + {"kp", 0.71}, + {"ki", 1.02}, + {"kc", 0.033}, + {"kf", 0.054}, + {"tf", 0.75}, + {"ke", 1.06}, + {"te", 0.57}, + {"efdmax", 9.8}}; + for (const auto& [name, value] : expected) { + EXPECT_DOUBLE_EQ(exciter->get(name), value) << name; + } + ASSERT_EQ(simulation->dynInitialize(), 0); + EXPECT_EQ(runResidualCheck(simulation, griddyn::cDaeSolverMode, false), 0); + EXPECT_EQ(runJacobianCheck(simulation, griddyn::cDaeSolverMode, false), 0); + } +} + TEST(AndesDyrReaderTests, MapsExpic1ParametersInPsseDyrOrder) { auto simulation = std::make_unique(); diff --git a/test/componentTests/testExciters.cpp b/test/componentTests/testExciters.cpp index 33f6c694..010e6409 100644 --- a/test/componentTests/testExciters.cpp +++ b/test/componentTests/testExciters.cpp @@ -11,16 +11,20 @@ #include "fileInput/fileInput.h" #include "gmlc/utilities/vectorOps.hpp" #include "griddyn/Generator.h" +#include "griddyn/exciters/ExciterAC7B.h" +#include "griddyn/exciters/ExciterAC8B.h" #include "griddyn/exciters/ExciterDC1A.h" #include "griddyn/exciters/ExciterDC2A.h" #include "griddyn/exciters/ExciterESAC6A.h" #include "griddyn/exciters/ExciterESST1A.h" +#include "griddyn/exciters/ExciterESST2A.h" #include "griddyn/exciters/ExciterESST3A.h" #include "griddyn/exciters/ExciterESST4B.h" #include "griddyn/exciters/ExciterEXAC1.h" #include "griddyn/exciters/ExciterEXAC4.h" #include "griddyn/exciters/ExciterEXPIC1.h" #include "griddyn/exciters/ExciterEXST1.h" +#include "griddyn/exciters/ExciterIEEET3.h" #include "griddyn/exciters/ExciterIEEEX1.h" #include "griddyn/exciters/ExciterIEEEtype1.h" #include "griddyn/exciters/ExciterSCRX.h" @@ -29,6 +33,7 @@ #include "griddyn/genmodels/GenModelGENSAL.h" #include "solvers/SolverMode.hpp" #include "utilities/MatrixDataSparse.hpp" +#include #include #include #include @@ -107,7 +112,8 @@ void verifyStabilityCase(ExciterTests& fixture, (exciterName == "esst3a") || (exciterName == "esst4b") || (exciterName == "exst1") || (exciterName == "exac1") || (exciterName == "esac1a") || (exciterName == "exac2") || (exciterName == "exac4") || (exciterName == "expic1") || (exciterName == "scrx") || - (exciterName == "esac6a")) { + (exciterName == "esac6a") || (exciterName == "ac7b") || (exciterName == "ac8b") || + (exciterName == "ieeet3") || (exciterName == "esst2a")) { continue; } if (std::find(skippedExcters.begin(), skippedExcters.end(), exciterName) != @@ -173,6 +179,62 @@ void verifyDefaultPsseSaturation(Exciter& exciter, } } +void expectExciterJacobian(Exciter& exciter, + const IOdata& inputs, + const std::vector& state, + double tolerance = 2e-5) +{ + constexpr double step = 1e-6; + ASSERT_EQ(exciter.stateSize(cDaeSolverMode), state.size()); + exciter.setOffset(0, cDaeSolverMode); + std::vector stateDerivative(state.size(), 0.0); + StateData stateData(0.0, state.data(), stateDerivative.data()); + stateData.stateSize = static_cast(state.size()); + stateData.cj = 0.0; + MatrixDataSparse jacobian; + IOlocs inputLocs(exciterInputCount, kNullLocation); + constexpr index_t inputColumnBase = 20; + for (index_t index = 0; index < exciterInputCount; ++index) { + inputLocs[index] = inputColumnBase + index; + } + exciter.jacobianElements(inputs, stateData, jacobian, inputLocs, cDaeSolverMode); + + const auto residualAt = [&exciter, &stateDerivative](const IOdata& trialInputs, + const std::vector& trialState) { + StateData trialData(0.0, trialState.data(), stateDerivative.data()); + trialData.stateSize = static_cast(trialState.size()); + std::vector residual(trialState.size(), 0.0); + exciter.residual(trialInputs, trialData, residual.data(), cDaeSolverMode); + return residual; + }; + for (index_t column = 0; column < state.size(); ++column) { + auto plus = state; + auto minus = state; + plus[column] += step; + minus[column] -= step; + const auto plusResidual = residualAt(inputs, plus); + const auto minusResidual = residualAt(inputs, minus); + for (index_t row = 0; row < state.size(); ++row) { + const double numerical = (plusResidual[row] - minusResidual[row]) / (2.0 * step); + EXPECT_NEAR(jacobian.at(row, column), numerical, tolerance) + << "state row " << row << " column " << column; + } + } + for (index_t column = 0; column < exciterInputCount; ++column) { + auto plus = inputs; + auto minus = inputs; + plus[column] += step; + minus[column] -= step; + const auto plusResidual = residualAt(plus, state); + const auto minusResidual = residualAt(minus, state); + for (index_t row = 0; row < state.size(); ++row) { + const double numerical = (plusResidual[row] - minusResidual[row]) / (2.0 * step); + EXPECT_NEAR(jacobian.at(row, inputColumnBase + column), numerical, tolerance) + << "input row " << row << " column " << column; + } + } +} + } // namespace TEST_F(ExciterTests, RootExciterTest) @@ -1637,6 +1699,337 @@ TEST_F(ExciterTests, Exst1CouplesToSynchronousGeneratorFamiliesAndHasAnalyticJac } } +TEST(ExciterModelTests, Ac8bMatchesPidAndRotatingExciterEquations) +{ + auto factory = CoreObjectFactory::instance(); + std::unique_ptr object(factory->createObject("exciter", "ac8b")); + auto* exciter = dynamic_cast(object.get()); + ASSERT_NE(exciter, nullptr); + exciter->set("tr", 0.1); + exciter->set("kpr", 2.0); + exciter->set("kir", 0.5); + exciter->set("kdr", 0.1); + exciter->set("tdr", 0.2); + exciter->set("vpmax", 10.0); + exciter->set("vpmin", -10.0); + exciter->set("vrmax", 10.0); + exciter->set("vrmin", -10.0); + exciter->set("vfemax", 10.0); + exciter->set("vemin", -1.0); + exciter->set("ka", 3.0); + exciter->set("ta", 0.4); + exciter->set("te", 0.5); + exciter->set("kc", 0.0); + exciter->set("kd", 0.0); + exciter->set("ke", 1.0); + exciter->set("e1", 0.0); + exciter->set("se1", 0.0); + exciter->set("e2", 1.0); + exciter->set("se2", 0.0); + exciter->dynInitializeA(0.0, 0); + EXPECT_EQ(exciter->localStateNames(), + (stringVec{"efd", "vmeas", "pidint", "pidder", "vr", "ve"})); + + IOdata inputs(exciterInputCount, 0.0); + inputs[exciterVoltageInLocation] = 1.0; + inputs[exciterVsetInLocation] = 1.02; + inputs[exciterVssInLocation] = 0.01; + IOdata initializationInputs = inputs; + initializationInputs[exciterVoltageInLocation] = 1.03; + IOdata fieldSet(4, 0.0); + exciter->dynInitializeB(initializationInputs, {0.18}, fieldSet); + std::vector state{0.75, 0.95, 0.1, 0.8, 0.7, 0.8}; + std::vector stateDerivative(state.size(), 0.0); + exciter->setState(0.0, state.data(), stateDerivative.data(), cLocalSolverMode); + std::vector derivative(state.size(), 0.0); + exciter->derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode); + EXPECT_NEAR(derivative[1], 0.5, 1e-12); + EXPECT_NEAR(derivative[2], 0.04, 1e-12); + EXPECT_NEAR(derivative[3], -3.6, 1e-12); + EXPECT_NEAR(derivative[4], -2.5, 1e-12); + EXPECT_NEAR(derivative[5], -0.2, 1e-12); + + std::vector residual(state.size(), 0.0); + exciter->residual(inputs, emptyStateData, residual.data(), cLocalSolverMode); + EXPECT_NEAR(residual[0], 0.05, 1e-12); +} + +TEST(ExciterModelTests, Ac7bMatchesNestedRegulatorEquations) +{ + exciters::ExciterAC7B exciter; + exciter.set("tr", 0.1); + exciter.set("kpr", 2.0); + exciter.set("kir", 0.5); + exciter.set("kdr", 0.0); + exciter.set("tdr", 0.0); + exciter.set("vrmax", 10.0); + exciter.set("vrmin", -10.0); + exciter.set("kpa", 3.0); + exciter.set("kia", 0.5); + exciter.set("vamax", 10.0); + exciter.set("vamin", -10.0); + exciter.set("kp", 2.0); + exciter.set("kl", 10.0); + exciter.set("te", 0.5); + exciter.set("kc", 0.0); + exciter.set("kd", 0.0); + exciter.set("ke", 1.0); + exciter.set("kf1", 0.2); + exciter.set("kf2", 0.3); + exciter.set("kf3", 0.0); + exciter.set("tf3", 0.0); + exciter.set("vfemax", 1.0); + exciter.set("vemin", -10.0); + exciter.set("e1", 0.0); + exciter.set("se1", 0.0); + exciter.set("e2", 1.0); + exciter.set("se2", 0.0); + exciter.dynInitializeA(0.0, 0); + + IOdata inputs(exciterInputCount, 0.0); + inputs[exciterVoltageInLocation] = 1.0; + inputs[exciterVsetInLocation] = 1.02; + inputs[exciterVssInLocation] = 0.01; + IOdata initializationInputs = inputs; + initializationInputs[exciterVoltageInLocation] = 1.03; + IOdata fieldSet(4, 0.0); + exciter.dynInitializeB(initializationInputs, {0.8}, fieldSet); + std::vector state{0.8, 0.95, 0.1, 0.08, 0.2, 0.8, 0.8}; + std::vector stateDerivative(state.size(), 0.0); + exciter.setState(0.0, state.data(), stateDerivative.data(), cLocalSolverMode); + std::vector derivative(state.size(), 0.0); + exciter.derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode); + EXPECT_NEAR(derivative[1], 0.5, 1e-12); + EXPECT_NEAR(derivative[2], 0.04, 1e-12); + EXPECT_NEAR(derivative[3], 0.0, 1e-12); + EXPECT_NEAR(derivative[4], -0.07, 1e-12); + EXPECT_NEAR(derivative[5], 0.0, 1e-12); + EXPECT_NEAR(derivative[6], -2.48, 1e-12); + + // VEMIN/VFEMAX constrain the VE integrator state, not I_C. A large + // inner-PI output therefore remains present in the differential drive. + state[4] = 2.0; + exciter.setState(0.0, state.data(), stateDerivative.data(), cLocalSolverMode); + exciter.derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode); + EXPECT_NEAR(derivative[6], 4.72, 1e-12); +} + +TEST(ExciterModelTests, Ac8bUsesScaledSaturationAndFullPidAntiWindup) +{ + exciters::ExciterAC8B exciter; + exciter.set("tr", 0.1); + exciter.set("kpr", 2.0); + exciter.set("kir", 0.5); + exciter.set("kdr", 0.0); + exciter.set("tdr", 0.0); + exciter.set("vpmax", 0.2); + exciter.set("vpmin", -0.2); + exciter.set("vrmax", 10.0); + exciter.set("vrmin", -10.0); + exciter.set("vfemax", 10.0); + exciter.set("vemin", -10.0); + exciter.set("ka", 1.0); + exciter.set("ta", 0.4); + exciter.set("te", 0.5); + exciter.set("kc", 0.0); + exciter.set("kd", 0.0); + exciter.set("ke", 1.0); + exciter.set("e1", 2.0); + exciter.set("se1", 0.1); + exciter.set("e2", 3.0); + exciter.set("se2", 0.3); + exciter.dynInitializeA(0.0, 0); + + IOdata inputs(exciterInputCount, 0.0); + inputs[exciterVoltageInLocation] = 1.0; + inputs[exciterVsetInLocation] = 1.0; + IOdata fieldSet(4, 0.0); + exciter.dynInitializeB(inputs, {0.1}, fieldSet); + + // The nonlinear equation uses S_E(V_E), not E*S_E(V_E), inside the gain. + std::vector state{0.1, 1.0, 0.0, 0.0, 0.7, 2.5}; + std::vector stateDerivative(state.size(), 0.0); + exciter.setState(0.0, state.data(), stateDerivative.data(), cLocalSolverMode); + std::vector derivative(state.size(), 0.0); + exciter.derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode); + const double ratio = std::sqrt((0.1 * 2.0) / (0.3 * 3.0)); + const double cutoff = (2.0 - (ratio * 3.0)) / (1.0 - ratio); + const double coefficient = 0.1 * 2.0 / std::pow(2.0 - cutoff, 2.0); + const double saturation = coefficient * std::pow(2.5 - cutoff, 2.0) / 2.5; + EXPECT_NEAR(derivative[5], (0.7 - ((1.0 + saturation) * 2.5)) / 0.5, 1e-12); + + // The integral state is inside +/-0.2, but the complete PID output is + // saturated. The outward integral drive must be held. + inputs[exciterVsetInLocation] = 2.0; + state = {0.1, 1.0, 0.0, 0.0, 0.0, 0.1}; + exciter.setState(0.0, state.data(), stateDerivative.data(), cLocalSolverMode); + EXPECT_EQ( + exciter.rootCheck(inputs, emptyStateData, cLocalSolverMode, CheckLevel::REVERSABLE_ONLY), + ChangeCode::JACOBIAN_CHANGE); + exciter.derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode); + EXPECT_DOUBLE_EQ(derivative[2], 0.0); + inputs[exciterVsetInLocation] = 0.0; + EXPECT_EQ( + exciter.rootCheck(inputs, emptyStateData, cLocalSolverMode, CheckLevel::REVERSABLE_ONLY), + ChangeCode::JACOBIAN_CHANGE); + exciter.derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode); + EXPECT_LT(derivative[2], 0.0); +} + +TEST(ExciterModelTests, Ieeet3AndEsst2aMatchStaticSourceEquations) +{ + IOdata inputs(exciterInputCount, 0.0); + inputs[exciterVoltageInLocation] = 1.0; + inputs[exciterVsetInLocation] = 1.02; + inputs[exciterVssInLocation] = 0.01; + IOdata fieldSet(4, 0.0); + + exciters::ExciterIEEET3 ieeet3; + ieeet3.set("tr", 0.1); + ieeet3.set("ka", 5.0); + ieeet3.set("ta", 0.2); + ieeet3.set("vrmax", 10.0); + ieeet3.set("vrmin", -10.0); + ieeet3.set("vbmax", 10.0); + ieeet3.set("ke", 1.0); + ieeet3.set("te", 0.5); + ieeet3.set("kf", 0.1); + ieeet3.set("tf", 0.5); + ieeet3.set("kp", 0.0); + ieeet3.set("ki", 0.0); + ieeet3.dynInitializeA(0.0, 0); + ieeet3.dynInitializeB(inputs, {0.15}, fieldSet); + std::vector state{0.8, 0.95, 0.1, 0.7, 0.8}; + std::vector stateDerivative(state.size(), 0.0); + ieeet3.setState(0.0, state.data(), stateDerivative.data(), cLocalSolverMode); + std::vector derivative(state.size(), 0.0); + ieeet3.derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode); + EXPECT_NEAR(derivative[1], 0.5, 1e-12); + EXPECT_NEAR(derivative[2], 1.0, 1e-12); + EXPECT_NEAR(derivative[3], 0.2, 1e-12); + EXPECT_NEAR(derivative[4], -1.4, 1e-12); + + exciters::ExciterESST2A esst2a; + esst2a.set("tr", 0.1); + esst2a.set("ka", 5.0); + esst2a.set("ta", 0.2); + esst2a.set("vrmax", 10.0); + esst2a.set("vrmin", -10.0); + esst2a.set("kp", 0.0); + esst2a.set("ki", 0.0); + esst2a.set("kf", 0.1); + esst2a.set("tf", 0.5); + esst2a.set("ke", 1.0); + esst2a.set("te", 0.5); + esst2a.set("efdmax", 10.0); + esst2a.dynInitializeA(0.0, 0); + esst2a.dynInitializeB(inputs, {0.15}, fieldSet); + state = {0.8, 0.95, 0.8, 0.7, 0.8}; + esst2a.setState(0.0, state.data(), stateDerivative.data(), cLocalSolverMode); + derivative.assign(state.size(), 0.0); + esst2a.derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode); + EXPECT_NEAR(derivative[1], 0.5, 1e-12); + EXPECT_NEAR(derivative[2], -2.5, 1e-12); + EXPECT_NEAR(derivative[3], 0.2, 1e-12); + EXPECT_NEAR(derivative[4], 0.0, 1e-12); +} + +TEST(ExciterModelTests, Ieeet3AndEsst2aProjectAndReleaseDynamicLimits) +{ + IOdata inputs(exciterInputCount, 0.0); + inputs[exciterVoltageInLocation] = 1.0; + inputs[exciterVsetInLocation] = 1.0; + IOdata fieldSet(4, 0.0); + std::vector stateDerivative(5, 0.0); + std::vector derivative(5, 0.0); + + exciters::ExciterIEEET3 ieeet3; + ieeet3.set("vrmax", 1.0); + ieeet3.set("vrmin", -1.0); + ieeet3.set("kp", 0.0); + ieeet3.set("ki", 0.0); + ieeet3.dynInitializeA(0.0, 0); + ieeet3.dynInitializeB(inputs, {0.1}, fieldSet); + inputs[exciterVsetInLocation] = 2.0; + std::vector state{0.1, 1.0, 1.1, 0.1, 0.1}; + ieeet3.setState(0.0, state.data(), stateDerivative.data(), cLocalSolverMode); + EXPECT_EQ( + ieeet3.rootCheck(inputs, emptyStateData, cLocalSolverMode, CheckLevel::REVERSABLE_ONLY), + ChangeCode::JACOBIAN_CHANGE); + EXPECT_DOUBLE_EQ(ieeet3.getStates()[2], 1.0); + ieeet3.derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode); + EXPECT_DOUBLE_EQ(derivative[2], 0.0); + inputs[exciterVsetInLocation] = 0.0; + EXPECT_EQ( + ieeet3.rootCheck(inputs, emptyStateData, cLocalSolverMode, CheckLevel::REVERSABLE_ONLY), + ChangeCode::JACOBIAN_CHANGE); + ieeet3.derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode); + EXPECT_LT(derivative[2], 0.0); + + inputs[exciterVsetInLocation] = 1.0; + exciters::ExciterESST2A esst2a; + esst2a.set("vrmax", 10.0); + esst2a.set("vrmin", -10.0); + esst2a.set("kp", 0.0); + esst2a.set("ki", 0.0); + esst2a.set("efdmax", 1.0); + esst2a.dynInitializeA(0.0, 0); + esst2a.dynInitializeB(inputs, {0.1}, fieldSet); + state = {1.1, 1.0, 2.0, 1.1, 1.1}; + esst2a.setState(0.0, state.data(), stateDerivative.data(), cLocalSolverMode); + EXPECT_EQ( + esst2a.rootCheck(inputs, emptyStateData, cLocalSolverMode, CheckLevel::REVERSABLE_ONLY), + ChangeCode::JACOBIAN_CHANGE); + EXPECT_DOUBLE_EQ(esst2a.getStates()[4], 1.0); + esst2a.derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode); + EXPECT_DOUBLE_EQ(derivative[4], 0.0); + state = {1.0, 1.0, 0.0, 1.0, 1.0}; + esst2a.setState(0.0, state.data(), stateDerivative.data(), cLocalSolverMode); + EXPECT_EQ( + esst2a.rootCheck(inputs, emptyStateData, cLocalSolverMode, CheckLevel::REVERSABLE_ONLY), + ChangeCode::JACOBIAN_CHANGE); + esst2a.derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode); + EXPECT_LT(derivative[4], 0.0); +} + +TEST(ExciterModelTests, NewExciterJacobiansMatchOffEquilibriumFiniteDifferences) +{ + IOdata inputs(exciterInputCount, 0.0); + inputs[exciterVoltageInLocation] = 1.01; + inputs[exciterVsetInLocation] = 1.0; + inputs[exciterVssInLocation] = 0.02; + inputs[exciterXadIfdInLocation] = 0.4; + inputs[exciterIdInLocation] = 0.2; + inputs[exciterIqInLocation] = 0.1; + inputs[exciterVdInLocation] = 0.1; + inputs[exciterVqInLocation] = 1.0; + IOdata fieldSet(4, 0.0); + + exciters::ExciterAC7B ac7b; + ac7b.dynInitializeA(0.0, 0); + ac7b.dynInitializeB(inputs, {1.0}, fieldSet); + expectExciterJacobian(ac7b, inputs, {5.0, 0.98, 0.2, 0.03, 0.4, 5.0, 5.2}); + + exciters::ExciterAC8B ac8b; + ac8b.set("e1", 3.0); + ac8b.set("se1", 0.2); + ac8b.set("e2", 2.0); + ac8b.set("se2", 0.05); + ac8b.dynInitializeA(0.0, 0); + ac8b.dynInitializeB(inputs, {1.0}, fieldSet); + expectExciterJacobian(ac8b, inputs, {2.4, 0.98, 0.2, 0.03, 2.0, 2.5}); + + exciters::ExciterIEEET3 ieeet3; + ieeet3.dynInitializeA(0.0, 0); + ieeet3.dynInitializeB(inputs, {1.0}, fieldSet); + expectExciterJacobian(ieeet3, inputs, {0.9, 0.98, 0.4, 0.8, 1.1}); + + exciters::ExciterESST2A esst2a; + esst2a.dynInitializeA(0.0, 0); + esst2a.dynInitializeB(inputs, {1.0}, fieldSet); + expectExciterJacobian(esst2a, inputs, {0.9, 0.98, 0.4, 0.8, 1.1}); +} + TEST_F(ExciterTests, BasicStabilityTest1) { static const exciter_parameter_map parameters{ @@ -1721,7 +2114,8 @@ TEST_F(ExciterTests, ExciterTest2AlgDiffTests) if (excname.starts_with("fmi") || (excname == "esst3a") || (excname == "exst1") || (excname == "esst4b") || (excname == "exac1") || (excname == "exac2") || (excname == "esac1a") || (excname == "exac4") || (excname == "scrx") || - (excname == "esac6a")) { + (excname == "esac6a") || (excname == "ac7b") || (excname == "ac8b") || + (excname == "ieeet3") || (excname == "esst2a")) { continue; } gds = readSimXMLFile(fileName); @@ -1771,7 +2165,8 @@ TEST_F(ExciterTests, ExciterAlgDiffJacobianTests) if (excname.starts_with("fmi") || (excname == "esst3a") || (excname == "exst1") || (excname == "esst4b") || (excname == "exac1") || (excname == "exac2") || (excname == "esac1a") || (excname == "exac4") || (excname == "scrx") || - (excname == "esac6a")) { + (excname == "esac6a") || (excname == "ac7b") || (excname == "ac8b") || + (excname == "ieeet3") || (excname == "esst2a")) { continue; } gds = readSimXMLFile(fileName); diff --git a/test/libraryTests/testSaturation.cpp b/test/libraryTests/testSaturation.cpp index b97ad8ef..1247b497 100644 --- a/test/libraryTests/testSaturation.cpp +++ b/test/libraryTests/testSaturation.cpp @@ -51,6 +51,19 @@ TEST(SaturationTests, CutoffQuadraticMatchesAndesExciterSaturation) EXPECT_NEAR(saturation.compute(2.0), 0.4, 1e-14); } +TEST(SaturationTests, CutoffScaledQuadraticMatchesArbitraryReferencePoints) +{ + Saturation saturation(Saturation::SaturationType::CUTOFF_SCALED_QUADRATIC); + saturation.setParam(6.3, 0.44, 4.725, 0.075); + + EXPECT_NEAR(saturation.compute(6.3), 0.44, 1e-13); + EXPECT_NEAR(saturation.compute(4.725), 0.075, 1e-13); + const double step = 1e-6; + const double finiteDifference = + (saturation.compute(5.4 + step) - saturation.compute(5.4 - step)) / (2.0 * step); + EXPECT_NEAR(saturation.deriv(5.4), finiteDifference, 1e-9); +} + TEST(SaturationTests, DisabledCharacteristicIsFinite) { Saturation saturation(Saturation::SaturationType::CUTOFF_SCALED_QUADRATIC); diff --git a/test/test_files/andes_tests/ieee14_ac7b.dyr b/test/test_files/andes_tests/ieee14_ac7b.dyr new file mode 100644 index 00000000..da9c3342 --- /dev/null +++ b/test/test_files/andes_tests/ieee14_ac7b.dyr @@ -0,0 +1 @@ +1 'AC7B' '1' 0.011 4.1 4.2 0.13 0.14 20.1 -20.2 2.1 1.1 20.3 -20.4 1.2 0.3 1.3 0.01 0.02 1.01 0.2 0.03 0.04 0.5 -5.1 20.5 3.1 0.1 2.1 0.02 / diff --git a/test/test_files/andes_tests/ieee14_ac8b.dyr b/test/test_files/andes_tests/ieee14_ac8b.dyr new file mode 100644 index 00000000..92519912 --- /dev/null +++ b/test/test_files/andes_tests/ieee14_ac8b.dyr @@ -0,0 +1 @@ +1 'AC8B' '1' 0.011 10.1 11.2 12.3 0.24 99.5 -99.6 99.7 -99.8 99.9 -9.1 0.42 40.3 0.84 0.015 0.016 1.07 3.8 0.18 2.6 0.05 / diff --git a/test/test_files/andes_tests/ieee14_esst2a.dyr b/test/test_files/andes_tests/ieee14_esst2a.dyr new file mode 100644 index 00000000..a15693e9 --- /dev/null +++ b/test/test_files/andes_tests/ieee14_esst2a.dyr @@ -0,0 +1 @@ +1 'ESST2A' '1' 0.011 40.2 0.053 99.4 -99.5 0.71 1.02 0.033 0.054 0.75 1.06 0.57 9.8 / diff --git a/test/test_files/andes_tests/ieee14_ieeet3.dyr b/test/test_files/andes_tests/ieee14_ieeet3.dyr new file mode 100644 index 00000000..5946b563 --- /dev/null +++ b/test/test_files/andes_tests/ieee14_ieeet3.dyr @@ -0,0 +1 @@ +1 'IEEET3' '1' 0.011 5.1 0.041 20.1 -20.2 18.3 1.04 1.05 0.106 1.07 4.08 0.109 / From 76dd7aee6d6a201bd26dad96aa90ea60121e8998 Mon Sep 17 00:00:00 2001 From: Philip Top Date: Mon, 7 Sep 2026 06:07:26 -0700 Subject: [PATCH 2/3] clang-tidy fixes --- .clang-tidy | 2 ++ src/griddyn/exciters/ExciterAC7B.cpp | 11 ++++++++--- src/griddyn/exciters/ExciterAC8B.cpp | 11 ++++++++--- src/griddyn/exciters/ExciterIEEET3.cpp | 5 +++-- test/componentTests/testExciters.cpp | 7 ++++--- 5 files changed, 25 insertions(+), 11 deletions(-) diff --git a/.clang-tidy b/.clang-tidy index b1b4154f..54715881 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -12,6 +12,7 @@ modernize*, -modernize-use-trailing-return-type, -modernize-use-override, -modernize-use-ranges, +-modernize-use-std-numbers, -modernize-avoid-c-arrays, -modernize-pass-by-value, -modernize-use-equals-default, @@ -55,6 +56,7 @@ modernize*, -modernize-use-trailing-return-type, -modernize-use-override, -modernize-use-ranges, +-modernize-use-std-numbers, -modernize-use-using, -modernize-avoid-c-arrays, -modernize-pass-by-value, diff --git a/src/griddyn/exciters/ExciterAC7B.cpp b/src/griddyn/exciters/ExciterAC7B.cpp index e2e2a1c8..65a779ac 100644 --- a/src/griddyn/exciters/ExciterAC7B.cpp +++ b/src/griddyn/exciters/ExciterAC7B.cpp @@ -28,6 +28,11 @@ namespace { constexpr index_t feedbackState = 4; constexpr index_t exciterVoltageState = 5; constexpr double initializationTolerance = 1e-7; + // PSS/E's legacy rectifier initialization uses rounded reciprocal-sqrt(3) + // and sqrt(3) constants. These intentionally remain model constants rather + // than std::numbers::egamma or std::numbers::sqrt3. + constexpr double rectifierReciprocalSqrtThree = 0.577; + constexpr double rectifierSqrtThree = 1.732; bool finiteMachineSignal(double value) { @@ -112,13 +117,13 @@ double ExciterAC7B::solveExciterVoltage(double fieldVoltage, double fieldCurrent std::sqrt(((fieldVoltage * fieldVoltage) + (loading * loading)) / 0.75); const double middleSign = (fieldVoltage != 0.0) ? fieldVoltage : loading; const std::array candidates{fieldVoltage, - fieldVoltage + (0.577 * loading), + fieldVoltage + (rectifierReciprocalSqrtThree * loading), std::copysign(middleMagnitude, middleSign), - (fieldVoltage / 1.732) + loading, + (fieldVoltage / rectifierSqrtThree) + loading, loading}; double bestVoltage = candidates.front(); double bestMismatch = std::numeric_limits::infinity(); - for (double candidate : candidates) { + for (const double candidate : candidates) { const double normalizedCurrent = (std::abs(candidate) > 1e-14) ? loading / candidate : 0.0; const double output = candidate * detail::computeRectifierFactor(normalizedCurrent).factor; const double mismatch = std::abs(output - fieldVoltage); diff --git a/src/griddyn/exciters/ExciterAC8B.cpp b/src/griddyn/exciters/ExciterAC8B.cpp index 33e169c2..f27f16c4 100644 --- a/src/griddyn/exciters/ExciterAC8B.cpp +++ b/src/griddyn/exciters/ExciterAC8B.cpp @@ -27,6 +27,11 @@ namespace { constexpr index_t regulatorState = 3; constexpr index_t exciterVoltageState = 4; constexpr double initializationTolerance = 1e-7; + // PSS/E's legacy rectifier initialization uses rounded reciprocal-sqrt(3) + // and sqrt(3) constants. These intentionally remain model constants rather + // than std::numbers::egamma or std::numbers::sqrt3. + constexpr double rectifierReciprocalSqrtThree = 0.577; + constexpr double rectifierSqrtThree = 1.732; bool finiteMachineSignal(double value) { @@ -104,13 +109,13 @@ double ExciterAC8B::solveExciterVoltage(double fieldVoltage, double fieldCurrent std::sqrt(((fieldVoltage * fieldVoltage) + (loading * loading)) / 0.75); const double middleSign = (fieldVoltage != 0.0) ? fieldVoltage : loading; const std::array candidates{fieldVoltage, - fieldVoltage + (0.577 * loading), + fieldVoltage + (rectifierReciprocalSqrtThree * loading), std::copysign(middleMagnitude, middleSign), - (fieldVoltage / 1.732) + loading, + (fieldVoltage / rectifierSqrtThree) + loading, loading}; double bestVoltage = candidates.front(); double bestMismatch = std::numeric_limits::infinity(); - for (double candidate : candidates) { + for (const double candidate : candidates) { const double normalizedCurrent = (std::abs(candidate) > 1e-14) ? loading / candidate : 0.0; const double output = candidate * detail::computeRectifierFactor(normalizedCurrent).factor; const double mismatch = std::abs(output - fieldVoltage); diff --git a/src/griddyn/exciters/ExciterIEEET3.cpp b/src/griddyn/exciters/ExciterIEEET3.cpp index c3550139..c8133b47 100644 --- a/src/griddyn/exciters/ExciterIEEET3.cpp +++ b/src/griddyn/exciters/ExciterIEEET3.cpp @@ -148,8 +148,9 @@ ExciterIEEET3::Evaluation rates[regulatorState] = scaleSignal(regulatorDrive, 1.0 / Ta); } const Signal vbRaw = addSignals(regulator, v40); - const Signal vb = clampSignal(vbRaw, 0.0, Vbmax); - const Signal fieldDrive = subtractSignals(vb, scaleSignal(fieldStateSignal, Ke)); + const Signal boundedSourceVoltage = clampSignal(vbRaw, 0.0, Vbmax); + const Signal fieldDrive = + subtractSignals(boundedSourceVoltage, scaleSignal(fieldStateSignal, Ke)); rates[fieldState] = scaleSignal(fieldDrive, 1.0 / Te); const Signal algebraicResidual = subtractSignals(fieldStateSignal, algebraicSignal(fieldVoltage)); diff --git a/test/componentTests/testExciters.cpp b/test/componentTests/testExciters.cpp index 010e6409..442ee14f 100644 --- a/test/componentTests/testExciters.cpp +++ b/test/componentTests/testExciters.cpp @@ -185,6 +185,7 @@ void expectExciterJacobian(Exciter& exciter, double tolerance = 2e-5) { constexpr double step = 1e-6; + const index_t stateCount = static_cast(state.size()); ASSERT_EQ(exciter.stateSize(cDaeSolverMode), state.size()); exciter.setOffset(0, cDaeSolverMode); std::vector stateDerivative(state.size(), 0.0); @@ -207,14 +208,14 @@ void expectExciterJacobian(Exciter& exciter, exciter.residual(trialInputs, trialData, residual.data(), cDaeSolverMode); return residual; }; - for (index_t column = 0; column < state.size(); ++column) { + for (index_t column = 0; column < stateCount; ++column) { auto plus = state; auto minus = state; plus[column] += step; minus[column] -= step; const auto plusResidual = residualAt(inputs, plus); const auto minusResidual = residualAt(inputs, minus); - for (index_t row = 0; row < state.size(); ++row) { + for (index_t row = 0; row < stateCount; ++row) { const double numerical = (plusResidual[row] - minusResidual[row]) / (2.0 * step); EXPECT_NEAR(jacobian.at(row, column), numerical, tolerance) << "state row " << row << " column " << column; @@ -227,7 +228,7 @@ void expectExciterJacobian(Exciter& exciter, minus[column] -= step; const auto plusResidual = residualAt(plus, state); const auto minusResidual = residualAt(minus, state); - for (index_t row = 0; row < state.size(); ++row) { + for (index_t row = 0; row < stateCount; ++row) { const double numerical = (plusResidual[row] - minusResidual[row]) / (2.0 * step); EXPECT_NEAR(jacobian.at(row, inputColumnBase + column), numerical, tolerance) << "input row " << row << " column " << column; From aac1eefe54b7cccc87fa0ce95717e8b5a04b158a Mon Sep 17 00:00:00 2001 From: Philip Top Date: Mon, 7 Sep 2026 06:20:31 -0700 Subject: [PATCH 3/3] more clang-tidy --- .gitignore | 1 + test/componentTests/testExciters.cpp | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 6de724a1..fa412658 100644 --- a/.gitignore +++ b/.gitignore @@ -63,4 +63,5 @@ fmus /.pre-commit-cache-2/ .pre-commit-cache .vs +.julia output diff --git a/test/componentTests/testExciters.cpp b/test/componentTests/testExciters.cpp index 442ee14f..feeca9d8 100644 --- a/test/componentTests/testExciters.cpp +++ b/test/componentTests/testExciters.cpp @@ -185,7 +185,7 @@ void expectExciterJacobian(Exciter& exciter, double tolerance = 2e-5) { constexpr double step = 1e-6; - const index_t stateCount = static_cast(state.size()); + const auto stateCount = static_cast(state.size()); ASSERT_EQ(exciter.stateSize(cDaeSolverMode), state.size()); exciter.setOffset(0, cDaeSolverMode); std::vector stateDerivative(state.size(), 0.0);