From 2ac924d9ce13d355a0104d02f6d491008bb01c45 Mon Sep 17 00:00:00 2001 From: Slaven Peles Date: Tue, 3 Jun 2025 13:00:37 -0400 Subject: [PATCH 1/8] Copy of Abdou's model --- .../SynchronousMachine/CMakeLists.txt | 1 + .../GenClassical/CMakeLists.txt | 11 + .../GenClassical/GenClassical.cpp | 266 ++++++++++++++++++ .../GenClassical/GenClassical.hpp | 124 ++++++++ tests/UnitTests/PhasorDynamics/CMakeLists.txt | 8 +- .../PhasorDynamics/GenClassicalTests.hpp | 229 +++++++++++++++ .../PhasorDynamics/runGenClassicalTests.cpp | 15 + 7 files changed, 653 insertions(+), 1 deletion(-) create mode 100644 src/Model/PhasorDynamics/SynchronousMachine/GenClassical/CMakeLists.txt create mode 100644 src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp create mode 100644 src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp create mode 100644 tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp create mode 100644 tests/UnitTests/PhasorDynamics/runGenClassicalTests.cpp diff --git a/src/Model/PhasorDynamics/SynchronousMachine/CMakeLists.txt b/src/Model/PhasorDynamics/SynchronousMachine/CMakeLists.txt index 212813aae..96c3b2f94 100644 --- a/src/Model/PhasorDynamics/SynchronousMachine/CMakeLists.txt +++ b/src/Model/PhasorDynamics/SynchronousMachine/CMakeLists.txt @@ -6,3 +6,4 @@ # ]] add_subdirectory(GENROUwS) +add_subdirectory(GenClassical) diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/CMakeLists.txt b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/CMakeLists.txt new file mode 100644 index 000000000..67e36e73f --- /dev/null +++ b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/CMakeLists.txt @@ -0,0 +1,11 @@ +# [[ +# Author(s): +# - Cameron Rutherford +# - Slaven Peles +# ]] + +gridkit_add_library(phasor_dynamics_gen_classical + SOURCES + GenClassical.cpp + OUTPUT_NAME + gridkit_phasor_dynamics_gen_classical) diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp new file mode 100644 index 000000000..17676edb7 --- /dev/null +++ b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp @@ -0,0 +1,266 @@ +/** + * @file GenClassical.cpp + * @author Abdourahman Barry (abdourahman@vt.edu) + * @brief Definition of a Classical generator model. + * + * + */ + +#include "GenClassical.hpp" + +#include +#include + +#include + +#define _USE_MATH_DEFINES + +namespace GridKit +{ + namespace PhasorDynamics + { + /*! + * @brief Constructor for a pi-model branch + * + * Arguments passed to ModelEvaluatorImpl: + * - Number of equations = 0 + * - Number of independent variables = 0 + * - Number of quadratures = 0 + * - Number of optimization parameters = 0 + */ + template + GenClassical::GenClassical(bus_type* bus, int unit_id) + : bus_(bus), + busID_(0), + unit_id_(unit_id), + p0_(0.0), + q0_(0.0), + H_(3.0), + D_(0.0), + Ra_(0.0), + Xdp_(0.5) + { + size_ = 7; + setDerivedParams(); + + // Temporary, to eliminate compiler warnings + (void) busID_; + (void) unit_id_; + } + + /*! + * @brief Constructor for a pi-model branch + * + * Arguments passed to ModelEvaluatorImpl: + * - Number of equations = 0 + * - Number of independent variables = 0 + * - Number of quadratures = 0 + * - Number of optimization parameters = 0 + */ + template + GenClassical::GenClassical(bus_type* bus, + int unit_id, + ScalarT p0, + ScalarT q0, + real_type H, + real_type D, + real_type Ra, + real_type Xdp) + + : bus_(bus), + busID_(0), + unit_id_(unit_id), + p0_(p0), + q0_(q0), + H_(H), + D_(D), + Ra_(Ra), + Xdp_(Xdp) + { + size_ = 7; + setDerivedParams(); + } + + /*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ + template + int GenClassical::allocate() + { + f_.resize(size_); + y_.resize(size_); + yp_.resize(size_); + tag_.resize(size_); + fB_.resize(size_); + yB_.resize(size_); + ypB_.resize(size_); + return 0; + } + + /** + * Initialization of the branch model + * + */ + template + int GenClassical::initialize() + { + ScalarT vr = Vr(); + ScalarT vi = Vi(); + ScalarT p = p0_; + ScalarT q = q0_; + ScalarT vm2 = vr * vr + vi * vi; + ScalarT ir = (p * vr + q * vi) / vm2; + ScalarT ii = (p * vi - q * vr) / vm2; + ScalarT Er = (G_ * (ir + G_ * vr - B_ * vi) + B_ * (ii + B_ * vr + G_ * vi)) / (G_ * G_ + B_ * B_); + ScalarT Ei = (-B_ * (ir + G_ * vr - B_ * vi) + G_ * (ii + B_ * vr + G_ * vi)) / (G_ * G_ + B_ * B_); + ScalarT delta = atan2(Ei, Er); + ScalarT omega = 0; + ScalarT Ep = sqrt(Er * Er + Ei * Ei); + ScalarT Te = G_ * Ep * Ep - Ep * ((G_ * vr - B_ * vi) * cos(delta) + (B_ * vr + G_ * vi) * sin(delta)); + + y_[0] = delta; + y_[1] = omega; + y_[2] = Te; + y_[3] = ir; + y_[4] = ii; + y_[5] = pmech_set_ = Te; + y_[6] = ep_set_ = Ep; + + for (IdxT i = 0; i < size_; ++i) + yp_[i] = 0.0; + + return 0; + } + + /** + * \brief Identify differential variables. + */ + template + int GenClassical::tagDifferentiable() + { + + return 0; + } + + /** + * \brief Residual contribution of the branch is pushed to the + * two terminal buses. + * + */ + template + int GenClassical::evaluateResidual() + { + /* Read variables */ + ScalarT delta = y_[0]; + ScalarT omega = y_[1]; + ScalarT telec = y_[2]; + ScalarT ir = y_[3]; + ScalarT ii = y_[4]; + ScalarT pmech = y_[5]; + ScalarT ep = y_[6]; + + /* Read derivatives */ + ScalarT delta_dot = yp_[0]; + ScalarT omega_dot = yp_[1]; + + /* 6 GenClassical differential equations */ + f_[0] = delta_dot - omega * (2.0 * M_PI * 60.0); + f_[1] = omega_dot - (1.0 / (2 * H_)) * ((pmech - D_ * omega) / (1 + omega) - telec); + + /* 11 GenClassical algebraic equations */ + f_[2] = telec - (1.0 / (1.0 + omega)) * (G_ * ep * ep - ep * (cos(delta) * (G_ * Vr() - B_ * Vi()) + sin(delta) * (B_ * Vr() + G_ * Vi()))); + + f_[3] = ir + G_ * Vr() - B_ * Vi() - ep * (G_ * cos(delta) - B_ * sin(delta)); + f_[4] = ii + B_ * Vr() + G_ * Vi() - ep * (B_ * cos(delta) + G_ * sin(delta)); + + f_[5] = pmech - pmech_set_; + f_[6] = ep - ep_set_; + + Ir() += -G_ * Vr() + B_ * Vi() + ep * (G_ * cos(delta) - B_ * sin(delta)); + Ii() += -B_ * Vr() - G_ * Vi() + ep * (B_ * cos(delta) + G_ * sin(delta)); + + return 0; + } + + /** + * @brief Jacobian evaluation not implemented yet + * + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + * @return int - error code, 0 = success + */ + template + int GenClassical::evaluateJacobian() + { + return 0; + } + + /** + * @brief Integrand (objective) evaluation not implemented yet + * + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + * @return int - error code, 0 = success + */ + template + int GenClassical::evaluateIntegrand() + { + // std::cout << "Evaluate Integrand for GenClassical..." << std::endl; + return 0; + } + + /** + * @brief Adjoint initialization not implemented yet + * + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + * @return int - error code, 0 = success + */ + template + int GenClassical::initializeAdjoint() + { + // std::cout << "Initialize adjoint for GenClassical..." << std::endl; + return 0; + } + + /** + * @brief Adjoint residual evaluation not implemented yet + * + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + * @return int - error code, 0 = success + */ + template + int GenClassical::evaluateAdjointResidual() + { + // std::cout << "Evaluate adjoint residual for GenClassical..." << std::endl; + return 0; + } + + /** + * @brief Adjoint integrand (objective) evaluation not implemented yet + * + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + * @return int - error code, 0 = success + */ + template + int GenClassical::evaluateAdjointIntegrand() + { + // std::cout << "Evaluate adjoint Integrand for GenClassical..." << std::endl; + return 0; + } + + template + void GenClassical::setDerivedParams() + { + G_ = Ra_ / (Ra_ * Ra_ + Xdp_ * Xdp_); + B_ = Xdp_ / (Ra_ * Ra_ + Xdp_ * Xdp_); + } + + // Available template instantiations + template class GenClassical; + template class GenClassical; + + } // namespace PhasorDynamics +} // namespace GridKit \ No newline at end of file diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp new file mode 100644 index 000000000..bd41828da --- /dev/null +++ b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp @@ -0,0 +1,124 @@ +/** + * @file GenClassical.cpp + * @author Abdourahman Barry (abdourahman@vt.edu) + * + */ + +#pragma once + +#include + +// Forward declarations. +namespace GridKit +{ + namespace PhasorDynamics + { + template + class BusBase; + } +} // namespace GridKit + +namespace GridKit +{ + namespace PhasorDynamics + { + + template + class GenClassical : public Component + { + using Component::alpha_; + using Component::f_; + using Component::fB_; + using Component::g_; + using Component::gB_; + using Component::nnz_; + using Component::param_; + using Component::size_; + using Component::tag_; + using Component::time_; + using Component::y_; + using Component::yB_; + using Component::yp_; + using Component::ypB_; + + using bus_type = BusBase; + using real_type = typename Component::real_type; + + public: + GenClassical(bus_type* bus, int unit_id); + GenClassical(bus_type* bus, + int unit_id, + ScalarT p0, + ScalarT q0, + real_type H, + real_type D, + real_type Ra, + real_type Xdp); + ~GenClassical() = default; + + int allocate() override; + int initialize() override; + int tagDifferentiable() override; + int evaluateResidual() override; + + // Still to be implemented + int evaluateJacobian() override; + int evaluateIntegrand() override; + int initializeAdjoint() override; + int evaluateAdjointResidual() override; + int evaluateAdjointIntegrand() override; + + void updateTime(real_type /* t */, real_type /* a */) override + { + } + + private: + void setDerivedParams(); + + ScalarT& Vr() + { + return bus_->Vr(); + } + + ScalarT& Vi() + { + return bus_->Vi(); + } + + ScalarT& Ir() + { + return bus_->Ir(); + } + + ScalarT& Ii() + { + return bus_->Ii(); + } + + private: + /* Identification */ + bus_type* bus_; + const int busID_; + int unit_id_; + + /* Initial terminal conditions */ + ScalarT p0_; + ScalarT q0_; + + /* Input parameters */ + real_type H_; + real_type D_; + real_type Ra_; + real_type Xdp_; + + /* Derivied parameters */ + real_type G_; + real_type B_; + + /* Setpoints for control variables (determined at initialization) */ + real_type pmech_set_; + real_type ep_set_; + }; + + } // namespace PhasorDynamics +} // namespace GridKit \ No newline at end of file diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index 18876a998..dfe3c93d4 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -17,7 +17,11 @@ target_link_libraries(test_phasor_load GRIDKIT::phasor_dynamics_load add_executable(test_phasor_genrou runGenrouTests.cpp) target_link_libraries(test_phasor_genrou GRIDKIT::phasor_dynamics_genrou - GRIDKIT::phasor_dynamics_bus) + GRIDKIT::phasor_dynamics_bus) + +add_executable(test_phasor_gen_classical runGenClassicalTests.cpp) +target_link_libraries(test_phasor_gen_classical GRIDKIT::phasor_dynamics_gen_classical + GRIDKIT::phasor_dynamics_bus) add_executable(test_phasor_system runSystemTests.cpp) @@ -28,6 +32,7 @@ target_link_libraries(test_phasor_system GRIDKIT::phasor_dynamics_load add_test(NAME PhasorDynamicsBusTest COMMAND $) add_test(NAME PhasorDynamicsBranchTest COMMAND $) add_test(NAME PhasorDynamicsGenrouTest COMMAND $) +add_test(NAME PhasorDynamicsGenClassicalTest COMMAND $) add_test(NAME PhasorDynamicsLoadTest COMMAND $) add_test(NAME PhasorDynamicsSystemTest COMMAND $) @@ -35,4 +40,5 @@ install(TARGETS test_phasor_bus test_phasor_branch test_phasor_load test_phasor_genrou + test_phasor_gen_classical test_phasor_system RUNTIME DESTINATION bin) diff --git a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp new file mode 100644 index 000000000..94175a676 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp @@ -0,0 +1,229 @@ +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + + template + class GenClassicalTests + { + private: + using real_type = typename PhasorDynamics::Component::real_type; + + public: + GenClassicalTests() = default; + ~GenClassicalTests() = default; + + TestOutcome constructor() + { + TestStatus success = true; + + auto* bus = new PhasorDynamics::Bus(1.0, 0.0); + + PhasorDynamics::Component* machine = + new PhasorDynamics::GenClassical(bus, 1); + + success *= (machine != nullptr); + + if (machine) + { + delete machine; + } + delete bus; + + return success.report(__func__); + } + + /** + * A test case to verify residual values + */ + TestOutcome residual() + { + TestStatus success = true; + + // classical generator parameters + real_type H{0.5}; + real_type D{-1.0}; + real_type Ra{0.5}; + real_type Xdp{0.5}; + real_type pmech{1.0}; + real_type ep{2.0}; + + ScalarT Vr1{1.0}; ///< Bus-1 real voltage + ScalarT Vi1{1.0}; ///< Bus-1 imaginary voltage + + const ScalarT res0{0.0}; /// first residual + const ScalarT res1{0.0}; /// second residual + const ScalarT res2{0.0}; /// third residual + const ScalarT res3{0.0}; /// fourth residual + const ScalarT res4{0.0}; /// fifth residual + const ScalarT res5{0.0}; /// fifth residual + const ScalarT res6{0.0}; /// fifth residual + const ScalarT tol = 7 * (std::numeric_limits::epsilon()); // tolerance for comparing results + + PhasorDynamics::Bus bus(Vr1, Vi1); + PhasorDynamics::GenClassical gen(&bus, 1, 1.0, 1.0, H, D, Ra, Xdp); + bus.allocate(); + bus.initialize(); + gen.allocate(); + + gen.y()[0] = M_PI; // delta + gen.y()[1] = 1.0; // omega + gen.y()[2] = 2.0; // telec + gen.y()[1] = 1.0; // omega + gen.y()[2] = 2.0; // telec + gen.y()[3] = -2.0; // ir + gen.y()[4] = -4.0; // ii + gen.y()[5] = 1.0; // pmech + gen.y()[6] = 2.0; // Ep + gen.y()[5] = 1.0; // pmech + gen.y()[6] = 2.0; // Ep + + gen.yp()[0] = 2 * M_PI * 60.0; // delta_dot + gen.yp()[1] = -1.0; // omega_dot + gen.yp()[2] = 0.0; // telec + gen.yp()[3] = 0.0; // ir + gen.yp()[4] = 0.0; // ii + gen.yp()[5] = 0.0; // pmech + gen.yp()[6] = 0.0; // Ep + gen.yp()[0] = 2 * M_PI * 60.0; // delta_dot + gen.yp()[1] = -1.0; // omega_dot + gen.yp()[2] = 0.0; // telec + gen.yp()[3] = 0.0; // ir + gen.yp()[4] = 0.0; // ii + gen.yp()[5] = 0.0; // pmech + gen.yp()[6] = 0.0; // Ep + + gen.evaluateResidual(); + + std::vector residual = gen.getResidual(); + + success *= isEqual(residual[0], res0, tol); + success *= isEqual(residual[1], res1, tol); + success *= isEqual(residual[2], res2, tol); + success *= isEqual(residual[3], res3, tol); + success *= isEqual(residual[4], res4, tol); + // success *= isEqual(residual[5], res5, tol); + // success *= isEqual(residual[6], res6, tol); + + return success.report(__func__); + } + + /** + * + * Verifies correctness of the system initialization + */ + TestOutcome initial() + { + TestStatus success = true; + + // classical generator parameters + real_type p0{3.0}; + real_type q0{-1.0}; + real_type H{1.0}; + real_type D{1.0}; + real_type Ra{0.6}; + real_type Xdp{0.2}; + + ScalarT Vr1{1.0}; ///< Bus-1 real voltage + ScalarT Vi1{1.0}; ///< Bus-1 imaginary voltage + + const ScalarT delta{M_PI / 4.0}; /// first residual + const ScalarT omega{0.0}; /// second residual + const ScalarT Te{6.0}; /// third residual + const ScalarT ir{1.0}; /// fourth residual + const ScalarT ii{2.0}; /// fifth residual + const ScalarT pmech{6.0}; /// fifth residual + const ScalarT Ep{2.0 * sqrt(2.0)}; /// fifth residual + + const ScalarT tol = 5.0 * (std::numeric_limits::epsilon()); // tolerance for comparing result + + PhasorDynamics::Bus bus(Vr1, Vi1); + PhasorDynamics::GenClassical gen(&bus, 1, p0, q0, H, D, Ra, Xdp); + bus.allocate(); + bus.initialize(); + gen.allocate(); + gen.initialize(); + + success *= isEqual(gen.y()[0], delta, tol); + success *= isEqual(gen.y()[1], omega, tol); + success *= isEqual(gen.y()[2], Te, tol); + success *= isEqual(gen.y()[3], ir, tol); + success *= isEqual(gen.y()[4], ii, tol); + success *= isEqual(gen.y()[5], pmech, tol); + success *= isEqual(gen.y()[6], Ep, tol); + + success *= isEqual(gen.yp()[0], 0.0, tol); + success *= isEqual(gen.yp()[1], 0.0, tol); + success *= isEqual(gen.yp()[2], 0.0, tol); + success *= isEqual(gen.yp()[3], 0.0, tol); + success *= isEqual(gen.yp()[4], 0.0, tol); + success *= isEqual(gen.yp()[5], 0.0, tol); + success *= isEqual(gen.yp()[6], 0.0, tol); + + return success.report(__func__); + } + + /* + *Verifies the residual evaluates to zero for the initial conditions + */ + TestOutcome zeroInitialResidual() + { + TestStatus success = true; + + // classical generator parameters + real_type p0{3.0}; + real_type q0{-1.0}; + real_type H{1.0}; + real_type D{1.0}; + real_type Ra{0.6}; + real_type Xdp{0.2}; + + ScalarT Vr1{1.0}; ///< Bus-1 real voltage + ScalarT Vi1{1.0}; ///< Bus-1 imaginary voltage + + const ScalarT delta{M_PI / 4.0}; /// first residual + const ScalarT omega{0.0}; /// second residual + const ScalarT Te{6.0}; /// third residual + const ScalarT ir{1.0}; /// fourth residual + const ScalarT ii{2.0}; /// fifth residual + const ScalarT pmech{6.0}; /// fifth residual + const ScalarT Ep{2.0 * sqrt(2.0)}; /// fifth residual + + const ScalarT tol = 5.0 * (std::numeric_limits::epsilon()); // tolerance for comparing result + + PhasorDynamics::Bus bus(Vr1, Vi1); + PhasorDynamics::GenClassical gen(&bus, 1, p0, q0, H, D, Ra, Xdp); + bus.allocate(); + bus.initialize(); + gen.allocate(); + gen.initialize(); + gen.evaluateResidual(); + std::vector res = gen.getResidual(); + + success *= isEqual(res[0], 0.0, tol); + success *= isEqual(res[1], 0.0, tol); + success *= isEqual(res[2], 0.0, tol); + success *= isEqual(res[3], 0.0, tol); + success *= isEqual(res[4], 0.0, tol); + success *= isEqual(res[5], 0.0, tol); + success *= isEqual(res[6], 0.0, tol); + + return success.report(__func__); + } + + }; // class BranchTest + + } // namespace Testing +} // namespace GridKit \ No newline at end of file diff --git a/tests/UnitTests/PhasorDynamics/runGenClassicalTests.cpp b/tests/UnitTests/PhasorDynamics/runGenClassicalTests.cpp new file mode 100644 index 000000000..b01febe49 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runGenClassicalTests.cpp @@ -0,0 +1,15 @@ +#include "GenClassicalTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::GenClassicalTests test; + + result += test.constructor(); + result += test.residual(); + result += test.initial(); + result += test.zeroInitialResidual(); + + return result.summary(); +} \ No newline at end of file From 8952eb4c74d00aec73529f713f06e18cde39350e Mon Sep 17 00:00:00 2001 From: Slaven Peles Date: Tue, 3 Jun 2025 14:11:19 -0400 Subject: [PATCH 2/8] Partial, ill thought of changes. --- .../SynchronousMachine/GenClassical/GenClassical.cpp | 6 +++--- .../SynchronousMachine/GenClassical/README.md | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp index 17676edb7..8135b3511 100644 --- a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp +++ b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp @@ -164,11 +164,11 @@ namespace GridKit ScalarT omega_dot = yp_[1]; /* 6 GenClassical differential equations */ - f_[0] = delta_dot - omega * (2.0 * M_PI * 60.0); - f_[1] = omega_dot - (1.0 / (2 * H_)) * ((pmech - D_ * omega) / (1 + omega) - telec); + f_[0] = delta_dot - (omega - 1.0) * (2.0 * M_PI * 60.0); + f_[1] = omega_dot - (1.0 / (2 * H_)) * ((pmech - D_ * (omega - 1.0)) / omega - telec); /* 11 GenClassical algebraic equations */ - f_[2] = telec - (1.0 / (1.0 + omega)) * (G_ * ep * ep - ep * (cos(delta) * (G_ * Vr() - B_ * Vi()) + sin(delta) * (B_ * Vr() + G_ * Vi()))); + f_[2] = telec - (1.0 / omega) * (G_ * ep * ep - ep * ((G_ * Vr() - B_ * Vi()) * cos(delta) + (B_ * Vr() + G_ * Vi()) * sin(delta))); f_[3] = ir + G_ * Vr() - B_ * Vi() - ep * (G_ * cos(delta) - B_ * sin(delta)); f_[4] = ii + B_ * Vr() + G_ * Vi() - ep * (B_ * cos(delta) + G_ * sin(delta)); diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md index aa1dbe490..89618ee5f 100644 --- a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md +++ b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md @@ -85,7 +85,7 @@ $E_p$ | [p.u.] | field winding voltage | owned by exciter, constant if \begin{aligned} 0 &= T_{e} - \frac{1}{\omega}\left( g E_p^2 - E_p \left[(gV_r - bV_i)\cos\delta + (bV_r + gV_i)\sin\delta \right]\right)\\ 0 &= I_r + gV_r - bV_i - E_p(g \cos\delta - b \sin\delta) \\ - 0 &= I_i + gV_r + bV_i - E_p(b \cos\delta + g \sin\delta) + 0 &= I_i + bV_r + gV_i - E_p(b \cos\delta + g \sin\delta) \end{aligned} ``` As noted earlier, all three algebraic equations can be expressed as functions From d5ff2c8e8470150ff39681068b4b1d22f52270d9 Mon Sep 17 00:00:00 2001 From: Slaven Peles Date: Tue, 3 Jun 2025 18:56:30 -0400 Subject: [PATCH 3/8] Gen classical unit tests pass. --- .../GenClassical/GenClassical.cpp | 6 +- .../SynchronousMachine/GenClassical/README.md | 56 +++++++++++++------ .../PhasorDynamics/GenClassicalTests.hpp | 35 +++--------- 3 files changed, 51 insertions(+), 46 deletions(-) diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp index 8135b3511..d8107143e 100644 --- a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp +++ b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp @@ -111,10 +111,10 @@ namespace GridKit ScalarT vm2 = vr * vr + vi * vi; ScalarT ir = (p * vr + q * vi) / vm2; ScalarT ii = (p * vi - q * vr) / vm2; - ScalarT Er = (G_ * (ir + G_ * vr - B_ * vi) + B_ * (ii + B_ * vr + G_ * vi)) / (G_ * G_ + B_ * B_); - ScalarT Ei = (-B_ * (ir + G_ * vr - B_ * vi) + G_ * (ii + B_ * vr + G_ * vi)) / (G_ * G_ + B_ * B_); + ScalarT Er = (G_ * ir + B_ * ii) / (G_ * G_ + B_ * B_) + vr; + ScalarT Ei = (-B_ * ir + G_ * ii) / (G_ * G_ + B_ * B_) + vi; ScalarT delta = atan2(Ei, Er); - ScalarT omega = 0; + ScalarT omega = 1.0; ScalarT Ep = sqrt(Er * Er + Ei * Ei); ScalarT Te = G_ * Ep * Ep - Ep * ((G_ * vr - B_ * vi) * cos(delta) + (B_ * vr + G_ * vi) * sin(delta)); diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md index 89618ee5f..9e6eac597 100644 --- a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md +++ b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md @@ -17,8 +17,8 @@ $X_{dp}$ | [p.u.] | machine reactance parameter | ### Model Derived Parameters -- $g = \dfrac{R_a}{R_a^2 + X_{dp}^2}$ -- $b = \dfrac{-X_{dp}}{R_a^2 + X_{dp}^2}$ +- $G = \dfrac{R_a}{R_a^2 + X_{dp}^2}$ +- $B = \dfrac{-X_{dp}}{R_a^2 + X_{dp}^2}$
@@ -83,9 +83,9 @@ $E_p$ | [p.u.] | field winding voltage | owned by exciter, constant if ```math \begin{aligned} - 0 &= T_{e} - \frac{1}{\omega}\left( g E_p^2 - E_p \left[(gV_r - bV_i)\cos\delta + (bV_r + gV_i)\sin\delta \right]\right)\\ - 0 &= I_r + gV_r - bV_i - E_p(g \cos\delta - b \sin\delta) \\ - 0 &= I_i + bV_r + gV_i - E_p(b \cos\delta + g \sin\delta) + 0 &= T_{e} - \frac{1}{\omega}\left( G E_p^2 - E_p \left[(G V_r - B V_i)\cos\delta + (B V_r + G V_i)\sin\delta \right]\right)\\ + 0 &= I_r + G V_r - B V_i - E_p(G \cos\delta - B \sin\delta) \\ + 0 &= I_i + B V_r + G V_i - E_p(B \cos\delta + G \sin\delta) \end{aligned} ``` As noted earlier, all three algebraic equations can be expressed as functions @@ -100,8 +100,16 @@ To initialize the model, given bus voltages $V_r$, $V_i$, and initial generator injection active and reactive power, $P$ and $Q$, we take following steps to initialize the system: -First compute injection currents from initial power injection power and bus -voltages: +Complex power is defined as +```math +S=VI^{*} +``` +or +```math +P + jQ = (V_r + j V_i)(I_r - j I_i). +``` +From here, we compute injection currents from the initial power injection and bus +voltages as ```math \begin{aligned} I_r &= \frac{PV_r + QV_i}{V_r^2 + V_i^2} \\ @@ -109,27 +117,43 @@ I_i &= \frac{PV_i - QV_r}{V_r^2 + V_i^2} \end{aligned} ``` -Next compute field winding voltage and machine angle: +We substitute expressions above into equations for current injections and +obtain +```math +\begin{align} +E_p \sin\delta &= \dfrac{-B I_r + G I_i}{G^2 + B^2} + V_i \\ +E_p \cos\delta &= \dfrac{ G I_r + B I_i}{G^2 + B^2} + V_r +\end{align} +``` +By dividing (1) by (2) we get expression for machine angle at the +steady state +```math +\delta = \arctan \dfrac{E_i}{E_r} \, , +``` +and by squaring and adding (1) and (2), we get expression for field +winding voltage at the steady state +```math +E_p = \sqrt{E_r^2 + E_i^2} \, , +``` +where ```math \begin{aligned} -E_r &= \frac{ g(I_r + gV_r - bV_i) + b (I_i + bV_r + gV_i) }{g^2 + b^2} \\ -E_i &= \frac{ -b(I_r + gV_r - bV_i) + g (I_i + bV_r + gV_i) }{g^2 + b^2} \\ -E_p &= \sqrt{E_r^2 + E_i^2} \\ -\delta &= \arctan \dfrac{E_i}{E_r} +E_r &= \dfrac{ G I_r + B I_i}{G^2 + B^2} + V_r \, ,\\ +E_i &= \dfrac{-B I_r + G I_i}{G^2 + B^2} + V_i \, . \end{aligned} ``` -Set machine speed to the synchronous speed: +Next, we set machine speed to the synchronous speed: ```math \omega = 1 ``` Now, we can compute electrical torque and set mechanical torque to be equal -to the electrical. +to the electrical: ```math \begin{aligned} -T_{elec} &= gE_p^2 - E_p \left[ (gV_r - bV_i ) \cos\delta + (bV_r + gV_i )\sin\delta \right] \\ -P_{mech} &= T_{elec} +T_{e} &= G E_p^2 - E_p \left[ (G V_r - B V_i ) \cos\delta + (B V_r + G V_i )\sin\delta \right] \\ +P_{m} &= T_{e} \end{aligned} ``` diff --git a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp index 94175a676..738819678 100644 --- a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp @@ -68,8 +68,8 @@ namespace GridKit const ScalarT res2{0.0}; /// third residual const ScalarT res3{0.0}; /// fourth residual const ScalarT res4{0.0}; /// fifth residual - const ScalarT res5{0.0}; /// fifth residual - const ScalarT res6{0.0}; /// fifth residual + const ScalarT res5{pmech}; /// fifth residual + const ScalarT res6{ep}; /// fifth residual const ScalarT tol = 7 * (std::numeric_limits::epsilon()); // tolerance for comparing results PhasorDynamics::Bus bus(Vr1, Vi1); @@ -79,24 +79,13 @@ namespace GridKit gen.allocate(); gen.y()[0] = M_PI; // delta - gen.y()[1] = 1.0; // omega - gen.y()[2] = 2.0; // telec - gen.y()[1] = 1.0; // omega + gen.y()[1] = 2.0; // omega gen.y()[2] = 2.0; // telec gen.y()[3] = -2.0; // ir gen.y()[4] = -4.0; // ii - gen.y()[5] = 1.0; // pmech - gen.y()[6] = 2.0; // Ep - gen.y()[5] = 1.0; // pmech - gen.y()[6] = 2.0; // Ep + gen.y()[5] = pmech; // pmech + gen.y()[6] = ep; // Ep - gen.yp()[0] = 2 * M_PI * 60.0; // delta_dot - gen.yp()[1] = -1.0; // omega_dot - gen.yp()[2] = 0.0; // telec - gen.yp()[3] = 0.0; // ir - gen.yp()[4] = 0.0; // ii - gen.yp()[5] = 0.0; // pmech - gen.yp()[6] = 0.0; // Ep gen.yp()[0] = 2 * M_PI * 60.0; // delta_dot gen.yp()[1] = -1.0; // omega_dot gen.yp()[2] = 0.0; // telec @@ -114,8 +103,8 @@ namespace GridKit success *= isEqual(residual[2], res2, tol); success *= isEqual(residual[3], res3, tol); success *= isEqual(residual[4], res4, tol); - // success *= isEqual(residual[5], res5, tol); - // success *= isEqual(residual[6], res6, tol); + success *= isEqual(residual[5], res5, tol); + success *= isEqual(residual[6], res6, tol); return success.report(__func__); } @@ -140,7 +129,7 @@ namespace GridKit ScalarT Vi1{1.0}; ///< Bus-1 imaginary voltage const ScalarT delta{M_PI / 4.0}; /// first residual - const ScalarT omega{0.0}; /// second residual + const ScalarT omega{1.0}; /// second residual const ScalarT Te{6.0}; /// third residual const ScalarT ir{1.0}; /// fourth residual const ScalarT ii{2.0}; /// fifth residual @@ -193,14 +182,6 @@ namespace GridKit ScalarT Vr1{1.0}; ///< Bus-1 real voltage ScalarT Vi1{1.0}; ///< Bus-1 imaginary voltage - const ScalarT delta{M_PI / 4.0}; /// first residual - const ScalarT omega{0.0}; /// second residual - const ScalarT Te{6.0}; /// third residual - const ScalarT ir{1.0}; /// fourth residual - const ScalarT ii{2.0}; /// fifth residual - const ScalarT pmech{6.0}; /// fifth residual - const ScalarT Ep{2.0 * sqrt(2.0)}; /// fifth residual - const ScalarT tol = 5.0 * (std::numeric_limits::epsilon()); // tolerance for comparing result PhasorDynamics::Bus bus(Vr1, Vi1); From ce59e7d5c33454eb9d0e387115efdbc8106d1176 Mon Sep 17 00:00:00 2001 From: Slaven Peles Date: Wed, 4 Jun 2025 11:51:57 -0400 Subject: [PATCH 4/8] Clean up GenClassical unit tests --- .../PhasorDynamics/GenClassicalTests.hpp | 149 ++++++++++-------- 1 file changed, 83 insertions(+), 66 deletions(-) diff --git a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp index 738819678..b02aeac2d 100644 --- a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp @@ -1,3 +1,10 @@ +/** + * @file GenClassicalTests.hpp + * @author Slaven Peles (peless@ornl.gov) + * @author Abdourahman Barry (abdourahman@vt.edu) + * @brief Tests for classical generator model. + * + */ #include #include #include @@ -20,6 +27,7 @@ namespace GridKit { private: using real_type = typename PhasorDynamics::Component::real_type; + static constexpr ScalarT tol_ = 10 * std::numeric_limits::epsilon(); public: GenClassicalTests() = default; @@ -63,48 +71,54 @@ namespace GridKit ScalarT Vr1{1.0}; ///< Bus-1 real voltage ScalarT Vi1{1.0}; ///< Bus-1 imaginary voltage - const ScalarT res0{0.0}; /// first residual - const ScalarT res1{0.0}; /// second residual - const ScalarT res2{0.0}; /// third residual - const ScalarT res3{0.0}; /// fourth residual - const ScalarT res4{0.0}; /// fifth residual - const ScalarT res5{pmech}; /// fifth residual - const ScalarT res6{ep}; /// fifth residual - const ScalarT tol = 7 * (std::numeric_limits::epsilon()); // tolerance for comparing results + // Test answer keys + const std::vector res_answer = {0.0, + 0.0, + 0.0, + 0.0, + 0.0, + pmech, + ep}; PhasorDynamics::Bus bus(Vr1, Vi1); PhasorDynamics::GenClassical gen(&bus, 1, 1.0, 1.0, H, D, Ra, Xdp); bus.allocate(); bus.initialize(); + + // Allocate but not initialize genrator model gen.allocate(); - gen.y()[0] = M_PI; // delta - gen.y()[1] = 2.0; // omega - gen.y()[2] = 2.0; // telec - gen.y()[3] = -2.0; // ir - gen.y()[4] = -4.0; // ii - gen.y()[5] = pmech; // pmech - gen.y()[6] = ep; // Ep + // Set variable values matching the answer key + gen.y()[0] = M_PI; // delta + gen.y()[1] = 2.0; // omega + gen.y()[2] = 2.0; // telec + gen.y()[3] = -2.0; // ir + gen.y()[4] = -4.0; // ii + gen.y()[5] = pmech; // pmech + gen.y()[6] = ep; // Ep + // Set derivative values matching the answer key gen.yp()[0] = 2 * M_PI * 60.0; // delta_dot gen.yp()[1] = -1.0; // omega_dot - gen.yp()[2] = 0.0; // telec - gen.yp()[3] = 0.0; // ir - gen.yp()[4] = 0.0; // ii - gen.yp()[5] = 0.0; // pmech - gen.yp()[6] = 0.0; // Ep + gen.yp()[2] = 0.0; + gen.yp()[3] = 0.0; + gen.yp()[4] = 0.0; + gen.yp()[5] = 0.0; + gen.yp()[6] = 0.0; gen.evaluateResidual(); + std::vector& residual = gen.getResidual(); - std::vector residual = gen.getResidual(); - - success *= isEqual(residual[0], res0, tol); - success *= isEqual(residual[1], res1, tol); - success *= isEqual(residual[2], res2, tol); - success *= isEqual(residual[3], res3, tol); - success *= isEqual(residual[4], res4, tol); - success *= isEqual(residual[5], res5, tol); - success *= isEqual(residual[6], res6, tol); + for (size_t i = 0; i < res_answer.size(); ++i) + { + if (!isEqual(residual[i], res_answer[i], tol_)) + { + std::cout << "Incorrect result: " + << residual[i] << " != " << res_answer[i] << "\n"; + success = false; + break; + } + } return success.report(__func__); } @@ -128,15 +142,14 @@ namespace GridKit ScalarT Vr1{1.0}; ///< Bus-1 real voltage ScalarT Vi1{1.0}; ///< Bus-1 imaginary voltage - const ScalarT delta{M_PI / 4.0}; /// first residual - const ScalarT omega{1.0}; /// second residual - const ScalarT Te{6.0}; /// third residual - const ScalarT ir{1.0}; /// fourth residual - const ScalarT ii{2.0}; /// fifth residual - const ScalarT pmech{6.0}; /// fifth residual - const ScalarT Ep{2.0 * sqrt(2.0)}; /// fifth residual - - const ScalarT tol = 5.0 * (std::numeric_limits::epsilon()); // tolerance for comparing result + // Test answer keys + const std::vector var_answer = {M_PI / 4.0, // delta + 1.0, // omega + 6.0, // Te + 1.0, // Ir + 2.0, // Ii + 6.0, // Pm + 2.0 * sqrt(2.0)}; // Ep PhasorDynamics::Bus bus(Vr1, Vi1); PhasorDynamics::GenClassical gen(&bus, 1, p0, q0, H, D, Ra, Xdp); @@ -145,21 +158,24 @@ namespace GridKit gen.allocate(); gen.initialize(); - success *= isEqual(gen.y()[0], delta, tol); - success *= isEqual(gen.y()[1], omega, tol); - success *= isEqual(gen.y()[2], Te, tol); - success *= isEqual(gen.y()[3], ir, tol); - success *= isEqual(gen.y()[4], ii, tol); - success *= isEqual(gen.y()[5], pmech, tol); - success *= isEqual(gen.y()[6], Ep, tol); - - success *= isEqual(gen.yp()[0], 0.0, tol); - success *= isEqual(gen.yp()[1], 0.0, tol); - success *= isEqual(gen.yp()[2], 0.0, tol); - success *= isEqual(gen.yp()[3], 0.0, tol); - success *= isEqual(gen.yp()[4], 0.0, tol); - success *= isEqual(gen.yp()[5], 0.0, tol); - success *= isEqual(gen.yp()[6], 0.0, tol); + for (size_t i = 0; i < var_answer.size(); ++i) + { + if (!isEqual(gen.y()[i], var_answer[i], tol_)) + { + std::cout << "Incorrect result: " + << gen.y()[i] << " != " << var_answer[i] << "\n"; + success = false; + break; + } + + if (!isEqual(gen.yp()[i], 0.0, tol_)) + { + std::cout << "Incorrect result: " + << gen.yp()[i] << " != 0\n"; + success = false; + break; + } + } return success.report(__func__); } @@ -179,10 +195,8 @@ namespace GridKit real_type Ra{0.6}; real_type Xdp{0.2}; - ScalarT Vr1{1.0}; ///< Bus-1 real voltage - ScalarT Vi1{1.0}; ///< Bus-1 imaginary voltage - - const ScalarT tol = 5.0 * (std::numeric_limits::epsilon()); // tolerance for comparing result + ScalarT Vr1{1.0}; ///< Bus real voltage + ScalarT Vi1{1.0}; ///< Bus imaginary voltage PhasorDynamics::Bus bus(Vr1, Vi1); PhasorDynamics::GenClassical gen(&bus, 1, p0, q0, H, D, Ra, Xdp); @@ -191,20 +205,23 @@ namespace GridKit gen.allocate(); gen.initialize(); gen.evaluateResidual(); - std::vector res = gen.getResidual(); + std::vector res = gen.getResidual(); - success *= isEqual(res[0], 0.0, tol); - success *= isEqual(res[1], 0.0, tol); - success *= isEqual(res[2], 0.0, tol); - success *= isEqual(res[3], 0.0, tol); - success *= isEqual(res[4], 0.0, tol); - success *= isEqual(res[5], 0.0, tol); - success *= isEqual(res[6], 0.0, tol); + for (size_t i = 0; i < res.size(); ++i) + { + if (!isEqual(res[i], 0.0, tol_)) + { + std::cout << "Incorrect result: " + << gen.yp()[i] << " != 0\n"; + success = false; + break; + } + } return success.report(__func__); } - }; // class BranchTest + }; // class GenClassicalTests } // namespace Testing } // namespace GridKit \ No newline at end of file From 3d14a53e1fd5f8ca5ed18d219fe6da67d0877346 Mon Sep 17 00:00:00 2001 From: Slaven Peles Date: Wed, 4 Jun 2025 13:01:40 -0400 Subject: [PATCH 5/8] Fix incorrect sign in GenClassical. --- .../GenClassical/GenClassical.cpp | 38 +++++++------------ .../SynchronousMachine/GenClassical/README.md | 24 ++++++------ 2 files changed, 26 insertions(+), 36 deletions(-) diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp index d8107143e..82384684c 100644 --- a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp +++ b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp @@ -1,6 +1,7 @@ /** * @file GenClassical.cpp * @author Abdourahman Barry (abdourahman@vt.edu) + * @author Slaven Peles (peless@ornl.gov) * @brief Definition of a Classical generator model. * * @@ -19,14 +20,9 @@ namespace GridKit { namespace PhasorDynamics { - /*! - * @brief Constructor for a pi-model branch + /** + * @brief Constructor for a classical generator model. * - * Arguments passed to ModelEvaluatorImpl: - * - Number of equations = 0 - * - Number of independent variables = 0 - * - Number of quadratures = 0 - * - Number of optimization parameters = 0 */ template GenClassical::GenClassical(bus_type* bus, int unit_id) @@ -51,11 +47,6 @@ namespace GridKit /*! * @brief Constructor for a pi-model branch * - * Arguments passed to ModelEvaluatorImpl: - * - Number of equations = 0 - * - Number of independent variables = 0 - * - Number of quadratures = 0 - * - Number of optimization parameters = 0 */ template GenClassical::GenClassical(bus_type* bus, @@ -98,7 +89,7 @@ namespace GridKit } /** - * Initialization of the branch model + * Initialization of the generator model * */ template @@ -111,12 +102,12 @@ namespace GridKit ScalarT vm2 = vr * vr + vi * vi; ScalarT ir = (p * vr + q * vi) / vm2; ScalarT ii = (p * vi - q * vr) / vm2; - ScalarT Er = (G_ * ir + B_ * ii) / (G_ * G_ + B_ * B_) + vr; - ScalarT Ei = (-B_ * ir + G_ * ii) / (G_ * G_ + B_ * B_) + vi; + ScalarT Er = (G_ * ir - B_ * ii) / (G_ * G_ + B_ * B_) + vr; + ScalarT Ei = (B_ * ir + G_ * ii) / (G_ * G_ + B_ * B_) + vi; ScalarT delta = atan2(Ei, Er); ScalarT omega = 1.0; ScalarT Ep = sqrt(Er * Er + Ei * Ei); - ScalarT Te = G_ * Ep * Ep - Ep * ((G_ * vr - B_ * vi) * cos(delta) + (B_ * vr + G_ * vi) * sin(delta)); + ScalarT Te = G_ * Ep * Ep - Ep * ((G_ * vr + B_ * vi) * cos(delta) + (-B_ * vr + G_ * vi) * sin(delta)); y_[0] = delta; y_[1] = omega; @@ -143,8 +134,7 @@ namespace GridKit } /** - * \brief Residual contribution of the branch is pushed to the - * two terminal buses. + * \brief Residual for the generator model. * */ template @@ -168,16 +158,16 @@ namespace GridKit f_[1] = omega_dot - (1.0 / (2 * H_)) * ((pmech - D_ * (omega - 1.0)) / omega - telec); /* 11 GenClassical algebraic equations */ - f_[2] = telec - (1.0 / omega) * (G_ * ep * ep - ep * ((G_ * Vr() - B_ * Vi()) * cos(delta) + (B_ * Vr() + G_ * Vi()) * sin(delta))); + f_[2] = telec - (1.0 / omega) * (G_ * ep * ep - ep * ((G_ * Vr() + B_ * Vi()) * cos(delta) + (-B_ * Vr() + G_ * Vi()) * sin(delta))); - f_[3] = ir + G_ * Vr() - B_ * Vi() - ep * (G_ * cos(delta) - B_ * sin(delta)); - f_[4] = ii + B_ * Vr() + G_ * Vi() - ep * (B_ * cos(delta) + G_ * sin(delta)); + f_[3] = ir + G_ * Vr() + B_ * Vi() - ep * (G_ * cos(delta) + B_ * sin(delta)); + f_[4] = ii - B_ * Vr() + G_ * Vi() - ep * (-B_ * cos(delta) + G_ * sin(delta)); f_[5] = pmech - pmech_set_; f_[6] = ep - ep_set_; - Ir() += -G_ * Vr() + B_ * Vi() + ep * (G_ * cos(delta) - B_ * sin(delta)); - Ii() += -B_ * Vr() - G_ * Vi() + ep * (B_ * cos(delta) + G_ * sin(delta)); + Ir() += ir; + Ii() += ii; return 0; } @@ -255,7 +245,7 @@ namespace GridKit void GenClassical::setDerivedParams() { G_ = Ra_ / (Ra_ * Ra_ + Xdp_ * Xdp_); - B_ = Xdp_ / (Ra_ * Ra_ + Xdp_ * Xdp_); + B_ = -Xdp_ / (Ra_ * Ra_ + Xdp_ * Xdp_); } // Available template instantiations diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md index 9e6eac597..3fdb710a3 100644 --- a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md +++ b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md @@ -83,9 +83,9 @@ $E_p$ | [p.u.] | field winding voltage | owned by exciter, constant if ```math \begin{aligned} - 0 &= T_{e} - \frac{1}{\omega}\left( G E_p^2 - E_p \left[(G V_r - B V_i)\cos\delta + (B V_r + G V_i)\sin\delta \right]\right)\\ - 0 &= I_r + G V_r - B V_i - E_p(G \cos\delta - B \sin\delta) \\ - 0 &= I_i + B V_r + G V_i - E_p(B \cos\delta + G \sin\delta) + 0 &= T_{e} - \frac{1}{\omega}\left( G E_p^2 - E_p \left[(G V_r + B V_i)\cos\delta + (-B V_r + G V_i)\sin\delta \right]\right) \\ + 0 &= I_r + G V_r + B V_i - E_p(G \cos\delta + B \sin\delta) \\ + 0 &= I_i - B V_r + G V_i - E_p(-B \cos\delta + G \sin\delta) \end{aligned} ``` As noted earlier, all three algebraic equations can be expressed as functions @@ -120,17 +120,17 @@ I_i &= \frac{PV_i - QV_r}{V_r^2 + V_i^2} We substitute expressions above into equations for current injections and obtain ```math -\begin{align} -E_p \sin\delta &= \dfrac{-B I_r + G I_i}{G^2 + B^2} + V_i \\ -E_p \cos\delta &= \dfrac{ G I_r + B I_i}{G^2 + B^2} + V_r -\end{align} +\begin{aligned} +E_p \sin\delta &= \dfrac{B I_r + G I_i}{G^2 + B^2} + V_i \\ +E_p \cos\delta &= \dfrac{G I_r - B I_i}{G^2 + B^2} + V_r +\end{aligned} ``` -By dividing (1) by (2) we get expression for machine angle at the +By dividing these two equations we get expression for machine angle at the steady state ```math \delta = \arctan \dfrac{E_i}{E_r} \, , ``` -and by squaring and adding (1) and (2), we get expression for field +and by squaring and adding them, we get expression for field winding voltage at the steady state ```math E_p = \sqrt{E_r^2 + E_i^2} \, , @@ -138,8 +138,8 @@ E_p = \sqrt{E_r^2 + E_i^2} \, , where ```math \begin{aligned} -E_r &= \dfrac{ G I_r + B I_i}{G^2 + B^2} + V_r \, ,\\ -E_i &= \dfrac{-B I_r + G I_i}{G^2 + B^2} + V_i \, . +E_r &= \dfrac{G I_r - B I_i}{G^2 + B^2} + V_r \, ,\\ +E_i &= \dfrac{B I_r + G I_i}{G^2 + B^2} + V_i \, . \end{aligned} ``` @@ -152,7 +152,7 @@ Now, we can compute electrical torque and set mechanical torque to be equal to the electrical: ```math \begin{aligned} -T_{e} &= G E_p^2 - E_p \left[ (G V_r - B V_i ) \cos\delta + (B V_r + G V_i )\sin\delta \right] \\ +T_{e} &= G E_p^2 - E_p \left[ (G V_r + B V_i ) \cos\delta + (-B V_r + G V_i )\sin\delta \right] \\ P_{m} &= T_{e} \end{aligned} ``` From 753843cf91a0623bf130f267f4dc4fdabceba380 Mon Sep 17 00:00:00 2001 From: pelesh Date: Wed, 4 Jun 2025 17:10:26 +0000 Subject: [PATCH 6/8] Apply pre-commmit fixes --- .../GenClassical/GenClassical.cpp | 2 +- .../GenClassical/GenClassical.hpp | 2 +- .../PhasorDynamics/GenClassicalTests.hpp | 20 +++++++++---------- .../PhasorDynamics/runGenClassicalTests.cpp | 2 +- 4 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp index 82384684c..12db533b7 100644 --- a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp +++ b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp @@ -253,4 +253,4 @@ namespace GridKit template class GenClassical; } // namespace PhasorDynamics -} // namespace GridKit \ No newline at end of file +} // namespace GridKit diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp index bd41828da..1e8a46550 100644 --- a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp +++ b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp @@ -121,4 +121,4 @@ namespace GridKit }; } // namespace PhasorDynamics -} // namespace GridKit \ No newline at end of file +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp index b02aeac2d..3c17c8ab1 100644 --- a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp @@ -3,7 +3,7 @@ * @author Slaven Peles (peless@ornl.gov) * @author Abdourahman Barry (abdourahman@vt.edu) * @brief Tests for classical generator model. - * + * */ #include #include @@ -26,8 +26,8 @@ namespace GridKit class GenClassicalTests { private: - using real_type = typename PhasorDynamics::Component::real_type; - static constexpr ScalarT tol_ = 10 * std::numeric_limits::epsilon(); + using real_type = typename PhasorDynamics::Component::real_type; + static constexpr ScalarT tol_ = 10 * std::numeric_limits::epsilon(); public: GenClassicalTests() = default; @@ -143,12 +143,12 @@ namespace GridKit ScalarT Vi1{1.0}; ///< Bus-1 imaginary voltage // Test answer keys - const std::vector var_answer = {M_PI / 4.0, // delta - 1.0, // omega - 6.0, // Te - 1.0, // Ir - 2.0, // Ii - 6.0, // Pm + const std::vector var_answer = {M_PI / 4.0, // delta + 1.0, // omega + 6.0, // Te + 1.0, // Ir + 2.0, // Ii + 6.0, // Pm 2.0 * sqrt(2.0)}; // Ep PhasorDynamics::Bus bus(Vr1, Vi1); @@ -224,4 +224,4 @@ namespace GridKit }; // class GenClassicalTests } // namespace Testing -} // namespace GridKit \ No newline at end of file +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/runGenClassicalTests.cpp b/tests/UnitTests/PhasorDynamics/runGenClassicalTests.cpp index b01febe49..21b7f45ba 100644 --- a/tests/UnitTests/PhasorDynamics/runGenClassicalTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runGenClassicalTests.cpp @@ -12,4 +12,4 @@ int main() result += test.zeroInitialResidual(); return result.summary(); -} \ No newline at end of file +} From 25693a80b7ab110b26ee5d001c188f7dc029864f Mon Sep 17 00:00:00 2001 From: Slaven Peles Date: Wed, 4 Jun 2025 17:30:33 -0400 Subject: [PATCH 7/8] Fix warnings in GenClassical model. --- .../GenClassical/GenClassical.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp index 12db533b7..9d4a9add3 100644 --- a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp +++ b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp @@ -78,13 +78,14 @@ namespace GridKit template int GenClassical::allocate() { - f_.resize(size_); - y_.resize(size_); - yp_.resize(size_); - tag_.resize(size_); - fB_.resize(size_); - yB_.resize(size_); - ypB_.resize(size_); + auto size = static_cast(size_); + f_.resize(size); + y_.resize(size); + yp_.resize(size); + tag_.resize(size); + fB_.resize(size); + yB_.resize(size); + ypB_.resize(size); return 0; } @@ -117,7 +118,7 @@ namespace GridKit y_[5] = pmech_set_ = Te; y_[6] = ep_set_ = Ep; - for (IdxT i = 0; i < size_; ++i) + for (size_t i = 0; i < static_cast(size_); ++i) yp_[i] = 0.0; return 0; From 433d27bd1fb8d60fe9ab5efdb92b090b65a68eed Mon Sep 17 00:00:00 2001 From: pelesh Date: Thu, 5 Jun 2025 11:02:00 -0400 Subject: [PATCH 8/8] Minor tweaks to comments and docs. --- .../GenClassical/GenClassical.cpp | 30 +++++++++---------- .../SynchronousMachine/GenClassical/README.md | 4 +-- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp index 9d4a9add3..b34180330 100644 --- a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp +++ b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp @@ -57,7 +57,6 @@ namespace GridKit real_type D, real_type Ra, real_type Xdp) - : bus_(bus), busID_(0), unit_id_(unit_id), @@ -141,24 +140,24 @@ namespace GridKit template int GenClassical::evaluateResidual() { - /* Read variables */ - ScalarT delta = y_[0]; - ScalarT omega = y_[1]; - ScalarT telec = y_[2]; - ScalarT ir = y_[3]; - ScalarT ii = y_[4]; - ScalarT pmech = y_[5]; - ScalarT ep = y_[6]; + // Set variable aliases for better reliability + const ScalarT delta = y_[0]; + const ScalarT omega = y_[1]; + const ScalarT telec = y_[2]; + const ScalarT ir = y_[3]; + const ScalarT ii = y_[4]; + const ScalarT pmech = y_[5]; + const ScalarT ep = y_[6]; - /* Read derivatives */ - ScalarT delta_dot = yp_[0]; - ScalarT omega_dot = yp_[1]; + // Set derivative aliases for better reliability + const ScalarT delta_dot = yp_[0]; + const ScalarT omega_dot = yp_[1]; - /* 6 GenClassical differential equations */ + // GenClassical differential equations f_[0] = delta_dot - (omega - 1.0) * (2.0 * M_PI * 60.0); - f_[1] = omega_dot - (1.0 / (2 * H_)) * ((pmech - D_ * (omega - 1.0)) / omega - telec); + f_[1] = omega_dot - (1.0 / (2.0 * H_)) * ((pmech - D_ * (omega - 1.0)) / omega - telec); - /* 11 GenClassical algebraic equations */ + // GenClassical algebraic equations f_[2] = telec - (1.0 / omega) * (G_ * ep * ep - ep * ((G_ * Vr() + B_ * Vi()) * cos(delta) + (-B_ * Vr() + G_ * Vi()) * sin(delta))); f_[3] = ir + G_ * Vr() + B_ * Vi() - ep * (G_ * cos(delta) + B_ * sin(delta)); @@ -167,6 +166,7 @@ namespace GridKit f_[5] = pmech - pmech_set_; f_[6] = ep - ep_set_; + // GenClassical contribution to bus algebraic equations Ir() += ir; Ii() += ii; diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md index 3fdb710a3..44ac0be75 100644 --- a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md +++ b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md @@ -17,8 +17,8 @@ $X_{dp}$ | [p.u.] | machine reactance parameter | ### Model Derived Parameters -- $G = \dfrac{R_a}{R_a^2 + X_{dp}^2}$ -- $B = \dfrac{-X_{dp}}{R_a^2 + X_{dp}^2}$ +- $G = \dfrac{R_a}{R_a^2 + X_{dp}^2} ~~~$ equivalent stator winding conductance +- $B = \dfrac{-X_{dp}}{R_a^2 + X_{dp}^2} ~~~$ equivalent stator winding susceptance