From d1179997bc637de12769d6258c7461753e152a36 Mon Sep 17 00:00:00 2001 From: Paulm991 Date: Wed, 26 Jun 2024 11:34:06 -0400 Subject: [PATCH 1/5] Test and files for GENROU generator --- ComponentLib/CMakeLists.txt | 1 + ComponentLib/DynamicPhasor/CMakeLists.txt | 1 + .../SynchronousMachine/CMakeLists.txt | 1 + .../GENROUwS/CMakeLists.txt | 5 + .../SynchronousMachine/GENROUwS/GENROU.cpp | 536 ++++++++++++++++ .../SynchronousMachine/GENROUwS/GENROU.hpp | 607 ++++++++++++++++++ Examples/CMakeLists.txt | 1 + Examples/GenrouTest/CMakeLists.txt | 5 + Examples/GenrouTest/GenrouTest.cpp | 99 +++ ModelEvaluatorImpl.hpp | 2 +- Solver/Dynamic/Ida.cpp | 9 + SparseMatrix/COO_Matrix.hpp | 2 +- SystemModel.hpp | 39 +- 13 files changed, 1304 insertions(+), 4 deletions(-) create mode 100644 ComponentLib/DynamicPhasor/CMakeLists.txt create mode 100644 ComponentLib/DynamicPhasor/SynchronousMachine/CMakeLists.txt create mode 100644 ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/CMakeLists.txt create mode 100644 ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp create mode 100644 ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp create mode 100644 Examples/GenrouTest/CMakeLists.txt create mode 100644 Examples/GenrouTest/GenrouTest.cpp diff --git a/ComponentLib/CMakeLists.txt b/ComponentLib/CMakeLists.txt index a97c973..16b5a24 100644 --- a/ComponentLib/CMakeLists.txt +++ b/ComponentLib/CMakeLists.txt @@ -62,4 +62,5 @@ add_subdirectory(PowerFlow) add_subdirectory(PowerElectronics) +add_subdirectory(DynamicPhasor) diff --git a/ComponentLib/DynamicPhasor/CMakeLists.txt b/ComponentLib/DynamicPhasor/CMakeLists.txt new file mode 100644 index 0000000..b484860 --- /dev/null +++ b/ComponentLib/DynamicPhasor/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(SynchronousMachine) \ No newline at end of file diff --git a/ComponentLib/DynamicPhasor/SynchronousMachine/CMakeLists.txt b/ComponentLib/DynamicPhasor/SynchronousMachine/CMakeLists.txt new file mode 100644 index 0000000..279e0eb --- /dev/null +++ b/ComponentLib/DynamicPhasor/SynchronousMachine/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(GENROUwS) \ No newline at end of file diff --git a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/CMakeLists.txt b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/CMakeLists.txt new file mode 100644 index 0000000..8fd87d9 --- /dev/null +++ b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/CMakeLists.txt @@ -0,0 +1,5 @@ +gridkit_add_library(genrou + SOURCES + GENROU.cpp + OUTPUT_NAME + gridkit_genrou) \ No newline at end of file diff --git a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp new file mode 100644 index 0000000..2aea397 --- /dev/null +++ b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp @@ -0,0 +1,536 @@ +// Made by Paul Moon 6/7/2024 + + +#include +#include +#include +#include +#include "GENROU.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a simple generator model + * + * Arguments passed to ModelEvaluatorImpl: + * - Number of equations = 4 differential + 2 algebraic = 6 + * - Number of quadratures = 1 + * - Number of optimization parameters = 2 + */ +template +GENROU::GENROU(bus_type* bus, ScalarT P0, ScalarT Q0) + : ModelEvaluatorImpl(20, 1, 2), + H_(3.0), + D_(0.0), + Xq_(0.5), + Xd_(2.1), + Xqp_(0.5), + Xdp_(0.2), + Xqpp_(0.18), + Xdpp_(0.18), + Rs_(0.0), + Tq0p_(0.75), + Td0p_(7.0), + Tq0pp_(0.05), + Td0pp_(0.035), + Ef_(2.6904), + Pm_(1.0), + X_(0.22), + Xl_(0.15), + omega_s_(0.0), + omega0_(2.0*60.0*M_PI), + S1_(0.05), + S12_(0.2), + omega_up_(omega_s_ + 0.0001), + omega_lo_(omega_s_ - 0.0001), + c_(10000.0), + beta_(2), + P0_(P0), + Q0_(Q0), + bus_(bus) +{ +} + +template +GENROU::~GENROU() +{ +} + +/*! + * @brief This function will be used to allocate sparse Jacobian matrices. + * + */ +template +int GENROU::allocate() +{ + //std::cout << "Allocate GENROU..." << std::endl; + tag_.resize(size_); + + return 0; +} + +/** + * @brief Initialization of the generator model + * + * Initialization equations are derived from example 9.2 in Power System + * Modeling and Scripting, Federico Milano, Chapter 9, p. 225: + * \f{eqnarray*}{ + * &~& \omega_0 = 0, \\ + * &~& \delta_0 = \tan^{-1} \left(\frac{X_q P_0 - R_s Q_0}{V_0^2 + R_s P_0 + X_q Q_0} \right) + \theta_0, \\ + * &~& \phi_0 = \delta_0 - \theta_0 + \tan^{-1} \left( \frac{Q_0}{P_0} \right), \\ + * &~& I_{d0} = \frac{\sqrt{P_0^2 + Q_0^2}}{V_0} \sin(\phi_0), \\ + * &~& I_{q0} = \frac{\sqrt{P_0^2 + Q_0^2}}{V_0} \cos(\phi_0), \\ + * &~& E_{d0}' = V_0 \sin(\delta_0 - \theta_0) + R_s I_{d0} - X_q' I_{q0}, \\ + * &~& E_{q0}' = V_0 \cos(\delta_0 - \theta_0) + R_s I_{q0} + X_d' I_{d0} + * \f} + * + * The input from exciter and governor is set to the steady state value: + * \f{eqnarray*}{ + * &~& E_{f0} = E_{q0}' + (X_d - X_d') I_{d0}, \\ + * &~& P_{m0} = E_{d0}' I_{d0} + E_{q0}' I_{q0} + ( X_q' - X_d') I_{d0} I_{q0} + * \f} + * + */ +template +int GENROU::initialize() +{ + // std::cout << "Initialize GENROU..." << std::endl; + //const ScalarT Vr = V()*cos(theta()); + //const ScalarT Vi = V()*sin(theta()); + //std::cout << "Initial Vr & Vi: " << Vr << " " << Vi << std::endl; + + //const ScalarT Ir = (P0_*Vr + Q0_*Vi)/((Vr*Vr)+(Vi*Vi)); + //const ScalarT Ii = -(Q0_*Vr - P0_*Vi)/((Vr*Vr)+(Vi*Vi)); + //std::cout << "Initial Ir & Ii: " << Ir << " " << Ii << std::endl; + + // Compute initial guess for the generator voltage phase + //const ScalarT delta = atan((Xq_*P0_ - Rs_*Q0_) / (V()*V() + Rs_*P0_ + Xq_*Q0_)) + theta(); + + // Compute initial guess for the generator current phase + //const ScalarT phi = theta() - delta - atan(Q0_/P0_); + + //std::cout << "Initial delta & phi: " << delta << " " << phi << std::endl; + + // Compute initial gueses for generator currents and potentials in d-q frame + //const ScalarT Id = std::sqrt(P0_*P0_ + Q0_*Q0_)/V() * sin(phi); + //const ScalarT Iq = std::sqrt(P0_*P0_ + Q0_*Q0_)/V() * cos(phi); + //const ScalarT Id = 0.9909; + //const ScalarT Iq = 0.3553; + + //std::cout << "Initial Id & Iq: " << Id << " " << Iq << std::endl; + + // const ScalarT Vd = V()*sin(theta() - delta) + Rs_*Id + Xqp_*Iq; + // const ScalarT Vq = V()*cos(theta() - delta) + Rs_*Iq - Xdp_*Id; + + //std::cout << "Initial Vd & Vq: " << Vd << " " << Vq << std::endl; + + /*const ScalarT Vr = 1.0723; + const ScalarT Vi = 0.22; + const ScalarT Ir = 1; + const ScalarT Ii = -0.3286; + const ScalarT delta = 0.909; + const ScalarT Vd = 0.7107; + const ScalarT Vq = 0.8326; + const ScalarT Eqp = 1.1298; + const ScalarT Psidp = 0.9614; + const ScalarT Psiqp = 0.6645; + const ScalarT Edp = 0.533; + const ScalarT Psiqpp = -0.6111; + const ScalarT Psidpp = 1.110; + const ScalarT Psipp = 1.2671; + const ScalarT ksat = 1; + const ScalarT Telec = 1; + const ScalarT Pmech = 0.9999;*/ + + + //y0 = [delta_init; omega_init; E_p_q_init; psi_p_d_init; psi_p_q_init; E_p_d_init; psi_pp_q_init; psi_pp_d_init; psi_pp_init; ksat_init; Vd_init; Vq_init ; Telec_init; Id_init; Iq_init; Vr_init; Vi_init; Ir_init; Ii_init;Pmech_init] + /*y_[0] = delta; + y_[1] = omega_s_; + y_[2] = Eqp; + y_[3] = Psidp; + y_[4] = Psiqp; + y_[5] = Edp; + y_[6] = Psiqpp; + y_[7] = Psidpp; + y_[8] = Psipp; + y_[9] = ksat; + y_[10] = Vd; + y_[11] = Vq; + y_[12] = Telec; + y_[13] = Id; + y_[14] = Iq; + y_[15] = Vr; + y_[16] = Vi; + y_[17] = Ir; + y_[18] = Ii; + y_[19] = Pmech;*/ + y_[0] = 0.5273; + y_[1] = 0.0; + y_[2] = 1.1948; + y_[3] = 1.1554; + y_[4] = 0.2446; + y_[5] = 0.0; + y_[6] = -0.2236; + y_[7] = 1.1790; + y_[8] = 1.2000; + y_[9] = 0.0; + y_[10] = 0.3494; + y_[11] = 1.0373; + y_[12] = 1.0; + y_[13] = 0.7872; + y_[14] = 0.6989; + y_[15] = 1.0723; + y_[16] = 0.22; + y_[17] = 1.0; + y_[18] = -0.3286; + y_[19] = 1.0; + /*y_[0] = 0.9224; + y_[1] = 0.1765; + y_[2] = 1.2163; + y_[3] = 1.1554; + y_[4] = -0.3907; + y_[5] = 0.0; + y_[6] = 0.3572; + y_[7] = 0.2676; + y_[8] = 0.4464; + y_[9] = 2.5613; + y_[10] = -0.4203; + y_[11] = 0.3149; + y_[12] = 0.0; + y_[13] = 0.7872; + y_[14] = 1.0507; + y_[15] = 0.0797; + y_[16] = 0.2776; + y_[17] = 1.2620; + y_[18] = 0.3622; + y_[19] = 1.0;*/ + yp_[0] = 0.0; + yp_[1] = 0.0; + yp_[2] = 0.0; + yp_[3] = 0.0; + yp_[4] = 0.0; + yp_[5] = 0.0; + yp_[6] = 0.0; + yp_[7] = 0.0; + yp_[8] = 0.0; + yp_[9] = 0.0; + yp_[10] = 0.0; + yp_[11] = 0.0; + yp_[12] = 0.0; + yp_[13] = 0.0; + yp_[14] = 0.0; + yp_[15] = 0.0; + yp_[16] = 0.0; + yp_[17] = 0.0; + yp_[18] = 0.0; + yp_[19] = 0.0; + + /*for (int i = 0; i < 20; i++) { + std::cout << y_[i] << ", "; + } + std::cout << std::endl;*/ + + return 0; +} + +/** + * \brief Identify differential variables. + */ +template +int GENROU::tagDifferentiable() +{ + for (IdxT i=6; i < size_; ++i) + { + tag_[i] = false; + } + tag_[0] = true; + tag_[1] = true; + tag_[2] = true; + tag_[3] = true; + tag_[4] = true; + tag_[5] = true; + + + + return 0; +} + +/** + * @brief Computes residual vector for the generator model. + * + * Residual equations are given per model in Power System Modeling and + * Scripting, Federico Milano, Chapter 15, p. 334: + * \f{eqnarray*}{ + * f_0: &~& \dot{\delta} -\omega_b (\omega - \omega_s), \\ + * f_1: &~& 2H/\omega_s \dot{\omega} - L_m(P_m) + E_q' I_q + E_d' I_d + (X_q' - X_d')I_d I_q + D (\omega - \omega_s), \\ + * f_2: &~& T_{q0}' \dot{E}_d' + E_d' - (X_q - X_q')I_q, \\ + * f_3: &~& T_{d0}' \dot{E}_q' + E_q' + (X_d - X_d')I_d - E_f, \\ + * f_4: &~& R_s I_d - X_q' I_q + V \sin(\delta - \theta) - E_d', \\ + * f_5: &~& R_s I_q + X_d' I_d + V \cos(\delta - \theta) - E_q', + * \f} + * where \f$ \Omega_b \f$ is the synchronous frequency in [rad/s], and + * overdot denotes time derivative. + * + * Generator injection active and reactive power are + * \f{eqnarray*}{ + * P_g &=& E_d' I_d + E_q' I_q + (X_q' - X_d') I_d I_q - R_s (I_d^2 + I_q^2), \\ + * Q_q &=& E_q' I_d - E_d' I_q - X_q' I_q^2 - X_d' I_d^2, \\ + * \f} + * respectively. + * + * State variables are: + * \f$ y_0 = \omega \f$, \f$ y_1 = \delta \f$, \f$ y_2 = E_d' \f$, \f$ y_3 = E_q' \f$, + * \f$ y_4 = I_d \f$, \f$ y_5 = I_q \f$. + * + */ +template +int GENROU::evaluateResidual() +{ + double SA=std::log(S12_/S1_)/std::log(1.2); + double SB=S1_; + std::complex numerator(y_[15], y_[16]); // numerator is y(16) + 1i * y(17) + std::complex denominator(0.0, X_); // denominator is 1i * X (1i is imaginary unit in MATLAB) + std::complex division_result = numerator / denominator; + + double real_part = division_result.real(); // real part of (numerator / denominator) + double imag_part = division_result.imag(); // imaginary part of (numerator / denominator) + std::complex i(0.0, 1.0); // imaginary unit + + + //std::cout << "SA: " << SA << " SB: " << SB << std::endl; + // std::cout << "Evaluate residual for GENROU..." << std::endl; + f_[0] = -yp_[0] + omega0_*y_[1]; //f1 = -yp(1)+y(2)*w0_1; + f_[1] = -yp_[1] + (1/(2*H_))*((Pm_ - D_*y_[1])/(1 + y_[1]) - y_[12]); //f2 = -yp(2) + (1/(2*H_1))*((Pmech_init - D_1*y(2))/(1 + y(2)) - y(13)); + f_[2] = -yp_[2]*Td0p_ + Ef_ - (y_[2] + (Xd_-Xdp_)*y_[13]+((Xdp_-Xdpp_)/((Xdp_-Xl_)*(Xdp_-Xl_)))*(y_[2]-y_[3]-(Xdp_-Xl_)*y_[13])) + y_[7]*y_[9]; //f3 = -yp(3)*T_p_d0_1 + Efd - (y(3) + (Xd_1-X_p_d_1)*(y(14) + ((X_p_d_1-X_pp_d_1)/(X_p_d_1 - Xl_1)^2)*(y(3)-y(4)-(X_p_d_1-Xl_1)*y(14))) + y(8)*y(10)); + f_[3] = -yp_[3]+(1/Td0pp_)*(y_[2]-y_[3]-(Xdp_-Xl_)*y_[13]); //f4 = -yp(4)+(1/T_pp_d0_1)*(y(3)-y(4)-(X_p_d_1-Xl_1)*y(14)); + f_[4] = -yp_[4]+(1/Tq0pp_)*(y_[5]-y_[4]+(Xqp_-Xl_)*y_[14]); //f5 = -yp(5)+(1/T_pp_q0_1)*(y(6)-y(5)+(X_p_q_1-Xl_1)*y(15)); + f_[5] = -yp_[5] + (1/Tq0p_)*(-y_[5]+(((Xq_-Xl_)/(Xd_-Xl_))*y_[6]*y_[9] + (Xq_-Xqp_)*(y_[14]-((Xqp_-Xqpp_)/((Xqp_-Xl_)*(Xqp_-Xl_))*y_[5]+y_[14])*(Xqp_-Xl_)-y_[4]))); //f6 = -yp(6) + (1/T_p_q0_1)*(-y(6)+((Xq_1-Xl_1)/(Xd_1-Xl_1))*y(7)*y(10) +(Xq_1-X_p_q_1)*(y(15)-((X_p_q_1-X_pp_q_1)/(X_p_q_1-Xl_1)^2)*(y(6)+y(15)*(X_p_q_1-Xl_1)-y(5)))); + + f_[6] = -y_[6] - y_[4]*((Xqp_-Xqpp_)/(Xqp_-Xl_))-y_[5]*((Xqpp_-Xl_)/(Xqp_-Xl_)); //f7 = -y(7) - y(5)*((X_p_q_1-X_pp_q_1)/(X_p_q_1 - Xl_1))- y(6)*((X_pp_q_1 - Xl_1)/(X_p_q_1 - Xl_1)); + f_[7] = -y_[7] - y_[3]*((Xdp_-Xdpp_)/(Xdp_-Xl_))+y_[2]*((Xdpp_-Xl_)/(Xdp_-Xl_)); //f8 = -y(8) - y(4)*((X_p_d_1-X_pp_d_1)/(X_p_d_1 - Xl_1)) + y(3)*((X_pp_d_1 - Xl_1)/(X_p_d_1 - Xl_1)); + f_[8] = -y_[8] + sqrt((y_[7]*y_[7])+(y_[6]*y_[6]));//f9 = -y(9) + sqrt(y(8)^2+y(7)^2); + f_[9] = -y_[9] + SB*((y_[8]-SA)*(y_[8]-SA)); //f10 = -y(10) + SB * (y(9) - SA)^2; + f_[10] = -y_[10] - y_[6] * (1+y_[1]); //f11 = -y(11) - y(7) * (1+y(2)); + f_[11] = -y_[11] + y_[7] * (1+y_[1]); //f12 = -y(12) + y(8) * (1+y(2)); + f_[12] = -y_[12] + y_[14] * (y_[7] - y_[13] * Xdpp_) - y_[13] * (y_[6] - y_[14] * Xdpp_); //f13 = -y(13) + y(15) * (y(8) - y(14) * X_pp_d_1) - y(14) * (y(7) - y(15) * X_pp_d_1); + f_[13] = -y_[13] + y_[17] * sin(y_[0]) - y_[18] * cos(y_[0]); //f14 = -y(14) + y(18) * sin(y(1)) - y(19) * cos(y(1)); % prev version: -y(14) + Ireal * sin(y(1)) - Iimag * cos(y(1)); %update Ireal and Iimag along with state variables each iteration + f_[14] = -y_[14] + y_[17] * cos(y_[0]) + y_[18] * sin(y_[0]); //f15 = -y(15) + y(18) * cos(y(1)) + y(19) * sin(y(1)); % prev: -y(15) + Ireal * cos(y(1)) + Iimag * sin(y(1)); + f_[15] = -y_[10] + y_[15] * sin(y_[0]) - y_[16]*cos(y_[0]) + y_[13]*Rs_ - y_[14]*Xqpp_; //f16 = -y(11) + y(16)*sin(y(1)) - y(17)*cos(y(1))+y(14)*Ra_1 - y(15)*X_pp_q_1; + f_[16] = -y_[11] + y_[15] * cos(y_[0]) + y_[16]*sin(y_[0]) + y_[13]*Xqpp_ + y_[14]*Rs_; //f17 = -y(12) + y(16)*cos(y(1)) + y(17)*sin(y(1)) + y(14)*X_pp_q_1 + y(15)*Ra_1; + f_[17] = -y_[17] + std::real((y_[15]+i*y_[16])/(i*X_));// y_[16]/X_; //f18 = -y(18)+real((y(16) + 1i * y(17))/(1i*X)); + f_[18] = -y_[18] + std::imag((y_[15]+i*y_[16])/(i*X_));;// y_[15]/X_; //f19 = -y(19)+imag((y(16) + 1i * y(17))/(1i*X)); + f_[19] = -Pm_ + y_[19]; //f20 = -Pmech_init + y(20); + + return 0; +} + +template +int GENROU::evaluateJacobian() +{ + //std::cout << "Evaluate Jacobian for GENROU..." << std::endl; + J_.zeroMatrix(); + double SA=log(S12_/S1_)/log(1.2); + double SB=S1_; + + // These variables are for ease of reading + double df2_domega= y_[19]/((1+(y_[1])*(y_[1]))*2*H_) - D_/((1+y_[1])*2*H_) + D_*y_[1]/((1+(y_[1])*(y_[1]))*2*H_); + + double df3_dId = -(Xd_ - Xdp_) * (1 - (Xdp_ - Xl_) * ((Xdp_ - Xdpp_) / ((Xdp_ - Xl_)*(Xdp_ - Xl_)))); + double df3_dksat = -y_[7]; + double df3_dpsippd = -y_[9]; + double df3_dpsipd = (Xd_ - Xdp_)*(Xdp_ - Xdpp_)/((Xdp_ - Xl_)*(Xdp_ - Xl_)); + double df3_dEpq = -1 - (Xd_ - Xdp_) * ((Xdp_ - Xdpp_)/((Xdp_ - Xl_)*(Xdp_ - Xl_))); + + double df6_dpsippq = y_[9]*((Xq_-Xl_)/(Xd_-Xl_))/Tq0p_; + double df6_dEpd = -(1/Tq0p_) * (1 + (Xq_-Xqp_)*((Xqp_ - Xqpp_)/((Xqp_ - Xl_)*(Xqp_ - Xl_)))); + double df6_dpsipq = (Xq_-Xqp_)*((Xqp_ - Xqpp_)/((Xqp_ - Xl_)*(Xqp_ - Xl_)))/Tq0p_; + double df6_dksat = y_[6]*((Xq_-Xl_)/(Xd_-Xl_))/Tq0p_; + double df6_dIq = (Xq_-Xqp_)/Tq0p_ - (Xqp_-Xl_)*(Xq_-Xqp_)*((Xqp_ - Xqpp_)/((Xqp_ - Xl_)*(Xqp_ - Xl_)))/Tq0p_; + + double df7_dpsipq = -((Xqp_ - Xqpp_) / (Xqp_ - Xl_)); + double df7_dEpd = - ((Xqpp_ - Xl_)/(Xqp_ - Xl_)); + + double df8_dpsipd = -(Xdp_-Xdpp_)/(Xdp_ - Xl_); + double df8_dEpq = (Xdpp_ - Xl_)/(Xdp_ - Xl_); + + double dTelec_dId = -Xdpp_*y_[14] - y_[6] + y_[14]*Xdpp_; + double dTelec_dIq = y_[7] - y_[13]*Xdpp_ + Xdpp_*y_[13]; + + //Create dF/dy + std::vector rcord{0, + 1, 1, + 2, 2, 2, 2, 2, + 3, 3, 3, + 4, 4, 4, + 5, 5, 5, 5, 5, + 6, 6, 6, + 7, 7, 7, + 8, 8, 8, + 9, 9, + 10, 10, 10, + 11, 11, 11, + 12, 12, 12, 12, 12, + 13, 13, 13, 13, + 14, 14, 14, 14, + 15, 15, 15, 15, 15, 15, + 16, 16, 16, 16, 16, 16, + 17, 17, + 18, 18, + 19}; + std::vector ccord{1, + 1, 12, + 2, 3, 7, 9, 13, + 2, 3, 13, + 4, 5, 14, + 4, 5, 6, 9, 14, + 4, 5, 6, + 2, 3, 7, + 6, 7, 8, + 8, 9, + 1, 6, 10, + 1, 7, 11, + 6, 7, 12, 13, 14, + 0, 13, 17, 18, + 0, 14, 17, 18, + 0, 10, 13, 14, 15, 16, + 0, 11, 13, 14, 15, 16, + 16, 17, + 15, 18, + 19}; + std::vector vals{omega0_, + -df2_domega, (-1/(2*H_)), + df3_dEpq, df3_dpsipd, df3_dpsippd, df3_dksat, df3_dId, + 1/Td0pp_, -1/Td0pp_, (Xl_-Xdp_)/Td0pp_, + -1/Tq0pp_, 1/Tq0pp_, (Xqp_-Xl_)/Tq0pp_, + df6_dpsipq, df6_dEpd, df6_dpsippq, df6_dksat, df6_dIq, + df7_dpsipq, df7_dEpd, -1, + df8_dEpq, df8_dpsipd, -1, + y_[6]/sqrt((y_[7]*y_[7])+(y_[6]*y_[6])), y_[7]/sqrt((y_[7]*y_[7])+(y_[6]*y_[6])), -1, + 2*SB*(y_[8] - SA), -1, + -y_[6], -(1+y_[1]), -1, + y_[7], 1+y_[1], -1, + -y_[13], y_[14], -1, dTelec_dId, dTelec_dIq, + y_[17]*cos(y_[0])+y_[18]*sin(y_[0]), -1, sin(y_[0]), -cos(y_[0]), + -y_[17]*sin(y_[0])+y_[18]*cos(y_[0]), -1, cos(y_[0]), sin(y_[0]), + y_[15]*cos(y_[0])+y_[16]*sin(y_[0]), -1, Rs_, -Xqpp_, sin(y_[0]), -cos(y_[0]), + -y_[15]*sin(y_[0])+y_[16]*cos(y_[0]), -1, Xqpp_, Rs_, cos(y_[0]), sin(y_[0]), + 1.0/X_, -1, + -1.0/X_, -1, + 1}; + J_.setValues(rcord, ccord, vals); + + //Create dF/dy' + std::vector rcordder{0,1,2,3,4,5}; + std::vector ccordder{0,1,2,3,4,5}; + std::vector valsder{-1,-1,-Td0p_,-1,-1,-1}; + COO_Matrix Jacder = COO_Matrix(rcordder, ccordder, valsder,20,20); + + //Perform dF/dy + \alpha dF/dy' + J_.axpy(alpha_, Jacder); + return 0; +} + +template +int GENROU::evaluateIntegrand() +{ + // std::cout << "Evaluate Integrand for GENROU..." << std::endl; + return 0; +} + +template +int GENROU::initializeAdjoint() +{ + //std::cout << "Initialize adjoint for GENROU..." << std::endl; + return 0; +} + + +/** + * @brief Computes adjoint residual vector for the generator model. + * + * Adjoint residual equations are given as: + * \f{eqnarray*}{ + * f_{B0}: &~& \dot{y}_{B0} - y_{B4} V \cos(\delta - \theta) + y_{B5} V \sin(\delta - \theta), \\ + * f_{B1}: &~& 2H/\omega_s \dot{y}_{B1} + y_{B0} \omega_b - y_{B1} D + y_{B9} (1 - T_2/T_1) - y_{B10} K T_2/T_1 + g_{\omega}(\omega), \\ + * f_{B2}: &~& T_{q0}' \dot{y}_{B2} - y_{B1} I_d - y_{B2} + y_{B4} + y_{B6} I_d - y_{B7} I_q, \\ + * f_{B3}: &~& T_{d0}' \dot{y}_{B3} - y_{B1} I_q - y_{B3} + y_{B5} + y_{B6} I_q + y_{B7} I_d, \\ + * f_{B4}: &~& -y_{B1} (E_d' + (-X_d'+X_q') I_q) - y_{B3} (X_d - X_d') - y_{B4} R_s - y_{B5} X_d' + y_{B6} (E_d' + (X_q' - X_d') I_q - 2 R_s I_d) + y_{B7} (E_q' - 2 X_d' I_d), \\ + * f_{B5}: &~& -y_{B1} (E_q' + (-X_d'+X_q') I_d) + y_{B2} (X_q - X_q') + y_{B4} X_q' - y_{B5} R_s + y_{B6} (E_q' + (X_q' - X_d') I_d - 2 R_s I_q) - y_{B7} (E_d' + 2 X_q' I_q). \\ + * \f} + * + */ +template +int GENROU::evaluateAdjointResidual() +{ + return 0; +} + +// template +// int GENROU::evaluateAdjointJacobian() +// { +// std::cout << "Evaluate adjoint Jacobian for GENROU..." << std::endl; +// std::cout << "Adjoint Jacobian evaluation not implemented!" << std::endl; +// return 0; +// } + +template +int GENROU::evaluateAdjointIntegrand() +{ + return 0; +} + + +// +// Private functions +// + +/** + * Generator active power Pg. + * + * \f[ P_g = E_q' I_q + E_d' I_d + (X_q' - X_d') I_q I_d - R_a (I_d^2 + I_q^2) \f] + * + */ +template +ScalarT GENROU::Pg() +{ + return 1; +} + +/** + * Generator reactive power Qg. + * + * \f[ Q_g = E_q' I_d - E_d' I_q - X_d' I_d^2 - X_q' I_q^2 \f] + */ +template +ScalarT GENROU::Qg() +{ + return 1; +} + +/** + * Frequency penalty is used as the objective function for the generator model. + */ +template +ScalarT GENROU::frequencyPenalty(ScalarT omega) +{ + return 0; +} + +/** + * Derivative of frequency penalty cannot be written in terms of min/max functions. + * Need to expand conditional statements instead. + */ +template +ScalarT GENROU::frequencyPenaltyDer(ScalarT omega) +{ + return 0.0; +} + + + +// Available template instantiations +template class GENROU; +template class GENROU; + + +} // namespace ModelLib \ No newline at end of file diff --git a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp new file mode 100644 index 0000000..6d266ee --- /dev/null +++ b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp @@ -0,0 +1,607 @@ +// Made by Paul Moon 6/7/2024 + +#ifndef _GENROU_H_ +#define _GENROU_H_ + +#include + +namespace ModelLib +{ + template class BaseBus; +} + +namespace ModelLib +{ + /*! + * @brief Implementation of a fourth order generator model. + * + */ + template + class GENROU : public ModelEvaluatorImpl + { + using ModelEvaluatorImpl::size_; + using ModelEvaluatorImpl::nnz_; + using ModelEvaluatorImpl::time_; + using ModelEvaluatorImpl::alpha_; + using ModelEvaluatorImpl::y_; + using ModelEvaluatorImpl::yp_; + using ModelEvaluatorImpl::tag_; + using ModelEvaluatorImpl::f_; + using ModelEvaluatorImpl::g_; + using ModelEvaluatorImpl::yB_; + using ModelEvaluatorImpl::ypB_; + using ModelEvaluatorImpl::fB_; + using ModelEvaluatorImpl::gB_; + using ModelEvaluatorImpl::param_; + using ModelEvaluatorImpl::param_up_; + using ModelEvaluatorImpl::param_lo_; + using ModelEvaluatorImpl::J_; + + typedef typename ModelEvaluatorImpl::real_type real_type; + typedef BaseBus bus_type; + + public: + GENROU(BaseBus* bus, ScalarT P0 = 0.999, ScalarT Q0 = 0.5723); + virtual ~GENROU(); + + int allocate(); + int initialize(); + int tagDifferentiable(); + int evaluateResidual(); + int evaluateJacobian(); + int evaluateIntegrand(); + + int initializeAdjoint(); + int evaluateAdjointResidual(); + //int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); + + void updateTime(real_type t, real_type a) + { + time_ = t; + alpha_ = a; + } + + // Inline accesor functions + ScalarT& V() + { + return bus_->V(); + } + + const ScalarT& V() const + { + return bus_->V(); + } + + ScalarT& theta() + { + return bus_->theta(); + } + + const ScalarT& theta() const + { + return bus_->theta(); + } + + ScalarT& P() + { + return bus_->P(); + } + + const ScalarT& P() const + { + return bus_->P(); + } + + ScalarT& Q() + { + return bus_->Q(); + } + + const ScalarT& Q() const + { + return bus_->Q(); + } + + + private: + const ScalarT& Pm() const + { + return param_[0]; + } + + const ScalarT& Ef() const + { + return param_[1]; + } + + ScalarT Pg(); + ScalarT Qg(); + ScalarT frequencyPenalty(ScalarT omega); + ScalarT frequencyPenaltyDer(ScalarT omega); + + private: + // + // Private inlined accessor methods + // + + const ScalarT dotDelta() const + { + return yp_[0]; + } + + const ScalarT dotOmega() const + { + return yp_[1]; + } + + const ScalarT dotEqp() const + { + return yp_[2]; + } + + const ScalarT dotPsidp() const + { + return yp_[3]; + } + + const ScalarT dotPsiqp() const + { + return yp_[4]; + } + + const ScalarT dotEdp() const + { + return yp_[5]; + } + + const ScalarT dotPsiqpp() const + { + return yp_[6]; + } + + const ScalarT dotPsidpp() const + { + return yp_[7]; + } + + const ScalarT dotPsipp() const + { + return yp_[8]; + } + + const ScalarT dotksat() const + { + return yp_[9]; + } + + const ScalarT dotVd() const + { + return yp_[10]; + } + + const ScalarT dotVq() const + { + return yp_[11]; + } + + const ScalarT dotTelec() const + { + return yp_[12]; + } + + const ScalarT dotId() const + { + return yp_[13]; + } + + const ScalarT dotIq() const + { + return yp_[14]; + } + + const ScalarT dotVr() const + { + return yp_[15]; + } + + const ScalarT dotVi() const + { + return yp_[16]; + } + + const ScalarT dotIr() const + { + return yp_[17]; + } + + const ScalarT dotIi() const + { + return yp_[18]; + } + + const ScalarT dotPmech() const + { + return yp_[19]; + } + + // + // y_ values + // + + const ScalarT delta() const + { + return y_[0]; + } + + const ScalarT omega() const + { + return y_[1]; + } + + const ScalarT Eqp() const + { + return y_[2]; + } + + const ScalarT Psidp() const + { + return y_[3]; + } + + const ScalarT Psiqp() const + { + return y_[4]; + } + + const ScalarT Edp() const + { + return y_[5]; + } + + const ScalarT Psiqpp() const + { + return y_[6]; + } + + const ScalarT Psidpp() const + { + return y_[7]; + } + + const ScalarT Psipp() const + { + return y_[8]; + } + + const ScalarT ksat() const + { + return y_[9]; + } + + const ScalarT Vd() const + { + return y_[10]; + } + + const ScalarT Vq() const + { + return y_[11]; + } + + const ScalarT Telec() const + { + return y_[12]; + } + + const ScalarT Id() const + { + return y_[13]; + } + + const ScalarT Iq() const + { + return y_[14]; + } + + const ScalarT Vr() const + { + return y_[15]; + } + + const ScalarT Vi() const + { + return y_[16]; + } + + const ScalarT Ir() const + { + return y_[17]; + } + + const ScalarT Ii() const + { + return y_[18]; + } + + const ScalarT Pmech() const + { + return y_[19]; + } + + private: + real_type omega0_;///< Nominal frequency (assume 2pi60 rad/s) + real_type H_; ///< Inertia constant [s] + real_type D_; ///< Damping constant [pu] + real_type Rs_; ///< stator resistance [pu] + real_type Xl_; ///< stator leakage reactance [pu] + real_type Xd_; ///< direct axis synchronous reactance [pu] + real_type Xdp_; ///< direct axis transient reactance [pu] + real_type Xdpp_; ///< direct axis sub-transient reactance [pu] + real_type Xq_; ///< quadrature axis synchronous reactance [pu] + real_type Xqp_; ///< quadrature axis transient reactance [pu] + real_type Xqpp_; ///< quadrature axis sub-transient reactance [pu] + real_type Td0p_; ///< open circuit direct axis transient time constant [s] + real_type Td0pp_; ///< open circuit direct axis sub-transient time constant [s] + real_type Tq0p_; ///< open circuit quadrature axis transient time constant [s] + real_type Tq0pp_; ///< open circuit quadrature axis sub-transient time constant [s] + real_type S1_; ///< saturation at 1.0 pu flux + real_type S12_; ///< saturation at 1.2 pu flux + real_type X_; + + real_type Ef_; + real_type Pm_; + real_type omega_s_; + real_type omega_up_; + real_type omega_lo_; + real_type c_; + real_type beta_; + + ScalarT P0_; + ScalarT Q0_; + + bus_type* bus_; + }; + +} // namespace ModelLib + + +#endif // _GENROU_H_ + + +/* Paul Moon 6/6/2024 */ + +/*#include "ModelEvaluatorImpl.hpp" + +namespace ModelLib +{ + template class BaseBus; +} + +namespace ModelLib +{ + /** + * @brief Generator + * + * @tparam ScalarT - Scalar type + * @tparam IdxT - Matrix and vector index type + */ + /*template + class GENROU : public ModelEvaluatorImpl + { + using ModelEvaluatorImpl::size_; + using ModelEvaluatorImpl::nnz_; + using ModelEvaluatorImpl::time_; + using ModelEvaluatorImpl::alpha_; + using ModelEvaluatorImpl::y_; + using ModelEvaluatorImpl::yp_; + using ModelEvaluatorImpl::tag_; + using ModelEvaluatorImpl::f_; + using ModelEvaluatorImpl::g_; + using ModelEvaluatorImpl::yB_; + using ModelEvaluatorImpl::ypB_; + using ModelEvaluatorImpl::fB_; + using ModelEvaluatorImpl::gB_; + using ModelEvaluatorImpl::param_; + using ModelEvaluatorImpl::param_up_; + using ModelEvaluatorImpl::param_lo_; + + typedef typename ModelEvaluatorImpl::real_type real_type; + typedef BaseBus bus_type; + + public: + GENROU(BaseBus* bus, ScalarT P0 = 1.0, ScalarT Q0 = 0.0); + virtual ~GENROU(); + + int allocate(); + int initialize(); + int tagDifferentiable(); + int evaluateResidual(); + int evaluateJacobian(); + int evaluateIntegrand(); + + int initializeAdjoint(); + int evaluateAdjointResidual(); + //int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); + + void updateTime(real_type t, real_type a) + { + time_ = t; + alpha_ = a; + } + + // Inline accesor functions + ScalarT& V() + { + return bus_->V(); + } + + const ScalarT& V() const + { + return bus_->V(); + } + + ScalarT& theta() + { + return bus_->theta(); + } + + const ScalarT& theta() const + { + return bus_->theta(); + } + + ScalarT& P() + { + return bus_->P(); + } + + const ScalarT& P() const + { + return bus_->P(); + } + + ScalarT& Q() + { + return bus_->Q(); + } + + const ScalarT& Q() const + { + return bus_->Q(); + } + + + private: + const ScalarT& Pm() const + { + return param_[0]; + } + + const ScalarT& Ef() const + { + return param_[1]; + } + + ScalarT Pg(); + ScalarT Qg(); + ScalarT frequencyPenalty(ScalarT omega); + ScalarT frequencyPenaltyDer(ScalarT omega); + + private: + // + // Private inlined accessor methods + // + + const ScalarT dotDelta() const + { + return yp_[0]; + } + + const ScalarT dotOmega() const + { + return yp_[1]; + } + + const ScalarT dotEdp() const + { + return yp_[2]; + } + + const ScalarT dotEqp() const + { + return yp_[3]; + } + + const ScalarT delta() const + { + return y_[0]; + } + + const ScalarT omega() const + { + return y_[1]; + } + + const ScalarT Edp() const + { + return y_[2]; + } + + const ScalarT Eqp() const + { + return y_[3]; + } + + const ScalarT Id() const + { + return y_[4]; + } + + const ScalarT Iq() const + { + return y_[5]; + } + + private: + inline ScalarT frequencyPenalty(ScalarT omega); + inline ScalarT frequencyPenaltyDer(ScalarT omega); + + private: + real_type omega0_;///< Nominal frequency (assume 2pi60 rad/s) + real_type H_; ///< Inertia constant [s] + real_type D_; ///< Damping constant [pu] + real_type Rs_; ///< stator resistance [pu] + real_type Xl_; ///< stator leakage reactance [pu] + real_type Xd_; ///< direct axis synchronous reactance [pu] + real_type Xdp_; ///< direct axis transient reactance [pu] + real_type Xdpp_; ///< direct axis sub-transient reactance [pu] + real_type Xq_; ///< quadrature axis synchronous reactance [pu] + real_type Xqp_; ///< quadrature axis transient reactance [pu] + real_type Xqpp_; ///< quadrature axis sub-transient reactance [pu] + real_type Td0p_; ///< open circuit direct axis transient time constant [s] + real_type Td0pp_; ///< open circuit direct axis sub-transient time constant [s] + real_type Tq0p_; ///< open circuit quadrature axis transient time constant [s] + real_type Tq0pp_; ///< open circuit quadrature axis sub-transient time constant [s] + real_type S1_; ///< saturation at 1.0 pu flux + real_type S12_; ///< saturation at 1.2 pu flux + + real_type delta_; ///< Machine internal angle + real_type omega_; ///< + real_type Eqp_; ///< + real_type psidp_; ///< + real_type psiqp_; ///< + real_type Edp_; ///< + real_type psi_; ///< + real_type ksat_; ///< Saturation coefficient + real_type Vd_; ///< + real_type Vq_; ///< + real_type Telec_; ///< + real_type Ir_; ///< + real_type Ii_; ///< + real_type Vr_; ///< + real_type Vi_; ///< + real_type Id_; ///< + real_type Iq_; ///< + real_type Pmech_; ///< + real_type Efd_; ///< + + + + + real_type Ef_; // + real_type Pm_; + real_type omega_s_; + real_type omega_b_; + real_type omega_up_; + real_type omega_lo_; + real_type c_; + real_type beta_; + + ScalarT P0_; + ScalarT Q0_; + + bus_type* bus_; + }; +}*/ \ No newline at end of file diff --git a/Examples/CMakeLists.txt b/Examples/CMakeLists.txt index 67e0da1..2283111 100644 --- a/Examples/CMakeLists.txt +++ b/Examples/CMakeLists.txt @@ -73,5 +73,6 @@ if(TARGET SUNDIALS::idas) add_subdirectory(GenConstLoad) add_subdirectory(GenInfiniteBus) add_subdirectory(ParameterEstimation) + add_subdirectory(GenrouTest) endif() endif() diff --git a/Examples/GenrouTest/CMakeLists.txt b/Examples/GenrouTest/CMakeLists.txt new file mode 100644 index 0000000..b6435e8 --- /dev/null +++ b/Examples/GenrouTest/CMakeLists.txt @@ -0,0 +1,5 @@ +add_executable(genroutest GenrouTest.cpp) +target_link_libraries(genroutest GRIDKIT::bus GRIDKIT::genrou GRIDKIT::solvers_opt GRIDKIT::solvers_dyn) +install(TARGETS genroutest DESTINATION bin) + +add_test(NAME GenrouTest COMMAND $) \ No newline at end of file diff --git a/Examples/GenrouTest/GenrouTest.cpp b/Examples/GenrouTest/GenrouTest.cpp new file mode 100644 index 0000000..62d7c3a --- /dev/null +++ b/Examples/GenrouTest/GenrouTest.cpp @@ -0,0 +1,99 @@ +// Made by Paul Moon 6/7/2024 + + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + + +int main() +{ + using namespace ModelLib; + using namespace AnalysisManager::Sundials; + using namespace AnalysisManager; + using namespace GridKit::Testing; + + + // Create an infinite bus + BaseBus* bus = new BusSlack(1, 0); // roughly 1.072+0.22j + + // Attach a generator to that bus + GENROU* gen = new GENROU(bus); + + // Create a system model + SystemModel* model = new SystemModel(); + model->addBus(bus); + model->addComponent(gen); + + // allocate model components + model->allocate(); + + std::cout << "Size: " << model->y().size() << std::endl; + + model->initialize(); + model->evaluateResidual(); + + std::cout << "Verify Intial Residual is Zero: {"; + for (double i : model->getResidual()) + { + std::cout << i << ", "; + } + std::cout << "}\n"; + + + model->updateTime(0.0, 1.0); + model->evaluateJacobian(); + std::cout << "Intial Jacobian with alpha = 1:\n"; + model->getJacobian().printMatrix(); + + + // Create numerical integrator and configure it for the generator model + AnalysisManager::Sundials::Ida* idas = new AnalysisManager::Sundials::Ida(model); + + double t_init = 0.0; + double t_final = 1.0; + double t_timestep = 0.0001; + + idas->configureSimulation(); + idas->getDefaultInitialCondition(); + idas->initializeSimulation(t_init); + idas->runSimulation(t_final, t_final/t_timestep); + + /*for (double i = t_init; i <= t_final; i += t_timestep) { + //idas->configureSimulation(); + //idas->getDefaultInitialCondition(); + //idas->initializeSimulation(i); + + idas->runSimulation(i+t_timestep); + + std::vector& yfinial = model->y(); + + outFile << i; + for (size_t j = 0; j < yfinial.size(); j++) { + outFile << " " << yfinial[j]; + } + outFile << "\n"; + + std::cout << "Final Vector y\n"; + for (size_t j = 0; j < yfinial.size(); j++) + { + std::cout << yfinial[j] << "\n"; + } + }*/ + + delete idas; + delete model; + return 0; +} diff --git a/ModelEvaluatorImpl.hpp b/ModelEvaluatorImpl.hpp index 0f8969f..e98c526 100644 --- a/ModelEvaluatorImpl.hpp +++ b/ModelEvaluatorImpl.hpp @@ -113,7 +113,7 @@ namespace ModelLib virtual bool hasJacobian() { - return false; + return true; } virtual IdxT size_quad() diff --git a/Solver/Dynamic/Ida.cpp b/Solver/Dynamic/Ida.cpp index 637d110..9263604 100644 --- a/Solver/Dynamic/Ida.cpp +++ b/Solver/Dynamic/Ida.cpp @@ -60,6 +60,7 @@ #include #include +#include #include /* access to IDADls interface */ #include @@ -269,6 +270,7 @@ namespace Sundials template int Ida::runSimulation(real_type tf, int nout) { + std::ofstream outFile("genrououtputjac.txt"); int retval = 0; int iout = 0; real_type tret; @@ -283,6 +285,12 @@ namespace Sundials retval = IDASolve(solver_, tout, &tret, yy_, yp_, IDA_NORMAL); checkOutput(retval, "IDASolve"); //printOutput(tout); + realtype *yval = N_VGetArrayPointer_Serial(yy_); + outFile << tout; + for (size_t j = 0; j < model_->size(); j++) { + outFile << " " << yval[j]; + } + outFile << "\n"; if (retval == IDA_SUCCESS) { @@ -290,6 +298,7 @@ namespace Sundials tout += dt; } } + outFile.close(); //std::cout << "\n"; return retval; } diff --git a/SparseMatrix/COO_Matrix.hpp b/SparseMatrix/COO_Matrix.hpp index 9036085..107c3be 100644 --- a/SparseMatrix/COO_Matrix.hpp +++ b/SparseMatrix/COO_Matrix.hpp @@ -563,7 +563,7 @@ inline void COO_Matrix::printMatrix() { std::cout << "(" << this->row_indexes_[i] << ", " << this->column_indexes_[i] - << ", " << this->values_[i] << ")\n"; + << ", " << this->values_[i] << "),\n"; } std::cout << std::flush; } diff --git a/SystemModel.hpp b/SystemModel.hpp index 14b7e4e..0b9fb09 100644 --- a/SystemModel.hpp +++ b/SystemModel.hpp @@ -107,6 +107,7 @@ class SystemModel : public ModelEvaluatorImpl using ModelEvaluatorImpl::param_; using ModelEvaluatorImpl::param_up_; using ModelEvaluatorImpl::param_lo_; + using ModelEvaluatorImpl::J_; public: /** @@ -192,7 +193,14 @@ class SystemModel : public ModelEvaluatorImpl */ bool hasJacobian() { - return false; + for (const auto &component : components_) + { + if (!component->hasJacobian()) + { + return false; + } + } + return false; //Change this to false to make other tests run } /** @@ -390,7 +398,31 @@ class SystemModel : public ModelEvaluatorImpl * only. * */ - int evaluateJacobian(){return 0;} + int evaluateJacobian() + { + IdxT varOffset = 0; + IdxT optOffset = 0; + for(const auto& component : components_) { + for(IdxT j=0; jsize(); ++j) + { + component->y()[j] = y_[varOffset + j]; + component->yp()[j] = yp_[varOffset + j]; + } + varOffset += component->size(); + + for(IdxT j=0; jsize_opt(); ++j) + { + component->param()[j] = param_[optOffset + j]; + } + optOffset += component->size_opt(); + + component->evaluateJacobian(); + } + for(const auto& component : components_) { + J_ = component->getJacobian(); + } + return 0; + } /** * @brief Evaluate integrands for the system quadratures. @@ -692,6 +724,8 @@ class SystemModel : public ModelEvaluatorImpl { component->updateTime(t, a); } + time_ = t; + alpha_ = a; } void addBus(bus_type* bus) @@ -707,6 +741,7 @@ class SystemModel : public ModelEvaluatorImpl private: std::vector buses_; std::vector components_; + bool usejac_; }; // class SystemModel From 30d4a44b69e259cbec37aa6c8e4df235fed3108b Mon Sep 17 00:00:00 2001 From: Paulm991 Date: Thu, 27 Jun 2024 12:50:42 -0400 Subject: [PATCH 2/5] GENROU test working --- ComponentLib/DynamicPhasor/Bus/CMakeLists.txt | 5 + ComponentLib/DynamicPhasor/Bus/PhasorBus.hpp | 100 ++++ ComponentLib/DynamicPhasor/Bus/SlackBus.cpp | 87 ++++ ComponentLib/DynamicPhasor/Bus/SlackBus.hpp | 146 ++++++ ComponentLib/DynamicPhasor/CMakeLists.txt | 3 +- .../SynchronousMachine/GENROUwS/GENROU.cpp | 150 +++--- .../SynchronousMachine/GENROUwS/GENROU.hpp | 457 ------------------ Examples/GenrouTest/GenrouTest.cpp | 34 +- Solver/Dynamic/Ida.cpp | 19 +- SystemModel.hpp | 6 +- 10 files changed, 446 insertions(+), 561 deletions(-) create mode 100644 ComponentLib/DynamicPhasor/Bus/CMakeLists.txt create mode 100644 ComponentLib/DynamicPhasor/Bus/PhasorBus.hpp create mode 100644 ComponentLib/DynamicPhasor/Bus/SlackBus.cpp create mode 100644 ComponentLib/DynamicPhasor/Bus/SlackBus.hpp diff --git a/ComponentLib/DynamicPhasor/Bus/CMakeLists.txt b/ComponentLib/DynamicPhasor/Bus/CMakeLists.txt new file mode 100644 index 0000000..8d713c9 --- /dev/null +++ b/ComponentLib/DynamicPhasor/Bus/CMakeLists.txt @@ -0,0 +1,5 @@ +gridkit_add_library(slackbus + SOURCES + SlackBus.cpp + OUTPUT_NAME + gridkit_slackbus) \ No newline at end of file diff --git a/ComponentLib/DynamicPhasor/Bus/PhasorBus.hpp b/ComponentLib/DynamicPhasor/Bus/PhasorBus.hpp new file mode 100644 index 0000000..3e8ac00 --- /dev/null +++ b/ComponentLib/DynamicPhasor/Bus/PhasorBus.hpp @@ -0,0 +1,100 @@ +/* Paul Moon 6/26/2024 */ + +#ifndef _PHASOR_BUS_HPP_ +#define _PHASOR_BUS_HPP_ + +#include + +namespace ModelLib +{ + /*! + * @brief Base class for all dynamic phasor buses. + * + * Derived bus types: + * 0 - swing bus (Vr and Vi are constants) + * 1 - PV bus (P and V are constants) + * 2 - PQ bus (P and Q are constants) + * + * @todo Consider static instead of dynamic polymorphism for + * bus types. Create Bus class that takes template parameter + * BusType. + */ + template + class PhasorBus : public ModelEvaluatorImpl + { + protected: + using ModelEvaluatorImpl::size_; + using ModelEvaluatorImpl::nnz_; + using ModelEvaluatorImpl::time_; + using ModelEvaluatorImpl::alpha_; + using ModelEvaluatorImpl::rtol_; + using ModelEvaluatorImpl::atol_; + using ModelEvaluatorImpl::y_; + using ModelEvaluatorImpl::yp_; + using ModelEvaluatorImpl::tag_; + using ModelEvaluatorImpl::f_; + using ModelEvaluatorImpl::g_; + using ModelEvaluatorImpl::yB_; + using ModelEvaluatorImpl::ypB_; + using ModelEvaluatorImpl::fB_; + using ModelEvaluatorImpl::gB_; + using ModelEvaluatorImpl::param_; + using ModelEvaluatorImpl::param_up_; + using ModelEvaluatorImpl::param_lo_; + + public: + typedef typename ModelEvaluatorImpl::real_type real_type; + + enum BusType{PQ=1, PV, Slack, Isolated}; + + PhasorBus(IdxT id) : busID_(id) {} + virtual ~PhasorBus(){} + + // Set defaults for ModelEvaluator methods + virtual int allocate() { return 0;} + virtual int initialize() { return 0;} + virtual int tagDifferentiable() { return 0;} + virtual int evaluateResidual() { return 0;} + virtual int evaluateJacobian() { return 0;} + virtual int evaluateIntegrand() { return 0;} + + virtual int initializeAdjoint() { return 0;} + virtual int evaluateAdjointResidual() { return 0;} + //virtual int evaluateAdjointJacobian() { return 0;} + virtual int evaluateAdjointIntegrand() { return 0;} + virtual void updateTime(real_type, real_type) {} // <- throw exception here + + // Pure virtual methods specific to Bus types + virtual ScalarT& Vr() = 0; + virtual const ScalarT& Vr() const = 0; + virtual ScalarT& Vi() = 0; + virtual const ScalarT& Vi() const = 0; + virtual ScalarT& P() = 0; + virtual const ScalarT& P() const = 0; + virtual ScalarT& Q() = 0; + virtual const ScalarT& Q() const = 0; + + virtual ScalarT& lambdaP() = 0; + virtual const ScalarT& lambdaP() const = 0; + virtual ScalarT& lambdaQ() = 0; + virtual const ScalarT& lambdaQ() const = 0; + virtual ScalarT& PB() = 0; + virtual const ScalarT& PB() const = 0; + virtual ScalarT& QB() = 0; + virtual const ScalarT& QB() const = 0; + + virtual const int BusType() const = 0; + + virtual const IdxT BusID() const + { + return busID_; + } + + protected: + const IdxT busID_; + }; // class PhasorBus + +} // namespace ModelLib + + +#endif // _PHASOR_BUS_HPP_ diff --git a/ComponentLib/DynamicPhasor/Bus/SlackBus.cpp b/ComponentLib/DynamicPhasor/Bus/SlackBus.cpp new file mode 100644 index 0000000..d4924aa --- /dev/null +++ b/ComponentLib/DynamicPhasor/Bus/SlackBus.cpp @@ -0,0 +1,87 @@ +/* Created by Paul Moon 6/26/2024 */ +#include +#include +#include "SlackBus.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a slack bus + * + * Arguments passed to ModelEvaluatorImpl: + * - Number of equations = 0 (size_) + * - Number of variables = 0 (size_) + * - Number of quadratures = 0 + * - Number of optimization parameters = 0 + */ +template +SlackBus::SlackBus() + : PhasorBus(0), Vr_(0.0), Vi_(0.0), P_(0.0), Q_(0.0), PB_(0.0), QB_(0.0) +{ + //std::cout << "Create SlackBus..." << std::endl; + //std::cout << "Number of equations is " << size_ << std::endl; + + size_ = 0; +} + +/*! + * @brief SlackBus constructor. + * + * Arguments passed to ModelEvaluatorImpl: + * - Number of equations = 0 (size_) + * - Number of variables = 0 (size_) + * - Number of quadratures = 0 + * - Number of optimization parameters = 0 + */ +template +SlackBus::SlackBus(ScalarT Vr, ScalarT Vi) + : PhasorBus(0), Vr_(Vr), Vi_(Vi), P_(0.0), Q_(0.0), PB_(0.0), QB_(0.0) +{ + //std::cout << "Create SlackBus..." << std::endl; + //std::cout << "Number of equations is " << size_ << std::endl; + P() = 0.0; + Q() = 0.0; + size_ = 0; +} + +template +SlackBus::SlackBus(BusData& data) + : PhasorBus(data.bus_i), Vr_(data.Vm*cos(data.Va)), Vi_(data.Vm*sin(data.Va)) +{ + //std::cout << "Create SlackBus..." << std::endl; + //std::cout << "Number of equations is " << size_ << std::endl; + P() = 0.0; + Q() = 0.0; + size_ = 0; +} + +template +SlackBus::~SlackBus() +{ +} + +template +int SlackBus::evaluateResidual() +{ + // std::cout << "Evaluating residual of a slack bus ...\n"; + P() = 0.0; + Q() = 0.0; + return 0; +} + +template +int SlackBus::evaluateAdjointResidual() +{ + PB() = 0.0; + QB() = 0.0; + return 0; +} + + +// Available template instantiations +template class SlackBus; +template class SlackBus; + + +} // namespace ModelLib + diff --git a/ComponentLib/DynamicPhasor/Bus/SlackBus.hpp b/ComponentLib/DynamicPhasor/Bus/SlackBus.hpp new file mode 100644 index 0000000..c6bb630 --- /dev/null +++ b/ComponentLib/DynamicPhasor/Bus/SlackBus.hpp @@ -0,0 +1,146 @@ +/* Created by Paul Moon 6/26/2024 */ +#ifndef _SLACK_BUS_HPP_ +#define _SLACK_BUS_HPP_ + +#include +#include "PhasorBus.hpp" + +namespace ModelLib +{ + /*! + * @brief Implementation of a slack bus. + * + * Slack bus sets voltage _V_ and phase _theta_ as constants. + * Active and reactive power, _P_ and _Q_, are component model outputs, + * but are computed outside the SlackBus class. + * + * + */ + template + class SlackBus : public PhasorBus + { + using PhasorBus::size_; + using PhasorBus::y_; + using PhasorBus::yp_; + using PhasorBus::f_; + using PhasorBus::g_; + using PhasorBus::atol_; + using PhasorBus::rtol_; + + public: + using real_type = typename ModelEvaluatorImpl::real_type; + using BusData = GridKit::PowerSystemData::BusData; + + SlackBus(); + SlackBus(ScalarT Vr, ScalarT Vi); + SlackBus(BusData& data); + virtual ~SlackBus(); + virtual int evaluateResidual(); + virtual int evaluateAdjointResidual(); + + /// @todo Should slack bus allow changing voltage? + virtual ScalarT& Vr() + { + return Vr_; + } + + virtual const ScalarT& Vr() const + { + return Vr_; + } + + /// @todo Should slack bus allow changing phase? + virtual ScalarT& Vi() + { + return Vi_; + } + + virtual const ScalarT& Vi() const + { + return Vi_; + } + + virtual ScalarT& P() + { + return P_; + } + + virtual const ScalarT& P() const + { + return P_; + } + + virtual ScalarT& Q() + { + return Q_; + } + + virtual const ScalarT& Q() const + { + return Q_; + } + + /// @todo Should slack bus allow changing voltage? + virtual ScalarT& lambdaP() + { + return thetaB_; + } + + virtual const ScalarT& lambdaP() const + { + return thetaB_; + } + + /// @todo Should slack bus allow changing phase? + virtual ScalarT& lambdaQ() + { + return VB_; + } + + virtual const ScalarT& lambdaQ() const + { + return VB_; + } + + virtual ScalarT& PB() + { + return PB_; + } + + virtual const ScalarT& PB() const + { + return PB_; + } + + virtual ScalarT& QB() + { + return QB_; + } + + virtual const ScalarT& QB() const + { + return QB_; + } + + virtual const int BusType() const + { + return PhasorBus::BusType::Slack; + } + + private: + ScalarT Vr_; + ScalarT Vi_; + ScalarT P_; + ScalarT Q_; + + ScalarT VB_; + ScalarT thetaB_; + ScalarT PB_; + ScalarT QB_; + + }; // class SlackBus + +} // namespace ModelLib + + +#endif // _SLACK_BUS_HPP_ \ No newline at end of file diff --git a/ComponentLib/DynamicPhasor/CMakeLists.txt b/ComponentLib/DynamicPhasor/CMakeLists.txt index b484860..ee3a45e 100644 --- a/ComponentLib/DynamicPhasor/CMakeLists.txt +++ b/ComponentLib/DynamicPhasor/CMakeLists.txt @@ -1 +1,2 @@ -add_subdirectory(SynchronousMachine) \ No newline at end of file +add_subdirectory(SynchronousMachine) +add_subdirectory(Bus) \ No newline at end of file diff --git a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp index 2aea397..6a84d51 100644 --- a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp +++ b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp @@ -19,7 +19,7 @@ namespace ModelLib { */ template GENROU::GENROU(bus_type* bus, ScalarT P0, ScalarT Q0) - : ModelEvaluatorImpl(20, 1, 2), + : ModelEvaluatorImpl(20, 0, 0), H_(3.0), D_(0.0), Xq_(0.5), @@ -34,7 +34,7 @@ GENROU::GENROU(bus_type* bus, ScalarT P0, ScalarT Q0) Tq0pp_(0.05), Td0pp_(0.035), Ef_(2.6904), - Pm_(1.0), + Pm_(0.998273), X_(0.22), Xl_(0.15), omega_s_(0.0), @@ -45,8 +45,8 @@ GENROU::GENROU(bus_type* bus, ScalarT P0, ScalarT Q0) omega_lo_(omega_s_ - 0.0001), c_(10000.0), beta_(2), - P0_(P0), - Q0_(Q0), + P0_(1), + Q0_(0.5723), bus_(bus) { } @@ -95,17 +95,17 @@ template int GENROU::initialize() { // std::cout << "Initialize GENROU..." << std::endl; - //const ScalarT Vr = V()*cos(theta()); - //const ScalarT Vi = V()*sin(theta()); - //std::cout << "Initial Vr & Vi: " << Vr << " " << Vi << std::endl; - - //const ScalarT Ir = (P0_*Vr + Q0_*Vi)/((Vr*Vr)+(Vi*Vi)); - //const ScalarT Ii = -(Q0_*Vr - P0_*Vi)/((Vr*Vr)+(Vi*Vi)); - //std::cout << "Initial Ir & Ii: " << Ir << " " << Ii << std::endl; + const ScalarT Vr = V()*cos(theta()); + const ScalarT Vi = V()*sin(theta()); + std::cout << "Initial Vr & Vi: " << Vr << " " << Vi << std::endl; + const ScalarT Ir = Vi/X_; + const ScalarT Ii = -0.3286; + std::cout << "Initial Ir & Ii: " << Ir << " " << Ii << std::endl; + //delta_init = atan((Xq_1*Ir_init + Vi_init)/(- Xq_1 * Ii_init + Vr_init)); // Compute initial guess for the generator voltage phase - //const ScalarT delta = atan((Xq_*P0_ - Rs_*Q0_) / (V()*V() + Rs_*P0_ + Xq_*Q0_)) + theta(); - + const ScalarT delta = atan((Xq_*Ir + Vi) / (-Xq_*Ii + Vr)); + std::cout << "Delta: " << delta << std::endl; // Compute initial guess for the generator current phase //const ScalarT phi = theta() - delta - atan(Q0_/P0_); @@ -114,37 +114,52 @@ int GENROU::initialize() // Compute initial gueses for generator currents and potentials in d-q frame //const ScalarT Id = std::sqrt(P0_*P0_ + Q0_*Q0_)/V() * sin(phi); //const ScalarT Iq = std::sqrt(P0_*P0_ + Q0_*Q0_)/V() * cos(phi); - //const ScalarT Id = 0.9909; - //const ScalarT Iq = 0.3553; + //Id_init = Ir_init * sin(delta_init) - Ii_init * cos(delta_init) + //Iq_init = Ir_init * cos(delta_init) + Ii_init * sin(delta_init) + const ScalarT Id = Ir * sin(delta) - Ii * cos(delta); + const ScalarT Iq = Ir * cos(delta) + Ii * sin(delta); - //std::cout << "Initial Id & Iq: " << Id << " " << Iq << std::endl; + std::cout << "Initial Id & Iq: " << Id << " " << Iq << std::endl; - // const ScalarT Vd = V()*sin(theta() - delta) + Rs_*Id + Xqp_*Iq; - // const ScalarT Vq = V()*cos(theta() - delta) + Rs_*Iq - Xdp_*Id; + //Vq_init = Vr_init * cos(delta_init) + Vi_init * sin(delta_init) + //Vd_init = Vr_init * sin(delta_init) - Vi_init * cos(delta_init) - //std::cout << "Initial Vd & Vq: " << Vd << " " << Vq << std::endl; + const ScalarT Vq = Vr * cos(delta) + Vi * sin(delta); + const ScalarT Vd = Vr * sin(delta) - Vi * cos(delta); - /*const ScalarT Vr = 1.0723; - const ScalarT Vi = 0.22; - const ScalarT Ir = 1; - const ScalarT Ii = -0.3286; - const ScalarT delta = 0.909; - const ScalarT Vd = 0.7107; - const ScalarT Vq = 0.8326; - const ScalarT Eqp = 1.1298; - const ScalarT Psidp = 0.9614; - const ScalarT Psiqp = 0.6645; - const ScalarT Edp = 0.533; - const ScalarT Psiqpp = -0.6111; - const ScalarT Psidpp = 1.110; - const ScalarT Psipp = 1.2671; - const ScalarT ksat = 1; - const ScalarT Telec = 1; - const ScalarT Pmech = 0.9999;*/ + std::cout << "Initial Vq & Vd: " << Vq << " " << Vd << std::endl; + + //psi_pp_q_init = -Vd_init - Id_init * Ra_1 + Iq_init * X_pp_q_1; + //psi_pp_d_init = Vq_init + Id_init * X_pp_q_1 + Iq_init * Ra_1; + const ScalarT Psiqpp = -Vd - Id*Rs_ + Iq*Xqpp_; + const ScalarT Psidpp = Vq + Id*Xqpp_ + Iq*Rs_; + + std::cout << "Initial Psiqpp & Psidpp: " << Psiqpp << " " << Psidpp << std::endl; + + //psi_pp_init = sqrt(psi_pp_q_init^2 + psi_pp_d_init^2); + const ScalarT Psipp = sqrt((Psiqpp*Psiqpp) + (Psidpp*Psidpp)); + //ksat_init = 0; + const ScalarT ksat = 1.0; + //Telec_term1 = psi_pp_d_init + Id_init * X_pp_d_1 + //Telec_term2 = psi_pp_q_init + Iq_init * X_pp_d_1 + //Telec_init = Telec_term1 * Iq_init - Telec_term2 * Id_init; + //Pmech_init = Telec_init %1.0337 + const ScalarT Telec = (Psidpp + Id*Xdpp_) * Iq - (Psiqpp + Iq*Xdpp_) * Id; + const ScalarT Pmech = Telec; + + //psi_p_q_init = (1/(1+(X_p_q_1 - X_pp_q_1)/(X_pp_q_1 - Xl_1)))*(-psi_pp_q_init * (X_p_q_1 - Xl_1)/(X_pp_q_1 - Xl_1) + Iq_init * (X_p_q_1 - Xl_1)) + const ScalarT Psiqp = (1.0/(1.0+(Xqp_ - Xqpp_)/(Xqpp_ - Xl_)))*(-Psiqpp * (Xqp_ - Xl_)/(Xqpp_ - Xl_) + Iq * (Xqp_ - Xl_)); + + //E_p_d_init = psi_p_q_init - Iq_init * (X_p_q_1 - Xl_1) + const ScalarT Edp = Psiqp - Iq*(Xqp_-Xl_); + //psi_p_d_init = -(1/(1+ (X_p_d_1 - X_pp_d_1)/(X_pp_d_1 - Xl_1))) * (Id_init*(X_p_d_1 - Xl_1)-psi_pp_d_init * (X_p_d_1 - Xl_1)/(X_pp_d_1 - Xl_1)) + const ScalarT Psidp = -(1.0/(1.0+ (Xdp_ - Xdpp_)/(Xdpp_ - Xl_))) * (Id*(Xdp_ - Xl_)-Psidpp * (Xdp_ - Xl_)/(Xdpp_ - Xl_)); + //E_p_q_init = psi_p_d_init + Id_init * (X_p_d_1 - Xl_1) + const ScalarT Eqp = Psidp + Id * (Xdp_-Xl_); //y0 = [delta_init; omega_init; E_p_q_init; psi_p_d_init; psi_p_q_init; E_p_d_init; psi_pp_q_init; psi_pp_d_init; psi_pp_init; ksat_init; Vd_init; Vq_init ; Telec_init; Id_init; Iq_init; Vr_init; Vi_init; Ir_init; Ii_init;Pmech_init] - /*y_[0] = delta; + y_[0] = delta; y_[1] = omega_s_; y_[2] = Eqp; y_[3] = Psidp; @@ -163,8 +178,8 @@ int GENROU::initialize() y_[16] = Vi; y_[17] = Ir; y_[18] = Ii; - y_[19] = Pmech;*/ - y_[0] = 0.5273; + y_[19] = Pmech; + /*y_[0] = 0.5273; y_[1] = 0.0; y_[2] = 1.1948; y_[3] = 1.1554; @@ -183,26 +198,6 @@ int GENROU::initialize() y_[16] = 0.22; y_[17] = 1.0; y_[18] = -0.3286; - y_[19] = 1.0; - /*y_[0] = 0.9224; - y_[1] = 0.1765; - y_[2] = 1.2163; - y_[3] = 1.1554; - y_[4] = -0.3907; - y_[5] = 0.0; - y_[6] = 0.3572; - y_[7] = 0.2676; - y_[8] = 0.4464; - y_[9] = 2.5613; - y_[10] = -0.4203; - y_[11] = 0.3149; - y_[12] = 0.0; - y_[13] = 0.7872; - y_[14] = 1.0507; - y_[15] = 0.0797; - y_[16] = 0.2776; - y_[17] = 1.2620; - y_[18] = 0.3622; y_[19] = 1.0;*/ yp_[0] = 0.0; yp_[1] = 0.0; @@ -225,10 +220,10 @@ int GENROU::initialize() yp_[18] = 0.0; yp_[19] = 0.0; - /*for (int i = 0; i < 20; i++) { + for (int i = 0; i < 20; i++) { std::cout << y_[i] << ", "; } - std::cout << std::endl;*/ + std::cout << std::endl; return 0; } @@ -295,12 +290,12 @@ int GENROU::evaluateResidual() double real_part = division_result.real(); // real part of (numerator / denominator) double imag_part = division_result.imag(); // imaginary part of (numerator / denominator) std::complex i(0.0, 1.0); // imaginary unit - + //std::cout << y_[17] << " " << y_[18] << std::endl; //std::cout << "SA: " << SA << " SB: " << SB << std::endl; // std::cout << "Evaluate residual for GENROU..." << std::endl; f_[0] = -yp_[0] + omega0_*y_[1]; //f1 = -yp(1)+y(2)*w0_1; - f_[1] = -yp_[1] + (1/(2*H_))*((Pm_ - D_*y_[1])/(1 + y_[1]) - y_[12]); //f2 = -yp(2) + (1/(2*H_1))*((Pmech_init - D_1*y(2))/(1 + y(2)) - y(13)); + f_[1] = -yp_[1] + (1/(2*H_))*((y_[19] - D_*y_[1])/(1 + y_[1]) - y_[12]); //f2 = -yp(2) + (1/(2*H_1))*((Pmech_init - D_1*y(2))/(1 + y(2)) - y(13)); f_[2] = -yp_[2]*Td0p_ + Ef_ - (y_[2] + (Xd_-Xdp_)*y_[13]+((Xdp_-Xdpp_)/((Xdp_-Xl_)*(Xdp_-Xl_)))*(y_[2]-y_[3]-(Xdp_-Xl_)*y_[13])) + y_[7]*y_[9]; //f3 = -yp(3)*T_p_d0_1 + Efd - (y(3) + (Xd_1-X_p_d_1)*(y(14) + ((X_p_d_1-X_pp_d_1)/(X_p_d_1 - Xl_1)^2)*(y(3)-y(4)-(X_p_d_1-Xl_1)*y(14))) + y(8)*y(10)); f_[3] = -yp_[3]+(1/Td0pp_)*(y_[2]-y_[3]-(Xdp_-Xl_)*y_[13]); //f4 = -yp(4)+(1/T_pp_d0_1)*(y(3)-y(4)-(X_p_d_1-Xl_1)*y(14)); f_[4] = -yp_[4]+(1/Tq0pp_)*(y_[5]-y_[4]+(Xqp_-Xl_)*y_[14]); //f5 = -yp(5)+(1/T_pp_q0_1)*(y(6)-y(5)+(X_p_q_1-Xl_1)*y(15)); @@ -308,7 +303,7 @@ int GENROU::evaluateResidual() f_[6] = -y_[6] - y_[4]*((Xqp_-Xqpp_)/(Xqp_-Xl_))-y_[5]*((Xqpp_-Xl_)/(Xqp_-Xl_)); //f7 = -y(7) - y(5)*((X_p_q_1-X_pp_q_1)/(X_p_q_1 - Xl_1))- y(6)*((X_pp_q_1 - Xl_1)/(X_p_q_1 - Xl_1)); f_[7] = -y_[7] - y_[3]*((Xdp_-Xdpp_)/(Xdp_-Xl_))+y_[2]*((Xdpp_-Xl_)/(Xdp_-Xl_)); //f8 = -y(8) - y(4)*((X_p_d_1-X_pp_d_1)/(X_p_d_1 - Xl_1)) + y(3)*((X_pp_d_1 - Xl_1)/(X_p_d_1 - Xl_1)); - f_[8] = -y_[8] + sqrt((y_[7]*y_[7])+(y_[6]*y_[6]));//f9 = -y(9) + sqrt(y(8)^2+y(7)^2); + f_[8] = -y_[8] + sqrt((y_[7]*y_[7])+(y_[6]*y_[6])); //f9 = -y(9) + sqrt(y(8)^2+y(7)^2); f_[9] = -y_[9] + SB*((y_[8]-SA)*(y_[8]-SA)); //f10 = -y(10) + SB * (y(9) - SA)^2; f_[10] = -y_[10] - y_[6] * (1+y_[1]); //f11 = -y(11) - y(7) * (1+y(2)); f_[11] = -y_[11] + y_[7] * (1+y_[1]); //f12 = -y(12) + y(8) * (1+y(2)); @@ -317,10 +312,14 @@ int GENROU::evaluateResidual() f_[14] = -y_[14] + y_[17] * cos(y_[0]) + y_[18] * sin(y_[0]); //f15 = -y(15) + y(18) * cos(y(1)) + y(19) * sin(y(1)); % prev: -y(15) + Ireal * cos(y(1)) + Iimag * sin(y(1)); f_[15] = -y_[10] + y_[15] * sin(y_[0]) - y_[16]*cos(y_[0]) + y_[13]*Rs_ - y_[14]*Xqpp_; //f16 = -y(11) + y(16)*sin(y(1)) - y(17)*cos(y(1))+y(14)*Ra_1 - y(15)*X_pp_q_1; f_[16] = -y_[11] + y_[15] * cos(y_[0]) + y_[16]*sin(y_[0]) + y_[13]*Xqpp_ + y_[14]*Rs_; //f17 = -y(12) + y(16)*cos(y(1)) + y(17)*sin(y(1)) + y(14)*X_pp_q_1 + y(15)*Ra_1; - f_[17] = -y_[17] + std::real((y_[15]+i*y_[16])/(i*X_));// y_[16]/X_; //f18 = -y(18)+real((y(16) + 1i * y(17))/(1i*X)); - f_[18] = -y_[18] + std::imag((y_[15]+i*y_[16])/(i*X_));;// y_[15]/X_; //f19 = -y(19)+imag((y(16) + 1i * y(17))/(1i*X)); + f_[17] = -y_[17] + std::real((y_[15]+i*y_[16])/(i*X_)); // y_[16]/X_; //f18 = -y(18)+real((y(16) + 1i * y(17))/(1i*X)); + //std::cout << "f_[17]: " << f_[17] << " y17: " << -y_[17] << " y15: " << y_[15] << " y16: " << y_[16] << std::endl; + f_[18] = -y_[18] + std::imag((y_[15]+i*y_[16])/(i*X_)); // y_[15]/X_; //f19 = -y(19)+imag((y(16) + 1i * y(17))/(1i*X)); f_[19] = -Pm_ + y_[19]; //f20 = -Pmech_init + y(20); + P() += Pg(); + Q() += Qg(); + return 0; } @@ -493,7 +492,7 @@ int GENROU::evaluateAdjointIntegrand() template ScalarT GENROU::Pg() { - return 1; + return y_[14]*V()*cos(theta() - y_[0]) + y_[13]*V()*sin(theta() - y_[0]); } /** @@ -504,7 +503,7 @@ ScalarT GENROU::Pg() template ScalarT GENROU::Qg() { - return 1; + return y_[14]*V()*sin(theta() - y_[0]) - y_[13]*V()*cos(theta() - y_[0]); } /** @@ -513,7 +512,7 @@ ScalarT GENROU::Qg() template ScalarT GENROU::frequencyPenalty(ScalarT omega) { - return 0; + return c_ * pow(std::max(0.0, std::max(omega - omega_up_, omega_lo_ - omega)), beta_); } /** @@ -523,7 +522,18 @@ ScalarT GENROU::frequencyPenalty(ScalarT omega) template ScalarT GENROU::frequencyPenaltyDer(ScalarT omega) { - return 0.0; + if (omega > omega_up_) + { + return beta_ * c_ * pow(omega - omega_up_, beta_ - 1.0); + } + else if (omega < omega_lo_) + { + return beta_ * c_ * pow(omega - omega_lo_, beta_ - 1.0); + } + else + { + return 0.0; + } } diff --git a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp index 6d266ee..eecb00f 100644 --- a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp +++ b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp @@ -105,230 +105,11 @@ namespace ModelLib private: - const ScalarT& Pm() const - { - return param_[0]; - } - - const ScalarT& Ef() const - { - return param_[1]; - } - ScalarT Pg(); ScalarT Qg(); ScalarT frequencyPenalty(ScalarT omega); ScalarT frequencyPenaltyDer(ScalarT omega); - private: - // - // Private inlined accessor methods - // - - const ScalarT dotDelta() const - { - return yp_[0]; - } - - const ScalarT dotOmega() const - { - return yp_[1]; - } - - const ScalarT dotEqp() const - { - return yp_[2]; - } - - const ScalarT dotPsidp() const - { - return yp_[3]; - } - - const ScalarT dotPsiqp() const - { - return yp_[4]; - } - - const ScalarT dotEdp() const - { - return yp_[5]; - } - - const ScalarT dotPsiqpp() const - { - return yp_[6]; - } - - const ScalarT dotPsidpp() const - { - return yp_[7]; - } - - const ScalarT dotPsipp() const - { - return yp_[8]; - } - - const ScalarT dotksat() const - { - return yp_[9]; - } - - const ScalarT dotVd() const - { - return yp_[10]; - } - - const ScalarT dotVq() const - { - return yp_[11]; - } - - const ScalarT dotTelec() const - { - return yp_[12]; - } - - const ScalarT dotId() const - { - return yp_[13]; - } - - const ScalarT dotIq() const - { - return yp_[14]; - } - - const ScalarT dotVr() const - { - return yp_[15]; - } - - const ScalarT dotVi() const - { - return yp_[16]; - } - - const ScalarT dotIr() const - { - return yp_[17]; - } - - const ScalarT dotIi() const - { - return yp_[18]; - } - - const ScalarT dotPmech() const - { - return yp_[19]; - } - - // - // y_ values - // - - const ScalarT delta() const - { - return y_[0]; - } - - const ScalarT omega() const - { - return y_[1]; - } - - const ScalarT Eqp() const - { - return y_[2]; - } - - const ScalarT Psidp() const - { - return y_[3]; - } - - const ScalarT Psiqp() const - { - return y_[4]; - } - - const ScalarT Edp() const - { - return y_[5]; - } - - const ScalarT Psiqpp() const - { - return y_[6]; - } - - const ScalarT Psidpp() const - { - return y_[7]; - } - - const ScalarT Psipp() const - { - return y_[8]; - } - - const ScalarT ksat() const - { - return y_[9]; - } - - const ScalarT Vd() const - { - return y_[10]; - } - - const ScalarT Vq() const - { - return y_[11]; - } - - const ScalarT Telec() const - { - return y_[12]; - } - - const ScalarT Id() const - { - return y_[13]; - } - - const ScalarT Iq() const - { - return y_[14]; - } - - const ScalarT Vr() const - { - return y_[15]; - } - - const ScalarT Vi() const - { - return y_[16]; - } - - const ScalarT Ir() const - { - return y_[17]; - } - - const ScalarT Ii() const - { - return y_[18]; - } - - const ScalarT Pmech() const - { - return y_[19]; - } - private: real_type omega0_;///< Nominal frequency (assume 2pi60 rad/s) real_type H_; ///< Inertia constant [s] @@ -367,241 +148,3 @@ namespace ModelLib #endif // _GENROU_H_ - - -/* Paul Moon 6/6/2024 */ - -/*#include "ModelEvaluatorImpl.hpp" - -namespace ModelLib -{ - template class BaseBus; -} - -namespace ModelLib -{ - /** - * @brief Generator - * - * @tparam ScalarT - Scalar type - * @tparam IdxT - Matrix and vector index type - */ - /*template - class GENROU : public ModelEvaluatorImpl - { - using ModelEvaluatorImpl::size_; - using ModelEvaluatorImpl::nnz_; - using ModelEvaluatorImpl::time_; - using ModelEvaluatorImpl::alpha_; - using ModelEvaluatorImpl::y_; - using ModelEvaluatorImpl::yp_; - using ModelEvaluatorImpl::tag_; - using ModelEvaluatorImpl::f_; - using ModelEvaluatorImpl::g_; - using ModelEvaluatorImpl::yB_; - using ModelEvaluatorImpl::ypB_; - using ModelEvaluatorImpl::fB_; - using ModelEvaluatorImpl::gB_; - using ModelEvaluatorImpl::param_; - using ModelEvaluatorImpl::param_up_; - using ModelEvaluatorImpl::param_lo_; - - typedef typename ModelEvaluatorImpl::real_type real_type; - typedef BaseBus bus_type; - - public: - GENROU(BaseBus* bus, ScalarT P0 = 1.0, ScalarT Q0 = 0.0); - virtual ~GENROU(); - - int allocate(); - int initialize(); - int tagDifferentiable(); - int evaluateResidual(); - int evaluateJacobian(); - int evaluateIntegrand(); - - int initializeAdjoint(); - int evaluateAdjointResidual(); - //int evaluateAdjointJacobian(); - int evaluateAdjointIntegrand(); - - void updateTime(real_type t, real_type a) - { - time_ = t; - alpha_ = a; - } - - // Inline accesor functions - ScalarT& V() - { - return bus_->V(); - } - - const ScalarT& V() const - { - return bus_->V(); - } - - ScalarT& theta() - { - return bus_->theta(); - } - - const ScalarT& theta() const - { - return bus_->theta(); - } - - ScalarT& P() - { - return bus_->P(); - } - - const ScalarT& P() const - { - return bus_->P(); - } - - ScalarT& Q() - { - return bus_->Q(); - } - - const ScalarT& Q() const - { - return bus_->Q(); - } - - - private: - const ScalarT& Pm() const - { - return param_[0]; - } - - const ScalarT& Ef() const - { - return param_[1]; - } - - ScalarT Pg(); - ScalarT Qg(); - ScalarT frequencyPenalty(ScalarT omega); - ScalarT frequencyPenaltyDer(ScalarT omega); - - private: - // - // Private inlined accessor methods - // - - const ScalarT dotDelta() const - { - return yp_[0]; - } - - const ScalarT dotOmega() const - { - return yp_[1]; - } - - const ScalarT dotEdp() const - { - return yp_[2]; - } - - const ScalarT dotEqp() const - { - return yp_[3]; - } - - const ScalarT delta() const - { - return y_[0]; - } - - const ScalarT omega() const - { - return y_[1]; - } - - const ScalarT Edp() const - { - return y_[2]; - } - - const ScalarT Eqp() const - { - return y_[3]; - } - - const ScalarT Id() const - { - return y_[4]; - } - - const ScalarT Iq() const - { - return y_[5]; - } - - private: - inline ScalarT frequencyPenalty(ScalarT omega); - inline ScalarT frequencyPenaltyDer(ScalarT omega); - - private: - real_type omega0_;///< Nominal frequency (assume 2pi60 rad/s) - real_type H_; ///< Inertia constant [s] - real_type D_; ///< Damping constant [pu] - real_type Rs_; ///< stator resistance [pu] - real_type Xl_; ///< stator leakage reactance [pu] - real_type Xd_; ///< direct axis synchronous reactance [pu] - real_type Xdp_; ///< direct axis transient reactance [pu] - real_type Xdpp_; ///< direct axis sub-transient reactance [pu] - real_type Xq_; ///< quadrature axis synchronous reactance [pu] - real_type Xqp_; ///< quadrature axis transient reactance [pu] - real_type Xqpp_; ///< quadrature axis sub-transient reactance [pu] - real_type Td0p_; ///< open circuit direct axis transient time constant [s] - real_type Td0pp_; ///< open circuit direct axis sub-transient time constant [s] - real_type Tq0p_; ///< open circuit quadrature axis transient time constant [s] - real_type Tq0pp_; ///< open circuit quadrature axis sub-transient time constant [s] - real_type S1_; ///< saturation at 1.0 pu flux - real_type S12_; ///< saturation at 1.2 pu flux - - real_type delta_; ///< Machine internal angle - real_type omega_; ///< - real_type Eqp_; ///< - real_type psidp_; ///< - real_type psiqp_; ///< - real_type Edp_; ///< - real_type psi_; ///< - real_type ksat_; ///< Saturation coefficient - real_type Vd_; ///< - real_type Vq_; ///< - real_type Telec_; ///< - real_type Ir_; ///< - real_type Ii_; ///< - real_type Vr_; ///< - real_type Vi_; ///< - real_type Id_; ///< - real_type Iq_; ///< - real_type Pmech_; ///< - real_type Efd_; ///< - - - - - real_type Ef_; // - real_type Pm_; - real_type omega_s_; - real_type omega_b_; - real_type omega_up_; - real_type omega_lo_; - real_type c_; - real_type beta_; - - ScalarT P0_; - ScalarT Q0_; - - bus_type* bus_; - }; -}*/ \ No newline at end of file diff --git a/Examples/GenrouTest/GenrouTest.cpp b/Examples/GenrouTest/GenrouTest.cpp index 62d7c3a..bf07466 100644 --- a/Examples/GenrouTest/GenrouTest.cpp +++ b/Examples/GenrouTest/GenrouTest.cpp @@ -27,7 +27,7 @@ int main() // Create an infinite bus - BaseBus* bus = new BusSlack(1, 0); // roughly 1.072+0.22j + BaseBus* bus = new BusSlack(1.0946, 0.202); // roughly 1.072+0.22j // Attach a generator to that bus GENROU* gen = new GENROU(bus); @@ -53,45 +53,23 @@ int main() std::cout << "}\n"; - model->updateTime(0.0, 1.0); - model->evaluateJacobian(); + //model->updateTime(0.0, 1.0); + //model->evaluateJacobian(); std::cout << "Intial Jacobian with alpha = 1:\n"; - model->getJacobian().printMatrix(); + //model->getJacobian().printMatrix(); // Create numerical integrator and configure it for the generator model AnalysisManager::Sundials::Ida* idas = new AnalysisManager::Sundials::Ida(model); double t_init = 0.0; - double t_final = 1.0; + double t_final = 2.0; double t_timestep = 0.0001; idas->configureSimulation(); idas->getDefaultInitialCondition(); idas->initializeSimulation(t_init); - idas->runSimulation(t_final, t_final/t_timestep); - - /*for (double i = t_init; i <= t_final; i += t_timestep) { - //idas->configureSimulation(); - //idas->getDefaultInitialCondition(); - //idas->initializeSimulation(i); - - idas->runSimulation(i+t_timestep); - - std::vector& yfinial = model->y(); - - outFile << i; - for (size_t j = 0; j < yfinial.size(); j++) { - outFile << " " << yfinial[j]; - } - outFile << "\n"; - - std::cout << "Final Vector y\n"; - for (size_t j = 0; j < yfinial.size(); j++) - { - std::cout << yfinial[j] << "\n"; - } - }*/ + idas->runSimulation(t_final, 2000); delete idas; delete model; diff --git a/Solver/Dynamic/Ida.cpp b/Solver/Dynamic/Ida.cpp index 9263604..5cee776 100644 --- a/Solver/Dynamic/Ida.cpp +++ b/Solver/Dynamic/Ida.cpp @@ -286,11 +286,26 @@ namespace Sundials checkOutput(retval, "IDASolve"); //printOutput(tout); realtype *yval = N_VGetArrayPointer_Serial(yy_); + realtype *ypval = N_VGetArrayPointer_Serial(yp_); + /*yval[0] = fmod(yval[0], 2 * M_PI); + if (yval[0] < 0) { + yval[0] += 2 * M_PI; + }*/ outFile << tout; - for (size_t j = 0; j < model_->size(); j++) { + for (size_t j = 1; j < model_->size(); j++) { outFile << " " << yval[j]; } outFile << "\n"; + for (size_t j = 1; j < model_->size(); j++) { + if (yval[j] > 10) { + std::cout << "yval " << j << " is too high at " << yval[j] << " at time " << tout << std::endl; + } + } + + /*for (int i = 0; i < 6; i++) { + std::cout << ypval[i] << " "; + } + std::cout << std::endl;*/ if (retval == IDA_SUCCESS) { @@ -457,7 +472,7 @@ namespace Sundials checkOutput(retval, "IDASetUserDataB"); /// \todo Need to set max number of steps based on user input! - retval = IDASetMaxNumStepsB(solver_, backwardID_, 2000); + retval = IDASetMaxNumStepsB(solver_, backwardID_, 20000); checkOutput(retval, "IDASetMaxNumSteps"); // Set up linear solver diff --git a/SystemModel.hpp b/SystemModel.hpp index 0b9fb09..b6a9f53 100644 --- a/SystemModel.hpp +++ b/SystemModel.hpp @@ -116,9 +116,9 @@ class SystemModel : public ModelEvaluatorImpl SystemModel() : ModelEvaluatorImpl(0, 0, 0) { // Set system model tolerances - rtol_ = 1e-7; - atol_ = 1e-9; - this->max_steps_=2000; + rtol_ = 1e-12; + atol_ = 1e-12; + this->max_steps_=20000; } /** From c0569cd9d9725bcef9e547aa8f3c37324f94f509 Mon Sep 17 00:00:00 2001 From: Paulm991 Date: Tue, 9 Jul 2024 15:45:32 -0400 Subject: [PATCH 3/5] Functional GENROU Test --- .../SynchronousMachine/GENROUwS/GENROU.cpp | 219 +++++++----------- .../SynchronousMachine/GENROUwS/GENROU.hpp | 2 +- Examples/GenrouTest/GenrouTest.cpp | 6 +- ModelEvaluatorImpl.hpp | 2 +- Solver/Dynamic/Ida.cpp | 30 +-- SystemModel.hpp | 8 +- 6 files changed, 103 insertions(+), 164 deletions(-) diff --git a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp index 6a84d51..e9572ce 100644 --- a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp +++ b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp @@ -19,7 +19,7 @@ namespace ModelLib { */ template GENROU::GENROU(bus_type* bus, ScalarT P0, ScalarT Q0) - : ModelEvaluatorImpl(20, 0, 0), + : ModelEvaluatorImpl(15, 0, 0), H_(3.0), D_(0.0), Xq_(0.5), @@ -33,20 +33,20 @@ GENROU::GENROU(bus_type* bus, ScalarT P0, ScalarT Q0) Td0p_(7.0), Tq0pp_(0.05), Td0pp_(0.035), - Ef_(2.6904), - Pm_(0.998273), + Ef_(1.8336), + Pm_(1.0), X_(0.22), Xl_(0.15), omega_s_(0.0), - omega0_(2.0*60.0*M_PI), + omega0_(2.0*60.0*3.14159), S1_(0.05), S12_(0.2), omega_up_(omega_s_ + 0.0001), omega_lo_(omega_s_ - 0.0001), c_(10000.0), beta_(2), - P0_(1), - Q0_(0.5723), + P0_(1.0), + Q0_(0.0), bus_(bus) { } @@ -94,68 +94,71 @@ int GENROU::allocate() template int GENROU::initialize() { - // std::cout << "Initialize GENROU..." << std::endl; - const ScalarT Vr = V()*cos(theta()); - const ScalarT Vi = V()*sin(theta()); - std::cout << "Initial Vr & Vi: " << Vr << " " << Vi << std::endl; - - const ScalarT Ir = Vi/X_; - const ScalarT Ii = -0.3286; - std::cout << "Initial Ir & Ii: " << Ir << " " << Ii << std::endl; - //delta_init = atan((Xq_1*Ir_init + Vi_init)/(- Xq_1 * Ii_init + Vr_init)); - // Compute initial guess for the generator voltage phase - const ScalarT delta = atan((Xq_*Ir + Vi) / (-Xq_*Ii + Vr)); - std::cout << "Delta: " << delta << std::endl; - // Compute initial guess for the generator current phase - //const ScalarT phi = theta() - delta - atan(Q0_/P0_); - - //std::cout << "Initial delta & phi: " << delta << " " << phi << std::endl; - - // Compute initial gueses for generator currents and potentials in d-q frame - //const ScalarT Id = std::sqrt(P0_*P0_ + Q0_*Q0_)/V() * sin(phi); - //const ScalarT Iq = std::sqrt(P0_*P0_ + Q0_*Q0_)/V() * cos(phi); - //Id_init = Ir_init * sin(delta_init) - Ii_init * cos(delta_init) - //Iq_init = Ir_init * cos(delta_init) + Ii_init * sin(delta_init) - const ScalarT Id = Ir * sin(delta) - Ii * cos(delta); - const ScalarT Iq = Ir * cos(delta) + Ii * sin(delta); - - std::cout << "Initial Id & Iq: " << Id << " " << Iq << std::endl; - - //Vq_init = Vr_init * cos(delta_init) + Vi_init * sin(delta_init) - //Vd_init = Vr_init * sin(delta_init) - Vi_init * cos(delta_init) - - const ScalarT Vq = Vr * cos(delta) + Vi * sin(delta); - const ScalarT Vd = Vr * sin(delta) - Vi * cos(delta); - - std::cout << "Initial Vq & Vd: " << Vq << " " << Vd << std::endl; - - //psi_pp_q_init = -Vd_init - Id_init * Ra_1 + Iq_init * X_pp_q_1; - //psi_pp_d_init = Vq_init + Id_init * X_pp_q_1 + Iq_init * Ra_1; - const ScalarT Psiqpp = -Vd - Id*Rs_ + Iq*Xqpp_; - const ScalarT Psidpp = Vq + Id*Xqpp_ + Iq*Rs_; - - std::cout << "Initial Psiqpp & Psidpp: " << Psiqpp << " " << Psidpp << std::endl; - - //psi_pp_init = sqrt(psi_pp_q_init^2 + psi_pp_d_init^2); - const ScalarT Psipp = sqrt((Psiqpp*Psiqpp) + (Psidpp*Psidpp)); - //ksat_init = 0; - const ScalarT ksat = 1.0; - //Telec_term1 = psi_pp_d_init + Id_init * X_pp_d_1 - //Telec_term2 = psi_pp_q_init + Iq_init * X_pp_d_1 - //Telec_init = Telec_term1 * Iq_init - Telec_term2 * Id_init; - //Pmech_init = Telec_init %1.0337 - const ScalarT Telec = (Psidpp + Id*Xdpp_) * Iq - (Psiqpp + Iq*Xdpp_) * Id; - const ScalarT Pmech = Telec; - - //psi_p_q_init = (1/(1+(X_p_q_1 - X_pp_q_1)/(X_pp_q_1 - Xl_1)))*(-psi_pp_q_init * (X_p_q_1 - Xl_1)/(X_pp_q_1 - Xl_1) + Iq_init * (X_p_q_1 - Xl_1)) - const ScalarT Psiqp = (1.0/(1.0+(Xqp_ - Xqpp_)/(Xqpp_ - Xl_)))*(-Psiqpp * (Xqp_ - Xl_)/(Xqpp_ - Xl_) + Iq * (Xqp_ - Xl_)); - - //E_p_d_init = psi_p_q_init - Iq_init * (X_p_q_1 - Xl_1) - const ScalarT Edp = Psiqp - Iq*(Xqp_-Xl_); - //psi_p_d_init = -(1/(1+ (X_p_d_1 - X_pp_d_1)/(X_pp_d_1 - Xl_1))) * (Id_init*(X_p_d_1 - Xl_1)-psi_pp_d_init * (X_p_d_1 - Xl_1)/(X_pp_d_1 - Xl_1)) - const ScalarT Psidp = -(1.0/(1.0+ (Xdp_ - Xdpp_)/(Xdpp_ - Xl_))) * (Id*(Xdp_ - Xl_)-Psidpp * (Xdp_ - Xl_)/(Xdpp_ - Xl_)); - //E_p_q_init = psi_p_d_init + Id_init * (X_p_d_1 - Xl_1) - const ScalarT Eqp = Psidp + Id * (Xdp_-Xl_); + double SB= ((1/0.2)*(1/0.2))*((sqrt(S1_) - sqrt(S12_))*(sqrt(S1_) - sqrt(S12_))); + double SA = 1 - sqrt(S1_/SB); + const double Vr = V() * cos(theta()); + const double Vi = V() * sin(theta()); + //std::cout << "Initial voltage & theta: " << V() << " " << theta() << std::endl; + + //const double Ir = Vi / X_; + //const double Ii = -Vr / X_; + const double Ir = P0_/V(); + const double Ii = 0; + // delta_init = atan((Xq_1*Ir_init + Vi_init)/(- Xq_1 * Ii_init + Vr_init)); + // Compute initial guess for the generator voltage phase + //const double delta = atan((Xq_*P0_ - Rs_*Q0_) / (V*V + Rs_*P0_ + Xq_*Q0_)) + theta; + const double delta = atan((Xq_ * Ir + Vi) / (-Xq_ * Ii + Vr)); + // Compute initial guess for the generator current phase + // const double phi = theta() - delta - atan(Q0_/P0_); + + // std::cout << "Initial delta & phi: " << delta << " " << phi << std::endl; + + // Compute initial gueses for generator currents and potentials in d-q frame + // const double Id = std::sqrt(P0_*P0_ + Q0_*Q0_)/V() * sin(phi); + // const double Iq = std::sqrt(P0_*P0_ + Q0_*Q0_)/V() * cos(phi); + // Id_init = Ir_init * sin(delta_init) - Ii_init * cos(delta_init) + // Iq_init = Ir_init * cos(delta_init) + Ii_init * sin(delta_init) + const double Id = Ir * sin(delta) - Ii * cos(delta); + const double Iq = Ir * cos(delta) + Ii * sin(delta); + + // Vq_init = Vr_init * cos(delta_init) + Vi_init * sin(delta_init) + // Vd_init = Vr_init * sin(delta_init) - Vi_init * cos(delta_init) + + const double Vq = Vr * cos(delta) + Vi * sin(delta) + Id*Xqpp_ + Iq*Rs_; + const double Vd = Vr * sin(delta) - Vi * cos(delta) + Id*Rs_ - Iq*Xqpp_; + + // psi_pp_q_init = -Vd_init - Id_init * Ra_1 + Iq_init * X_pp_q_1; + // psi_pp_d_init = Vq_init + Id_init * X_pp_q_1 + Iq_init * Ra_1; + const double Psiqpp = -Vd / (1 + omega_s_); + const double Psidpp = Vq / (1 + omega_s_); + + // psi_pp_init = sqrt(psi_pp_q_init^2 + psi_pp_d_init^2); + const double Psipp = sqrt((Psiqpp * Psiqpp) + (Psidpp * Psidpp)); + // ksat_init = 0; + const double ksat = SB*((Psipp-SA)*(Psipp-SA)); + // Telec_term1 = psi_pp_d_init + Id_init * X_pp_d_1 + // Telec_term2 = psi_pp_q_init + Iq_init * X_pp_d_1 + // Telec_init = Telec_term1 * Iq_init - Telec_term2 * Id_init; + // Pmech_init = Telec_init %1.0337 + const double Telec = (Psidpp - Id * Xdpp_) * Iq - (Psiqpp - Iq * Xdpp_) * Id; + const double Pmech = Telec; + + // psi_p_q_init = (1/(1+(X_p_q_1 - X_pp_q_1)/(X_pp_q_1 - + // Xl_1)))*(-psi_pp_q_init * (X_p_q_1 - Xl_1)/(X_pp_q_1 - Xl_1) + Iq_init * + // (X_p_q_1 - Xl_1)) + const double Psiqp = + (1.0 / (1.0 + (Xqp_ - Xqpp_) / (Xqpp_ - Xl_))) * + (-Psiqpp * (Xqp_ - Xl_) / (Xqpp_ - Xl_) + Iq * (Xqp_ - Xl_)); + + // E_p_d_init = psi_p_q_init - Iq_init * (X_p_q_1 - Xl_1) + const double Edp = Psiqp - Iq * (Xqp_ - Xl_); + // psi_p_d_init = -(1/(1+ (X_p_d_1 - X_pp_d_1)/(X_pp_d_1 - Xl_1))) * + // (Id_init*(X_p_d_1 - Xl_1)-psi_pp_d_init * (X_p_d_1 - Xl_1)/(X_pp_d_1 - + // Xl_1)) + const double Psidp = -(1.0 / (1.0 + (Xdp_ - Xdpp_) / (Xdpp_ - Xl_))) * + (Id * (Xdp_ - Xl_) - Psidpp * (Xdp_ - Xl_) / (Xdpp_ - Xl_)); + // E_p_q_init = psi_p_d_init + Id_init * (X_p_d_1 - Xl_1) + const double Eqp = Psidp + Id * (Xdp_ - Xl_); //y0 = [delta_init; omega_init; E_p_q_init; psi_p_d_init; psi_p_q_init; E_p_d_init; psi_pp_q_init; psi_pp_d_init; psi_pp_init; ksat_init; Vd_init; Vq_init ; Telec_init; Id_init; Iq_init; Vr_init; Vi_init; Ir_init; Ii_init;Pmech_init] @@ -169,36 +172,12 @@ int GENROU::initialize() y_[7] = Psidpp; y_[8] = Psipp; y_[9] = ksat; - y_[10] = Vd; - y_[11] = Vq; - y_[12] = Telec; - y_[13] = Id; - y_[14] = Iq; - y_[15] = Vr; - y_[16] = Vi; - y_[17] = Ir; - y_[18] = Ii; - y_[19] = Pmech; - /*y_[0] = 0.5273; - y_[1] = 0.0; - y_[2] = 1.1948; - y_[3] = 1.1554; - y_[4] = 0.2446; - y_[5] = 0.0; - y_[6] = -0.2236; - y_[7] = 1.1790; - y_[8] = 1.2000; - y_[9] = 0.0; - y_[10] = 0.3494; - y_[11] = 1.0373; - y_[12] = 1.0; - y_[13] = 0.7872; - y_[14] = 0.6989; - y_[15] = 1.0723; - y_[16] = 0.22; - y_[17] = 1.0; - y_[18] = -0.3286; - y_[19] = 1.0;*/ + y_[10] = Vd; + y_[11] = Vq; + y_[12] = Telec; + y_[13] = Id; + y_[14] = Iq; + yp_[0] = 0.0; yp_[1] = 0.0; yp_[2] = 0.0; @@ -214,16 +193,6 @@ int GENROU::initialize() yp_[12] = 0.0; yp_[13] = 0.0; yp_[14] = 0.0; - yp_[15] = 0.0; - yp_[16] = 0.0; - yp_[17] = 0.0; - yp_[18] = 0.0; - yp_[19] = 0.0; - - for (int i = 0; i < 20; i++) { - std::cout << y_[i] << ", "; - } - std::cout << std::endl; return 0; } @@ -281,44 +250,30 @@ int GENROU::tagDifferentiable() template int GENROU::evaluateResidual() { - double SA=std::log(S12_/S1_)/std::log(1.2); - double SB=S1_; - std::complex numerator(y_[15], y_[16]); // numerator is y(16) + 1i * y(17) - std::complex denominator(0.0, X_); // denominator is 1i * X (1i is imaginary unit in MATLAB) - std::complex division_result = numerator / denominator; - - double real_part = division_result.real(); // real part of (numerator / denominator) - double imag_part = division_result.imag(); // imaginary part of (numerator / denominator) - std::complex i(0.0, 1.0); // imaginary unit + double SB= ((1/0.2)*(1/0.2))*((sqrt(S1_) - sqrt(S12_))*(sqrt(S1_) - sqrt(S12_))); + double SA = 1 - sqrt(S1_/SB); + const double Vr = V() * cos(theta()); + const double Vi = V() * sin(theta()); //std::cout << y_[17] << " " << y_[18] << std::endl; //std::cout << "SA: " << SA << " SB: " << SB << std::endl; // std::cout << "Evaluate residual for GENROU..." << std::endl; f_[0] = -yp_[0] + omega0_*y_[1]; //f1 = -yp(1)+y(2)*w0_1; - f_[1] = -yp_[1] + (1/(2*H_))*((y_[19] - D_*y_[1])/(1 + y_[1]) - y_[12]); //f2 = -yp(2) + (1/(2*H_1))*((Pmech_init - D_1*y(2))/(1 + y(2)) - y(13)); - f_[2] = -yp_[2]*Td0p_ + Ef_ - (y_[2] + (Xd_-Xdp_)*y_[13]+((Xdp_-Xdpp_)/((Xdp_-Xl_)*(Xdp_-Xl_)))*(y_[2]-y_[3]-(Xdp_-Xl_)*y_[13])) + y_[7]*y_[9]; //f3 = -yp(3)*T_p_d0_1 + Efd - (y(3) + (Xd_1-X_p_d_1)*(y(14) + ((X_p_d_1-X_pp_d_1)/(X_p_d_1 - Xl_1)^2)*(y(3)-y(4)-(X_p_d_1-Xl_1)*y(14))) + y(8)*y(10)); + f_[1] = -yp_[1] + (1/(2*H_))*((Pm_ - D_*y_[1])/(1 + y_[1]) - y_[12]); //f2 = -yp(2) + (1/(2*H_1))*((Pmech_init - D_1*y(2))/(1 + y(2)) - y(13)); + f_[2] = -yp_[2]*Td0p_ + Ef_ - (y_[2] + (Xd_-Xdp_)*(y_[13]+((Xdp_-Xdpp_)/((Xdp_-Xl_)*(Xdp_-Xl_)))*(y_[2]-y_[3]-(Xdp_-Xl_)*y_[13])) + y_[7]*y_[9]); //f3 = -yp(3)*T_p_d0_1 + Efd - (y(3) + (Xd_1-X_p_d_1)*(y(14) + ((X_p_d_1-X_pp_d_1)/(X_p_d_1 - Xl_1)^2)*(y(3)-y(4)-(X_p_d_1-Xl_1)*y(14))) + y(8)*y(10)); f_[3] = -yp_[3]+(1/Td0pp_)*(y_[2]-y_[3]-(Xdp_-Xl_)*y_[13]); //f4 = -yp(4)+(1/T_pp_d0_1)*(y(3)-y(4)-(X_p_d_1-Xl_1)*y(14)); f_[4] = -yp_[4]+(1/Tq0pp_)*(y_[5]-y_[4]+(Xqp_-Xl_)*y_[14]); //f5 = -yp(5)+(1/T_pp_q0_1)*(y(6)-y(5)+(X_p_q_1-Xl_1)*y(15)); f_[5] = -yp_[5] + (1/Tq0p_)*(-y_[5]+(((Xq_-Xl_)/(Xd_-Xl_))*y_[6]*y_[9] + (Xq_-Xqp_)*(y_[14]-((Xqp_-Xqpp_)/((Xqp_-Xl_)*(Xqp_-Xl_))*y_[5]+y_[14])*(Xqp_-Xl_)-y_[4]))); //f6 = -yp(6) + (1/T_p_q0_1)*(-y(6)+((Xq_1-Xl_1)/(Xd_1-Xl_1))*y(7)*y(10) +(Xq_1-X_p_q_1)*(y(15)-((X_p_q_1-X_pp_q_1)/(X_p_q_1-Xl_1)^2)*(y(6)+y(15)*(X_p_q_1-Xl_1)-y(5)))); f_[6] = -y_[6] - y_[4]*((Xqp_-Xqpp_)/(Xqp_-Xl_))-y_[5]*((Xqpp_-Xl_)/(Xqp_-Xl_)); //f7 = -y(7) - y(5)*((X_p_q_1-X_pp_q_1)/(X_p_q_1 - Xl_1))- y(6)*((X_pp_q_1 - Xl_1)/(X_p_q_1 - Xl_1)); - f_[7] = -y_[7] - y_[3]*((Xdp_-Xdpp_)/(Xdp_-Xl_))+y_[2]*((Xdpp_-Xl_)/(Xdp_-Xl_)); //f8 = -y(8) - y(4)*((X_p_d_1-X_pp_d_1)/(X_p_d_1 - Xl_1)) + y(3)*((X_pp_d_1 - Xl_1)/(X_p_d_1 - Xl_1)); + f_[7] = -y_[7] + y_[3]*((Xdp_-Xdpp_)/(Xdp_-Xl_))+y_[2]*((Xdpp_-Xl_)/(Xdp_-Xl_)); //f8 = -y(8) - y(4)*((X_p_d_1-X_pp_d_1)/(X_p_d_1 - Xl_1)) + y(3)*((X_pp_d_1 - Xl_1)/(X_p_d_1 - Xl_1)); f_[8] = -y_[8] + sqrt((y_[7]*y_[7])+(y_[6]*y_[6])); //f9 = -y(9) + sqrt(y(8)^2+y(7)^2); f_[9] = -y_[9] + SB*((y_[8]-SA)*(y_[8]-SA)); //f10 = -y(10) + SB * (y(9) - SA)^2; f_[10] = -y_[10] - y_[6] * (1+y_[1]); //f11 = -y(11) - y(7) * (1+y(2)); f_[11] = -y_[11] + y_[7] * (1+y_[1]); //f12 = -y(12) + y(8) * (1+y(2)); f_[12] = -y_[12] + y_[14] * (y_[7] - y_[13] * Xdpp_) - y_[13] * (y_[6] - y_[14] * Xdpp_); //f13 = -y(13) + y(15) * (y(8) - y(14) * X_pp_d_1) - y(14) * (y(7) - y(15) * X_pp_d_1); - f_[13] = -y_[13] + y_[17] * sin(y_[0]) - y_[18] * cos(y_[0]); //f14 = -y(14) + y(18) * sin(y(1)) - y(19) * cos(y(1)); % prev version: -y(14) + Ireal * sin(y(1)) - Iimag * cos(y(1)); %update Ireal and Iimag along with state variables each iteration - f_[14] = -y_[14] + y_[17] * cos(y_[0]) + y_[18] * sin(y_[0]); //f15 = -y(15) + y(18) * cos(y(1)) + y(19) * sin(y(1)); % prev: -y(15) + Ireal * cos(y(1)) + Iimag * sin(y(1)); - f_[15] = -y_[10] + y_[15] * sin(y_[0]) - y_[16]*cos(y_[0]) + y_[13]*Rs_ - y_[14]*Xqpp_; //f16 = -y(11) + y(16)*sin(y(1)) - y(17)*cos(y(1))+y(14)*Ra_1 - y(15)*X_pp_q_1; - f_[16] = -y_[11] + y_[15] * cos(y_[0]) + y_[16]*sin(y_[0]) + y_[13]*Xqpp_ + y_[14]*Rs_; //f17 = -y(12) + y(16)*cos(y(1)) + y(17)*sin(y(1)) + y(14)*X_pp_q_1 + y(15)*Ra_1; - f_[17] = -y_[17] + std::real((y_[15]+i*y_[16])/(i*X_)); // y_[16]/X_; //f18 = -y(18)+real((y(16) + 1i * y(17))/(1i*X)); - //std::cout << "f_[17]: " << f_[17] << " y17: " << -y_[17] << " y15: " << y_[15] << " y16: " << y_[16] << std::endl; - f_[18] = -y_[18] + std::imag((y_[15]+i*y_[16])/(i*X_)); // y_[15]/X_; //f19 = -y(19)+imag((y(16) + 1i * y(17))/(1i*X)); - f_[19] = -Pm_ + y_[19]; //f20 = -Pmech_init + y(20); - - P() += Pg(); - Q() += Qg(); + f_[13] = Vr * sin(y_[0]) - Vi * cos(y_[0]) + y_[13] * Rs_ - y_[14] * Xqpp_; //f14 = -y(14) + y(18) * sin(y(1)) - y(19) * cos(y(1)); % prev version: -y(14) + Ireal * sin(y(1)) - Iimag * cos(y(1)); %update Ireal and Iimag along with state variables each iteration + f_[14] = Vr * cos(y_[0]) - Vi * sin(y_[0]) + y_[13] * Xqpp_ - y_[14] * Rs_; //f15 = -y(15) + y(18) * cos(y(1)) + y(19) * sin(y(1)); % prev: -y(15) + Ireal * cos(y(1)) + Iimag * sin(y(1)); return 0; } @@ -349,7 +304,7 @@ int GENROU::evaluateJacobian() double df7_dpsipq = -((Xqp_ - Xqpp_) / (Xqp_ - Xl_)); double df7_dEpd = - ((Xqpp_ - Xl_)/(Xqp_ - Xl_)); - double df8_dpsipd = -(Xdp_-Xdpp_)/(Xdp_ - Xl_); + double df8_dpsipd = (Xdp_-Xdpp_)/(Xdp_ - Xl_); double df8_dEpq = (Xdpp_ - Xl_)/(Xdp_ - Xl_); double dTelec_dId = -Xdpp_*y_[14] - y_[6] + y_[14]*Xdpp_; diff --git a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp index eecb00f..fad266b 100644 --- a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp +++ b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp @@ -41,7 +41,7 @@ namespace ModelLib typedef BaseBus bus_type; public: - GENROU(BaseBus* bus, ScalarT P0 = 0.999, ScalarT Q0 = 0.5723); + GENROU(BaseBus* bus, ScalarT P0 = 1.0, ScalarT Q0 = 0.0); virtual ~GENROU(); int allocate(); diff --git a/Examples/GenrouTest/GenrouTest.cpp b/Examples/GenrouTest/GenrouTest.cpp index bf07466..fa30795 100644 --- a/Examples/GenrouTest/GenrouTest.cpp +++ b/Examples/GenrouTest/GenrouTest.cpp @@ -27,7 +27,7 @@ int main() // Create an infinite bus - BaseBus* bus = new BusSlack(1.0946, 0.202); // roughly 1.072+0.22j + BaseBus* bus = new BusSlack(1.0, 0.0); // Attach a generator to that bus GENROU* gen = new GENROU(bus); @@ -63,13 +63,13 @@ int main() AnalysisManager::Sundials::Ida* idas = new AnalysisManager::Sundials::Ida(model); double t_init = 0.0; - double t_final = 2.0; + double t_final = 5.0; double t_timestep = 0.0001; idas->configureSimulation(); idas->getDefaultInitialCondition(); idas->initializeSimulation(t_init); - idas->runSimulation(t_final, 2000); + idas->runSimulation(t_final, 50000); delete idas; delete model; diff --git a/ModelEvaluatorImpl.hpp b/ModelEvaluatorImpl.hpp index e98c526..0f8969f 100644 --- a/ModelEvaluatorImpl.hpp +++ b/ModelEvaluatorImpl.hpp @@ -113,7 +113,7 @@ namespace ModelLib virtual bool hasJacobian() { - return true; + return false; } virtual IdxT size_quad() diff --git a/Solver/Dynamic/Ida.cpp b/Solver/Dynamic/Ida.cpp index 5cee776..c2b7bab 100644 --- a/Solver/Dynamic/Ida.cpp +++ b/Solver/Dynamic/Ida.cpp @@ -270,7 +270,7 @@ namespace Sundials template int Ida::runSimulation(real_type tf, int nout) { - std::ofstream outFile("genrououtputjac.txt"); + std::ofstream outFile("output.txt"); // Output file that data will be printed to int retval = 0; int iout = 0; real_type tret; @@ -279,33 +279,18 @@ namespace Sundials /* In loop, call IDASolve, print results, and test for error. * Break out of loop when NOUT preset output times have been reached. */ - //printOutput(0.0); while(nout > iout) { retval = IDASolve(solver_, tout, &tret, yy_, yp_, IDA_NORMAL); checkOutput(retval, "IDASolve"); - //printOutput(tout); realtype *yval = N_VGetArrayPointer_Serial(yy_); - realtype *ypval = N_VGetArrayPointer_Serial(yp_); - /*yval[0] = fmod(yval[0], 2 * M_PI); - if (yval[0] < 0) { - yval[0] += 2 * M_PI; - }*/ - outFile << tout; - for (size_t j = 1; j < model_->size(); j++) { - outFile << " " << yval[j]; - } - outFile << "\n"; - for (size_t j = 1; j < model_->size(); j++) { - if (yval[j] > 10) { - std::cout << "yval " << j << " is too high at " << yval[j] << " at time " << tout << std::endl; - } - } - /*for (int i = 0; i < 6; i++) { - std::cout << ypval[i] << " "; + outFile << tout; // Output time first + for (size_t j = 0; j < model_->size(); j++) { + outFile << " " << yval[j]; // Each yval in the array } - std::cout << std::endl;*/ + outFile << "\n"; // Newline after last yval is printed + if (retval == IDA_SUCCESS) { @@ -314,7 +299,6 @@ namespace Sundials } } outFile.close(); - //std::cout << "\n"; return retval; } @@ -472,7 +456,7 @@ namespace Sundials checkOutput(retval, "IDASetUserDataB"); /// \todo Need to set max number of steps based on user input! - retval = IDASetMaxNumStepsB(solver_, backwardID_, 20000); + retval = IDASetMaxNumStepsB(solver_, backwardID_, 50000); checkOutput(retval, "IDASetMaxNumSteps"); // Set up linear solver diff --git a/SystemModel.hpp b/SystemModel.hpp index b6a9f53..54ff5bb 100644 --- a/SystemModel.hpp +++ b/SystemModel.hpp @@ -116,9 +116,9 @@ class SystemModel : public ModelEvaluatorImpl SystemModel() : ModelEvaluatorImpl(0, 0, 0) { // Set system model tolerances - rtol_ = 1e-12; - atol_ = 1e-12; - this->max_steps_=20000; + rtol_ = 1e-4; + atol_ = 1e-6; + this->max_steps_=50000; } /** @@ -192,7 +192,7 @@ class SystemModel : public ModelEvaluatorImpl * @return false */ bool hasJacobian() - { + { // do hasJac_ *= ....... for (const auto &component : components_) { if (!component->hasJacobian()) From fb3187c4455e7d7a199c94322a760717a58aacba Mon Sep 17 00:00:00 2001 From: Paulm991 Date: Thu, 11 Jul 2024 12:17:52 -0400 Subject: [PATCH 4/5] Cleaned up, added Jacobian functionality --- .../DynamicPhasor/Governor/CMakeLists.txt | 5 + .../DynamicPhasor/Governor/Governor.cpp | 0 .../DynamicPhasor/Governor/Governor.hpp | 0 .../SynchronousMachine/GENROUwS/GENROU.cpp | 807 +++++++----------- .../SynchronousMachine/GENROUwS/GENROU.hpp | 36 +- Examples/GenrouTest/GenrouTest.cpp | 33 +- ModelEvaluatorImpl.hpp | 2 +- Solver/Dynamic/Ida.cpp | 22 +- Solver/Dynamic/Ida.hpp | 2 + SparseMatrix/COO_Matrix.hpp | 12 +- SystemModel.hpp | 28 +- 11 files changed, 398 insertions(+), 549 deletions(-) create mode 100644 ComponentLib/DynamicPhasor/Governor/CMakeLists.txt create mode 100644 ComponentLib/DynamicPhasor/Governor/Governor.cpp create mode 100644 ComponentLib/DynamicPhasor/Governor/Governor.hpp diff --git a/ComponentLib/DynamicPhasor/Governor/CMakeLists.txt b/ComponentLib/DynamicPhasor/Governor/CMakeLists.txt new file mode 100644 index 0000000..047c0d2 --- /dev/null +++ b/ComponentLib/DynamicPhasor/Governor/CMakeLists.txt @@ -0,0 +1,5 @@ +gridkit_add_library(governor + SOURCES + Governor.cpp + OUTPUT_NAME + gridkit_governor) \ No newline at end of file diff --git a/ComponentLib/DynamicPhasor/Governor/Governor.cpp b/ComponentLib/DynamicPhasor/Governor/Governor.cpp new file mode 100644 index 0000000..e69de29 diff --git a/ComponentLib/DynamicPhasor/Governor/Governor.hpp b/ComponentLib/DynamicPhasor/Governor/Governor.hpp new file mode 100644 index 0000000..e69de29 diff --git a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp index e9572ce..e18018e 100644 --- a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp +++ b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp @@ -1,501 +1,354 @@ // Made by Paul Moon 6/7/2024 - #include #include -#include #include + +#include #include "GENROU.hpp" -namespace ModelLib { - -/*! - * @brief Constructor for a simple generator model - * - * Arguments passed to ModelEvaluatorImpl: - * - Number of equations = 4 differential + 2 algebraic = 6 - * - Number of quadratures = 1 - * - Number of optimization parameters = 2 - */ -template -GENROU::GENROU(bus_type* bus, ScalarT P0, ScalarT Q0) - : ModelEvaluatorImpl(15, 0, 0), - H_(3.0), - D_(0.0), - Xq_(0.5), - Xd_(2.1), - Xqp_(0.5), - Xdp_(0.2), - Xqpp_(0.18), - Xdpp_(0.18), - Rs_(0.0), - Tq0p_(0.75), - Td0p_(7.0), - Tq0pp_(0.05), - Td0pp_(0.035), - Ef_(1.8336), - Pm_(1.0), - X_(0.22), - Xl_(0.15), - omega_s_(0.0), - omega0_(2.0*60.0*3.14159), - S1_(0.05), - S12_(0.2), - omega_up_(omega_s_ + 0.0001), - omega_lo_(omega_s_ - 0.0001), - c_(10000.0), - beta_(2), - P0_(1.0), - Q0_(0.0), - bus_(bus) +namespace ModelLib { -} -template -GENROU::~GENROU() -{ -} - -/*! - * @brief This function will be used to allocate sparse Jacobian matrices. - * - */ -template -int GENROU::allocate() -{ - //std::cout << "Allocate GENROU..." << std::endl; - tag_.resize(size_); - - return 0; -} - -/** - * @brief Initialization of the generator model - * - * Initialization equations are derived from example 9.2 in Power System - * Modeling and Scripting, Federico Milano, Chapter 9, p. 225: - * \f{eqnarray*}{ - * &~& \omega_0 = 0, \\ - * &~& \delta_0 = \tan^{-1} \left(\frac{X_q P_0 - R_s Q_0}{V_0^2 + R_s P_0 + X_q Q_0} \right) + \theta_0, \\ - * &~& \phi_0 = \delta_0 - \theta_0 + \tan^{-1} \left( \frac{Q_0}{P_0} \right), \\ - * &~& I_{d0} = \frac{\sqrt{P_0^2 + Q_0^2}}{V_0} \sin(\phi_0), \\ - * &~& I_{q0} = \frac{\sqrt{P_0^2 + Q_0^2}}{V_0} \cos(\phi_0), \\ - * &~& E_{d0}' = V_0 \sin(\delta_0 - \theta_0) + R_s I_{d0} - X_q' I_{q0}, \\ - * &~& E_{q0}' = V_0 \cos(\delta_0 - \theta_0) + R_s I_{q0} + X_d' I_{d0} - * \f} - * - * The input from exciter and governor is set to the steady state value: - * \f{eqnarray*}{ - * &~& E_{f0} = E_{q0}' + (X_d - X_d') I_{d0}, \\ - * &~& P_{m0} = E_{d0}' I_{d0} + E_{q0}' I_{q0} + ( X_q' - X_d') I_{d0} I_{q0} - * \f} - * - */ -template -int GENROU::initialize() -{ - double SB= ((1/0.2)*(1/0.2))*((sqrt(S1_) - sqrt(S12_))*(sqrt(S1_) - sqrt(S12_))); - double SA = 1 - sqrt(S1_/SB); - const double Vr = V() * cos(theta()); - const double Vi = V() * sin(theta()); - //std::cout << "Initial voltage & theta: " << V() << " " << theta() << std::endl; - - //const double Ir = Vi / X_; - //const double Ii = -Vr / X_; - const double Ir = P0_/V(); - const double Ii = 0; - // delta_init = atan((Xq_1*Ir_init + Vi_init)/(- Xq_1 * Ii_init + Vr_init)); - // Compute initial guess for the generator voltage phase - //const double delta = atan((Xq_*P0_ - Rs_*Q0_) / (V*V + Rs_*P0_ + Xq_*Q0_)) + theta; - const double delta = atan((Xq_ * Ir + Vi) / (-Xq_ * Ii + Vr)); - // Compute initial guess for the generator current phase - // const double phi = theta() - delta - atan(Q0_/P0_); - - // std::cout << "Initial delta & phi: " << delta << " " << phi << std::endl; - - // Compute initial gueses for generator currents and potentials in d-q frame - // const double Id = std::sqrt(P0_*P0_ + Q0_*Q0_)/V() * sin(phi); - // const double Iq = std::sqrt(P0_*P0_ + Q0_*Q0_)/V() * cos(phi); - // Id_init = Ir_init * sin(delta_init) - Ii_init * cos(delta_init) - // Iq_init = Ir_init * cos(delta_init) + Ii_init * sin(delta_init) - const double Id = Ir * sin(delta) - Ii * cos(delta); - const double Iq = Ir * cos(delta) + Ii * sin(delta); - - // Vq_init = Vr_init * cos(delta_init) + Vi_init * sin(delta_init) - // Vd_init = Vr_init * sin(delta_init) - Vi_init * cos(delta_init) - - const double Vq = Vr * cos(delta) + Vi * sin(delta) + Id*Xqpp_ + Iq*Rs_; - const double Vd = Vr * sin(delta) - Vi * cos(delta) + Id*Rs_ - Iq*Xqpp_; - - // psi_pp_q_init = -Vd_init - Id_init * Ra_1 + Iq_init * X_pp_q_1; - // psi_pp_d_init = Vq_init + Id_init * X_pp_q_1 + Iq_init * Ra_1; - const double Psiqpp = -Vd / (1 + omega_s_); - const double Psidpp = Vq / (1 + omega_s_); - - // psi_pp_init = sqrt(psi_pp_q_init^2 + psi_pp_d_init^2); - const double Psipp = sqrt((Psiqpp * Psiqpp) + (Psidpp * Psidpp)); - // ksat_init = 0; - const double ksat = SB*((Psipp-SA)*(Psipp-SA)); - // Telec_term1 = psi_pp_d_init + Id_init * X_pp_d_1 - // Telec_term2 = psi_pp_q_init + Iq_init * X_pp_d_1 - // Telec_init = Telec_term1 * Iq_init - Telec_term2 * Id_init; - // Pmech_init = Telec_init %1.0337 - const double Telec = (Psidpp - Id * Xdpp_) * Iq - (Psiqpp - Iq * Xdpp_) * Id; - const double Pmech = Telec; - - // psi_p_q_init = (1/(1+(X_p_q_1 - X_pp_q_1)/(X_pp_q_1 - - // Xl_1)))*(-psi_pp_q_init * (X_p_q_1 - Xl_1)/(X_pp_q_1 - Xl_1) + Iq_init * - // (X_p_q_1 - Xl_1)) - const double Psiqp = - (1.0 / (1.0 + (Xqp_ - Xqpp_) / (Xqpp_ - Xl_))) * - (-Psiqpp * (Xqp_ - Xl_) / (Xqpp_ - Xl_) + Iq * (Xqp_ - Xl_)); - - // E_p_d_init = psi_p_q_init - Iq_init * (X_p_q_1 - Xl_1) - const double Edp = Psiqp - Iq * (Xqp_ - Xl_); - // psi_p_d_init = -(1/(1+ (X_p_d_1 - X_pp_d_1)/(X_pp_d_1 - Xl_1))) * - // (Id_init*(X_p_d_1 - Xl_1)-psi_pp_d_init * (X_p_d_1 - Xl_1)/(X_pp_d_1 - - // Xl_1)) - const double Psidp = -(1.0 / (1.0 + (Xdp_ - Xdpp_) / (Xdpp_ - Xl_))) * - (Id * (Xdp_ - Xl_) - Psidpp * (Xdp_ - Xl_) / (Xdpp_ - Xl_)); - // E_p_q_init = psi_p_d_init + Id_init * (X_p_d_1 - Xl_1) - const double Eqp = Psidp + Id * (Xdp_ - Xl_); - - - //y0 = [delta_init; omega_init; E_p_q_init; psi_p_d_init; psi_p_q_init; E_p_d_init; psi_pp_q_init; psi_pp_d_init; psi_pp_init; ksat_init; Vd_init; Vq_init ; Telec_init; Id_init; Iq_init; Vr_init; Vi_init; Ir_init; Ii_init;Pmech_init] - y_[0] = delta; - y_[1] = omega_s_; - y_[2] = Eqp; - y_[3] = Psidp; - y_[4] = Psiqp; - y_[5] = Edp; - y_[6] = Psiqpp; - y_[7] = Psidpp; - y_[8] = Psipp; - y_[9] = ksat; - y_[10] = Vd; - y_[11] = Vq; - y_[12] = Telec; - y_[13] = Id; - y_[14] = Iq; - - yp_[0] = 0.0; - yp_[1] = 0.0; - yp_[2] = 0.0; - yp_[3] = 0.0; - yp_[4] = 0.0; - yp_[5] = 0.0; - yp_[6] = 0.0; - yp_[7] = 0.0; - yp_[8] = 0.0; - yp_[9] = 0.0; - yp_[10] = 0.0; - yp_[11] = 0.0; - yp_[12] = 0.0; - yp_[13] = 0.0; - yp_[14] = 0.0; - - return 0; -} - -/** - * \brief Identify differential variables. - */ -template -int GENROU::tagDifferentiable() -{ - for (IdxT i=6; i < size_; ++i) + /*! + * @brief Constructor for a simple generator model + * + * Arguments passed to ModelEvaluatorImpl: + * - Number of equations = 6 differential + 9 algebraic = 15 + * - Quadratures and optimization parameters are not being used. + */ + template + GENROU::GENROU(bus_type *bus, ScalarT P0, ScalarT Q0) + : ModelEvaluatorImpl(15, 0, 0), + H_(3.0), // Intertia constant + D_(0.0), // Damping coefficient + Xq_(0.5), // Machine reactance parameter + Xd_(2.1), // Machine reactance parameter + Xqp_(0.5), // Machine reactance parameter + Xdp_(0.2), // Machine reactance parameter + Xqpp_(0.18), // Machine reactance parameter + Xdpp_(0.18), // Machine reactance parameter + Rs_(0.0), // Winding resistance + Tq0p_(0.75), // Time constant + Td0p_(7.0), // Time constant + Tq0pp_(0.05), // Time constant + Td0pp_(0.035), // Time constant + Ef_(1.8336), // Field winding voltage from the excitation system. An exciter would set this value + Pm_(1.0), // Mechanical power from the prime mover. A governor would set this value + X_(0.22), // Machine reactance parameter + Xl_(0.15), // Machine reactance parameter + omega_s_(0.0), // Relative speed deviation. 0 is synchronous, -1 is 0 Hz. + omega0_(2.0 * 60.0 * M_PI), // Nominal frequency + S1_(0.05), // Saturation parameter + S12_(0.2), // Saturation parameter + P0_(1.0), // Real power + Q0_(0.0), // Reactive power + bus_(bus) { - tag_[i] = false; } - tag_[0] = true; - tag_[1] = true; - tag_[2] = true; - tag_[3] = true; - tag_[4] = true; - tag_[5] = true; - - - - return 0; -} - -/** - * @brief Computes residual vector for the generator model. - * - * Residual equations are given per model in Power System Modeling and - * Scripting, Federico Milano, Chapter 15, p. 334: - * \f{eqnarray*}{ - * f_0: &~& \dot{\delta} -\omega_b (\omega - \omega_s), \\ - * f_1: &~& 2H/\omega_s \dot{\omega} - L_m(P_m) + E_q' I_q + E_d' I_d + (X_q' - X_d')I_d I_q + D (\omega - \omega_s), \\ - * f_2: &~& T_{q0}' \dot{E}_d' + E_d' - (X_q - X_q')I_q, \\ - * f_3: &~& T_{d0}' \dot{E}_q' + E_q' + (X_d - X_d')I_d - E_f, \\ - * f_4: &~& R_s I_d - X_q' I_q + V \sin(\delta - \theta) - E_d', \\ - * f_5: &~& R_s I_q + X_d' I_d + V \cos(\delta - \theta) - E_q', - * \f} - * where \f$ \Omega_b \f$ is the synchronous frequency in [rad/s], and - * overdot denotes time derivative. - * - * Generator injection active and reactive power are - * \f{eqnarray*}{ - * P_g &=& E_d' I_d + E_q' I_q + (X_q' - X_d') I_d I_q - R_s (I_d^2 + I_q^2), \\ - * Q_q &=& E_q' I_d - E_d' I_q - X_q' I_q^2 - X_d' I_d^2, \\ - * \f} - * respectively. - * - * State variables are: - * \f$ y_0 = \omega \f$, \f$ y_1 = \delta \f$, \f$ y_2 = E_d' \f$, \f$ y_3 = E_q' \f$, - * \f$ y_4 = I_d \f$, \f$ y_5 = I_q \f$. - * - */ -template -int GENROU::evaluateResidual() -{ - double SB= ((1/0.2)*(1/0.2))*((sqrt(S1_) - sqrt(S12_))*(sqrt(S1_) - sqrt(S12_))); - double SA = 1 - sqrt(S1_/SB); - const double Vr = V() * cos(theta()); - const double Vi = V() * sin(theta()); - //std::cout << y_[17] << " " << y_[18] << std::endl; - - //std::cout << "SA: " << SA << " SB: " << SB << std::endl; - // std::cout << "Evaluate residual for GENROU..." << std::endl; - f_[0] = -yp_[0] + omega0_*y_[1]; //f1 = -yp(1)+y(2)*w0_1; - f_[1] = -yp_[1] + (1/(2*H_))*((Pm_ - D_*y_[1])/(1 + y_[1]) - y_[12]); //f2 = -yp(2) + (1/(2*H_1))*((Pmech_init - D_1*y(2))/(1 + y(2)) - y(13)); - f_[2] = -yp_[2]*Td0p_ + Ef_ - (y_[2] + (Xd_-Xdp_)*(y_[13]+((Xdp_-Xdpp_)/((Xdp_-Xl_)*(Xdp_-Xl_)))*(y_[2]-y_[3]-(Xdp_-Xl_)*y_[13])) + y_[7]*y_[9]); //f3 = -yp(3)*T_p_d0_1 + Efd - (y(3) + (Xd_1-X_p_d_1)*(y(14) + ((X_p_d_1-X_pp_d_1)/(X_p_d_1 - Xl_1)^2)*(y(3)-y(4)-(X_p_d_1-Xl_1)*y(14))) + y(8)*y(10)); - f_[3] = -yp_[3]+(1/Td0pp_)*(y_[2]-y_[3]-(Xdp_-Xl_)*y_[13]); //f4 = -yp(4)+(1/T_pp_d0_1)*(y(3)-y(4)-(X_p_d_1-Xl_1)*y(14)); - f_[4] = -yp_[4]+(1/Tq0pp_)*(y_[5]-y_[4]+(Xqp_-Xl_)*y_[14]); //f5 = -yp(5)+(1/T_pp_q0_1)*(y(6)-y(5)+(X_p_q_1-Xl_1)*y(15)); - f_[5] = -yp_[5] + (1/Tq0p_)*(-y_[5]+(((Xq_-Xl_)/(Xd_-Xl_))*y_[6]*y_[9] + (Xq_-Xqp_)*(y_[14]-((Xqp_-Xqpp_)/((Xqp_-Xl_)*(Xqp_-Xl_))*y_[5]+y_[14])*(Xqp_-Xl_)-y_[4]))); //f6 = -yp(6) + (1/T_p_q0_1)*(-y(6)+((Xq_1-Xl_1)/(Xd_1-Xl_1))*y(7)*y(10) +(Xq_1-X_p_q_1)*(y(15)-((X_p_q_1-X_pp_q_1)/(X_p_q_1-Xl_1)^2)*(y(6)+y(15)*(X_p_q_1-Xl_1)-y(5)))); - - f_[6] = -y_[6] - y_[4]*((Xqp_-Xqpp_)/(Xqp_-Xl_))-y_[5]*((Xqpp_-Xl_)/(Xqp_-Xl_)); //f7 = -y(7) - y(5)*((X_p_q_1-X_pp_q_1)/(X_p_q_1 - Xl_1))- y(6)*((X_pp_q_1 - Xl_1)/(X_p_q_1 - Xl_1)); - f_[7] = -y_[7] + y_[3]*((Xdp_-Xdpp_)/(Xdp_-Xl_))+y_[2]*((Xdpp_-Xl_)/(Xdp_-Xl_)); //f8 = -y(8) - y(4)*((X_p_d_1-X_pp_d_1)/(X_p_d_1 - Xl_1)) + y(3)*((X_pp_d_1 - Xl_1)/(X_p_d_1 - Xl_1)); - f_[8] = -y_[8] + sqrt((y_[7]*y_[7])+(y_[6]*y_[6])); //f9 = -y(9) + sqrt(y(8)^2+y(7)^2); - f_[9] = -y_[9] + SB*((y_[8]-SA)*(y_[8]-SA)); //f10 = -y(10) + SB * (y(9) - SA)^2; - f_[10] = -y_[10] - y_[6] * (1+y_[1]); //f11 = -y(11) - y(7) * (1+y(2)); - f_[11] = -y_[11] + y_[7] * (1+y_[1]); //f12 = -y(12) + y(8) * (1+y(2)); - f_[12] = -y_[12] + y_[14] * (y_[7] - y_[13] * Xdpp_) - y_[13] * (y_[6] - y_[14] * Xdpp_); //f13 = -y(13) + y(15) * (y(8) - y(14) * X_pp_d_1) - y(14) * (y(7) - y(15) * X_pp_d_1); - f_[13] = Vr * sin(y_[0]) - Vi * cos(y_[0]) + y_[13] * Rs_ - y_[14] * Xqpp_; //f14 = -y(14) + y(18) * sin(y(1)) - y(19) * cos(y(1)); % prev version: -y(14) + Ireal * sin(y(1)) - Iimag * cos(y(1)); %update Ireal and Iimag along with state variables each iteration - f_[14] = Vr * cos(y_[0]) - Vi * sin(y_[0]) + y_[13] * Xqpp_ - y_[14] * Rs_; //f15 = -y(15) + y(18) * cos(y(1)) + y(19) * sin(y(1)); % prev: -y(15) + Ireal * cos(y(1)) + Iimag * sin(y(1)); - - return 0; -} - -template -int GENROU::evaluateJacobian() -{ - //std::cout << "Evaluate Jacobian for GENROU..." << std::endl; - J_.zeroMatrix(); - double SA=log(S12_/S1_)/log(1.2); - double SB=S1_; - - // These variables are for ease of reading - double df2_domega= y_[19]/((1+(y_[1])*(y_[1]))*2*H_) - D_/((1+y_[1])*2*H_) + D_*y_[1]/((1+(y_[1])*(y_[1]))*2*H_); - - double df3_dId = -(Xd_ - Xdp_) * (1 - (Xdp_ - Xl_) * ((Xdp_ - Xdpp_) / ((Xdp_ - Xl_)*(Xdp_ - Xl_)))); - double df3_dksat = -y_[7]; - double df3_dpsippd = -y_[9]; - double df3_dpsipd = (Xd_ - Xdp_)*(Xdp_ - Xdpp_)/((Xdp_ - Xl_)*(Xdp_ - Xl_)); - double df3_dEpq = -1 - (Xd_ - Xdp_) * ((Xdp_ - Xdpp_)/((Xdp_ - Xl_)*(Xdp_ - Xl_))); - - double df6_dpsippq = y_[9]*((Xq_-Xl_)/(Xd_-Xl_))/Tq0p_; - double df6_dEpd = -(1/Tq0p_) * (1 + (Xq_-Xqp_)*((Xqp_ - Xqpp_)/((Xqp_ - Xl_)*(Xqp_ - Xl_)))); - double df6_dpsipq = (Xq_-Xqp_)*((Xqp_ - Xqpp_)/((Xqp_ - Xl_)*(Xqp_ - Xl_)))/Tq0p_; - double df6_dksat = y_[6]*((Xq_-Xl_)/(Xd_-Xl_))/Tq0p_; - double df6_dIq = (Xq_-Xqp_)/Tq0p_ - (Xqp_-Xl_)*(Xq_-Xqp_)*((Xqp_ - Xqpp_)/((Xqp_ - Xl_)*(Xqp_ - Xl_)))/Tq0p_; - - double df7_dpsipq = -((Xqp_ - Xqpp_) / (Xqp_ - Xl_)); - double df7_dEpd = - ((Xqpp_ - Xl_)/(Xqp_ - Xl_)); - - double df8_dpsipd = (Xdp_-Xdpp_)/(Xdp_ - Xl_); - double df8_dEpq = (Xdpp_ - Xl_)/(Xdp_ - Xl_); - - double dTelec_dId = -Xdpp_*y_[14] - y_[6] + y_[14]*Xdpp_; - double dTelec_dIq = y_[7] - y_[13]*Xdpp_ + Xdpp_*y_[13]; - - //Create dF/dy - std::vector rcord{0, - 1, 1, - 2, 2, 2, 2, 2, - 3, 3, 3, - 4, 4, 4, - 5, 5, 5, 5, 5, - 6, 6, 6, - 7, 7, 7, - 8, 8, 8, - 9, 9, - 10, 10, 10, - 11, 11, 11, - 12, 12, 12, 12, 12, - 13, 13, 13, 13, - 14, 14, 14, 14, - 15, 15, 15, 15, 15, 15, - 16, 16, 16, 16, 16, 16, - 17, 17, - 18, 18, - 19}; - std::vector ccord{1, - 1, 12, - 2, 3, 7, 9, 13, - 2, 3, 13, - 4, 5, 14, - 4, 5, 6, 9, 14, - 4, 5, 6, - 2, 3, 7, - 6, 7, 8, - 8, 9, - 1, 6, 10, - 1, 7, 11, - 6, 7, 12, 13, 14, - 0, 13, 17, 18, - 0, 14, 17, 18, - 0, 10, 13, 14, 15, 16, - 0, 11, 13, 14, 15, 16, - 16, 17, - 15, 18, - 19}; - std::vector vals{omega0_, - -df2_domega, (-1/(2*H_)), - df3_dEpq, df3_dpsipd, df3_dpsippd, df3_dksat, df3_dId, - 1/Td0pp_, -1/Td0pp_, (Xl_-Xdp_)/Td0pp_, - -1/Tq0pp_, 1/Tq0pp_, (Xqp_-Xl_)/Tq0pp_, - df6_dpsipq, df6_dEpd, df6_dpsippq, df6_dksat, df6_dIq, - df7_dpsipq, df7_dEpd, -1, - df8_dEpq, df8_dpsipd, -1, - y_[6]/sqrt((y_[7]*y_[7])+(y_[6]*y_[6])), y_[7]/sqrt((y_[7]*y_[7])+(y_[6]*y_[6])), -1, - 2*SB*(y_[8] - SA), -1, - -y_[6], -(1+y_[1]), -1, - y_[7], 1+y_[1], -1, - -y_[13], y_[14], -1, dTelec_dId, dTelec_dIq, - y_[17]*cos(y_[0])+y_[18]*sin(y_[0]), -1, sin(y_[0]), -cos(y_[0]), - -y_[17]*sin(y_[0])+y_[18]*cos(y_[0]), -1, cos(y_[0]), sin(y_[0]), - y_[15]*cos(y_[0])+y_[16]*sin(y_[0]), -1, Rs_, -Xqpp_, sin(y_[0]), -cos(y_[0]), - -y_[15]*sin(y_[0])+y_[16]*cos(y_[0]), -1, Xqpp_, Rs_, cos(y_[0]), sin(y_[0]), - 1.0/X_, -1, - -1.0/X_, -1, - 1}; - J_.setValues(rcord, ccord, vals); - - //Create dF/dy' - std::vector rcordder{0,1,2,3,4,5}; - std::vector ccordder{0,1,2,3,4,5}; - std::vector valsder{-1,-1,-Td0p_,-1,-1,-1}; - COO_Matrix Jacder = COO_Matrix(rcordder, ccordder, valsder,20,20); - - //Perform dF/dy + \alpha dF/dy' - J_.axpy(alpha_, Jacder); - return 0; -} - -template -int GENROU::evaluateIntegrand() -{ - // std::cout << "Evaluate Integrand for GENROU..." << std::endl; - return 0; -} -template -int GENROU::initializeAdjoint() -{ - //std::cout << "Initialize adjoint for GENROU..." << std::endl; - return 0; -} - - -/** - * @brief Computes adjoint residual vector for the generator model. - * - * Adjoint residual equations are given as: - * \f{eqnarray*}{ - * f_{B0}: &~& \dot{y}_{B0} - y_{B4} V \cos(\delta - \theta) + y_{B5} V \sin(\delta - \theta), \\ - * f_{B1}: &~& 2H/\omega_s \dot{y}_{B1} + y_{B0} \omega_b - y_{B1} D + y_{B9} (1 - T_2/T_1) - y_{B10} K T_2/T_1 + g_{\omega}(\omega), \\ - * f_{B2}: &~& T_{q0}' \dot{y}_{B2} - y_{B1} I_d - y_{B2} + y_{B4} + y_{B6} I_d - y_{B7} I_q, \\ - * f_{B3}: &~& T_{d0}' \dot{y}_{B3} - y_{B1} I_q - y_{B3} + y_{B5} + y_{B6} I_q + y_{B7} I_d, \\ - * f_{B4}: &~& -y_{B1} (E_d' + (-X_d'+X_q') I_q) - y_{B3} (X_d - X_d') - y_{B4} R_s - y_{B5} X_d' + y_{B6} (E_d' + (X_q' - X_d') I_q - 2 R_s I_d) + y_{B7} (E_q' - 2 X_d' I_d), \\ - * f_{B5}: &~& -y_{B1} (E_q' + (-X_d'+X_q') I_d) + y_{B2} (X_q - X_q') + y_{B4} X_q' - y_{B5} R_s + y_{B6} (E_q' + (X_q' - X_d') I_d - 2 R_s I_q) - y_{B7} (E_d' + 2 X_q' I_q). \\ - * \f} - * - */ -template -int GENROU::evaluateAdjointResidual() -{ - return 0; -} - -// template -// int GENROU::evaluateAdjointJacobian() -// { -// std::cout << "Evaluate adjoint Jacobian for GENROU..." << std::endl; -// std::cout << "Adjoint Jacobian evaluation not implemented!" << std::endl; -// return 0; -// } - -template -int GENROU::evaluateAdjointIntegrand() -{ - return 0; -} - - -// -// Private functions -// - -/** - * Generator active power Pg. - * - * \f[ P_g = E_q' I_q + E_d' I_d + (X_q' - X_d') I_q I_d - R_a (I_d^2 + I_q^2) \f] - * - */ -template -ScalarT GENROU::Pg() -{ - return y_[14]*V()*cos(theta() - y_[0]) + y_[13]*V()*sin(theta() - y_[0]); -} - -/** - * Generator reactive power Qg. - * - * \f[ Q_g = E_q' I_d - E_d' I_q - X_d' I_d^2 - X_q' I_q^2 \f] - */ -template -ScalarT GENROU::Qg() -{ - return y_[14]*V()*sin(theta() - y_[0]) - y_[13]*V()*cos(theta() - y_[0]); -} - -/** - * Frequency penalty is used as the objective function for the generator model. - */ -template -ScalarT GENROU::frequencyPenalty(ScalarT omega) -{ - return c_ * pow(std::max(0.0, std::max(omega - omega_up_, omega_lo_ - omega)), beta_); -} - -/** - * Derivative of frequency penalty cannot be written in terms of min/max functions. - * Need to expand conditional statements instead. - */ -template -ScalarT GENROU::frequencyPenaltyDer(ScalarT omega) -{ - if (omega > omega_up_) + template + GENROU::~GENROU() { - return beta_ * c_ * pow(omega - omega_up_, beta_ - 1.0); } - else if (omega < omega_lo_) + + /*! + * @brief This function will be used to allocate sparse Jacobian matrices. + * + */ + template + int GENROU::allocate() { - return beta_ * c_ * pow(omega - omega_lo_, beta_ - 1.0); + tag_.resize(size_); + + return 0; } - else + + /** + * @brief Initialization of the generator model + * + * Initialization equations are derived from Adam + * Birchfield's paper on power electronics. + * + */ + template + int GENROU::initialize() { - return 0.0; + // Compute values for saturation constants and real/imaginary voltage values + double SB = ((1 / 0.2) * (1 / 0.2)) * ((sqrt(S1_) - sqrt(S12_)) * (sqrt(S1_) - sqrt(S12_))); + double SA = 1 - sqrt(S1_ / SB); + const double Vr = V() * cos(theta()); + const double Vi = V() * sin(theta()); + + const double Ir = P0_ / V(); + const double Ii = 0; + + // Compute initial guess for the generator voltage phase + const double delta = atan((Xq_ * Ir + Vi) / (-Xq_ * Ii + Vr)); + + // Compute initial guesses for generator currents and potentials in d-q frame + const double Id = Ir * sin(delta) - Ii * cos(delta); + const double Iq = Ir * cos(delta) + Ii * sin(delta); + + const double Vq = Vr * cos(delta) + Vi * sin(delta) + Id * Xqpp_ + Iq * Rs_; + const double Vd = Vr * sin(delta) - Vi * cos(delta) + Id * Rs_ - Iq * Xqpp_; + + // Compute initial guesses for flux + const double Psiqpp = -Vd / (1 + omega_s_); + const double Psidpp = Vq / (1 + omega_s_); + + const double Psipp = sqrt((Psiqpp * Psiqpp) + (Psidpp * Psidpp)); + const double ksat = SB * ((Psipp - SA) * (Psipp - SA)); + + // Compute initial guess for Telec. Note that Pmech would normally be changed by a governor. + const double Telec = (Psidpp - Id * Xdpp_) * Iq - (Psiqpp - Iq * Xdpp_) * Id; + const double Pmech = Telec; + + // Compute remaining initial guesses. + const double Psiqp = (1.0 / (1.0 + (Xqp_ - Xqpp_) / (Xqpp_ - Xl_))) * (-Psiqpp * (Xqp_ - Xl_) / (Xqpp_ - Xl_) + Iq * (Xqp_ - Xl_)); + + const double Edp = Psiqp - Iq * (Xqp_ - Xl_); + + const double Psidp = -(1.0 / (1.0 + (Xdp_ - Xdpp_) / (Xdpp_ - Xl_))) * (Id * (Xdp_ - Xl_) - Psidpp * (Xdp_ - Xl_) / (Xdpp_ - Xl_)); + + const double Eqp = Psidp + Id * (Xdp_ - Xl_); + + y_[0] = delta; + y_[1] = omega_s_; + y_[2] = Eqp; + y_[3] = Psidp; + y_[4] = Psiqp; + y_[5] = Edp; + y_[6] = Psiqpp; + y_[7] = Psidpp; + y_[8] = Psipp; + y_[9] = ksat; + y_[10] = Vd; + y_[11] = Vq; + y_[12] = Telec; + y_[13] = Id; + y_[14] = Iq; + + yp_[0] = 0.0; + yp_[1] = 0.0; + yp_[2] = 0.0; + yp_[3] = 0.0; + yp_[4] = 0.0; + yp_[5] = 0.0; + yp_[6] = 0.0; + yp_[7] = 0.0; + yp_[8] = 0.0; + yp_[9] = 0.0; + yp_[10] = 0.0; + yp_[11] = 0.0; + yp_[12] = 0.0; + yp_[13] = 0.0; + yp_[14] = 0.0; + + return 0; } -} + /** + * \brief Identify differential variables. + */ + template + int GENROU::tagDifferentiable() + { + tag_[0] = true; + tag_[1] = true; + tag_[2] = true; + tag_[3] = true; + tag_[4] = true; + tag_[5] = true; + for (IdxT i = 6; i < size_; ++i) + { + tag_[i] = false; + } + + return 0; + } + + /** + * @brief Computes residual vector for the generator model. + * + * Residual equations are given by Adam Birchfield's paper + * on power electronics. + * + */ + template + int GENROU::evaluateResidual() + { + // These constants are the same as in the initialization state. Maybe move them to .hpp file? + double SB = ((1 / 0.2) * (1 / 0.2)) * ((sqrt(S1_) - sqrt(S12_)) * (sqrt(S1_) - sqrt(S12_))); + double SA = 1 - sqrt(S1_ / SB); + const double Vr = V() * cos(theta()); + const double Vi = V() * sin(theta()); + + // The first 6 functions are differential equations (involving yp_) + f_[0] = -yp_[0] + omega0_ * y_[1]; + f_[1] = -yp_[1] + (1 / (2 * H_)) * ((Pm_ - D_ * y_[1]) / (1 + y_[1]) - y_[12]); + f_[2] = -yp_[2] * Td0p_ + Ef_ - (y_[2] + (Xd_ - Xdp_) * (y_[13] + ((Xdp_ - Xdpp_) / ((Xdp_ - Xl_) * (Xdp_ - Xl_))) * (y_[2] - y_[3] - (Xdp_ - Xl_) * y_[13])) + y_[7] * y_[9]); + f_[3] = -yp_[3] + (1 / Td0pp_) * (y_[2] - y_[3] - (Xdp_ - Xl_) * y_[13]); + f_[4] = -yp_[4] + (1 / Tq0pp_) * (y_[5] - y_[4] + (Xqp_ - Xl_) * y_[14]); + f_[5] = -yp_[5] + (1 / Tq0p_) * (-y_[5] + (((Xq_ - Xl_) / (Xd_ - Xl_)) * y_[6] * y_[9] + (Xq_ - Xqp_) * (y_[14] - ((Xqp_ - Xqpp_) / ((Xqp_ - Xl_) * (Xqp_ - Xl_)) * y_[5] + y_[14]) * (Xqp_ - Xl_) - y_[4]))); + + // Remaining equations are all algebraic + f_[6] = -y_[6] - y_[4] * ((Xqp_ - Xqpp_) / (Xqp_ - Xl_)) - y_[5] * ((Xqpp_ - Xl_) / (Xqp_ - Xl_)); + f_[7] = -y_[7] + y_[3] * ((Xdp_ - Xdpp_) / (Xdp_ - Xl_)) + y_[2] * ((Xdpp_ - Xl_) / (Xdp_ - Xl_)); + f_[8] = -y_[8] + sqrt((y_[7] * y_[7]) + (y_[6] * y_[6])); + f_[9] = -y_[9] + SB * ((y_[8] - SA) * (y_[8] - SA)); + f_[10] = -y_[10] - y_[6] * (1 + y_[1]); + f_[11] = -y_[11] + y_[7] * (1 + y_[1]); + f_[12] = -y_[12] + y_[14] * (y_[7] - y_[13] * Xdpp_) - y_[13] * (y_[6] - y_[14] * Xdpp_); + f_[13] = Vr * sin(y_[0]) - Vi * cos(y_[0]) + y_[13] * Rs_ - y_[14] * Xqpp_; + f_[14] = Vr * cos(y_[0]) - Vi * sin(y_[0]) + y_[13] * Xqpp_ - y_[14] * Rs_; + + return 0; + } + /** + * @brief Creates the Jacobian using J and Jp. + * + * Values were calculated by hand and verified with MATLAB. + * + */ + template + int GENROU::evaluateJacobian() + { + J_.zeroMatrix(); + double SB = ((1 / 0.2) * (1 / 0.2)) * ((sqrt(S1_) - sqrt(S12_)) * (sqrt(S1_) - sqrt(S12_))); + double SA = 1 - sqrt(S1_ / SB); + const double Vr = V() * cos(theta()); + const double Vi = V() * sin(theta()); + + // These variables are for ease of reading + double df2_domega = Pm_ / ((1 + (y_[1]) * (y_[1])) * 2 * H_) - D_ / ((1 + y_[1]) * 2 * H_) + D_ * y_[1] / ((1 + (y_[1]) * (y_[1])) * 2 * H_); + double df2_dy12 = -1 / (2 * H_); + + double df3_dId = -(Xd_ - Xdp_) * (1 - (Xdp_ - Xl_) * ((Xdp_ - Xdpp_) / ((Xdp_ - Xl_) * (Xdp_ - Xl_)))); + double df3_dksat = -y_[7]; + double df3_dpsippd = -y_[9]; + double df3_dpsipd = (Xd_ - Xdp_) * (Xdp_ - Xdpp_) / ((Xdp_ - Xl_) * (Xdp_ - Xl_)); + double df3_dEpq = -1 - (Xd_ - Xdp_) * ((Xdp_ - Xdpp_) / ((Xdp_ - Xl_) * (Xdp_ - Xl_))); + + double df6_dpsippq = y_[9] * ((Xq_ - Xl_) / (Xd_ - Xl_)) / Tq0p_; + double df6_dEpd = -(1 / Tq0p_) * (1 + (Xq_ - Xqp_) * ((Xqp_ - Xqpp_) / ((Xqp_ - Xl_) * (Xqp_ - Xl_)))); + double df6_dpsipq = (Xq_ - Xqp_) * ((Xqp_ - Xqpp_) / ((Xqp_ - Xl_) * (Xqp_ - Xl_))) / Tq0p_; + double df6_dksat = y_[6] * ((Xq_ - Xl_) / (Xd_ - Xl_)) / Tq0p_; + double df6_dIq = (Xq_ - Xqp_) / Tq0p_ - (Xqp_ - Xl_) * (Xq_ - Xqp_) * ((Xqp_ - Xqpp_) / ((Xqp_ - Xl_) * (Xqp_ - Xl_))) / Tq0p_; + + double df7_dpsipq = -((Xqp_ - Xqpp_) / (Xqp_ - Xl_)); + double df7_dEpd = -((Xqpp_ - Xl_) / (Xqp_ - Xl_)); + + double df8_dpsipd = (Xdp_ - Xdpp_) / (Xdp_ - Xl_); + double df8_dEpq = (Xdpp_ - Xl_) / (Xdp_ - Xl_); + + double dTelec_dId = -Xdpp_ * y_[14] - y_[6] + y_[14] * Xdpp_; + double dTelec_dIq = y_[7] - y_[13] * Xdpp_ + Xdpp_ * y_[13]; + + // Create dF/dy + std::vector rcord{0, + 1, 1, + 2, 2, 2, 2, 2, + 3, 3, 3, + 4, 4, 4, + 5, 5, 5, 5, 5, + 6, 6, 6, + 7, 7, 7, + 8, 8, 8, + 9, 9, + 10, 10, 10, + 11, 11, 11, + 12, 12, 12, 12, 12, + 13, 13, 13, + 14, 14, 14}; + + std::vector ccord{1, + 1, 12, + 2, 3, 7, 9, 13, + 2, 3, 13, + 4, 5, 14, + 4, 5, 6, 9, 14, + 4, 5, 6, + 2, 3, 7, + 6, 7, 8, + 8, 9, + 1, 6, 10, + 1, 7, 11, + 6, 7, 12, 13, 14, + 0, 13, 14, + 0, 13, 14}; + + std::vector vals{omega0_, + -df2_domega, df2_dy12, + df3_dEpq, df3_dpsipd, df3_dpsippd, df3_dksat, df3_dId, + 1 / Td0pp_, -1 / Td0pp_, (Xl_ - Xdp_) / Td0pp_, + -1 / Tq0pp_, 1 / Tq0pp_, (Xqp_ - Xl_) / Tq0pp_, + df6_dpsipq, df6_dEpd, df6_dpsippq, df6_dksat, df6_dIq, + df7_dpsipq, df7_dEpd, -1, + df8_dEpq, df8_dpsipd, -1, + y_[6] / sqrt((y_[7] * y_[7]) + (y_[6] * y_[6])), y_[7] / sqrt((y_[7] * y_[7]) + (y_[6] * y_[6])), -1, + 2.0 * SB * (y_[8] - SA), -1, + -y_[6], -(1 + y_[1]), -1, + y_[7], 1 + y_[1], -1, + -y_[13], y_[14], -1, dTelec_dId, dTelec_dIq, + Vr * cos(y_[0]) + Vi * sin(y_[0]), Rs_, -Xqpp_, + -Vr * sin(y_[0]) + Vi * cos(y_[0]), Xqpp_, Rs_}; + J_.setValues(rcord, ccord, vals); + + // Create dF/dy' + std::vector rcordder{0, 1, 2, 3, 4, 5}; + std::vector ccordder{0, 1, 2, 3, 4, 5}; + std::vector valsder{-1, -1, -Td0p_, -1, -1, -1}; + COO_Matrix Jacder = COO_Matrix(rcordder, ccordder, valsder, 15, 15); + + // Perform dF/dy + \alpha dF/dy' + J_.axpy(alpha_, Jacder); + return 0; + } -// Available template instantiations -template class GENROU; -template class GENROU; + /** + * @brief This is not used, but is implemented to keep SystemModel.hpp happy. + * + */ + template + int GENROU::evaluateIntegrand() + { + return 0; + } + + /** + * @brief This is not used, but is implemented to keep SystemModel.hpp happy. + * + */ + template + int GENROU::initializeAdjoint() + { + return 0; + } + + /** + * @brief This is not used, but is implemented to keep SystemModel.hpp happy. + * + */ + template + int GENROU::evaluateAdjointResidual() + { + return 0; + } + + /** + * @brief This is not used, but is implemented to keep SystemModel.hpp happy. + * + */ + template + int GENROU::evaluateAdjointIntegrand() + { + return 0; + } + // Available template instantiations + template class GENROU; + template class GENROU; } // namespace ModelLib \ No newline at end of file diff --git a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp index fad266b..16ac561 100644 --- a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp +++ b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp @@ -13,7 +13,7 @@ namespace ModelLib namespace ModelLib { /*! - * @brief Implementation of a fourth order generator model. + * @brief Implementation of a GENROU synchronous machine model. * */ template @@ -49,11 +49,10 @@ namespace ModelLib int tagDifferentiable(); int evaluateResidual(); int evaluateJacobian(); - int evaluateIntegrand(); + int evaluateIntegrand(); int initializeAdjoint(); int evaluateAdjointResidual(); - //int evaluateAdjointJacobian(); int evaluateAdjointIntegrand(); void updateTime(real_type t, real_type a) @@ -83,33 +82,6 @@ namespace ModelLib return bus_->theta(); } - ScalarT& P() - { - return bus_->P(); - } - - const ScalarT& P() const - { - return bus_->P(); - } - - ScalarT& Q() - { - return bus_->Q(); - } - - const ScalarT& Q() const - { - return bus_->Q(); - } - - - private: - ScalarT Pg(); - ScalarT Qg(); - ScalarT frequencyPenalty(ScalarT omega); - ScalarT frequencyPenaltyDer(ScalarT omega); - private: real_type omega0_;///< Nominal frequency (assume 2pi60 rad/s) real_type H_; ///< Inertia constant [s] @@ -133,10 +105,6 @@ namespace ModelLib real_type Ef_; real_type Pm_; real_type omega_s_; - real_type omega_up_; - real_type omega_lo_; - real_type c_; - real_type beta_; ScalarT P0_; ScalarT Q0_; diff --git a/Examples/GenrouTest/GenrouTest.cpp b/Examples/GenrouTest/GenrouTest.cpp index fa30795..12de7b9 100644 --- a/Examples/GenrouTest/GenrouTest.cpp +++ b/Examples/GenrouTest/GenrouTest.cpp @@ -10,13 +10,26 @@ #include #include #include +#include +#include #include #include -#include -#include #include +void printOutput(realtype time, const N_Vector yy) +{ + static std::ofstream outFile("output.txt"); + + realtype* yval = N_VGetArrayPointer_Serial(yy); + + outFile << time; + for (size_t j = 0; j < N_VGetLength_Serial(yy); j++) + { + outFile << " " << yval[j]; + } + outFile << "\n"; +} int main() { @@ -45,17 +58,17 @@ int main() model->initialize(); model->evaluateResidual(); - std::cout << "Verify Intial Residual is Zero: {"; + std::cout << "Verify Initial Residual is Zero: {"; for (double i : model->getResidual()) { std::cout << i << ", "; } - std::cout << "}\n"; + std::cout << "}" << std::endl; - //model->updateTime(0.0, 1.0); - //model->evaluateJacobian(); - std::cout << "Intial Jacobian with alpha = 1:\n"; + model->updateTime(0.0, 1.0); + model->evaluateJacobian(); + //std::cout << "Initial Jacobian with alpha = 1:\n"; //model->getJacobian().printMatrix(); @@ -64,12 +77,14 @@ int main() double t_init = 0.0; double t_final = 5.0; - double t_timestep = 0.0001; + double t_timestep = 0.001; + + idas->setOutputCallback(printOutput); idas->configureSimulation(); idas->getDefaultInitialCondition(); idas->initializeSimulation(t_init); - idas->runSimulation(t_final, 50000); + idas->runSimulation(t_final, 5000); delete idas; delete model; diff --git a/ModelEvaluatorImpl.hpp b/ModelEvaluatorImpl.hpp index 0f8969f..28b4c81 100644 --- a/ModelEvaluatorImpl.hpp +++ b/ModelEvaluatorImpl.hpp @@ -113,7 +113,7 @@ namespace ModelLib virtual bool hasJacobian() { - return false; + return !J_.isEmpty(); } virtual IdxT size_quad() diff --git a/Solver/Dynamic/Ida.cpp b/Solver/Dynamic/Ida.cpp index c2b7bab..d857dd3 100644 --- a/Solver/Dynamic/Ida.cpp +++ b/Solver/Dynamic/Ida.cpp @@ -60,7 +60,6 @@ #include #include -#include #include /* access to IDADls interface */ #include @@ -270,7 +269,6 @@ namespace Sundials template int Ida::runSimulation(real_type tf, int nout) { - std::ofstream outFile("output.txt"); // Output file that data will be printed to int retval = 0; int iout = 0; real_type tret; @@ -279,18 +277,17 @@ namespace Sundials /* In loop, call IDASolve, print results, and test for error. * Break out of loop when NOUT preset output times have been reached. */ + //printOutput(0.0); while(nout > iout) { retval = IDASolve(solver_, tout, &tret, yy_, yp_, IDA_NORMAL); checkOutput(retval, "IDASolve"); - realtype *yval = N_VGetArrayPointer_Serial(yy_); + //printOutput(tout); - outFile << tout; // Output time first - for (size_t j = 0; j < model_->size(); j++) { - outFile << " " << yval[j]; // Each yval in the array + if (outputCallback_) + { + outputCallback_(tret, yy_); } - outFile << "\n"; // Newline after last yval is printed - if (retval == IDA_SUCCESS) { @@ -298,7 +295,7 @@ namespace Sundials tout += dt; } } - outFile.close(); + //std::cout << "\n"; return retval; } @@ -799,6 +796,13 @@ namespace Sundials } } + template + void Ida::setOutputCallback(std::function callback) + { + outputCallback_ = callback; + } + + // Compiler will prevent building modules with data type incompatible with realtype template class Ida; template class Ida; diff --git a/Solver/Dynamic/Ida.hpp b/Solver/Dynamic/Ida.hpp index a4cea95..bf12b85 100644 --- a/Solver/Dynamic/Ida.hpp +++ b/Solver/Dynamic/Ida.hpp @@ -160,6 +160,7 @@ namespace AnalysisManager void printOutput(realtype t); void printSpecial(realtype t, N_Vector x); void printFinalStats(); + void setOutputCallback(std::function callback); private: static int Residual(realtype t, @@ -192,6 +193,7 @@ namespace AnalysisManager SUNMatrix JacobianMatB_; SUNLinearSolver linearSolver_; SUNLinearSolver linearSolverB_; + std::function outputCallback_; real_type t_init_; real_type t_final_; diff --git a/SparseMatrix/COO_Matrix.hpp b/SparseMatrix/COO_Matrix.hpp index 107c3be..c729642 100644 --- a/SparseMatrix/COO_Matrix.hpp +++ b/SparseMatrix/COO_Matrix.hpp @@ -69,7 +69,7 @@ class COO_Matrix std::tuple getDimensions(); void printMatrix(); - + bool isEmpty(); static void sortSparseCOO(std::vector &rows, std::vector &columns, std::vector &vals); @@ -726,6 +726,16 @@ inline void COO_Matrix::sortSparseCOO(std::vector &rows, } } +/** + * @brief Return if there are no values, meaning that the matrix is empty. + * + */ +template +inline bool COO_Matrix::isEmpty() +{ + return values_.empty(); +} + template inline COO_Matrix::COO_Matrix(std::vector r, std::vector c, std::vector v, Intdx m, Intdx n) { diff --git a/SystemModel.hpp b/SystemModel.hpp index 54ff5bb..de38ce2 100644 --- a/SystemModel.hpp +++ b/SystemModel.hpp @@ -116,9 +116,9 @@ class SystemModel : public ModelEvaluatorImpl SystemModel() : ModelEvaluatorImpl(0, 0, 0) { // Set system model tolerances - rtol_ = 1e-4; - atol_ = 1e-6; - this->max_steps_=50000; + rtol_ = 1e-7; + atol_ = 1e-9; + this->max_steps_=5000; } /** @@ -192,15 +192,13 @@ class SystemModel : public ModelEvaluatorImpl * @return false */ bool hasJacobian() - { // do hasJac_ *= ....... - for (const auto &component : components_) - { - if (!component->hasJacobian()) - { - return false; - } - } - return false; //Change this to false to make other tests run + { + bool hasJac; + for (const auto &component : components_) + { + hasJac *= component->hasJacobian(); + } + return hasJac; } /** @@ -393,9 +391,6 @@ class SystemModel : public ModelEvaluatorImpl /** * @brief Evaluate system Jacobian. * - * @todo Need to implement Jacobian. For now, using finite difference - * approximation provided by IDA. This works for dense Jacobian matrix - * only. * */ int evaluateJacobian() @@ -724,8 +719,6 @@ class SystemModel : public ModelEvaluatorImpl { component->updateTime(t, a); } - time_ = t; - alpha_ = a; } void addBus(bus_type* bus) @@ -741,7 +734,6 @@ class SystemModel : public ModelEvaluatorImpl private: std::vector buses_; std::vector components_; - bool usejac_; }; // class SystemModel From 4d92e333d5f930739dec81c7fb9fe599a6a385fc Mon Sep 17 00:00:00 2001 From: Paulm991 Date: Thu, 8 Aug 2024 10:16:59 -0400 Subject: [PATCH 5/5] Added all files for reference purposes --- ComponentLib/DynamicPhasor/CMakeLists.txt | 1 + .../DynamicPhasor/Governor/Governor.cpp | 360 ++++++++++++++++++ .../DynamicPhasor/Governor/Governor.hpp | 149 ++++++++ .../Generator4Governor/Generator4Governor.cpp | 3 +- Examples/GenConstLoad/GenConstLoad.cpp | 55 ++- Examples/GenrouTest/CMakeLists.txt | 6 +- Examples/GenrouTest/GenrouGovTest.cpp | 102 +++++ Examples/GenrouTest/GenrouTest.cpp | 2 +- Solver/Dynamic/Ida.cpp | 4 +- SystemModel.hpp | 2 +- 10 files changed, 677 insertions(+), 7 deletions(-) create mode 100644 Examples/GenrouTest/GenrouGovTest.cpp diff --git a/ComponentLib/DynamicPhasor/CMakeLists.txt b/ComponentLib/DynamicPhasor/CMakeLists.txt index ee3a45e..00dc18e 100644 --- a/ComponentLib/DynamicPhasor/CMakeLists.txt +++ b/ComponentLib/DynamicPhasor/CMakeLists.txt @@ -1,2 +1,3 @@ add_subdirectory(SynchronousMachine) +add_subdirectory(Governor) add_subdirectory(Bus) \ No newline at end of file diff --git a/ComponentLib/DynamicPhasor/Governor/Governor.cpp b/ComponentLib/DynamicPhasor/Governor/Governor.cpp index e69de29..42e13d8 100644 --- a/ComponentLib/DynamicPhasor/Governor/Governor.cpp +++ b/ComponentLib/DynamicPhasor/Governor/Governor.cpp @@ -0,0 +1,360 @@ +// Created by Paul Moon 7/11/2024 + +#include +#include +#include + +#include +#include "Governor.hpp" + +namespace ModelLib +{ + + /*! + * @brief Constructor for a simple governor model + * + * Arguments passed to ModelEvaluatorImpl: + * - Number of equations = 6 differential + 9 algebraic = 15 + * - Quadratures and optimization parameters are not being used. + */ + template + Governor::Governor(bus_type *bus, ScalarT P0, ScalarT Q0) + : ModelEvaluatorImpl(18, 0, 0), + H_(3.0), // Intertia constant + D_(0.0), // Damping coefficient + Xq_(0.5), // Machine reactance parameter + Xd_(2.1), // Machine reactance parameter + Xqp_(0.5), // Machine reactance parameter + Xdp_(0.2), // Machine reactance parameter + Xqpp_(0.18), // Machine reactance parameter + Xdpp_(0.18), // Machine reactance parameter + Rs_(0.0), // Winding resistance + Tq0p_(0.75), // Time constant + Td0p_(7.0), // Time constant + Tq0pp_(0.05), // Time constant + Td0pp_(0.035), // Time constant + Ef_(1.8336), // Field winding voltage from the excitation system. An exciter would set this value + X_(0.22), // Machine reactance parameter + Xl_(0.15), // Machine reactance parameter + omega_s_(0.0), // Relative speed deviation. 0 is synchronous, -1 is 0 Hz. + omega0_(2.0 * 60.0 * M_PI), // Nominal frequency + S1_(0.05), // Saturation parameter + S12_(0.2), // Saturation parameter + R_(0.05), // Droop constant + T1_(0.5), // Valve time delay + T2_(2.5), // Turbine numerator time constant + T3_(7.5), // Turbine delay + Vmax_(1), // Max valve position + Vmin_(-1), // Min valve position + Dt_(0), // Turbine damping + P0_(P0), // Real power + Q0_(Q0), // Reactive power + offsetGov_(15), + bus_(bus) + { + } + + template + Governor::~Governor() + { + } + + /*! + * @brief This function will be used to allocate sparse Jacobian matrices. + * + */ + template + int Governor::allocate() + { + tag_.resize(size_); + + return 0; + } + + /** + * @brief Initialization of the governor model + * + * Initialization equations are derived from Adam + * Birchfield's paper on power electronics. + * + */ + template + int Governor::initialize() + { + // Compute values for saturation constants and real/imaginary voltage values + double SB = ((1 / 0.2) * (1 / 0.2)) * ((sqrt(S1_) - sqrt(S12_)) * (sqrt(S1_) - sqrt(S12_))); + double SA = 1 - sqrt(S1_ / SB); + const double Vr = V() * cos(theta()); + const double Vi = V() * sin(theta()); + + const double Ir = P0_ / V(); + const double Ii = 0.0; + + // Compute initial guess for the generator voltage phase + const double delta = atan((Xq_ * Ir + Vi) / (-Xq_ * Ii + Vr)); + + // Compute initial guesses for generator currents and potentials in d-q frame + const double Id = Ir * sin(delta) - Ii * cos(delta); + const double Iq = Ir * cos(delta) + Ii * sin(delta); + + const double Vq = Vr * cos(delta) + Vi * sin(delta) + Id * Xqpp_ + Iq * Rs_; + const double Vd = Vr * sin(delta) - Vi * cos(delta) + Id * Rs_ - Iq * Xqpp_; + + // Compute initial guesses for flux + const double Psiqpp = -Vd / (1 + omega_s_); + const double Psidpp = Vq / (1 + omega_s_); + + const double Psipp = sqrt((Psiqpp * Psiqpp) + (Psidpp * Psidpp)); + const double ksat = SB * ((Psipp - SA) * (Psipp - SA)); + + // Compute initial guess for Telec. Note that Pmech would normally be changed by a governor. + const double Telec = (Psidpp - Id * Xdpp_) * Iq - (Psiqpp - Iq * Xdpp_) * Id; + const double Pmech = Telec; + + // Compute remaining initial guesses. + const double Psiqp = (1.0 / (1.0 + (Xqp_ - Xqpp_) / (Xqpp_ - Xl_))) * (-Psiqpp * (Xqp_ - Xl_) / (Xqpp_ - Xl_) + Iq * (Xqp_ - Xl_)); + + const double Edp = Psiqp - Iq * (Xqp_ - Xl_); + + const double Psidp = -(1.0 / (1.0 + (Xdp_ - Xdpp_) / (Xdpp_ - Xl_))) * (Id * (Xdp_ - Xl_) - Psidpp * (Xdp_ - Xl_) / (Xdpp_ - Xl_)); + + const double Eqp = Psidp + Id * (Xdp_ - Xl_); + + const double Pv_init = Vmin_; + const double Ptx_init = Pmech * T3_ - T2_ * Pv_init-T3_*Dt_*omega_s_; + const double Pref_init = P0_; + + y_[0] = delta; + y_[1] = omega_s_; + y_[2] = Eqp; + y_[3] = Psidp; + y_[4] = Psiqp; + y_[5] = Edp; + y_[6] = Psiqpp; + y_[7] = Psidpp; + y_[8] = Psipp; + y_[9] = ksat; + y_[10] = Vd; + y_[11] = Vq; + y_[12] = Telec; + y_[13] = Id; + y_[14] = Iq; + y_[15] = Pmech; + + yp_[0] = 0.0; + yp_[1] = 0.0; + yp_[2] = 0.0; + yp_[3] = 0.0; + yp_[4] = 0.0; + yp_[5] = 0.0; + yp_[6] = 0.0; + yp_[7] = 0.0; + yp_[8] = 0.0; + yp_[9] = 0.0; + yp_[10] = 0.0; + yp_[11] = 0.0; + yp_[12] = 0.0; + yp_[13] = 0.0; + yp_[14] = 0.0; + yp_[15] = 0.0; + + y_[offsetGov_ + 1] = Ptx_init; + y_[offsetGov_ + 2] = Pv_init; + + yp_[offsetGov_ + 1] = 0.0; + yp_[offsetGov_ + 2] = 0.0; + + for (int i = 0; i < size_; ++i) { + std::cout << y_[i] << ", "; + } + std::cout << std::endl; + + return 0; + } + + /** + * \brief Identify differential variables. + */ + template + int Governor::tagDifferentiable() + { + tag_[0] = true; + tag_[1] = true; + tag_[2] = true; + tag_[3] = true; + tag_[4] = true; + tag_[5] = true; + for (IdxT i = 6; i < size_; ++i) + { + tag_[i] = false; + } + tag_[offsetGov_ + 1] = true; + tag_[offsetGov_ + 2] = true; + + return 0; + } + + /** + * @brief Computes residual vector for the generator model. + * + * Residual equations are given by Adam Birchfield's paper + * on power electronics. + * + */ + template + int Governor::evaluateResidual() + { + // These constants are the same as in the initialization state. Maybe move them to .hpp file? + double SB = ((1 / 0.2) * (1 / 0.2)) * ((sqrt(S1_) - sqrt(S12_)) * (sqrt(S1_) - sqrt(S12_))); + double SA = 1 - sqrt(S1_ / SB); + const double Vr = V() * cos(theta()); + const double Vi = V() * sin(theta()); + double Pref = 1.0; + + // The first 6 functions are differential equations (involving yp_) + f_[0] = -yp_[0] + omega0_ * y_[1]; + f_[1] = -yp_[1] + (1 / (2 * H_)) * ((y_[15] - D_ * y_[1]) / (1 + y_[1]) - y_[12]); + f_[2] = -yp_[2] * Td0p_ + Ef_ - (y_[2] + (Xd_ - Xdp_) * (y_[13] + ((Xdp_ - Xdpp_) / ((Xdp_ - Xl_) * (Xdp_ - Xl_))) * (y_[2] - y_[3] - (Xdp_ - Xl_) * y_[13])) + y_[7] * y_[9]); + f_[3] = -yp_[3] + (1 / Td0pp_) * (y_[2] - y_[3] - (Xdp_ - Xl_) * y_[13]); + f_[4] = -yp_[4] + (1 / Tq0pp_) * (y_[5] - y_[4] + (Xqp_ - Xl_) * y_[14]); + f_[5] = -yp_[5] + (1 / Tq0p_) * (-y_[5] + (((Xq_ - Xl_) / (Xd_ - Xl_)) * y_[6] * y_[9] + (Xq_ - Xqp_) * (y_[14] - ((Xqp_ - Xqpp_) / ((Xqp_ - Xl_) * (Xqp_ - Xl_)) * y_[5] + y_[14]) * (Xqp_ - Xl_) - y_[4]))); + + // Remaining equations are all algebraic + f_[6] = -y_[6] - y_[4] * ((Xqp_ - Xqpp_) / (Xqp_ - Xl_)) - y_[5] * ((Xqpp_ - Xl_) / (Xqp_ - Xl_)); + f_[7] = -y_[7] + y_[3] * ((Xdp_ - Xdpp_) / (Xdp_ - Xl_)) + y_[2] * ((Xdpp_ - Xl_) / (Xdp_ - Xl_)); + f_[8] = -y_[8] + sqrt((y_[7] * y_[7]) + (y_[6] * y_[6])); + f_[9] = -y_[9] + SB * ((y_[8] - SA) * (y_[8] - SA)); + f_[10] = -y_[10] - y_[6] * (1 + y_[1]); + f_[11] = -y_[11] + y_[7] * (1 + y_[1]); + f_[12] = -y_[12] + y_[14] * (y_[7] - y_[13] * Xdpp_) - y_[13] * (y_[6] - y_[14] * Xdpp_); + f_[13] = Vr * sin(y_[0]) - Vi * cos(y_[0]) + y_[13] * Rs_ - y_[14] * Xqpp_; + f_[14] = Vr * cos(y_[0]) - Vi * sin(y_[0]) + y_[13] * Xqpp_ - y_[14] * Rs_; + f_[15] = -y_[offsetGov_] + (1.0 / T3_) * (y_[offsetGov_ + 1] + T2_ * y_[offsetGov_ + 2]) - Dt_ * y_[1]; + + P() += Pg(); + Q() += Qg(); + + if (y_[offsetGov_ + 2] >= Vmax_) + { + y_[offsetGov_ + 2] = Vmax_; + if (((1.0/T1_) * (-y_[offsetGov_ + 2]) + (1.0/R_) * (Pref - y_[1])) > 0) + { + f_[offsetGov_+1] = -yp_[offsetGov_ + 2]; + } + else + { + f_[offsetGov_+1] = -yp_[offsetGov_ + 2] + (1.0/T1_) * (-y_[offsetGov_ + 2] + (1.0/R_) * (Pref - y_[1])); + } + } + else if (y_[offsetGov_ + 2] <= Vmin_) + { + y_[offsetGov_ + 2] = Vmin_; + if (((1.0/T1_) * (-y_[offsetGov_ + 2]) + (1.0/R_) * (Pref - y_[1])) > 0) + { + f_[offsetGov_+1] = -yp_[offsetGov_ + 2]; + } + else + { + f_[offsetGov_+1] = -yp_[offsetGov_ + 2] + (1.0/T1_) * (-y_[offsetGov_ + 2] + (1.0/R_) * (Pref - y_[1])); + } + } + else + { + f_[offsetGov_+1] = -yp_[offsetGov_ + 2] + (1.0/T1_) * (-y_[offsetGov_ + 2] + (1.0/R_) * (Pref - y_[1])); + } + f_[offsetGov_+2] = -yp_[offsetGov_ + 1] + y_[offsetGov_ + 2] - (1.0 / T3_) * (y_[offsetGov_ + 1] + T2_ * y_[offsetGov_ + 2]); + //f_[offsetGov_+3] = -y_[offsetGov_] + (1.0 / T3_) * (y_[offsetGov_ + 1] + T2_ * y_[offsetGov_ + 2]) - Dt_ * y_[1]; + return 0; + } + + /** + * @brief Creates the Jacobian using J and Jp. + * + * Values were calculated by hand and verified with MATLAB. + * + */ + template + int Governor::evaluateJacobian() + { + return 0; + } + + template + ScalarT Governor::Pg() + { + return y_[14]*y_[2] + y_[13]*y_[5] + (Xqp_ - Xdp_)*y_[13]*y_[14] - Rs_*(y_[13]*y_[13] + y_[14]*y_[14]); + } + + /*y_[0] = delta; + y_[1] = omega_s_; + y_[2] = Eqp; + y_[3] = Psidp; + y_[4] = Psiqp; + y_[5] = Edp; + y_[6] = Psiqpp; + y_[7] = Psidpp; + y_[8] = Psipp; + y_[9] = ksat; + y_[10] = Vd; + y_[11] = Vq; + y_[12] = Telec; + y_[13] = Id; + y_[14] = Iq; + y_[15] = Pmech;*/ + +/** + * Generator reactive power Qg. + * + * \f[ Q_g = E_q' I_d - E_d' I_q - X_d' I_d^2 - X_q' I_q^2 \f] + */ + template + ScalarT Governor::Qg() + { + return -y_[14]*y_[5] + y_[13]*y_[2] - Xdp_*y_[13]*y_[13] - Xqp_*y_[14]*y_[14]; + } + + /** + * @brief This is not used, but is implemented to keep SystemModel.hpp happy. + * + */ + template + int Governor::evaluateIntegrand() + { + return 0; + } + + /** + * @brief This is not used, but is implemented to keep SystemModel.hpp happy. + * + */ + template + int Governor::initializeAdjoint() + { + return 0; + } + + /** + * @brief This is not used, but is implemented to keep SystemModel.hpp happy. + * + */ + template + int Governor::evaluateAdjointResidual() + { + return 0; + } + + /** + * @brief This is not used, but is implemented to keep SystemModel.hpp happy. + * + */ + template + int Governor::evaluateAdjointIntegrand() + { + return 0; + } + + // Available template instantiations + template class Governor; + template class Governor; + +} // namespace ModelLib \ No newline at end of file diff --git a/ComponentLib/DynamicPhasor/Governor/Governor.hpp b/ComponentLib/DynamicPhasor/Governor/Governor.hpp index e69de29..98c0919 100644 --- a/ComponentLib/DynamicPhasor/Governor/Governor.hpp +++ b/ComponentLib/DynamicPhasor/Governor/Governor.hpp @@ -0,0 +1,149 @@ +// Created by Paul Moon 7/11/2024 + +#ifndef _GOVERNOR_H_ +#define _GOVERNOR_H_ + +#include + +namespace ModelLib +{ + template class BaseBus; +} + +namespace ModelLib +{ + /*! + * @brief Implementation of a TGOV1 Turbine-Governor model. + * + */ + template + class Governor : public ModelEvaluatorImpl + { + using ModelEvaluatorImpl::size_; + using ModelEvaluatorImpl::nnz_; + using ModelEvaluatorImpl::time_; + using ModelEvaluatorImpl::alpha_; + using ModelEvaluatorImpl::y_; + using ModelEvaluatorImpl::yp_; + using ModelEvaluatorImpl::tag_; + using ModelEvaluatorImpl::f_; + using ModelEvaluatorImpl::g_; + using ModelEvaluatorImpl::yB_; + using ModelEvaluatorImpl::ypB_; + using ModelEvaluatorImpl::fB_; + using ModelEvaluatorImpl::gB_; + using ModelEvaluatorImpl::param_; + using ModelEvaluatorImpl::param_up_; + using ModelEvaluatorImpl::param_lo_; + using ModelEvaluatorImpl::J_; + + typedef typename ModelEvaluatorImpl::real_type real_type; + typedef BaseBus bus_type; + + public: + Governor(BaseBus* bus, ScalarT P0, ScalarT Q0); + virtual ~Governor(); + + int allocate(); + int initialize(); + int tagDifferentiable(); + int evaluateResidual(); + int evaluateJacobian(); + + int evaluateIntegrand(); + int initializeAdjoint(); + int evaluateAdjointResidual(); + int evaluateAdjointIntegrand(); + + ScalarT Pg(); + ScalarT Qg(); + + + void updateTime(real_type t, real_type a) + { + time_ = t; + alpha_ = a; + } + + // Inline accesor functions + ScalarT& V() + { + return bus_->V(); + } + + const ScalarT& V() const + { + return bus_->V(); + } + + ScalarT& theta() + { + return bus_->theta(); + } + + const ScalarT& theta() const + { + return bus_->theta(); + } + ScalarT& P() + { + return bus_->P(); + } + + const ScalarT& P() const + { + return bus_->P(); + } + + ScalarT& Q() + { + return bus_->Q(); + } + + const ScalarT& Q() const + { + return bus_->Q(); + } + + + private: + real_type omega0_;///< Nominal frequency (assume 2pi60 rad/s) + real_type H_; ///< Inertia constant [s] + real_type D_; ///< Damping constant [pu] + real_type Rs_; ///< stator resistance [pu] + real_type Xl_; ///< stator leakage reactance [pu] + real_type Xd_; ///< direct axis synchronous reactance [pu] + real_type Xdp_; ///< direct axis transient reactance [pu] + real_type Xdpp_; ///< direct axis sub-transient reactance [pu] + real_type Xq_; ///< quadrature axis synchronous reactance [pu] + real_type Xqp_; ///< quadrature axis transient reactance [pu] + real_type Xqpp_; ///< quadrature axis sub-transient reactance [pu] + real_type Td0p_; ///< open circuit direct axis transient time constant [s] + real_type Td0pp_; ///< open circuit direct axis sub-transient time constant [s] + real_type Tq0p_; ///< open circuit quadrature axis transient time constant [s] + real_type Tq0pp_; ///< open circuit quadrature axis sub-transient time constant [s] + real_type S1_; ///< saturation at 1.0 pu flux + real_type S12_; ///< saturation at 1.2 pu flux + real_type X_; + real_type Ef_; + real_type omega_s_; + + real_type R_; + real_type T1_; + real_type T2_; + real_type T3_; + real_type Vmax_; + real_type Vmin_; + real_type Dt_; + real_type offsetGov_; + + ScalarT P0_; + ScalarT Q0_; + + bus_type* bus_; + }; + +} // namespace ModelLib + + +#endif // _GOVERNOR_H \ No newline at end of file diff --git a/ComponentLib/PowerFlow/Generator4Governor/Generator4Governor.cpp b/ComponentLib/PowerFlow/Generator4Governor/Generator4Governor.cpp index db26213..11dba5c 100644 --- a/ComponentLib/PowerFlow/Generator4Governor/Generator4Governor.cpp +++ b/ComponentLib/PowerFlow/Generator4Governor/Generator4Governor.cpp @@ -157,6 +157,7 @@ int Generator4Governor::initialize() // Compute generator voltage phase const ScalarT delta = atan((Xq_*P0_ - Rs_*Q0_) / (V()*V() + Rs_*P0_ + Xq_*Q0_)) + theta(); + std::cout << delta << std::endl; // Compute generator current phase const ScalarT phi = delta - theta() + atan(Q0_/P0_); @@ -169,7 +170,7 @@ int Generator4Governor::initialize() // Initialize generator y_[offsetGen_ + 0] = delta; - y_[offsetGen_ + 1] = omega_s_ + 0.2; // <~ this is hack to perturb omega + y_[offsetGen_ + 1] = omega_s_; // <~ this is hack to perturb omega y_[offsetGen_ + 2] = Edp; y_[offsetGen_ + 3] = Eqp; y_[offsetGen_ + 4] = Id; diff --git a/Examples/GenConstLoad/GenConstLoad.cpp b/Examples/GenConstLoad/GenConstLoad.cpp index fbbc271..f5486a2 100644 --- a/Examples/GenConstLoad/GenConstLoad.cpp +++ b/Examples/GenConstLoad/GenConstLoad.cpp @@ -60,6 +60,8 @@ #include #include +#include +#include #include #include @@ -73,6 +75,20 @@ #include #include +void printOutput(realtype time, const N_Vector yy) +{ + static std::ofstream outFile("output.txt"); + + realtype* yval = N_VGetArrayPointer_Serial(yy); + + outFile << time; + for (size_t j = 0; j < N_VGetLength_Serial(yy); j++) + { + outFile << " " << yval[j]; + } + outFile << "\n"; +} + int main() { using namespace ModelLib; @@ -98,8 +114,45 @@ int main() // Create numerical integrator and configure it for the generator model Ida* idas = new Ida(model); + // allocate model components model->allocate(); + std::cout << "Size: " << model->y().size() << std::endl; + + model->initialize(); + model->evaluateResidual(); + + std::cout << "Verify Initial Residual is Zero: {"; + for (double i : model->getResidual()) + { + std::cout << i << ", "; + } + std::cout << "}" << std::endl; + + + //model->updateTime(0.0, 1.0); + //model->evaluateJacobian(); + //std::cout << "Initial Jacobian with alpha = 1:\n"; + //model->getJacobian().printMatrix(); + + + // Create numerical integrator and configure it for the generator model + //AnalysisManager::Sundials::Ida* idas = new AnalysisManager::Sundials::Ida(model); + + double t_init = 0.0; + double t_final = 5.0; + double t_timestep = 0.001; + + idas->setOutputCallback(printOutput); + + idas->configureSimulation(); + idas->getDefaultInitialCondition(); + idas->initializeSimulation(t_init, true); + idas->runSimulation(t_final, 5000); + + + /*model->allocate(); + double t_init = 0.0; double t_final = 15.0; @@ -171,7 +224,7 @@ int main() if(retval != 0) { std::cout << "The two results differ beyond solver tolerance!\n"; - } + }*/ delete idas; diff --git a/Examples/GenrouTest/CMakeLists.txt b/Examples/GenrouTest/CMakeLists.txt index b6435e8..c10f723 100644 --- a/Examples/GenrouTest/CMakeLists.txt +++ b/Examples/GenrouTest/CMakeLists.txt @@ -1,5 +1,9 @@ add_executable(genroutest GenrouTest.cpp) +add_executable(genrougovtest GenrouGovTest.cpp) target_link_libraries(genroutest GRIDKIT::bus GRIDKIT::genrou GRIDKIT::solvers_opt GRIDKIT::solvers_dyn) +target_link_libraries(genrougovtest GRIDKIT::bus GRIDKIT::genrou GRIDKIT::solvers_opt GRIDKIT::solvers_dyn GRIDKIT::governor GRIDKIT::load) install(TARGETS genroutest DESTINATION bin) +install(TARGETS genrougovtest DESTINATION bin) -add_test(NAME GenrouTest COMMAND $) \ No newline at end of file +add_test(NAME GenrouTest COMMAND $) +add_test(NAME GenrouGovTest COMMAND $) \ No newline at end of file diff --git a/Examples/GenrouTest/GenrouGovTest.cpp b/Examples/GenrouTest/GenrouGovTest.cpp new file mode 100644 index 0000000..df723e4 --- /dev/null +++ b/Examples/GenrouTest/GenrouGovTest.cpp @@ -0,0 +1,102 @@ +// Made by Paul Moon 7/11/2024 + + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +void printOutput(realtype time, const N_Vector yy) +{ + static std::ofstream outFile("output.txt"); + + realtype* yval = N_VGetArrayPointer_Serial(yy); + + outFile << time; + for (size_t j = 0; j < N_VGetLength_Serial(yy); j++) + { + outFile << " " << yval[j]; + } + outFile << "\n"; +} + +int main() +{ + using namespace ModelLib; + using namespace AnalysisManager::Sundials; + using namespace AnalysisManager; + using namespace GridKit::Testing; + + + // Create a bus + BaseBus* bus = new BusPQ(1.0, 0.0); + + // In GenrouGovTest.cpp + ModelEvaluatorImpl* gen = new Governor(bus, 1.0, 0.0); + + // Attach load to the bus + ModelEvaluatorImpl* load = new Load(bus, 1.0, 0.0); + + // Create a system model + SystemModel* model = new SystemModel(); + model->addBus(bus); + model->addComponent(gen); + model->addComponent(load); + + Ida* idas = new Ida(model); + + // allocate model components + model->allocate(); + + std::cout << "Size: " << model->y().size() << std::endl; + + model->initialize(); + model->evaluateResidual(); + + std::cout << "Verify Initial Residual is Zero: {"; + for (double i : model->getResidual()) + { + std::cout << i << ", "; + } + std::cout << "}" << std::endl; + + + //model->updateTime(0.0, 1.0); + //model->evaluateJacobian(); + //std::cout << "Initial Jacobian with alpha = 1:\n"; + //model->getJacobian().printMatrix(); + + + // Create numerical integrator and configure it for the generator model + //AnalysisManager::Sundials::Ida* idas = new AnalysisManager::Sundials::Ida(model); + + double t_init = 0.0; + double t_final = 1.0; + double t_timestep = 0.001; + + idas->setOutputCallback(printOutput); + + idas->configureSimulation(); + idas->getDefaultInitialCondition(); + idas->initializeSimulation(t_init, true); + idas->runSimulation(t_final, 5000); + + delete idas; + delete gen; + delete load; + delete bus; + delete model; + return 0; +} diff --git a/Examples/GenrouTest/GenrouTest.cpp b/Examples/GenrouTest/GenrouTest.cpp index 12de7b9..8d410a0 100644 --- a/Examples/GenrouTest/GenrouTest.cpp +++ b/Examples/GenrouTest/GenrouTest.cpp @@ -83,7 +83,7 @@ int main() idas->configureSimulation(); idas->getDefaultInitialCondition(); - idas->initializeSimulation(t_init); + idas->initializeSimulation(t_init, true); idas->runSimulation(t_final, 5000); delete idas; diff --git a/Solver/Dynamic/Ida.cpp b/Solver/Dynamic/Ida.cpp index d857dd3..5f0c091 100644 --- a/Solver/Dynamic/Ida.cpp +++ b/Solver/Dynamic/Ida.cpp @@ -256,8 +256,8 @@ namespace Sundials { int initType = IDA_Y_INIT; - if (tag_) - initType = IDA_YA_YDP_INIT; + /*if (tag_) + initType = IDA_YA_YDP_INIT;*/ retval = IDACalcIC(solver_, initType, 0.1); checkOutput(retval, "IDACalcIC"); diff --git a/SystemModel.hpp b/SystemModel.hpp index de38ce2..f2c4587 100644 --- a/SystemModel.hpp +++ b/SystemModel.hpp @@ -198,7 +198,7 @@ class SystemModel : public ModelEvaluatorImpl { hasJac *= component->hasJacobian(); } - return hasJac; + return false; //return hasJac; } /**