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..b34180330 --- /dev/null +++ b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp @@ -0,0 +1,257 @@ +/** + * @file GenClassical.cpp + * @author Abdourahman Barry (abdourahman@vt.edu) + * @author Slaven Peles (peless@ornl.gov) + * @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 classical generator model. + * + */ + 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 + * + */ + 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() + { + 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; + } + + /** + * Initialization of the generator 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 - 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)); + + 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 (size_t i = 0; i < static_cast(size_); ++i) + yp_[i] = 0.0; + + return 0; + } + + /** + * \brief Identify differential variables. + */ + template + int GenClassical::tagDifferentiable() + { + + return 0; + } + + /** + * \brief Residual for the generator model. + * + */ + template + int GenClassical::evaluateResidual() + { + // 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]; + + // Set derivative aliases for better reliability + const ScalarT delta_dot = yp_[0]; + const ScalarT omega_dot = yp_[1]; + + // GenClassical differential equations + f_[0] = delta_dot - (omega - 1.0) * (2.0 * M_PI * 60.0); + f_[1] = omega_dot - (1.0 / (2.0 * H_)) * ((pmech - D_ * (omega - 1.0)) / omega - telec); + + // 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)); + f_[4] = ii - B_ * Vr() + G_ * Vi() - ep * (-B_ * cos(delta) + G_ * sin(delta)); + + f_[5] = pmech - pmech_set_; + f_[6] = ep - ep_set_; + + // GenClassical contribution to bus algebraic equations + Ir() += ir; + Ii() += ii; + + 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 diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp new file mode 100644 index 000000000..1e8a46550 --- /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 diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md index aa1dbe490..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
@@ -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 + gV_r + bV_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{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 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 them, 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/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..3c17c8ab1 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp @@ -0,0 +1,227 @@ +/** + * @file GenClassicalTests.hpp + * @author Slaven Peles (peless@ornl.gov) + * @author Abdourahman Barry (abdourahman@vt.edu) + * @brief Tests for classical generator model. + * + */ +#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; + static constexpr ScalarT tol_ = 10 * std::numeric_limits::epsilon(); + + 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 + + // 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(); + + // 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; + 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(); + + 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__); + } + + /** + * + * 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 + + // 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); + bus.allocate(); + bus.initialize(); + gen.allocate(); + gen.initialize(); + + 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__); + } + + /* + *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 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); + bus.allocate(); + bus.initialize(); + gen.allocate(); + gen.initialize(); + gen.evaluateResidual(); + std::vector res = gen.getResidual(); + + 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 GenClassicalTests + + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/runGenClassicalTests.cpp b/tests/UnitTests/PhasorDynamics/runGenClassicalTests.cpp new file mode 100644 index 000000000..21b7f45ba --- /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(); +}