diff --git a/CMakeLists.txt b/CMakeLists.txt index 770855b9..a8de08e3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -20,6 +20,11 @@ set(PACKAGE_VERSION_PATCH "7") set(PACKAGE_VERSION "${PACKAGE_VERSION_MAJOR}.${PACKAGE_VERSION_MINOR}.${PACKAGE_VERSION_PATCH}") +# TODO: Probably beter to set a debug interface target +# set(CMAKE_CXX_FLAGS_DEBUG "-Wall -O0 -g -DDEBUG") + +set(CMAKE_CXX_STANDARD 17) + # Ipopt support is disabled by default option(GRIDKIT_ENABLE_IPOPT "Enable Ipopt support" OFF) @@ -33,6 +38,7 @@ option(GRIDKIT_ENABLE_SUNDIALS_SPARSE "Enable SUNDIALS sparse linear solvers" ON option(GRIDKIT_ENABLE_ENZYME "Enable automatic differentiation with Enzyme" OFF) set(CMAKE_MACOSX_RPATH 1) +set(CMAKE_WINDOWS_EXPORT_ALL_SYMBOLS 1) list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake) @@ -57,11 +63,6 @@ if("${isSystemDir}" STREQUAL "-1") endif("${isSystemDir}" STREQUAL "-1") -# TODO: Probably beter to set a debug interface target -set(CMAKE_CXX_FLAGS_DEBUG "-Wall -O0 -g -DDEBUG") - -set(CMAKE_CXX_STANDARD 17) - include_directories(${CMAKE_CURRENT_SOURCE_DIR}/src) option(GRIDKIT_BUILD_SHARED "Build shared libraries" ON) diff --git a/examples/Bus3Phasor/Bus3Example/CMakeLists.txt b/examples/Bus3Phasor/Bus3Example/CMakeLists.txt new file mode 100644 index 00000000..afd9f579 --- /dev/null +++ b/examples/Bus3Phasor/Bus3Example/CMakeLists.txt @@ -0,0 +1,9 @@ +add_executable(bus3grid bus3phasor.cpp) +target_link_libraries(bus3grid + GRIDKIT::bus + SUNDIALS::sunlinsolklu + SUNDIALS::core + SUNDIALS::ida + SUNDIALS::idas + SUNDIALS::sunmatrixdense) +install(TARGETS bus3grid RUNTIME DESTINATION bin) \ No newline at end of file diff --git a/examples/Bus3Phasor/CMakeLists.txt b/examples/Bus3Phasor/CMakeLists.txt new file mode 100644 index 00000000..a2b06cde --- /dev/null +++ b/examples/Bus3Phasor/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(Bus3Example) \ No newline at end of file diff --git a/examples/Bus3Phasor/README.md b/examples/Bus3Phasor/README.md new file mode 100644 index 00000000..e9d6e5aa --- /dev/null +++ b/examples/Bus3Phasor/README.md @@ -0,0 +1,189 @@ + + + +Notation: $I_{a,b}$ represents the current flowing into Bus b from the branch connecting Bus a and Bus b, $I_{G,a}$ the current from generator a, $I_{L,a}$ the current from load a, and $V_a$ the voltage at Bus a. +Each component is sub-index by the components ID. The 3Bus model is composed of 3 buses, 3 branches, 2 loads and 2 generators. For simplicity in some equation, $I = \left[\begin{array}{c} I_r \\ I_i \end{array}\right]$ to represent complex variables. + +## Bus equations: +```math +I_{G,1}+I_{2,1}+I_{3,1}+I_{L,1} = 0 +``` +```math +I_{L,2}+I_{1,2}+I_{3,2} = 0 +``` +```math +I_{G,3}+I_{1,3}+I_{2,3} = 0 +``` + +## Branch equations: + +(branch connecting Bus 1 and Bus 2) +```math +\left[ \begin{array}{c} + I_{2,1} \\ + I_{1,2} + \end{array}\right] = + Y^{[1]}\left[ \begin{array}{c} + V_1 \\ + V_2 +\end{array}\right] +``` + +(branch connecting Bus 2 and Bus 3) +```math +\left[ \begin{array}{c} + I_{3,2} \\ + I_{2,3} + \end{array}\right] = + Y^{[2]}\left[ \begin{array}{c} + V_2 \\ + V_3 + \end{array}\right] +``` + +(branch connecting Bus 1 and Bus 3) +```math +\left[ \begin{array}{c} + I_{3,1} \\ + I_{1,3} + \end{array}\right] = + Y^{[3]}\left[ \begin{array}{c} + V_1 \\ + V_3 + \end{array}\right] +``` + + ## Load equations: + Loads connected to bus 1: +```math + I_{L,1} = \frac{V_1}{R_1} +``` +Loads connected to bus 2: +```math + I_{L,2} = \frac{V_2}{R_2} +``` + +## Generator 1 +Generator connected to bus 1:
+Differential equations: +```math +\dot{\delta}_1 = \omega_1 \cdot \omega_0 +``` +```math +\dot{\omega}_1 = \frac{1}{2H_1}\left( \frac{P_{mech,1} - D_1\omega_1}{1+\omega_1} - T_{elec,1}\right) +``` +```math +\dot{E}_{qp,1} = \frac{1}{T_{dop,1}}\left( E_{fd,1} - \left( E_{qp,1} + X_{d1,1}\left( I_{d,1} + X_{d3,1}\left( E_{qp,1} - \psi_{dp,1} - X_{d2,1}I_{d,1} \right) \right) + \psi_{dpp,1}\cdot k_{sat,1} \right) \right) +``` +```math +\dot{\psi}_{dp,1} = \frac{1}{T_{dopp,1}} \left( E_{qp,1} - \psi_{dp,1} - X_{d2,1} I_{d, 1} \right) +``` +```math +\dot{\psi}_{qp,1} = \frac{1}{T_{qopp,1}} \left( E_{dp,1} - \psi_{qp,1} + X_{q2,1}I_{q,1}\right) +``` +```math +\dot{E}_{dp,1} = \frac{1}{T_{qop,1}}\left(-E_{dp,1} + X_{qd,1}\psi_{qpp,1}\cdot k_{sat,1} + X_{q1,1} \left(I_{q,1} - X_{q3,1} \left(E_{dp,1} + I_{q,1}X_{q2,1} - \psi_{qp,1}\right)\right)\right) +``` +Algebraic equations: +```math +V_{d,1} = -\psi_{qpp,1}(\psi_{qp,1},E_{dp,1})(1+\omega_1) +``` +```math +V_{q,1} = \psi_{dpp,1}(\psi_{dp,1},E_{qp,1})(1+\omega_1) +``` +```math +I_{d,1} = I_{G,r,1}\sin(\delta_1) - I_{G,i,1}\cos(\delta_1) +``` +```math +I_{q,1} = I_{G,r,1}\cos(\delta_1) + I_{G,i,1}\sin(\delta_1) +``` + +Network interfaces: +```math +V_{d,1} = V_{r,1}\sin(\delta_1) - V_{i,1}\cos(\delta_1) + I_{d,1} R_{a,1} -I_{q,1} X_{qpp,1} +``` +```math +V_{q,1} = V_{r, 1}\cos(\delta_1) + V_{i,1}\sin(\delta_1) + I_{d,1}X_{qpp,1} + I_{q,1}R_{a,1} +``` + +We can lump together the algebraic and network interface equations to get the following equivalent constraint equations: +```math +\left[\begin{array}{c} +I_{d} \\ +I_{q} +\end{array}\right] = +\left[\begin{array}{cr} +\sin(\delta_1) & -\cos(\delta_1)\\ +\cos(\delta_1) & \sin(\delta_1) +\end{array}\right] +\left[\begin{array}{c} +I_{G,r,1} \\ +I_{G,i,1} +\end{array}\right] +``` +```math +\left[\begin{array}{c} +-\psi_{qpp,1}(\psi_{qp,1},E_{dp,1})(1+\omega_1) \\ +\psi_{dpp,1}(\psi_{dp,1},E_{qp,1})(1+\omega_1) +\end{array}\right] +- \left[\begin{array}{cr} +\sin(\delta_1) & -\cos(\delta_1)\\ +\cos(\delta_1) & \sin(\delta_1) +\end{array}\right] +\left[\begin{array}{c} +V_{r,1} \\ +V_{i,1} +\end{array}\right] += +\left[\begin{array}{cr} +R_a & -X_{qpp}\\ +X_{qpp} & R_a +\end{array}\right] +\left[\begin{array}{c} +I_{d} \\ +I_{q} +\end{array}\right] +``` +## Generator 3 +Generator connected to bus 3:
+Differential equations: +```math +\dot{\delta}_3 = \omega_3 \cdot \omega_0 +``` +```math +\dot{\omega}_3 = \frac{1}{2H_3}\left( \frac{P_{mech,3} - D_3\omega_3}{1+\omega_3} - T_{elec,3}\right) +``` +```math +\dot{E_{qp,3}} = \frac{1}{T_{dop,3}}\left( E_{fd,3} - \left( E_{qp,3} + X_{d1,3}\left( I_{d,3} + X_{d3,3}\left( E_{qp,3} - \psi_{dp,3} - X_{d2,3}I_{d,3} \right) \right) + \psi_{dpp,3}\cdot k_{sat,3} \right) \right) +``` +```math +\dot{\psi}_{dp,3} = \frac{1}{T_{dopp,3}} \left( E_{qp,3} - \psi_{dp,3} - X_{d2,3} I_{d,3} \right) +``` +```math +\dot{\psi}_{qp,3} = \frac{1}{T_{qopp,3}} \left( E_{dp,3} - \psi_{qp,3} + X_{q2,3}I_{q,3}\right) +``` +```math +\dot{E}_{dp,3} = \frac{1}{T_{qop,3}}\left(-E_{dp,3} + X_{qd,3}\psi_{qpp,3}\cdot k_{sat,3} + X_{q1,3} \left(I_{q,3} - X_{q3,3} \left(E_{dp,3} + I_{q,3}X_{q2,3} - \psi_{qp,3}\right)\right)\right) +``` +Algebraic equations: +```math +V_{d,3} = -\psi_{qpp,3}(\psi_{qp^{\prime},3}E_{dp,3})(1+\omega_3) +``` +```math +V_{q,3} = \psi_{dpp,3}(\psi_{dp^{\prime},3}E_{qp,3})(1+\omega_3) +``` +```math +I_{d,3} = I_{G,r,3}\sin(\delta_3) - I_{G,i,3}\cos(\delta_3) +``` +```math +I_{q,3} = I_{G,r, 3}\cos(\delta_3) + I_{G,i,3}\sin(\delta_3) +``` + +Network interfaces: +```math +V_{d,3} = V_{r,3}\sin(\delta_3) - V_{i,3}\cos(\delta_3) + I_{d,3} R_{a,3} -I_{q,3} X_{qpp,3} +``` +```math +V_{q,3} = V_{r,3}\cos(\delta_3) + V_{i,3} \sin(\delta_3) + I_{d,3}X_{qpp,3} + I_{q,3}R_{a,3} +``` + diff --git a/examples/Bus3Phasor/bus3phasor.cpp b/examples/Bus3Phasor/bus3phasor.cpp new file mode 100644 index 00000000..5394796d --- /dev/null +++ b/examples/Bus3Phasor/bus3phasor.cpp @@ -0,0 +1,545 @@ +#include +#include +#include +#include +#include +#include +using namespace std; + +#define _DEL1 0 +#define _W1 1 +#define _E_qp1 2 +#define _PSI_dp1 3 +#define _PSI_qp1 4 +#define _E_dp1 5 +#define _I_Gr1 6 +#define _I_Gi1 7 +#define _I_r21 8 +#define _I_i21 9 +#define _I_r31 10 +#define _I_i31 11 +#define _I_Lr1 12 +#define _I_Li1 13 +#define _I_Lr2 14 +#define _I_Li2 15 +#define _I_r12 16 +#define _I_i12 17 +#define _I_r32 18 +#define _I_i32 19 +#define _I_Gr3 20 +#define _I_Gi3 21 +#define _I_r13 22 +#define _I_i13 23 +#define _I_r23 24 +#define _I_i23 25 +#define _V_r1 26 +#define _V_i1 27 +#define _V_r2 28 +#define _V_i2 29 +#define _V_r3 30 +#define _V_i3 31 +#define _I_d1 32 +#define _I_q1 33 +#define _DEL2 34 +#define _W2 35 +#define _E_qp2 36 +#define _PSI_dp2 37 +#define _PSI_qp2 38 +#define _E_dp2 39 +#define _I_d2 40 +#define _I_q2 41 + +struct generator_params { + double w_0; + double H; + double D; + double T_dop; + double T_dopp; + double T_qopp; + double T_qop; + double X_d1; + double X_d2; + double X_d3; + double X_d4; + double X_d5; + double X_q1; + double X_q2; + double X_q3; + double X_q4; + double X_q5; + double X_d; + double X_dp; + double X_dpp; + double X_q; + double X_qp; + double X_qpp; + double X_L; + double X_qd; + double R_a; + double S_A; + double S_B; +}; + +/** + * @brief function for calculating the helper variables + */ +int get_psi_qpp(double &psi_qpp, double psi_qp, double X_q4, double E_dp, + double X_q5) { + psi_qpp = -psi_qp * X_q4 - E_dp * X_q5; + return 0; +} + +/** + * @brief function for calculating the helper variables + */ +int get_psi_dpp(double &psi_dpp, double psi_dp, double X_d4, double E_qp, + double X_d5) { + psi_dpp = psi_dp * X_d4 + E_qp * X_d5; + return 0; +} + +int diff_gen(std::vector Y, std::vector &Yp, + std::vector Z, generator_params P, double P_mech, + double E_fd) { + double delta = Y[_DEL1]; + double w = Y[_W1]; + double E_qp = Y[_E_qp1]; + double psi_dp = Y[_PSI_dp1]; + double psi_qp = Y[_PSI_qp1]; + double E_dp = Y[_E_dp1]; + + double I_d = Z[0]; + double I_q = Z[1]; + + double psi_qpp; + double psi_dpp; + double psi_pp; + double k_sat; + double T_elec; + + get_psi_qpp(psi_qpp, psi_qp, P.X_q4, E_dp, P.X_q5); + get_psi_dpp(psi_dpp, psi_dp, P.X_d4, E_qp, P.X_d5); + psi_pp = sqrt(psi_qpp * psi_qpp + psi_dpp * psi_dpp); + k_sat = P.S_B * (psi_pp - P.S_A) * (psi_pp - P.S_A); + T_elec = (psi_dpp - I_d * P.X_dpp) * I_q - (psi_qpp - I_q * P.X_dpp) * I_d; + + Yp[0] = w * P.w_0; + Yp[1] = (1 / (2 * P.H)) * ((P_mech - P.D * w) / (1 + w) - T_elec); + Yp[2] = + (1 / P.T_dop) * + (E_fd - (E_qp + P.X_d1 * (I_d + P.X_d3 * (E_qp - psi_dp - P.X_d2 * I_d)) + + psi_dpp * k_sat)); + Yp[3] = (1 / P.T_dopp) * (E_qp - psi_dp - P.X_d2 * I_d); + Yp[4] = (1 / P.T_qopp) * (E_dp - psi_qp + P.X_q2 * I_q); + Yp[5] = (1 / P.T_qop) * + (-E_dp + P.X_qd * psi_qpp * k_sat + + P.X_q1 * (I_q - P.X_q3 * (E_dp + I_q * P.X_q2 - psi_qp))); + + return 0; +} + +int get_algebraic_res(std::vector &res, std::vector Y, + std::vector Z, std::vector V, + std::vector I, generator_params P, double R_a) { + double delta = Y[_DEL1]; + double w = Y[_W1]; + double E_qp = Y[_E_qp1]; + double psi_dp = Y[_PSI_dp1]; + double psi_qp = Y[_PSI_qp1]; + double E_dp = Y[_E_dp1]; + + double I_d = Z[0]; + double I_q = Z[1]; + + double I_r = I[0]; + double I_i = I[1]; + + double V_r = V[0]; + double V_i = V[1]; + + double psi_qpp; + double psi_dpp; + + get_psi_qpp(psi_qpp, psi_qp, P.X_q4, E_dp, P.X_q5); + get_psi_dpp(psi_dpp, psi_dp, P.X_d4, E_qp, P.X_d5); + + res[0] = I_d - I_r * sin(delta) + I_i * cos(delta); + res[1] = I_q - I_r * cos(delta) - I_i * sin(delta); + res[2] = -psi_qpp * (1 + w) - V_r * sin(delta) + V_i * cos(delta) - + I_d * R_a + I_q * P.X_qpp; + res[3] = psi_dpp * (1 + w) - V_r * cos(delta) - V_i * sin(delta) - + I_d * P.X_qpp - I_q * R_a; + + return 0; +} + +/** + * @brief Takes a vector of currents I and computes the residuals using + * kirchhoff's law + * @param I contains the currents flowing to and from the Bus + * @param r residual output + */ +int buseq(std::vector &I, std::vector &r) { + r[0] = 0.0; + r[1] = 0.0; + + for (size_t i = 0; i < I.size(); i++) { + if (i % 2 == 0) { + r[0] += I[i]; + } else { + r[1] += I[i]; + } + } + + return 0; +} + +/** + * @brief Calculates the current given the voltages + * I[0]: I_r1 + * I[1]: I_i1 + * I[2]: I_r2 + * I[3]: I_i2 + * + * V[0]: V_r1 + * V[1]: V_i1 + * V[2]: V_r2 + * V[3]: V_i2 + */ +int brancheq(std::vector &res, std::vector &I, + std::vector &V, double g, double b, double G, double B) { + res[0] = + -I[0] - (g + G / 2) * V[0] + (b + B / 2) * V[1] + g * V[2] - b * V[3]; + res[1] = + -I[1] - (b + B / 2) * V[0] - (g + G / 2) * V[1] + b * V[2] + g * V[3]; + res[2] = + -I[2] + g * V[0] - b * V[1] - (g + G / 2) * V[2] + (b + B / 2) * V[3]; + res[3] = + -I[3] + b * V[0] + g * V[1] - (b + B / 2) * V[2] - (g + G / 2) * V[3]; + + return 0; +} + +/** + * @brief calculates the current for the load given the voltages and Resistances + * I[0]: I_r1 + * I[1]: I_i1 + * + * V[0]: V_r1 + * V[1]: V_i1 + */ +int loadeq(std::vector &I, std::vector &V, + std::vector &res, double R1, double R2) { + res[0] = I[0] - (V[0] * R1 + V[1] * R2) / (R1 * R1 + R2 * R2); + res[1] = I[1] - (V[1] * R1 - V[0] * R2) / (R1 * R1 + R2 * R2); + + return 0; +} + +int copy_to_res_vector(std::vector &gen_res, + std::vector &towrite, int &at) { + for (int i = 0; i < towrite.size(); i++) { + gen_res[at] = towrite[i]; + at++; + } + return 0; +} + +int diff_res(double t, std::vector &Y, std::vector &Yp, + generator_params params, std::vector &gen_res) { + + int residual_index = 0; + + std::vector sub_Y(Y.begin(), Y.begin() + 6); + std::vector sub_Yp(Yp.begin(), Yp.begin() + 6); + std::vector Z(2); + std::vector diff_Yp(6); + Z[0] = Y[_I_d1]; + Z[1] = Y[_I_q1]; + + double P_mech = 3.0; + double E_fd = 1.0; + + diff_gen(sub_Y, diff_Yp, Z, params, P_mech, E_fd); + + for (int i = 0; i < diff_Yp.size(); i++) { + gen_res[i] = diff_Yp[i] - Yp[i]; + residual_index++; + } + + std::vector algeb_res(4); + std::vector V = {Y[_V_r1], Y[_V_i1]}; + std::vector I = {Y[_I_Gr1], Y[_I_Gi1]}; + double R_a = 2.0; + + get_algebraic_res(algeb_res, sub_Y, Z, V, I, params, R_a); + copy_to_res_vector(gen_res, algeb_res, residual_index); + + // Calculate residuals from bus equations and + // append to the global residual vector gen_res + + std::vector bus_res(2); + + std::vector I_bus_eqn(Y.begin() + 6, Y.begin() + 14); + buseq(I_bus_eqn, bus_res); + copy_to_res_vector(gen_res, bus_res, residual_index); + + I_bus_eqn.assign(Y.begin() + 14, Y.begin() + 20); + buseq(I_bus_eqn, bus_res); + copy_to_res_vector(gen_res, bus_res, residual_index); + + I_bus_eqn.assign(Y.begin() + 20, Y.begin() + 26); + buseq(I_bus_eqn, bus_res); + copy_to_res_vector(gen_res, bus_res, residual_index); + + // Calculate residuals from branch equations and + // append to the global residual vector gen_res + std::vector branch_res(4); + double b = 1.0; + double g = 1.0; + double B = 2.0; + double G = 2.0; + + std::vector I_branch = {Y[_I_r21], Y[_I_i21], Y[_I_r12], Y[_I_i12]}; + std::vector V_branch(Y.begin() + 26, Y.begin() + 30); + brancheq(branch_res, I_branch, V_branch, g, b, G, B); + copy_to_res_vector(gen_res, branch_res, residual_index); + + I_branch = {Y[_I_r32], Y[_I_i32], Y[_I_r23], Y[_I_i23]}; + V_branch.assign(Y.begin() + 28, Y.begin() + 32); + brancheq(branch_res, I_branch, V_branch, g, b, G, B); + copy_to_res_vector(gen_res, branch_res, residual_index); + + I_branch = {Y[_I_r31], Y[_I_i31], Y[_I_r13], Y[_I_i13]}; + V_branch = {Y[_V_r1], Y[_V_i1], Y[_V_r3], Y[_V_i3]}; + brancheq(branch_res, I_branch, V_branch, g, b, G, B); + copy_to_res_vector(gen_res, branch_res, residual_index); + + // Calculate residuals from load equations and + // append to the global residual vector gen_res + std::vector load_res(2); + double R1 = 1.0; + double R2 = 1.0; + + std::vector I_load = {Y[_I_Lr1], Y[_I_Li1]}; + std::vector V_load = {Y[_V_r1], Y[_V_i1]}; + loadeq(I_load, V_load, load_res, R1, R2); + copy_to_res_vector(gen_res, load_res, residual_index); + + I_load = {Y[_I_Lr2], Y[_I_Li2]}; + V_load = {Y[_V_r2], Y[_V_i2]}; + loadeq(I_load, V_load, load_res, R1, R2); + copy_to_res_vector(gen_res, load_res, residual_index); + + sub_Y.assign(Y.begin() + 34, Y.begin() + 40); + sub_Yp.assign(Yp.begin() + 34, Yp.begin() + 40); + Z[0] = Y[_I_d2]; + Z[1] = Y[_I_q2]; + + P_mech = 3.0; + E_fd = 1.0; + + diff_gen(sub_Y, diff_Yp, Z, params, P_mech, E_fd); + + for (int i = 0; i < diff_Yp.size(); i++) { + gen_res[i] = diff_Yp[i] - Yp[i]; + residual_index++; + } + + V = {Y[_V_r3], Y[_V_i3]}; + I = {Y[_I_Gr3], Y[_I_Gi3]}; + R_a = 2.0; + + get_algebraic_res(algeb_res, sub_Y, Z, V, I, params, R_a); + copy_to_res_vector(gen_res, algeb_res, residual_index); + + return 0; +} + +int main(int argc, char const *argv[]) { + // These are for branch 1 + double b1 = 0.0; + double g1 = 0.0; + double B1 = 0.0; + double G1 = 0.0; + + // These are for branch 2 + double b2 = 0.0; + double g2 = 0.0; + double B2 = 0.0; + double G2 = 0.0; + + // These are for branch 3 + double b3 = 0.0; + double g3 = 0.0; + double B3 = 0.0; + double G3 = 0.0; + + generator_params params; + + params.w_0 = 1; + params.H = 1; + params.D = 1; + params.T_dop = 1; + params.T_dopp = 1; + params.T_qopp = 1; + params.T_qop = 1; + params.X_d1 = 1; + params.X_d2 = 1; + params.X_d3 = 1; + params.X_d4 = 1; + params.X_d5 = 1; + params.X_q1 = 1; + params.X_q2 = 1; + params.X_q3 = 1; + params.X_q4 = 1; + params.X_q5 = 1; + params.X_d = 1; + params.X_dp = 1; + params.X_dpp = 1; + params.X_q = 1; + params.X_qp = 1; + params.X_qpp = 1; + params.X_L = 1; + params.X_qd = 1; + params.R_a = 1; + params.S_A = 1; + params.S_B = 1; + + std::vector gen_res(42); + std::vector Y = {1, 1, 2, 2, 2, 1, 1, 2, -6, -5, 4, + 1, -2, -2, 0.5, 4.5, -1, 15, 0, -18, 1, 1, + 2, 2, -3, -2, 2, 1, 1, 2, 3, 4, 1, + 2, 1, 1, 2, 2, 2, 1, 1, 2}; + + std::vector Yp(42); + Yp[0] = 1, Yp[1] = -5, Yp[2] = -65, Yp[3] = -1, Yp[4] = 1, Yp[5] = -48; + Yp[34] = 1, Yp[35] = -5, Yp[36] = -65, Yp[37] = -1, Yp[38] = 1, Yp[39] = -48; + + diff_res(1, Y, Yp, params, gen_res); + + for (int i = 0; i < gen_res.size(); i++) { + std::cout << "index: " << i << " - value: " << gen_res[i] << std::endl; + } + + // vector Y = {1, 1, 2, 2, 2, 1}; + // vector Z = {1, 2}; + // vector Yp(6); + // double P_mech = 3.0; + // double E_fd = 1.0; + + // diff_gen(Y, Yp, Z, params, P_mech, E_fd); + + // for(auto val: Yp) + // { + // std::cout << val << " "; + // } + + // std::cout << "\n--------------------Algebraics----------------------" << + // std::endl; + + // std::vector algeb_res(4); + // std::vector V = {2, 1}; + // std::vector I = {1, 2}; + // double R_a = 2; + // get_algebraic_res(algeb_res, Y, Z, V, I, params, R_a); + + // for(auto val: algeb_res) + // { + // std::cout << val << " "; + // } + + // std::cout << "\n--------------------Buses----------------------" << + // std::endl; + + // std::vector bus_res(2); + // I = {3, 4, -6, -5, 4, 1, -2, -2}; + // buseq(I, bus_res); + // for(auto val: bus_res) + // { + // std::cout << val << " "; + // } + // std::cout < branch_res(4); + // I = {-6, -5, -1, 15}; + // V = {2, 1, 1, 2}; + // brancheq(branch_res, I, V, g, b, B, G); + + // for(auto val: branch_res) + // { + // std::cout << val << " "; + // } + // std::cout < &I, std::vector &V, + // std::vector &res, double R1, double R2) + + // std::vector load_res(2); + // I = {-2, -2}; + // V = {2, 1}; + // double R1 = 1; + // double R2 = 1; + + // loadeq(I, V, load_res, R1, R2); + + // for(auto val: load_res) + // { + // std::cout << val << " "; + // } + // std::cout < +#define _USE_MATH_DEFINES +#include +#include + +// #include +#include +#include +#include +#include +#include + +#include "Model/PhasorDynamics/Branch/Branch.cpp" +#include "Model/PhasorDynamics/Branch/Branch.hpp" +#include "Model/PhasorDynamics/Bus/Bus.cpp" +#include "Model/PhasorDynamics/Bus/Bus.hpp" +#include "Model/PhasorDynamics/Load/Load.cpp" +#include "Model/PhasorDynamics/Load/Load.hpp" +#include "Model/PhasorDynamics/Bus/BusInfinite.cpp" +#include "Model/PhasorDynamics/Bus/BusInfinite.hpp" +#include "Model/PhasorDynamics/BusFault/BusFault.hpp" +#include "Model/PhasorDynamics/SynchronousMachine/GENROUwS/GEN2.hpp" +#include "Model/PhasorDynamics/SystemModel.hpp" +#include "Solver/Dynamic/Ida.cpp" +#include "Solver/Dynamic/Ida.hpp" + +#define _CRT_SECURE_NO_WARNINGS + +int main() +{ + using namespace GridKit::PhasorDynamics; + using namespace AnalysisManager::Sundials; + + printf("Example 1 version GENERATION 2\n"); + + SystemModel sys; + Bus bus1(0.9949877346411762, 0.09999703952427966); + BusInfinite bus2(1.0, 0.0); + Branch branch(&bus1, &bus2, 0.0, 0.1, 0, 0); + GEN2 gen(&bus1, 1,3, 0, 0, 0.2, 0.1, 0.2); + sys.allocate(); + + + /* Connect everything together */ + sys.addBus(&bus1); + sys.addBus(&bus2); + sys.addComponent(&branch); + sys.addComponent(&gen); + + sys.allocate(); + + + double dt = 1.0 / 4.0 / 60.0; + + /* Output file header */ + FILE* f = fopen("example1_v2_results.csv", "w"); + if (!f) + printf("ERROR writing to output file!\n"); + fprintf(f, "%s,%s", "t", "IDA Return Value"); + for (int i = 0; i < sys.size(); ++i) + fprintf(f, ",Y[%d]", i); + for (int i = 0; i < sys.size(); ++i) + fprintf(f, ",Yp[%d]", i); + fprintf(f, "\n"); + + /* Set up simulation */ + Ida ida(&sys); + ida.configureSimulation(); + + /* Run simulation */ + double start = (double) clock(); + ida.printOutputF(0, 0, f); + ida.initializeSimulation(0.0, false); + ida.runSimulationFixed(0.0, dt, 1.0, f); + + + printf("Complete in %.4g seconds\n", (clock() - start) / CLOCKS_PER_SEC); + fclose(f); + + return 0; +} \ No newline at end of file diff --git a/examples/PhasorDynamics/CMakeLists.txt b/examples/PhasorDynamics/CMakeLists.txt new file mode 100644 index 00000000..c758e4e5 --- /dev/null +++ b/examples/PhasorDynamics/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(Example1) diff --git a/examples/PhasorDynamics/Example1/CMakeLists.txt b/examples/PhasorDynamics/Example1/CMakeLists.txt new file mode 100644 index 00000000..78d3b018 --- /dev/null +++ b/examples/PhasorDynamics/Example1/CMakeLists.txt @@ -0,0 +1,9 @@ +add_executable(phasordynamics_example1 example1.cpp) +target_link_libraries(phasordynamics_example1 + GRIDKIT::bus + SUNDIALS::sunlinsolklu + SUNDIALS::core + SUNDIALS::ida + SUNDIALS::idas + SUNDIALS::sunmatrixdense) +install(TARGETS phasordynamics_example1 RUNTIME DESTINATION bin) diff --git a/examples/PhasorDynamics/Example1/example1.cpp b/examples/PhasorDynamics/Example1/example1.cpp new file mode 100644 index 00000000..8e498c83 --- /dev/null +++ b/examples/PhasorDynamics/Example1/example1.cpp @@ -0,0 +1,84 @@ +#include +#define _USE_MATH_DEFINES +#include +#include + +// #include +#include +#include +#include +#include +#include + +#include "Model/PhasorDynamics/Branch/Branch.cpp" +#include "Model/PhasorDynamics/Branch/Branch.hpp" +#include "Model/PhasorDynamics/Bus/Bus.cpp" +#include "Model/PhasorDynamics/Bus/Bus.hpp" +#include "Model/PhasorDynamics/Bus/BusInfinite.cpp" +#include "Model/PhasorDynamics/Bus/BusInfinite.hpp" +#include "Model/PhasorDynamics/BusFault/BusFault.hpp" +#include "Model/PhasorDynamics/SynchronousMachine/GENROUwS/GENROU.hpp" +#include "Model/PhasorDynamics/SystemModel.hpp" +#include "Solver/Dynamic/Ida.cpp" +#include "Solver/Dynamic/Ida.hpp" + +#define _CRT_SECURE_NO_WARNINGS + +int main() +{ + using namespace GridKit::PhasorDynamics; + using namespace AnalysisManager::Sundials; + + printf("Example 1 version 2\n"); + + /* Create model parts */ + SystemModel sys; + Bus bus1(0.9949877346411762, 0.09999703952427966); + BusInfinite bus2(1.0, 0.0); + Branch branch(&bus1, &bus2, 0, 0.1, 0, 0); + BusFault fault(&bus1, 0, 1e-3, 0); + + GENROU gen(&bus1, 1, 1, 0.05013, 3, 0, 0, 7, .04, .05, .75, 2.1, 0.2, 0.18, 0.5, 0.5, 0.18, 0.15, 0, 0); + + /* Connect everything together */ + sys.addBus(&bus1); + sys.addBus(&bus2); + sys.addComponent(&branch); + sys.addComponent(&fault); + sys.addComponent(&gen); + sys.allocate(); + + double dt = 1.0 / 4.0 / 60.0; + + /* Output file header */ + FILE* f = fopen("example1_v2_results.csv", "w"); + if (!f) + printf("ERROR writing to output file!\n"); + fprintf(f, "%s,%s", "t", "IDA Return Value"); + for (int i = 0; i < sys.size(); ++i) + fprintf(f, ",Y[%d]", i); + for (int i = 0; i < sys.size(); ++i) + fprintf(f, ",Yp[%d]", i); + fprintf(f, "\n"); + + /* Set up simulation */ + Ida ida(&sys); + ida.configureSimulation(); + + /* Run simulation */ + double start = (double) clock(); + ida.printOutputF(0, 0, f); + ida.initializeSimulation(0.0, false); + ida.runSimulationFixed(0.0, dt, 1.0, f); + fault.setStatus(1); + ida.initializeSimulation(1.0, false); + ida.runSimulationFixed(1.0, dt, 1.1, f); + fault.setStatus(0); + ida.initializeSimulation(1.1, false); + ida.runSimulationFixed(1.1, dt, 10.0, f); + + printf("Complete in %.4g seconds\n", (clock() - start) / CLOCKS_PER_SEC); + fclose(f); + + return 0; +} diff --git a/src/Model/PhasorDynamics/BusFault/BusFault.cpp b/src/Model/PhasorDynamics/BusFault/BusFault.cpp new file mode 100644 index 00000000..277e80e1 --- /dev/null +++ b/src/Model/PhasorDynamics/BusFault/BusFault.cpp @@ -0,0 +1 @@ +#include "BusFault.hpp" diff --git a/src/Model/PhasorDynamics/BusFault/BusFault.hpp b/src/Model/PhasorDynamics/BusFault/BusFault.hpp new file mode 100644 index 00000000..4de69fbd --- /dev/null +++ b/src/Model/PhasorDynamics/BusFault/BusFault.hpp @@ -0,0 +1,151 @@ +/* Bus Fault Component - Adam Birchfield */ +#pragma once + +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + using ComponentT = Component; + using BaseBusT = BusBase; + + class BusFault : public ComponentT + { + using ComponentT::alpha_; + using ComponentT::f_; + using ComponentT::fB_; + using ComponentT::g_; + using ComponentT::gB_; + using ComponentT::nnz_; + using ComponentT::param_; + using ComponentT::size_; + using ComponentT::tag_; + using ComponentT::time_; + using ComponentT::y_; + using ComponentT::yB_; + using ComponentT::yp_; + using ComponentT::ypB_; + + public: + BusFault(BaseBusT* bus) + : bus_(bus), R_(0), X_(0.01), status_(0), busID_(0) + { + size_ = 0; + } + + BusFault(BaseBusT* bus, double R, double X, int status) + : bus_(bus), R_(R), X_(X), status_(status), busID_(0) + { + size_ = 0; + } + + ~BusFault() + { + } + + int allocate() override + { + return 0; + } + + int initialize() override + { + return 0; + } + + int tagDifferentiable() override + { + return 0; + } + + int evaluateResidual() override + { + if (status_) + { + double B = -X_ / (X_ * X_ + R_ * R_); + double G = R_ / (X_ * X_ + R_ * R_); + + Ir() += -Vr() * G + Vi() * B; + Ii() += -Vr() * B - Vi() * G; + } + return 0; + } + + int evaluateJacobian() override + { + return 0; + } + + int evaluateIntegrand() override + { + return 0; + } + + int initializeAdjoint() override + { + return 0; + } + + int evaluateAdjointResidual() override + { + return 0; + } + + int evaluateAdjointIntegrand() override + { + return 0; + } + + void updateTime(double t, double a) override + { + } + + public: + void setR(double R) + { + R_ = R; + } + + void setX(double X) + { + X_ = X; + } + + void setStatus(int status) + { + status_ = status; + } + + private: + double& Vr() + { + return bus_->Vr(); + } + + double& Vi() + { + return bus_->Vi(); + } + + double& Ir() + { + return bus_->Ir(); + } + + double& Ii() + { + return bus_->Ii(); + } + + private: + BaseBusT* bus_; + double R_; + double X_; + int status_; + const int busID_; + }; + + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/src/Model/PhasorDynamics/BusFault/CMakeLists.txt b/src/Model/PhasorDynamics/BusFault/CMakeLists.txt new file mode 100644 index 00000000..57279fe4 --- /dev/null +++ b/src/Model/PhasorDynamics/BusFault/CMakeLists.txt @@ -0,0 +1,5 @@ +gridkit_add_library(phasor_dynamics_BusFault + SOURCES + BusFault.cpp + OUTPUT_NAME + gridkit_phasor_dynamics_BusFault) diff --git a/src/Model/PhasorDynamics/SynchronousMachine/CMakeLists.txt b/src/Model/PhasorDynamics/SynchronousMachine/CMakeLists.txt index f8239be4..07667d32 100644 --- a/src/Model/PhasorDynamics/SynchronousMachine/CMakeLists.txt +++ b/src/Model/PhasorDynamics/SynchronousMachine/CMakeLists.txt @@ -10,3 +10,5 @@ gridkit_add_library(phasor_dynamics_synchronous_machine SynchronousMachine.cpp OUTPUT_NAME gridkit_phasor_dynamics_synchronous_machine) + +add_subdirectory(GENROUwS) diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/CMakeLists.txt b/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/CMakeLists.txt new file mode 100644 index 00000000..fbe3da17 --- /dev/null +++ b/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/CMakeLists.txt @@ -0,0 +1,5 @@ +gridkit_add_library(phasor_dynamics_GENROU + SOURCES + GENROU.cpp + OUTPUT_NAME + gridkit_phasor_dynamics_GENROU) diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GEN2.cpp b/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GEN2.cpp new file mode 100644 index 00000000..a58c4d30 --- /dev/null +++ b/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GEN2.cpp @@ -0,0 +1 @@ +#include "GEN2.hpp" \ No newline at end of file diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GEN2.hpp b/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GEN2.hpp new file mode 100644 index 00000000..612b598f --- /dev/null +++ b/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GEN2.hpp @@ -0,0 +1,211 @@ +/* GENROU Component - Adam Birchfield */ +#pragma once + +#define _USE_MATH_DEFINES +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + using ComponentT = Component; + using BaseBusT = BusBase; + + class GEN2 : public ComponentT + { + using ComponentT::alpha_; + using ComponentT::f_; + using ComponentT::fB_; + using ComponentT::g_; + using ComponentT::gB_; + using ComponentT::nnz_; + using ComponentT::param_; + using ComponentT::size_; + using ComponentT::tag_; + using ComponentT::time_; + using ComponentT::y_; + using ComponentT::yB_; + using ComponentT::yp_; + using ComponentT::ypB_; + + public: + GEN2(BaseBusT* bus, int unit_id) + : bus_(bus), + unit_id_(unit_id), + busID_(0), + H_(3), + D_(0), + Ra_(0), + Xdp_(0.2) + { + size_ = 5; + set_derived_params(); + } + + GEN2(BaseBusT* bus, + int unit_id, + double H, + double D, + double Ra, + double Xdp, + double pmech_, + double ep_) + : bus_(bus), + unit_id_(unit_id), + busID_(0), + H_(H), + D_(D), + Ra_(Ra), + Xdp_(Xdp), + pmech(pmech_), + ep(ep_) + { + size_ = 5; + set_derived_params(); + } + + void set_derived_params() + { + g = Ra_ / (Ra_ * Ra_ + Xdp_ * Xdp_); + b = -Xdp_ / (Ra_ * Ra_ + Xdp_ * Xdp_); + } + + ~GEN2() + { + } + + int allocate() override + { + f_.resize(size_); + y_.resize(size_); + yp_.resize(size_); + tag_.resize(size_); + fB_.resize(size_); + yB_.resize(size_); + ypB_.resize(size_); + return 0; + } + + int initialize() override + { + return 0; + } + + int tagDifferentiable() override + { + return 0; + } + + int evaluateResidual() override + { + /* Read variables */ + double delta = y_[0]; + double omega = y_[1]; + double telec = y_[2]; + double ir = y_[3]; + double ii = y_[4]; + + /* Read derivatives */ + double delta_dot = yp_[0]; + double omega_dot = yp_[1]; + + /* 6 GENROU differential equations */ + f_[0] = delta_dot - omega * (2 * M_PI * 60); + f_[1] = omega_dot - (1.0 / (2 * H_)) * ((pmech - D_ * omega) / (1 + omega) - telec); + + /* 11 GENROU algebraic equations */ + f_[2] = telec - (1.0/(1.0 + omega))*(g*ep*ep - ep*(cos(delta)*(g*Vr() - b*Vi()) + sin(delta)*(b*Vr() + g*Vi()))); + + f_[3] = ir + g*Vr() - b * Vi() - ep*(g*cos(delta) -b*sin(delta)); + f_[4] = ii + b*Vr() + g * Vi() - ep*(b*cos(delta) + g*sin(delta)); + + + /* Current balance */ + Ir() += - (g*Vr() - b * Vi() - ep*(g*cos(delta) -b*sin(delta))); + Ii() += - (b*Vr() + g * Vi() - ep*(b*cos(delta) + g*sin(delta))); + + // printf("GENROU residual\n"); + // for (int i = 0 ; i < 21; ++i) printf("%d: %g\n", i, f_[i]); + + // printf("GENROU inr %g Vr %g B %g Vi %g G %g\n", inr, Vr(), B_, Vi(), G_); + // printf("GENROU Ii = %g\n", inr - Vr()*B_ - Vi()*G_); + + return 0; + } + + int evaluateJacobian() override + { + /* TODO */ + return 0; + } + + /* Don't know what to do with any of these */ + int evaluateIntegrand() override + { + return 0; + } + + int initializeAdjoint() override + { + return 0; + } + + int evaluateAdjointResidual() override + { + return 0; + } + + int evaluateAdjointIntegrand() override + { + return 0; + } + + void updateTime(double t, double a) override + { + } + + private: + double& Vr() + { + return bus_->Vr(); + } + + double& Vi() + { + return bus_->Vi(); + } + + double& Ir() + { + return bus_->Ir(); + } + + double& Ii() + { + return bus_->Ii(); + } + + private: + /* Identification */ + BaseBusT* bus_; + const int busID_; + int unit_id_; + + /* Input parameters */ + double H_; + double D_; + double Ra_; + double Xdp_; + + /* Derivied parameters */ + double g; + double b; + + /* Setpoints for control variables (determined at initialization) */ + double pmech; + double ep; + }; + + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GENROU.cpp b/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GENROU.cpp new file mode 100644 index 00000000..3540acfd --- /dev/null +++ b/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GENROU.cpp @@ -0,0 +1 @@ +#include "GENROU.hpp" diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GENROU.hpp b/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GENROU.hpp new file mode 100644 index 00000000..2b1a6692 --- /dev/null +++ b/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GENROU.hpp @@ -0,0 +1,402 @@ +/* GENROU Component - Adam Birchfield */ +#pragma once + +#define _USE_MATH_DEFINES +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + using ComponentT = Component; + using BaseBusT = BusBase; + + class GENROU : public ComponentT + { + using ComponentT::alpha_; + using ComponentT::f_; + using ComponentT::fB_; + using ComponentT::g_; + using ComponentT::gB_; + using ComponentT::nnz_; + using ComponentT::param_; + using ComponentT::size_; + using ComponentT::tag_; + using ComponentT::time_; + using ComponentT::y_; + using ComponentT::yB_; + using ComponentT::yp_; + using ComponentT::ypB_; + + public: + GENROU(BaseBusT* bus, int unit_id) + : bus_(bus), + unit_id_(unit_id), + p0_(0), + q0_(0), + busID_(0), + H_(3), + D_(0), + Ra_(0), + Tdop_(7), + Tdopp_(.04), + Tqopp_(.05), + Tqop_(.75), + Xd_(2.1), + Xdp_(0.2), + Xdpp_(0.18), + Xq_(.5), + Xqp_(.5), + Xqpp_(.18), + Xl_(.15), + S10_(0), + S12_(0) + { + size_ = 21; + set_derived_params(); + } + + GENROU(BaseBusT* bus, + int unit_id, + double p0, + double q0, + double H, + double D, + double Ra, + double Tdop, + double Tdopp, + double Tqopp, + double Tqop, + double Xd, + double Xdp, + double Xdpp, + double Xq, + double Xqp, + double Xqpp, + double Xl, + double S10, + double S12) + : bus_(bus), + unit_id_(unit_id), + p0_(p0), + q0_(q0), + busID_(0), + H_(H), + D_(D), + Ra_(Ra), + Tdop_(Tdop), + Tdopp_(Tdopp), + Tqopp_(Tqopp), + Tqop_(Tqop), + Xd_(Xd), + Xdp_(Xdp), + Xdpp_(Xdpp), + Xq_(Xq), + Xqp_(Xqp), + Xqpp_(Xqpp), + Xl_(Xl), + S10_(S10), + S12_(S12) + { + size_ = 21; + set_derived_params(); + } + + void set_derived_params() + { + SA_ = 0; + SB_ = 0; + if (S12_ != 0) + { + double s112 = sqrt(S10_ / S12_); + + SA_ = (1.2 * s112 + 1) / (s112 + 1); + SB_ = (1.2 * s112 - 1) / (s112 - 1); + if (SB_ < SA_) + SA_ = SB_; + SB_ = S12_ / pow(SA_ - 1.2, 2); + } + Xd1_ = Xd_ - Xdp_; + Xd2_ = Xdp_ - Xl_; + Xd3_ = (Xdp_ - Xdpp_) / (Xd2_ * Xd2_); + Xd4_ = (Xdp_ - Xdpp_) / Xd2_; + Xd5_ = (Xdpp_ - Xl_) / Xd2_; + Xq1_ = Xq_ - Xqp_; + Xq2_ = Xqp_ - Xl_; + Xq3_ = (Xqp_ - Xqpp_) / (Xq2_ * Xq2_); + Xq4_ = (Xqp_ - Xqpp_) / Xq2_; + Xq5_ = (Xqpp_ - Xl_) / Xq2_; + Xqd_ = (Xq_ - Xl_) / (Xd_ - Xl_); + G_ = Ra_ / (Ra_ * Ra_ + Xqpp_ * Xqpp_); + B_ = -Xqpp_ / (Ra_ * Ra_ + Xqpp_ * Xqpp_); + } + + ~GENROU() + { + } + + int allocate() override + { + f_.resize(size_); + y_.resize(size_); + yp_.resize(size_); + tag_.resize(size_); + fB_.resize(size_); + yB_.resize(size_); + ypB_.resize(size_); + return 0; + } + + int initialize() override + { + /* Initialization tricks -- assuming NO saturation */ + double vr = Vr(); + double vi = Vi(); + double p = p0_; + double q = q0_; + double vm2 = vr * vr + vi * vi; + double Er = vr + (Ra_ * p * vr + Ra_ * q * vi - Xq_ * p * vi + Xq_ * q * vr) / vm2; + double Ei = vi + (Ra_ * p * vi - Ra_ * q * vr + Xq_ * p * vr + Xq_ * q * vi) / vm2; + double delta = atan2(Ei, Er); + double omega = 0; + double ir = (p * vr + q * vi) / vm2; + double ii = (p * vi - q * vr) / vm2; + double id = ir * sin(delta) - ii * cos(delta); + double iq = ir * cos(delta) + ii * sin(delta); + double vd = vr * sin(delta) - vi * cos(delta) + id * Ra_ - iq * Xqpp_; + double vq = vr * cos(delta) + vi * sin(delta) + id * Xqpp_ - iq * Ra_; + double psiqpp = -vd / (1 + omega); + double psidpp = vq / (1 + omega); + double Te = (psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id; + double psiqp = -(-(Xqp_ - Xl_) * iq + psiqpp * (Xqp_ - Xl_) / (Xqpp_ - Xl_)) + / (1 + (Xqp_ - Xqpp_) / (Xqpp_ - Xl_)); + double Edp = psiqp - (Xqp_ - Xl_) * iq; + double psidp = -((Xdp_ - Xl_) * id - psidpp * (Xdp_ - Xl_) / (Xdpp_ - Xl_)) + / (1 + (Xdp_ - Xdpp_) / (Xdpp_ - Xl_)); + double Eqp = psidp + (Xdp_ - Xl_) * id; + + /* Now we have the state variables, solve for alg. variables */ + double ksat; + double psipp; + + y_[0] = delta; //= 0.55399038; + y_[1] = omega; // = 0; + y_[2] = Eqp; // = 0.995472581; + y_[3] = psidp; // = 0.971299567; + y_[4] = psiqp; // = 0.306880069; + y_[5] = Edp; // = 0; + + y_[6] = psiqpp = -psiqp * Xq4_ - Edp * Xq5_; + y_[7] = psidpp = psidp * Xd4_ + Eqp * Xd5_; + y_[8] = psipp = sqrt(psiqpp * psiqpp + psidpp * psidpp); + y_[9] = ksat = SB_ * pow(psipp - SA_, 2); + y_[10] = vd = -psiqpp * (1 + omega); + y_[11] = vq = psidpp * (1 + omega); + y_[12] = Te = (psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id; + y_[13] = id; + y_[14] = iq; + y_[15] = ir; + y_[16] = ii; + y_[17] = pmech_set_ = Te; + y_[18] = efd_set_ = Eqp + Xd1_ * (id + Xd3_ * (Eqp - psidp - Xd2_ * id)) + psidpp * ksat; + y_[19] = G_ * (vd * sin(delta) + vq * cos(delta)) + - B_ * (vd * -cos(delta) + vq * sin(delta)); /* inort, real */ + y_[20] = B_ * (vd * sin(delta) + vq * cos(delta)) + + G_ * (vd * -cos(delta) + vq * sin(delta)); /* inort, imag */ + + for (int i = 0; i < size_; ++i) + yp_[i] = 0.0; + return 0; + } + + int tagDifferentiable() override + { + for (int i = 0; i < size_; ++i) + { + tag_[i] = i < 6; + } + return 0; + } + + int evaluateResidual() override + { + /* Read variables */ + double delta = y_[0]; + double omega = y_[1]; + double Eqp = y_[2]; + double psidp = y_[3]; + double psiqp = y_[4]; + double Edp = y_[5]; + double psiqpp = y_[6]; + double psidpp = y_[7]; + double psipp = y_[8]; + double ksat = y_[9]; + double vd = y_[10]; + double vq = y_[11]; + double telec = y_[12]; + double id = y_[13]; + double iq = y_[14]; + double ir = y_[15]; + double ii = y_[16]; + double pmech = y_[17]; + double efd = y_[18]; + double inr = y_[19]; + double ini = y_[20]; + double vr = Vr(); + double vi = Vi(); + + /* Read derivatives */ + double delta_dot = yp_[0]; + double omega_dot = yp_[1]; + double Eqp_dot = yp_[2]; + double psidp_dot = yp_[3]; + double psiqp_dot = yp_[4]; + double Edp_dot = yp_[5]; + + /* 6 GENROU differential equations */ + f_[0] = delta_dot - omega * (2 * M_PI * 60); + f_[1] = omega_dot - (1 / (2 * H_)) * ((pmech - D_ * omega) / (1 + omega) - telec); + f_[2] = Eqp_dot - (1 / Tdop_) * (efd - (Eqp + Xd1_ * (id + Xd3_ * (Eqp - psidp - Xd2_ * id)) + psidpp * ksat)); + f_[3] = psidp_dot - (1 / Tdopp_) * (Eqp - psidp - Xd2_ * id); + f_[4] = psiqp_dot - (1 / Tqopp_) * (Edp - psiqp + Xq2_ * iq); + f_[5] = Edp_dot - (1 / Tqop_) * (-Edp + Xqd_ * psiqpp * ksat + Xq1_ * (iq - Xq3_ * (Edp + iq * Xq2_ - psiqp))); + + /* 11 GENROU algebraic equations */ + f_[6] = psiqpp - (-psiqp * Xq4_ - Edp * Xq5_); + f_[7] = psidpp - (psidp * Xd4_ + Eqp * Xd5_); + f_[8] = psipp - sqrt(pow(psidpp, 2.0) + pow(psiqpp, 2.0)); + f_[9] = ksat - SB_ * pow(psipp - SA_, 2.0); + f_[10] = vd + psiqpp * (1 + omega); + f_[11] = vq - psidpp * (1 + omega); + f_[12] = telec - ((psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id); + f_[13] = id - (ir * sin(delta) - ii * cos(delta)); + f_[14] = iq - (ir * cos(delta) + ii * sin(delta)); + f_[15] = ir + G_ * vr - B_ * vi - inr; + f_[16] = ii + B_ * vr + G_ * vi - ini; + + /* 2 GENROU control inputs are set to constant for this example */ + f_[17] = pmech - pmech_set_; + f_[18] = efd - efd_set_; + + /* 2 GENROU current source definitions */ + f_[19] = inr - (G_ * (sin(delta) * vd + cos(delta) * vq) - B_ * (-cos(delta) * vd + sin(delta) * vq)); + f_[20] = ini - (B_ * (sin(delta) * vd + cos(delta) * vq) + G_ * (-cos(delta) * vd + sin(delta) * vq)); + + /* Current balance */ + Ir() += inr - Vr() * G_ + Vi() * B_; + Ii() += ini - Vr() * B_ - Vi() * G_; + + // printf("GENROU residual\n"); + // for (int i = 0 ; i < 21; ++i) printf("%d: %g\n", i, f_[i]); + + // printf("GENROU inr %g Vr %g B %g Vi %g G %g\n", inr, Vr(), B_, Vi(), G_); + // printf("GENROU Ii = %g\n", inr - Vr()*B_ - Vi()*G_); + + return 0; + } + + int evaluateJacobian() override + { + /* TODO */ + return 0; + } + + /* Don't know what to do with any of these */ + int evaluateIntegrand() override + { + return 0; + } + + int initializeAdjoint() override + { + return 0; + } + + int evaluateAdjointResidual() override + { + return 0; + } + + int evaluateAdjointIntegrand() override + { + return 0; + } + + void updateTime(double t, double a) override + { + } + + private: + double& Vr() + { + return bus_->Vr(); + } + + double& Vi() + { + return bus_->Vi(); + } + + double& Ir() + { + return bus_->Ir(); + } + + double& Ii() + { + return bus_->Ii(); + } + + private: + /* Identification */ + BaseBusT* bus_; + const int busID_; + int unit_id_; + + /* Initial terminal conditions */ + double p0_; + double q0_; + + /* Input parameters */ + double H_; + double D_; + double Ra_; + double Tdop_; + double Tdopp_; + double Tqopp_; + double Tqop_; + double Xd_; + double Xdp_; + double Xdpp_; + double Xq_; + double Xqp_; + double Xqpp_; + double Xl_; + double S10_; + double S12_; + + /* Derivied parameters */ + double SA_; + double SB_; + double Xd1_; + double Xd2_; + double Xd3_; + double Xd4_; + double Xd5_; + double Xq1_; + double Xq2_; + double Xq3_; + double Xq4_; + double Xq5_; + double Xqd_; + double G_; + double B_; + + /* Setpoints for control variables (determined at initialization) */ + double pmech_set_; + double efd_set_; + }; + + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/README.md b/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/README.md index 029c89b6..fec362fe 100644 --- a/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/README.md +++ b/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/README.md @@ -93,9 +93,14 @@ then Sat(\psi'')=Sat(\vert V_{r}+jV_{i} \vert) ``` -It is important to point out that finding the initial value of $`\delta`$ for the model without saturation direct method can be used. In case when saturation is considered some "claver" math is needed. Key insight for determining initial $`\delta`$ is that the magnitude of the saturation depends upon the magnitude of $`\psi''`$, which is independent of $`\delta`$. +It is important to point out that finding the initial value of $`\delta`$ for +the model without saturation direct method can be used. In case when saturation +is considered some "claver" math is needed. Key insight for determining initial +$`\delta`$ is that the magnitude of the saturation depends upon the magnitude +of $`\psi''`$, which is independent of $`\delta`$. ```math -\delta=tan^{-1}(\dfrac{K_{sat}V_{iterm}+K_{sat}R_{a}I_{i}+(K_{sat}X''_{d}+X_{q}-X''_{q})I_{r}}{K_{sat}V_{rterm}+K_{sat}R_{a}I_{r}-(K_{sat}X''_{d}+X_{q}-X''_{q})I_{i}}) +\delta=\tan^{-1}\left(\dfrac{K_{sat}V_{iterm}+K_{sat}R_{a}I_{i}+(K_{sat}X''_{d}+X_{q}-X''_{q})I_{r}} + {K_{sat}V_{rterm}+K_{sat}R_{a}I_{r}-(K_{sat}X''_{d}+X_{q}-X''_{q})I_{i}} \right) ``` where ```math diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp index 456c4dcb..acb9dd2f 100644 --- a/src/Solver/Dynamic/Ida.cpp +++ b/src/Solver/Dynamic/Ida.cpp @@ -75,7 +75,7 @@ namespace AnalysisManager // Create vectors to store restart initial condition yy0_ = N_VClone(yy_); checkAllocation((void*) yy0_, "N_VClone"); - yp0_ = N_VClone(yy_); + yp0_ = N_VClone(yp_); checkAllocation((void*) yp0_, "N_VClone"); // Dummy initial time; will be overridden. @@ -205,6 +205,34 @@ namespace AnalysisManager return retval; } + template + int Ida::runSimulationFixed(real_type t0, real_type dt, real_type tmax, FILE* f) + { + int retval = 0; + int iout = 0; + real_type t, tret; + + // printOutputF(t0, 0, f); + for (t = t0 + dt; t < tmax; t += dt) + { + retval = IDASolve(solver_, t, &tret, yy_, yp_, IDA_NORMAL); + checkOutput(retval, "IDASolve"); + printOutputF(t, retval, f); + + if (retval != IDA_SUCCESS) + { + printf("IDA Failure! %dn", retval); + break; + } + } + + model_->updateTime(t, 0.0); + copyVec(yy_, model_->y()); + copyVec(yp_, model_->yp()); + + return retval; + } + template int Ida::runSimulation(real_type tf, int nout) { @@ -651,6 +679,24 @@ namespace AnalysisManager } } + template + void Ida::printOutputF(sunrealtype t, int res, FILE* f) + { + sunrealtype* yval = N_VGetArrayPointer_Serial(yy_); + sunrealtype* ypval = N_VGetArrayPointer_Serial(yp_); + + fprintf(f, "%g,%d", t, res); + for (IdxT i = 0; i < model_->size(); ++i) + { + fprintf(f, ",%g", yval[i]); + } + for (IdxT i = 0; i < model_->size(); ++i) + { + fprintf(f, ",%g", ypval[i]); + } + fprintf(f, "\n"); + } + template void Ida::printOutput(sunrealtype t) { diff --git a/src/Solver/Dynamic/Ida.hpp b/src/Solver/Dynamic/Ida.hpp index 71cc6a40..8ad99b51 100644 --- a/src/Solver/Dynamic/Ida.hpp +++ b/src/Solver/Dynamic/Ida.hpp @@ -36,6 +36,9 @@ namespace AnalysisManager int runSimulation(real_type tf, int nout = 1); int deleteSimulation(); + // TODO: Temporary + int runSimulationFixed(real_type t0, real_type dt, real_type tmax, FILE* f); + int configureQuadrature(); int initializeQuadrature(); int runSimulationQuadrature(real_type tf, int nout = 1); @@ -101,6 +104,7 @@ namespace AnalysisManager void printOutput(sunrealtype t); void printSpecial(sunrealtype t, N_Vector x); void printFinalStats(); + void printOutputF(sunrealtype t, int res, FILE* f); private: static int Residual(sunrealtype t, @@ -109,7 +113,16 @@ namespace AnalysisManager N_Vector rr, void* user_data); - static int Jac(sunrealtype t, sunrealtype cj, N_Vector yy, N_Vector yp, N_Vector resvec, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); + static int Jac(sunrealtype t, + sunrealtype cj, + N_Vector yy, + N_Vector yp, + N_Vector resvec, + SUNMatrix J, + void* user_data, + N_Vector tmp1, + N_Vector tmp2, + N_Vector tmp3); static int Integrand(sunrealtype t, N_Vector yy,