diff --git a/CMakeLists.txt b/CMakeLists.txt index 24f15839..0cbf1aaf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -21,7 +21,7 @@ set(PACKAGE_VERSION_PATCH "0") set(PACKAGE_VERSION "${PACKAGE_VERSION_MAJOR}.${PACKAGE_VERSION_MINOR}.${PACKAGE_VERSION_PATCH}") # TODO: Probably beter to set a debug interface target -# set(CMAKE_CXX_FLAGS_DEBUG "-Wall -O0 -g -DDEBUG") +set(CMAKE_CXX_FLAGS_DEBUG "-Wall -Wpedantic -Wconversion -Wextra -O0 -g -DDEBUG") set(CMAKE_CXX_STANDARD 17) diff --git a/examples/PhasorDynamics/CMakeLists.txt b/examples/PhasorDynamics/CMakeLists.txt index 58ded3a1..12d04958 100644 --- a/examples/PhasorDynamics/CMakeLists.txt +++ b/examples/PhasorDynamics/CMakeLists.txt @@ -1,2 +1,3 @@ add_subdirectory(Example1) add_subdirectory(Example2) +add_subdirectory(Gen2Example) diff --git a/examples/PhasorDynamics/Gen2Example/CMakeLists.txt b/examples/PhasorDynamics/Gen2Example/CMakeLists.txt new file mode 100644 index 00000000..c7252a7e --- /dev/null +++ b/examples/PhasorDynamics/Gen2Example/CMakeLists.txt @@ -0,0 +1,10 @@ +add_executable(gen2_example example.cpp) +target_link_libraries(gen2_example + GRIDKIT::phasor_dynamics_bus + GRIDKIT::phasor_dynamics_bus_fault + GRIDKIT::phasor_dynamics_branch + GRIDKIT::phasor_dynamics_classical_gen + GRIDKIT::solvers_dyn) +install(TARGETS gen2_example RUNTIME DESTINATION bin) + +add_test(NAME GenClassicalTest1 COMMAND $) diff --git a/examples/PhasorDynamics/Gen2Example/example.cpp b/examples/PhasorDynamics/Gen2Example/example.cpp new file mode 100644 index 00000000..dd87c82d --- /dev/null +++ b/examples/PhasorDynamics/Gen2Example/example.cpp @@ -0,0 +1,124 @@ +#include +#define _USE_MATH_DEFINES +#include +#include +#include +#include +#include + +#include "Model/PhasorDynamics/Branch/Branch.hpp" +#include "Model/PhasorDynamics/Bus/Bus.hpp" +#include "Model/PhasorDynamics/Bus/BusInfinite.hpp" +#include "Model/PhasorDynamics/BusFault/BusFault.hpp" +#include "Model/PhasorDynamics/Load/Load.hpp" +#include "Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp" +#include "Model/PhasorDynamics/SystemModel.hpp" +#include "Solver/Dynamic/Ida.hpp" + +#define _CRT_SECURE_NO_WARNINGS + +int main(int argc, char* argv[]) +{ + using namespace GridKit::PhasorDynamics; + using namespace AnalysisManager::Sundials; + + std::cout << "Example 1 version GENERATION 2" << std::endl; + + // bus voltages + double vr = 1.0; + double vi = 0.0; + + // branch parameters + double R = 0.0; // line series resistance + double X = 0.1; // line series reactance + double G = 0.0; // line shunt conductance + double B = 0.0; // line shunt charging + + // Generator parameters + double p0 = 1.0; // real power output + double q0 = 0.05013; // reactive power output + double H = 3.0; // Initia constant + double D = 0.0; // Damping coefficient + double Ra = 0.0; // Winding resistance + double Xdp = 0.2; // Machine reactance parameter + + SystemModel sys; + Bus bus1(vr, vi); + BusInfinite bus2(vr, vi); + Branch branch(&bus1, &bus2, R, X, G, B); + GenClassical gen(&bus1, 1, p0, q0, H, D, Ra, Xdp); + + /* Connect everything together */ + sys.addBus(&bus1); + sys.addBus(&bus2); + sys.addComponent(&branch); + sys.addComponent(&gen); + sys.allocate(); + sys.initialize(); + + std::vector> outputData; + + // callback for outputting solution + auto output_cb = [&](double t) + { + std::vector yval; + + yval.push_back(t); + for (auto val : sys.y()) + { + yval.push_back(val); + } + + outputData.push_back(yval); + }; + + output_cb(0); + + /* Set up simulation */ + Ida ida(&sys); + ida.configureSimulation(); + + /* Run simulation */ + double start = static_cast(clock()); + ida.initializeSimulation(0.0, false); + size_t nout = 200; + ida.runSimulation(10.0, nout, output_cb); + + // write solution to file if the user passes in a file name + if (argc >= 2) + { + std::ofstream outfile(argv[1]); + + if (!outfile) + { + std::cout << "ERROR writing to output file!" << std::endl; + } + + outfile << "t"; + for (int i = 0; i < sys.size(); ++i) + { + outfile << ",Y[" + std::to_string(i) + "]"; + } + outfile << "\n"; + + for (auto y_tstep : outputData) + { + outfile << y_tstep[0]; + for (int i = 1; i < y_tstep.size(); i++) + { + outfile << "," << y_tstep[i]; + } + outfile << "\n"; + } + + outfile.close(); + } + else + { + std::cout << "To print the solution, add a file name(eg:./gen2_example output.csv)" << std::endl; + } + + printf("Complete in %.4g seconds\n", (clock() - start) / CLOCKS_PER_SEC); + + return 0; +} diff --git a/src/Model/PhasorDynamics/SynchronousMachine/CMakeLists.txt b/src/Model/PhasorDynamics/SynchronousMachine/CMakeLists.txt index 212813aa..96c3b2f9 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 00000000..dea88aaa --- /dev/null +++ b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/CMakeLists.txt @@ -0,0 +1,5 @@ +gridkit_add_library(phasor_dynamics_classical_gen + SOURCES + GenClassical.cpp + OUTPUT_NAME + gridkit_phasor_dynamics_classical_gen) diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp new file mode 100644 index 00000000..b64104ae --- /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 * M_PI * 60); + 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 diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp new file mode 100644 index 00000000..651e3995 --- /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 new file mode 100644 index 00000000..d271eaf5 --- /dev/null +++ b/src/Model/PhasorDynamics/SynchronousMachine/GenClassical/README.md @@ -0,0 +1,40 @@ +Differential equations: + +```math +\begin{aligned} +\dot{\delta} &= \omega \cdot \omega _0 \\ +\dot{\omega} &= \frac{1}{2H}\bigg( \frac{P_{mech} - D\omega _0}{1 + \omega} - T_{elec}\bigg) +\end{aligned} +``` + +Algebraic Equations: + +```math + T_{elec} = \frac{1}{1+\omega}\bigg( g E_p^2 - E_p \bigg((gV_r - bV_i)cos\,\delta + (bV_r + gV_i)sin\,\delta \bigg)\bigg) +``` + +Network Interface Equations: + +```math +\begin{aligned} +I_r &= -gV_r + bV_i + E_p(g \cos \delta - b \sin \delta)\\ +I_i &= -gV_r - bV_i + E_p(b \cos \delta + g \sin \delta) +\end{aligned} +``` + +Intialization notes:
+To initialize the model, given $V_r$, $V_i$, $P$ and $Q$, we use the following equations: +

+```math +\begin{aligned} +I_r &= \frac{PV_r + QV_i}{V_r^2 + V_i^2} \\ +I_i &= \frac{PV_i - QV_r}{V_r^2 + V_i^2} \\ +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 &= atan2(E_i, E_r) \\ +\omega &= 0 \\ +T_{elec} &= gE_p^2 - E_p \bigg( \bigg(gV_r - bV_i \bigg) \cos \delta + \bigg(bV_r + gV_i \bigg)\sin \delta\bigg) \\ +P_{mech} &= T_{elec} +\end{aligned} +``` diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index 18876a99..8d33bec7 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -13,13 +13,16 @@ target_link_libraries(test_phasor_branch GRIDKIT::phasor_dynamics_branch add_executable(test_phasor_load runLoadTests.cpp) target_link_libraries(test_phasor_load GRIDKIT::phasor_dynamics_load - GRIDKIT::phasor_dynamics_bus) + GRIDKIT::phasor_dynamics_bus) 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_classical_gen runGenClassicalTests.cpp) +target_link_libraries(test_phasor_classical_gen GRIDKIT::phasor_dynamics_classical_gen + GRIDKIT::phasor_dynamics_bus) + add_executable(test_phasor_system runSystemTests.cpp) target_link_libraries(test_phasor_system GRIDKIT::phasor_dynamics_load GRIDKIT::phasor_dynamics_branch @@ -28,11 +31,13 @@ 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 $) install(TARGETS test_phasor_bus test_phasor_branch test_phasor_load + test_phasor_classical_gen test_phasor_genrou 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 00000000..6dafc7c7 --- /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 diff --git a/tests/UnitTests/PhasorDynamics/runGenClassicalTests.cpp b/tests/UnitTests/PhasorDynamics/runGenClassicalTests.cpp new file mode 100644 index 00000000..f50bc825 --- /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(); +}