diff --git a/docs/user_manual/tex/dae_model_equations.tex b/docs/user_manual/tex/dae_model_equations.tex index ff8439f3f..9d396bd2b 100644 --- a/docs/user_manual/tex/dae_model_equations.tex +++ b/docs/user_manual/tex/dae_model_equations.tex @@ -43,6 +43,8 @@ \section{Introduction} The document assumes the bus voltage is expressed in cartesian form, i.e. $\bar{V} = V_R + iV_I$, where $V_R$ and $V_I$ are the real and imaginary components of the complex bus voltage $V$. +\textbf{Note:} This document covers the GENCLS, GENROU, and GENSAL generator models, the EXDC1 and ESST1A exciter models, and the WSIEG1 turbine governor model. Additional models implemented in GridPACK (including REGCA1/REGCB1/REGCC1, GDFORM, SEXS, ESST4B, IEEET1, REECA1, REPCA1, TGOV1, GAST, HYGOV, WSHYGP, PSSSIM, and dynamic load models) are documented in their respective source files and PSS/E model references. + \section{Generator models} The inputs and outputs for the generator models are given below. These inputs/outputs representing the coupling with the bus and the generator controls such as exciters, turbine-governors. @@ -105,7 +107,7 @@ \subsubsection{Equations} &\frac{1}{\omega_s}\frac{d{\delta}}{dt}={\Delta\omega}&\\ &\frac{1}{2H}\frac{d{\Delta\omega}}{dt} = \frac{\Pm-D*{\Delta\omega}}{1+{\Delta\omega}}-T_{elec}&\\ &T^{'}_{do}\frac{dE^{'}_q}{dt}= \Efd - L_{ad}\Ifd&\\ -\nonumber &T^{'}_{qo}\frac{dE^{'}_d}{dt}= (X_q - X^{'}_q)\left( I_q -\frac{X^{'}_q - X^{''}_q}{(X^{'}_q - X_l)^2}(E^{'}_d +(X^{'}_q - X_l)I_q - \psi^{'}_q)\right) -&\\ &\qquad E^{'}_d + \frac{X_q - X_l}{X_d - X_l}\psi^{''}_q\frac{Sat(\sqrt{\psi^{''}_d + \psi^{''2}_q})}{\sqrt{\psi^{''}_d + \psi^{''2}_q}} +\nonumber &T^{'}_{qo}\frac{dE^{'}_d}{dt}= (X_q - X^{'}_q)\left( I_q -\frac{X^{'}_q - X^{''}_q}{(X^{'}_q - X_l)^2}(E^{'}_d +(X^{'}_q - X_l)I_q - \psi^{'}_q)\right) -&\\ &\qquad E^{'}_d + \frac{X_q - X_l}{X_d - X_l}\psi^{''}_q\frac{Sat(\sqrt{\psi^{''2}_d + \psi^{''2}_q})}{\sqrt{\psi^{''2}_d + \psi^{''2}_q}} &\\ &T^{''}_{do}\frac{d\psi^{'}_d}{dt} = E^{'}_q - \psi^{'}_d - (X^{'}_{d}-X_{l})I_d&\\ &{T^{''}_{qo}}\frac{d\psi^{'}_q}{dt} = E^{'}_d - \psi^{'}_q + (X^{'}_q-X_{l})I_q& @@ -121,7 +123,7 @@ \subsubsection{Outputs} \item Field current \begin{flalign*} L_{ad}\Ifd = &E^{'}_q + (X_d-X^{'}_{d})\left(I_d + \dfrac{X^{'}_d - X^{"}_d}{(X^{'}_d-X_{l})^2}\left(E^{'}_q - \psi^{'}_d - (X^{'}_d - X_l)I_d\right)\right) + \\ -&\psi^{''}_d\frac{Sat(\sqrt{\psi^{''}_d + \psi^{''2}_q})}{\sqrt{\psi^{''}_d + \psi^{''2}_q}} +&\psi^{''}_d\frac{Sat(\sqrt{\psi^{''2}_d + \psi^{''2}_q})}{\sqrt{\psi^{''2}_d + \psi^{''2}_q}} %L_{ad}\Ifd = E^{'}_q + f_{sat}(E^{'}_q) +(X_d-X^{'}_{d})\left[ \dfrac{X_d - X^{"}_d}{X_d-X_{l}^2}\left(E^{'}_q - \psi^{'}_d - I_d(X^{'}_d - X_l)\right) + I_d\right] % L_{ad}\Ifd=E^{'}_q+(X_d-X_{d})*(I_d+T_{empD})&\\ \end{flalign*} @@ -360,7 +362,7 @@ \subsubsection{Equations} Here, the intermediate variables used in the calculations are as follows: \begin{align*} &\begin{cases} - y_{LL} = x_{LL} + \dfrac{T_2}{T_1}K\dw,&\text{if}~~T_1 > 0 \& T_2 > \\ + y_{LL} = x_{LL} + \dfrac{T_2}{T_1}K\dw,&\text{if}~~T_1 > 0 \& T_2 > 0 \\ y_{LL} = x_{LL},&\text{if}~~T_1 = T_2 = 0 \end{cases}& \\ &u_{GV} = \max(\min(\dfrac{1}{T_3}\left(P_{ref} - x_{GV} - y_{LL}\right),U_o),U_c)& diff --git a/docs/user_manual/tex/ds_mod.tex b/docs/user_manual/tex/ds_mod.tex index f78d0df94..9084e394a 100644 --- a/docs/user_manual/tex/ds_mod.tex +++ b/docs/user_manual/tex/ds_mod.tex @@ -9,6 +9,9 @@ \section{Dynamic Simulation Module using Full Y-Matrix} GENCLS GENSAL GENROU + REGCA1, REGCB1, REGCC1 (grid-following IBR) + GDFORM (grid-forming IBR) + EPRIA1 \end{Verbatim} } @@ -19,6 +22,10 @@ \section{Dynamic Simulation Module using Full Y-Matrix} \begin{Verbatim}[fontseries=b] EXDC1 ESST1A + ESST4B + IEEET1 + SEXS + REECA1 (IBR electrical control) \end{Verbatim} } @@ -29,6 +36,38 @@ \section{Dynamic Simulation Module using Full Y-Matrix} \begin{Verbatim}[fontseries=b] WSIEG1 WSHYGP + TGOV1 + GAST + HYGOV +\end{Verbatim} +} + +Power System Stabilizers: + +{ +\color{red} +\begin{Verbatim}[fontseries=b] + PSSSIM +\end{Verbatim} +} + +Plant Controllers (IBR): + +{ +\color{red} +\begin{Verbatim}[fontseries=b] + REPCA1 (REPCTA1) +\end{Verbatim} +} + +IBR Mechanical Models: + +{ +\color{red} +\begin{Verbatim}[fontseries=b] + WTARA1 + WTDTA1 + WTTQA1 \end{Verbatim} } @@ -49,6 +88,7 @@ \section{Dynamic Simulation Module using Full Y-Matrix} \color{red} \begin{Verbatim}[fontseries=b] ACMTBLU1 + CMLDBLU1 IEEL MOTORW CIM6BL @@ -67,6 +107,7 @@ \section{Dynamic Simulation Module using Full Y-Matrix} IEEE_145.dyr 30 0.005 + false 2.00 @@ -101,7 +142,7 @@ \section{Dynamic Simulation Module using Full Y-Matrix} \end{Verbatim} } -The input for dynamic simulation module is contained in the \texttt{\textbf{Dynamic\_simulation}} block. Two features are important, the blocks describing faults and the blocks describing monitored generators. Faults are described in the \texttt{\textbf{Event}}s block. The code currently only handles faults on branches. Inside the \texttt{\textbf{Events}} block are individual faults, described by a \texttt{\textbf{faultEvent}} block. Multiple \texttt{\textbf{faultEvent}} blocks can be contained within the \texttt{\textbf{Events }}block. As will be described below, it is possible for the faults to be listed in a separate file. This can be convenient for describing a task-based calculation that may contain a lot of faults. The parameters describing the fault include the time (in seconds) that the fault is initiated, the time that it is terminated, the timestep used while integrating the fault and the indices of the two buses at either end of the fault branch. +The input for dynamic simulation module is contained in the \texttt{\textbf{Dynamic\_simulation}} block. The optional \texttt{\textbf{equilibriumInit}} parameter, when set to true, runs the simulation to a steady-state equilibrium before applying faults. This can be useful for verifying that the initial conditions are stable. Two features are important, the blocks describing faults and the blocks describing monitored generators. Faults are described in the \texttt{\textbf{Event}}s block. The code currently only handles faults on branches. Inside the \texttt{\textbf{Events}} block are individual faults, described by a \texttt{\textbf{faultEvent}} block. Multiple \texttt{\textbf{faultEvent}} blocks can be contained within the \texttt{\textbf{Events }}block. As will be described below, it is possible for the faults to be listed in a separate file. This can be convenient for describing a task-based calculation that may contain a lot of faults. The parameters describing the fault include the time (in seconds) that the fault is initiated, the time that it is terminated, the timestep used while integrating the fault and the indices of the two buses at either end of the fault branch. When running a dynamic simulation, it is generally desirable to monitor the behavior of a few generators in the system and this can be done by setting generator watch parameters. The \texttt{\textbf{generatorWatch}} block specifies which generators are to be monitored. Each generator is described within a \texttt{\textbf{generator}} block that contains the index of the bus that the generator is located on and the character string ID of the generator. The results of monitoring the generator are written to the file listed in the \texttt{\textbf{generatorWatchFileName}} field and the frequency for storing generator parameters in this file is set in the \texttt{\textbf{generatorWatchFrequency}} field. This parameter describes the time step interval for writing results (an integer), not the actual time interval. diff --git a/docs/user_manual/tex/se_mod.tex b/docs/user_manual/tex/se_mod.tex index 1f98a8666..eae9514c7 100644 --- a/docs/user_manual/tex/se_mod.tex +++ b/docs/user_manual/tex/se_mod.tex @@ -25,7 +25,7 @@ \section{State Estimation Module} 1 2 BL - -0.2040 + -0.2040 0.0100 @@ -44,7 +44,7 @@ \section{State Estimation Module} \end{Verbatim} } -for the five types of measurements \texttt{\textbf{VM}}, \texttt{\textbf{PIJ}}, \texttt{\textbf{QIJ}}, \texttt{\textbf{PI}}, and \texttt{\textbf{PJ}}. Measurements can appear on any element of the network and multiple measurements are allowed on each element. The state estimation module does not have any error checking ability to determine if there are sufficient measurements to guarantee solvability, if not enough measurements are available then the calculation will simply crash or fail to converge. +The supported measurement types are: \texttt{\textbf{VM}} (voltage magnitude), \texttt{\textbf{VA}} (voltage angle), \texttt{\textbf{PI}} (bus real power injection), \texttt{\textbf{QI}} (bus reactive power injection), \texttt{\textbf{PJ}} (bus real power injection at receiving end), \texttt{\textbf{QJ}} (bus reactive power injection at receiving end), \texttt{\textbf{PIJ}} (branch real power flow from-to), \texttt{\textbf{QIJ}} (branch reactive power flow from-to), \texttt{\textbf{PJI}} (branch real power flow to-from), and \texttt{\textbf{QJI}} (branch reactive power flow to-from). Measurements can appear on any element of the network and multiple measurements are allowed on each element. The state estimation module does not have any error checking ability to determine if there are sufficient measurements to guarantee solvability, if not enough measurements are available then the calculation will simply crash or fail to converge. The state estimation module is represented by the \texttt{\textbf{SEAppModule}} class which is in the \texttt{\textbf{gridpack::state\_estimation}} namespace. The \texttt{\textbf{gridpack.hpp}} file contains a definition for the state estimation network \texttt{\textbf{SENetwork}}. After instantiating an \texttt{\textbf{SEAppModule}} object and a shared pointer to an \texttt{\textbf{SENetwork}}, the state estimation calculation can read in an external network configuration file using the command { @@ -79,6 +79,26 @@ \section{State Estimation Module} The name of the measurement file is in the input deck and a pointer to this file has already been internally cached in the \texttt{\textbf{SEAppModule}} when the network was assigned. The measurement file name is stored in the \texttt{\textbf{measurementList}} field within the \texttt{\textbf{State\_estimation}} block. +The \texttt{\textbf{State\_estimation}} block supports the following optional configuration parameters: + +{ +\color{blue} +\begin{Verbatim}[fontseries=b] + + ieee14.raw + measurements.xml + 1.0e-3 + 20 + 1.0 + 3.0 + 5 + basic + +\end{Verbatim} +} + +The \texttt{\textbf{tolerance}} parameter sets the convergence threshold (default 1.0e-3). The \texttt{\textbf{maxIteration}} parameter limits the number of Newton-Raphson iterations (default 20). The \texttt{\textbf{dampingFactor}} scales the Newton update step (default 1.0). The \texttt{\textbf{badDataThreshold}} sets the normalized residual threshold for bad data detection (default 3.0), and \texttt{\textbf{maxBadDataIterations}} limits the number of bad data elimination passes (default 5). The \texttt{\textbf{diagnosticOutputLevel}} controls the verbosity of diagnostic output (default ``basic''). + The network object can be initialized and the exchange buffers set up by calling the { diff --git a/src/CHANGELOG.md b/src/CHANGELOG.md index bf2985398..f54163753 100644 --- a/src/CHANGELOG.md +++ b/src/CHANGELOG.md @@ -34,7 +34,7 @@ functionality appears in the develop branch. - Q-limit (QLIM) handling with PV-to-PQ bus switching and island detection - RMPCT-based reactive power distribution for multi-generator buses - Fixed shunt status check, Qg distribution when QLIM is reached, and - flat/warm start (VM/VA) + flat/warm start - Switched shunt support - ZIP load model with voltage-dependent load representation - Added Pinj, Qinj to power flow screen outputs @@ -49,7 +49,7 @@ functionality appears in the develop branch. - Generator ID display in contingency results - JSON and CSV export of contingency analysis results - Dynamic Simulation Enhancements - - GENROU PSS (power system stabilizer) integration + - GENROU PSS integration - Iterative equilibrium initialization (XML: equilibriumInit, default false) - Added build trigger when releasing a new version of GridPACK - Changed @@ -58,16 +58,29 @@ functionality appears in the develop branch. - Migrated QLIM parameter to bool type in power flow and CA for consistency - Optimized CMake builds - Fixed - - Fixed PSSSIM washout output divided by Tw, reducing PSS gain by 10x + - Fixed GENROU saturation: use scaled quadratic Se=B*(x-A)^2/x per PowerWorld + standard, evaluate at air-gap flux magnitude, iterative q-axis saturation + in initialization, and corrected LadIfd Se*Psidpp term + - Fixed ESST1A regulator limits from Vrmin/Vrmax to Vamin/Vamax (IEEE 421.5) + - Fixed ESST1A second lead-lag stage bypassed when TA != 0 + - Fixed WSIEG1 parser T3 parameter read from wrong field + - Fixed WSIEG1 NGV lookup table yin[] index + - Fixed REECA1 CurrentLimitLogic VDL1/VDL2 table assignment swapped + - Fixed REPCA1 parser FEMIN read from wrong dyr column + - Fixed REGCA1 Iq_olim sign inconsistency in display functions + - Fixed GENSAL predictor using stale speed deviation for governor + - Fixed GENSAL saturation initialization + - Fixed ForceSerial off-by-one in math library setElementRange half-open range, + causing stale values on multi-process runs with composite load models + - Fixed PSSSIM washout output divided by Tw - Fixed HYGOV fallback setting R instead of r when parameter missing - Fixed Cblock crash when model time constant is zero (e.g., REECA1 with TP=0) - Fixed SEXS exciter EMAX default from 0.0 to 999.0 - Fixed REGCA1 P and Q swapped in output reporting - Fixed null pointer crash on unrecognized generator model in dsf_components - Fixed REECA1 parser reading TIQ parameter from wrong dyr column - - Fixed incomplete REGCA1 data for bus 3433 in 240-bus WECC dyr file - Fixed contingency analysis multi-process deadlock in GA one-sided operations - during result export (getData must be collective for GA progress) + during result export - Fixed 3-winding transformer parsing in PSS/E parsers - Fixed warm start with qlim=false scenario - Fixed incomplete printing of many parallel lines in power flow output diff --git a/src/applications/contingency_analysis/ca_driver.cpp b/src/applications/contingency_analysis/ca_driver.cpp index 642f3cf76..6b60f0d78 100644 --- a/src/applications/contingency_analysis/ca_driver.cpp +++ b/src/applications/contingency_analysis/ca_driver.cpp @@ -1177,46 +1177,12 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) } #endif - // Export CA results to JSON or CSV - if (outputFormat == "csv") { - // Rank 0 writes base case (creates files with headers) - if (world.rank() == 0) { - gridpack::utility::ResultsExporter::writePFCSV(outputFile, - baseCaseResults, "base_case"); - } - world.barrier(); - // Each process writes its contingencies in order (files opened in append mode) - for (int p = 0; p < world.size(); p++) { - if (world.rank() == p) { - for (size_t ci = 0; ci < localContingencies.size(); ci++) { - gridpack::utility::ResultsExporter::writePFCSV( - outputFile, localContingencies[ci].solution, - localContingencies[ci].name); - } - } - world.barrier(); - } - // Rank 0 writes summary CSV using gathered data -#ifdef USE_SUCCESS - if (world.rank() == 0) { - std::string summaryFile = outputFile + "_summary.csv"; - std::ofstream sout(summaryFile.c_str()); - sout << "contingency,type,converged,has_voltage_violation,has_branch_violation\n"; - for (int ci = 0; ci < ntasks; ci++) { - bool converged = (contingency_violation[ci] > 0); - sout << events[ci].p_name << "," - << (events[ci].p_type == Branch ? "branch" : "generator") << "," - << (converged ? "true" : "false") << "," - << ((contingency_violation[ci] == 2 || contingency_violation[ci] == 4) - ? "true" : "false") << "," - << ((contingency_violation[ci] == 3 || contingency_violation[ci] == 4) - ? "true" : "false") << "\n"; - } - sout.close(); - } -#endif - } + // Sync GA before MPI collectives to flush any pending one-sided operations + world.sync(); + // Export CA results to JSON or CSV. + // Only rank 0 writes output files. Non-zero ranks send their serialized + // data to rank 0 using point-to-point MPI send/recv. if (outputFormat == "json") { // Each process serializes its contingency results as JSON text std::ostringstream localJsonStream; @@ -1229,41 +1195,40 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) } std::string localJson = localJsonStream.str(); - // Gather JSON fragment sizes - int localLen = static_cast(localJson.size()); - std::vector allLens(world.size()); - MPI_Allgather(&localLen, 1, MPI_INT, &allLens[0], 1, MPI_INT, - static_cast(world)); - - // Compute displacements - std::vector displs(world.size()); - int totalLen = 0; - for (int p = 0; p < world.size(); p++) { - displs[p] = totalLen; - totalLen += allLens[p]; + // Gather all JSON fragments on rank 0 using point-to-point send/recv + MPI_Comm mpi_comm = static_cast(world); + std::vector allFragments(world.size()); + allFragments[0] = localJson; // rank 0's own data + if (world.rank() == 0) { + for (int p = 1; p < world.size(); p++) { + int len; + MPI_Recv(&len, 1, MPI_INT, p, 0, mpi_comm, MPI_STATUS_IGNORE); + allFragments[p].resize(len); + if (len > 0) { + MPI_Recv(&allFragments[p][0], len, MPI_CHAR, p, 1, mpi_comm, + MPI_STATUS_IGNORE); + } + } + } else { + int len = static_cast(localJson.size()); + MPI_Send(&len, 1, MPI_INT, 0, 0, mpi_comm); + if (len > 0) { + MPI_Send(const_cast(localJson.c_str()), len, MPI_CHAR, 0, 1, + mpi_comm); + } } - // Gather all JSON fragments on rank 0 - std::vector allJsonBuf; - if (world.rank() == 0) allJsonBuf.resize(totalLen); - MPI_Gatherv(localJson.c_str(), localLen, MPI_CHAR, - world.rank() == 0 ? &allJsonBuf[0] : NULL, - &allLens[0], &displs[0], MPI_CHAR, - 0, static_cast(world)); - // Rank 0 writes the final JSON file if (world.rank() == 0) { std::string jsonFile = outputFile + ".json"; std::ofstream jout(jsonFile.c_str()); - // Write base case and open contingencies array gridpack::utility::ResultsExporter::writeCAJSONHeader(jout, baseCaseResults); - // Write contingencies from gathered fragments bool firstFragment = true; for (int p = 0; p < world.size(); p++) { - if (allLens[p] > 0) { + if (!allFragments[p].empty()) { if (!firstFragment) jout << ",\n"; - jout.write(&allJsonBuf[displs[p]], allLens[p]); + jout << allFragments[p]; firstFragment = false; } } @@ -1273,6 +1238,161 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) } } + if (outputFormat == "csv") { + // Each process serializes its contingency CSV data into 4 strings + // (buses, branches, generators, convergence rows without headers) + std::ostringstream localBus, localBranch, localGen, localConv; + localBus << std::fixed; + localBranch << std::fixed; + localGen << std::fixed; + for (size_t ci = 0; ci < localContingencies.size(); ci++) { + const gridpack::utility::ContingencyResult& ct = localContingencies[ci]; + const gridpack::utility::PowerFlowResults& r = ct.solution; + for (size_t bi = 0; bi < r.buses.size(); bi++) { + const gridpack::utility::BusResult& b = r.buses[bi]; + localBus << ct.name << "," + << b.busId << "," << b.type << "," + << b.area << "," << b.zone << "," + << std::setprecision(2) << b.baseKV << "," + << std::setprecision(6) << b.voltage << "," + << std::setprecision(6) << b.angle << "," + << std::setprecision(4) << b.pInjection << "," + << std::setprecision(4) << b.qInjection << "," + << std::setprecision(4) << b.pLoad << "," + << std::setprecision(4) << b.qLoad << "," + << std::setprecision(4) << b.pGen << "," + << std::setprecision(4) << b.qGen << "," + << std::setprecision(4) << b.shuntMvar << "\n"; + } + for (size_t bi = 0; bi < r.branches.size(); bi++) { + const gridpack::utility::BranchResult& br = r.branches[bi]; + localBranch << ct.name << "," + << br.fromBus << "," << br.toBus << "," + << br.circuitId << "," + << std::setprecision(4) << br.pFrom << "," + << std::setprecision(4) << br.qFrom << "," + << std::setprecision(4) << br.pTo << "," + << std::setprecision(4) << br.qTo << "," + << std::setprecision(4) << br.pLoss << "," + << std::setprecision(4) << br.qLoss << "," + << std::setprecision(4) << br.mvaFrom << "," + << std::setprecision(4) << br.mvaTo << "," + << std::setprecision(4) << br.rateA << "," + << std::setprecision(2) << br.loadingPercent << "\n"; + } + for (size_t gi = 0; gi < r.generators.size(); gi++) { + const gridpack::utility::GeneratorResult& g = r.generators[gi]; + localGen << ct.name << "," + << g.busId << "," << g.genId << "," + << std::setprecision(4) << g.pGen << "," + << std::setprecision(4) << g.qGen << "," + << std::setprecision(4) << g.qMax << "," + << std::setprecision(4) << g.qMin << "," + << std::setprecision(6) << g.voltageSetpoint << "," + << g.status << "\n"; + } + localConv << ct.name << "," + << (r.convergence.converged ? "true" : "false") << "," + << r.convergence.iterations << "," + << std::scientific << r.convergence.finalTolerance << "," + << std::fixed + << r.convergence.finalMismatch.maxPBus << "," + << std::setprecision(4) << r.convergence.finalMismatch.maxPMismatch << "," + << r.convergence.finalMismatch.maxQBus << "," + << std::setprecision(4) << r.convergence.finalMismatch.maxQMismatch << "\n"; + } + + // Gather all CSV fragments on rank 0 using point-to-point send/recv + MPI_Comm mpi_comm = static_cast(world); + std::vector allBus(world.size()), allBranch(world.size()); + std::vector allGen(world.size()), allConv(world.size()); + allBus[0] = localBus.str(); + allBranch[0] = localBranch.str(); + allGen[0] = localGen.str(); + allConv[0] = localConv.str(); + if (world.rank() == 0) { + for (int p = 1; p < world.size(); p++) { + int lens[4]; + MPI_Recv(lens, 4, MPI_INT, p, 0, mpi_comm, MPI_STATUS_IGNORE); + allBus[p].resize(lens[0]); + allBranch[p].resize(lens[1]); + allGen[p].resize(lens[2]); + allConv[p].resize(lens[3]); + if (lens[0] > 0) + MPI_Recv(&allBus[p][0], lens[0], MPI_CHAR, p, 1, mpi_comm, + MPI_STATUS_IGNORE); + if (lens[1] > 0) + MPI_Recv(&allBranch[p][0], lens[1], MPI_CHAR, p, 2, mpi_comm, + MPI_STATUS_IGNORE); + if (lens[2] > 0) + MPI_Recv(&allGen[p][0], lens[2], MPI_CHAR, p, 3, mpi_comm, + MPI_STATUS_IGNORE); + if (lens[3] > 0) + MPI_Recv(&allConv[p][0], lens[3], MPI_CHAR, p, 4, mpi_comm, + MPI_STATUS_IGNORE); + } + } else { + std::string sBus = localBus.str(), sBranch = localBranch.str(); + std::string sGen = localGen.str(), sConv = localConv.str(); + int lens[4] = {(int)sBus.size(), (int)sBranch.size(), + (int)sGen.size(), (int)sConv.size()}; + MPI_Send(lens, 4, MPI_INT, 0, 0, mpi_comm); + if (lens[0] > 0) + MPI_Send(const_cast(sBus.c_str()), lens[0], MPI_CHAR, 0, 1, + mpi_comm); + if (lens[1] > 0) + MPI_Send(const_cast(sBranch.c_str()), lens[1], MPI_CHAR, 0, 2, + mpi_comm); + if (lens[2] > 0) + MPI_Send(const_cast(sGen.c_str()), lens[2], MPI_CHAR, 0, 3, + mpi_comm); + if (lens[3] > 0) + MPI_Send(const_cast(sConv.c_str()), lens[3], MPI_CHAR, 0, 4, + mpi_comm); + } + + // Rank 0 writes the CSV files + if (world.rank() == 0) { + // Write base case first (creates files with headers) + gridpack::utility::ResultsExporter::writePFCSV(outputFile, + baseCaseResults, "base_case"); + // Append contingency data from all processes + { + std::ofstream out((outputFile + "_buses.csv").c_str(), std::ios::app); + for (size_t p = 0; p < allBus.size(); p++) out << allBus[p]; + } + { + std::ofstream out((outputFile + "_branches.csv").c_str(), std::ios::app); + for (size_t p = 0; p < allBranch.size(); p++) out << allBranch[p]; + } + { + std::ofstream out((outputFile + "_generators.csv").c_str(), std::ios::app); + for (size_t p = 0; p < allGen.size(); p++) out << allGen[p]; + } + { + std::ofstream out((outputFile + "_convergence.csv").c_str(), std::ios::app); + for (size_t p = 0; p < allConv.size(); p++) out << allConv[p]; + } +#ifdef USE_SUCCESS + // Write summary CSV + std::string summaryFile = outputFile + "_summary.csv"; + std::ofstream sout(summaryFile.c_str()); + sout << "contingency,type,converged,has_voltage_violation,has_branch_violation\n"; + for (int ci = 0; ci < ntasks; ci++) { + bool converged = (contingency_violation[ci] > 0); + sout << events[ci].p_name << "," + << (events[ci].p_type == Branch ? "branch" : "generator") << "," + << (converged ? "true" : "false") << "," + << ((contingency_violation[ci] == 2 || contingency_violation[ci] == 4) + ? "true" : "false") << "," + << ((contingency_violation[ci] == 3 || contingency_violation[ci] == 4) + ? "true" : "false") << "\n"; + } + sout.close(); +#endif + } + } + // Print out statistics on contingencies #ifdef USE_STATBLOCK int t_stats = timer->createCategory("Write Statistics"); diff --git a/src/applications/data_sets/dyr/case9_GENSAL.dyr b/src/applications/data_sets/dyr/case9_GENSAL.dyr index 3f63c5ec3..63f6970d5 100644 --- a/src/applications/data_sets/dyr/case9_GENSAL.dyr +++ b/src/applications/data_sets/dyr/case9_GENSAL.dyr @@ -1,3 +1,3 @@ 1,'GENSAL',1 , 8.96, 0.05, 0.05, 23.64, 0.0125, 0.1460, 0.0969, 0.0608, 0.05, 0.026, 0.000, 0.000,/USRWHT -2,'GENSAL',1 , 6.0, 0.037, 0.074, 6.470, 0.0068, 0.8958, 0.8645, 0.1198, 0.12, 0.06, 0.000, 0.000,/USRWHT +2,'GENSAL',1 , 6.0, 0.037, 0.074, 6.470, 0.0068, 0.8958, 0.8645, 0.18, 0.12, 0.06, 0.000, 0.000,/USRWHT 3,'GENSAL',1 , 5.89, 0.032, 0.079, 3.010, 0.0048, 1.3125, 1.2578, 0.1813, 0.1, 0.08, 0.000, 0.000,/USRWHT diff --git a/src/applications/data_sets/dyr/case9_GENSAL_ESST1A.dyr b/src/applications/data_sets/dyr/case9_GENSAL_ESST1A.dyr index c5d10dae8..ed7658a05 100644 --- a/src/applications/data_sets/dyr/case9_GENSAL_ESST1A.dyr +++ b/src/applications/data_sets/dyr/case9_GENSAL_ESST1A.dyr @@ -1,5 +1,5 @@ 1,'GENSAL',1 , 8.96, 0.05, 0.05, 23.64, 0.0125, 0.1460, 0.0969, 0.0608, 0.05, 0.026, 0.000, 0.000,/USRWHT -2,'GENSAL',1 , 6.0, 0.037, 0.074, 6.470, 0.0068, 0.8958, 0.8645, 0.1198, 0.12, 0.06, 0.000, 0.000,/USRWHT +2,'GENSAL',1 , 6.0, 0.037, 0.074, 6.470, 0.0068, 0.8958, 0.8645, 0.18, 0.12, 0.06, 0.000, 0.000,/USRWHT 3,'GENSAL',1 , 5.89, 0.032, 0.079, 3.010, 0.0048, 1.3125, 1.2578, 0.1813, 0.1, 0.08, 0.000, 0.000,/USRWHT 1, 'ESST1A', 1, 1, 1, 0.0, 999.000, -999.000, 0.510000, 2.01000, 0.00000, 0.00000, 178.900, 0.290000E-01, 999.000, -999.000, 24.48000, -21.79000, 0.110000, 0.00000, 1.00000, 0.00000, 2.80000 / 2, 'ESST1A', 1, 1, 1, 0.0, 999.000, -999.000, 0.510000, 2.01000, 0.00000, 0.00000, 178.900, 0.290000E-01, 999.000, -999.000, 24.48000, -21.79000, 0.110000, 0.00000, 1.00000, 0.00000, 2.80000 / diff --git a/src/applications/data_sets/dyr/case9_GENSAL_ESST1A_WSIEG1.dyr b/src/applications/data_sets/dyr/case9_GENSAL_ESST1A_WSIEG1.dyr index e93ebb9a4..c1271f0ac 100644 --- a/src/applications/data_sets/dyr/case9_GENSAL_ESST1A_WSIEG1.dyr +++ b/src/applications/data_sets/dyr/case9_GENSAL_ESST1A_WSIEG1.dyr @@ -1,5 +1,5 @@ 1,'GENSAL',1 , 8.96, 0.05, 0.05, 23.64, 0.0125, 0.1460, 0.0969, 0.0608, 0.05, 0.026, 0.000, 0.000,/USRWHT -2,'GENSAL',1 , 6.0, 0.037, 0.074, 6.470, 0.0068, 0.8958, 0.8645, 0.1198, 0.12, 0.06, 0.000, 0.000,/USRWHT +2,'GENSAL',1 , 6.0, 0.037, 0.074, 6.470, 0.0068, 0.8958, 0.8645, 0.18, 0.12, 0.06, 0.000, 0.000,/USRWHT 3,'GENSAL',1 , 5.89, 0.032, 0.079, 3.010, 0.0048, 1.3125, 1.2578, 0.1813, 0.1, 0.08, 0.000, 0.000,/USRWHT 1, 'ESST1A', 1, 1, 1, 0.0, 999.000, -999.000, 0.510000, 2.01000, 0.00000, 0.00000, 178.900, 0.290000E-01, 999.000, -999.000, 24.48000, -21.79000, 0.110000, 0.00000, 1.00000, 0.00000, 2.80000 / 2, 'ESST1A', 1, 1, 1, 0.0, 999.000, -999.000, 0.510000, 2.01000, 0.00000, 0.00000, 178.900, 0.290000E-01, 999.000, -999.000, 24.48000, -21.79000, 0.110000, 0.00000, 1.00000, 0.00000, 2.80000 / diff --git a/src/applications/data_sets/dyr/kundur-twoarea_4renewable_mech.dyr b/src/applications/data_sets/dyr/kundur-twoarea_4renewable_mech.dyr index eb8147c66..31aeb5de0 100644 --- a/src/applications/data_sets/dyr/kundur-twoarea_4renewable_mech.dyr +++ b/src/applications/data_sets/dyr/kundur-twoarea_4renewable_mech.dyr @@ -42,7 +42,7 @@ 1.0000 0.50000E-01 0.25000 -1.0000 1.0000 99.000 -99.000 1.0000 0.0000 0.10000 0.0000 0.0000 / -4 'WTARA1' '1 ' 1 0.007 0 / +4 'WTARA1' '1 ' 1 0.007 / 4 'WTPTA1' '1 ' 25 150 30 0 0 0.3 27 0 10 -10 / 4 'WTTQA1' '1 ' 1 0.5 2.5 0.066667 60 1.1 0 0.2 0.58 0.4 0.72 0.72 0.866 0.86 1 0 / -4 'WTDTA1' '1 ' 5 0.5 0.0 0.0 / +4 'WTDTA1' '1 ' 5 0.5 0.0 0.0 0.0 / diff --git a/src/applications/data_sets/input/ds/input_240bus.xml b/src/applications/data_sets/input/ds/input_240bus.xml index 5dee69aac..3597485d0 100644 --- a/src/applications/data_sets/input/ds/input_240bus.xml +++ b/src/applications/data_sets/input/ds/input_240bus.xml @@ -57,6 +57,7 @@ 10 0.005 + true 1.0 @@ -76,7 +77,7 @@ 1 - gen_watch.csv + gen_watch_240bus.csv diff --git a/src/applications/data_sets/input/ds/input_3000.xml b/src/applications/data_sets/input/ds/input_3000.xml index b96c21f2e..965431b71 100644 --- a/src/applications/data_sets/input/ds/input_3000.xml +++ b/src/applications/data_sets/input/ds/input_3000.xml @@ -56,7 +56,7 @@ classical_model_3000bus.dyr - 20 + 5 0.001 @@ -83,9 +83,10 @@ - -ksp_type preonly + -ksp_type richardson -pc_type lu - -pc_factor_mat_ordering_type amd + -pc_factor_mat_solver_type superlu_dist + -ksp_max_it 1 diff --git a/src/applications/data_sets/input/ds/input_300_cmpld.xml b/src/applications/data_sets/input/ds/input_300_cmpld.xml index ef2d4e9c8..9f8df182b 100644 --- a/src/applications/data_sets/input/ds/input_300_cmpld.xml +++ b/src/applications/data_sets/input/ds/input_300_cmpld.xml @@ -59,6 +59,7 @@ 300bus_detail_model_cmpld_combine.dyr 10 0.001 + true 2.0 @@ -75,7 +76,7 @@ 1 - 300bus_cmpld_gridpack_test_gen10063.csv + gen_watch_300_cmpld.csv 1.0E-12 true @@ -85,9 +86,10 @@ - -ksp_type preonly + -ksp_type richardson -pc_type lu - -pc_factor_mat_ordering_type amd + -pc_factor_mat_solver_type superlu_dist + -ksp_max_it 1 diff --git a/src/applications/data_sets/input/ds/input_9b3g.xml b/src/applications/data_sets/input/ds/input_9b3g.xml index 0f3577201..4b6cf93cb 100644 --- a/src/applications/data_sets/input/ds/input_9b3g.xml +++ b/src/applications/data_sets/input/ds/input_9b3g.xml @@ -70,7 +70,7 @@ 1 - gen_watch.csv + gen_watch_9b3g.csv diff --git a/src/applications/data_sets/input/ds/input_twoarea.xml b/src/applications/data_sets/input/ds/input_twoarea.xml index 298366a5c..651e978ef 100644 --- a/src/applications/data_sets/input/ds/input_twoarea.xml +++ b/src/applications/data_sets/input/ds/input_twoarea.xml @@ -54,7 +54,7 @@ kundur-twoarea.dyr - 30 + 10 0.005 @@ -77,9 +77,13 @@ 3 1 + + 4 + 1 + 1 - gen_watch.csv + gen_watch_twoarea.csv diff --git a/src/applications/data_sets/input/ds/input_unit_GAST.xml b/src/applications/data_sets/input/ds/input_unit_GAST.xml index 507654044..ddc24260a 100644 --- a/src/applications/data_sets/input/ds/input_unit_GAST.xml +++ b/src/applications/data_sets/input/ds/input_unit_GAST.xml @@ -48,7 +48,7 @@ case9_GAST.dyr 1.0 - 0.01 + 0.001 0.03 @@ -62,9 +62,17 @@ 1 1 + + 2 + 1 + + + 3 + 1 + 1 - gen_watch.csv + gen_watch_GAST.csv -ksp_type richardson diff --git a/src/applications/data_sets/input/ds/input_unit_GENROU.xml b/src/applications/data_sets/input/ds/input_unit_GENROU.xml index 870be5feb..adf80c032 100644 --- a/src/applications/data_sets/input/ds/input_unit_GENROU.xml +++ b/src/applications/data_sets/input/ds/input_unit_GENROU.xml @@ -48,7 +48,8 @@ case9_GENROU.dyr 1.0 - 0.01 + 0.001 + true 0.03 @@ -62,9 +63,17 @@ 1 1 + + 2 + 1 + + + 3 + 1 + 1 - gen_watch.csv + gen_watch_GENROU.csv -ksp_type richardson diff --git a/src/applications/data_sets/input/ds/input_unit_GENROU_ESST1A.xml b/src/applications/data_sets/input/ds/input_unit_GENROU_ESST1A.xml index cd0582475..4e72ce558 100644 --- a/src/applications/data_sets/input/ds/input_unit_GENROU_ESST1A.xml +++ b/src/applications/data_sets/input/ds/input_unit_GENROU_ESST1A.xml @@ -48,7 +48,7 @@ case9_GENROU_ESST1A.dyr 1.0 - 0.01 + 0.001 0.03 @@ -62,9 +62,17 @@ 1 1 + + 2 + 1 + + + 3 + 1 + 1 - gen_watch.csv + gen_watch_GENROU_ESST1A.csv -ksp_type richardson diff --git a/src/applications/data_sets/input/ds/input_unit_GENROU_ESST1A_WSIEG1.xml b/src/applications/data_sets/input/ds/input_unit_GENROU_ESST1A_WSIEG1.xml index 950813863..6b97133fb 100644 --- a/src/applications/data_sets/input/ds/input_unit_GENROU_ESST1A_WSIEG1.xml +++ b/src/applications/data_sets/input/ds/input_unit_GENROU_ESST1A_WSIEG1.xml @@ -48,7 +48,7 @@ case9_GENROU_ESST1A_WSIEG1.dyr 1.0 - 0.01 + 0.001 0.03 @@ -62,9 +62,17 @@ 1 1 + + 2 + 1 + + + 3 + 1 + 1 - gen_watch.csv + gen_watch_GENROU_ESST1A_WSIEG1.csv -ksp_type richardson diff --git a/src/applications/data_sets/input/ds/input_unit_GENROU_EXDC1.xml b/src/applications/data_sets/input/ds/input_unit_GENROU_EXDC1.xml index cad85f4c1..b7ebd0eeb 100644 --- a/src/applications/data_sets/input/ds/input_unit_GENROU_EXDC1.xml +++ b/src/applications/data_sets/input/ds/input_unit_GENROU_EXDC1.xml @@ -48,7 +48,7 @@ case9_GENROU_EXDC1.dyr 1.0 - 0.01 + 0.001 0.03 @@ -62,9 +62,17 @@ 1 1 + + 2 + 1 + + + 3 + 1 + 1 - gen_watch.csv + gen_watch_GENROU_EXDC1.csv -ksp_type richardson diff --git a/src/applications/data_sets/input/ds/input_unit_GENROU_EXDC1_WSIEG1.xml b/src/applications/data_sets/input/ds/input_unit_GENROU_EXDC1_WSIEG1.xml index 6d5efe97b..93ef216d2 100644 --- a/src/applications/data_sets/input/ds/input_unit_GENROU_EXDC1_WSIEG1.xml +++ b/src/applications/data_sets/input/ds/input_unit_GENROU_EXDC1_WSIEG1.xml @@ -48,7 +48,7 @@ case9_GENROU_EXDC1_WSIEG1.dyr 1.0 - 0.01 + 0.001 0.03 @@ -62,9 +62,17 @@ 1 1 + + 2 + 1 + + + 3 + 1 + 1 - gen_watch.csv + gen_watch_GENROU_EXDC1_WSIEG1.csv -ksp_type richardson diff --git a/src/applications/data_sets/input/ds/input_unit_GENSAL.xml b/src/applications/data_sets/input/ds/input_unit_GENSAL.xml index 73f53f861..a6fc0137f 100644 --- a/src/applications/data_sets/input/ds/input_unit_GENSAL.xml +++ b/src/applications/data_sets/input/ds/input_unit_GENSAL.xml @@ -48,7 +48,7 @@ case9_GENSAL.dyr 1.0 - 0.01 + 0.001 0.03 @@ -62,9 +62,17 @@ 1 1 + + 2 + 1 + + + 3 + 1 + 1 - gen_watch.csv + gen_watch_GENSAL.csv -ksp_type richardson diff --git a/src/applications/data_sets/input/ds/input_unit_GENSAL_ESST1A.xml b/src/applications/data_sets/input/ds/input_unit_GENSAL_ESST1A.xml index f9f3dbe10..719fb2aaf 100644 --- a/src/applications/data_sets/input/ds/input_unit_GENSAL_ESST1A.xml +++ b/src/applications/data_sets/input/ds/input_unit_GENSAL_ESST1A.xml @@ -48,7 +48,7 @@ case9_GENSAL_ESST1A.dyr 1.0 - 0.01 + 0.001 0.03 @@ -62,9 +62,17 @@ 1 1 + + 2 + 1 + + + 3 + 1 + 1 - gen_watch.csv + gen_watch_GENSAL_ESST1A.csv -ksp_type richardson diff --git a/src/applications/data_sets/input/ds/input_unit_GENSAL_ESST1A_WSIEG1.xml b/src/applications/data_sets/input/ds/input_unit_GENSAL_ESST1A_WSIEG1.xml index cf99d9300..5072963d7 100644 --- a/src/applications/data_sets/input/ds/input_unit_GENSAL_ESST1A_WSIEG1.xml +++ b/src/applications/data_sets/input/ds/input_unit_GENSAL_ESST1A_WSIEG1.xml @@ -48,7 +48,7 @@ case9_GENSAL_ESST1A_WSIEG1.dyr 1.0 - 0.01 + 0.001 0.03 @@ -62,9 +62,17 @@ 1 1 + + 2 + 1 + + + 3 + 1 + 1 - gen_watch.csv + gen_watch_GENSAL_ESST1A_WSIEG1.csv -ksp_type richardson diff --git a/src/applications/data_sets/input/ds/input_unit_HYGOV.xml b/src/applications/data_sets/input/ds/input_unit_HYGOV.xml index 57891d70f..03b9ac7f3 100644 --- a/src/applications/data_sets/input/ds/input_unit_HYGOV.xml +++ b/src/applications/data_sets/input/ds/input_unit_HYGOV.xml @@ -48,7 +48,7 @@ case9_HYGOV.dyr 1.0 - 0.01 + 0.001 0.03 @@ -62,9 +62,17 @@ 1 1 + + 2 + 1 + + + 3 + 1 + 1 - gen_watch.csv + gen_watch_HYGOV.csv -ksp_type richardson diff --git a/src/applications/data_sets/input/ds/input_unit_REGCA1_REECA1_REPCA1.xml b/src/applications/data_sets/input/ds/input_unit_REGCA1_REECA1_REPCA1.xml index cc0f05d73..0a72af67b 100644 --- a/src/applications/data_sets/input/ds/input_unit_REGCA1_REECA1_REPCA1.xml +++ b/src/applications/data_sets/input/ds/input_unit_REGCA1_REECA1_REPCA1.xml @@ -48,7 +48,7 @@ case9_REGCA1_REECA1_REPCA1.dyr 1.0 - 0.01 + 0.001 0.03 @@ -62,9 +62,17 @@ 1 1 + + 2 + 1 + + + 3 + 1 + 1 - gen_watch.csv + gen_watch_REGCA1_REECA1_REPCA1.csv -ksp_type richardson diff --git a/src/applications/data_sets/input/ds/input_unit_SEXS.xml b/src/applications/data_sets/input/ds/input_unit_SEXS.xml index 9d4d7b6c9..b60b1ba3c 100644 --- a/src/applications/data_sets/input/ds/input_unit_SEXS.xml +++ b/src/applications/data_sets/input/ds/input_unit_SEXS.xml @@ -48,7 +48,7 @@ case9_SEXS.dyr 1.0 - 0.01 + 0.001 0.03 @@ -62,9 +62,17 @@ 1 1 + + 2 + 1 + + + 3 + 1 + 1 - gen_watch.csv + gen_watch_SEXS.csv -ksp_type richardson diff --git a/src/applications/dynamic_simulation_full_y/CMakeLists.txt b/src/applications/dynamic_simulation_full_y/CMakeLists.txt index 91cf30041..abdc99a8b 100644 --- a/src/applications/dynamic_simulation_full_y/CMakeLists.txt +++ b/src/applications/dynamic_simulation_full_y/CMakeLists.txt @@ -88,26 +88,11 @@ endif() target_link_libraries(dsf.x ${target_libraries}) target_link_libraries(dsf2.x ${target_libraries}) -gridpack_set_lu_solver( - "${GRIDPACK_DATA_DIR}/input/ds/input_145.xml" - "${CMAKE_CURRENT_BINARY_DIR}/input_145.xml" -) - gridpack_set_lu_solver( "${GRIDPACK_DATA_DIR}/input/ds/input_9b3g.xml" "${CMAKE_CURRENT_BINARY_DIR}/input_9b3g.xml" ) -gridpack_set_lu_solver( - "${GRIDPACK_DATA_DIR}/input/ds/input_300_cmpld.xml" - "${CMAKE_CURRENT_BINARY_DIR}/input_300_cmpld.xml" -) - -gridpack_set_lu_solver( - "${GRIDPACK_DATA_DIR}/input/ds/input_3000.xml" - "${CMAKE_CURRENT_BINARY_DIR}/input_3000.xml" -) - gridpack_set_lu_solver( "${GRIDPACK_DATA_DIR}/input/ds/input_240bus.xml" "${CMAKE_CURRENT_BINARY_DIR}/input_240bus.xml" @@ -118,223 +103,42 @@ gridpack_set_lu_solver( "${CMAKE_CURRENT_BINARY_DIR}/input_twoarea.xml" ) -gridpack_set_lu_solver( - "${GRIDPACK_DATA_DIR}/input/ds/input_twoarea_renewable_mech.xml" - "${CMAKE_CURRENT_BINARY_DIR}/input_twoarea_renewable_mech.xml" -) - -gridpack_set_lu_solver( - "${GRIDPACK_DATA_DIR}/input/ds/input_unit_GAST.xml" - "${CMAKE_CURRENT_BINARY_DIR}/input_unit_GAST.xml" -) - -gridpack_set_lu_solver( - "${GRIDPACK_DATA_DIR}/input/ds/input_unit_HYGOV.xml" - "${CMAKE_CURRENT_BINARY_DIR}/input_unit_HYGOV.xml" -) - -gridpack_set_lu_solver( - "${GRIDPACK_DATA_DIR}/input/ds/input_unit_REGCA1_REECA1_REPCA1.xml" - "${CMAKE_CURRENT_BINARY_DIR}/input_unit_REGCA1_REECA1_REPCA1.xml" -) - -gridpack_set_lu_solver( - "${GRIDPACK_DATA_DIR}/input/ds/input_unit_SEXS.xml" - "${CMAKE_CURRENT_BINARY_DIR}/input_unit_SEXS.xml" -) - -gridpack_set_lu_solver( - "${GRIDPACK_DATA_DIR}/input/ds/input_unit_GENROU.xml" - "${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENROU.xml" -) - -gridpack_set_lu_solver( - "${GRIDPACK_DATA_DIR}/input/ds/input_unit_GENROU_ESST1A.xml" - "${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENROU_ESST1A.xml" -) - -gridpack_set_lu_solver( - "${GRIDPACK_DATA_DIR}/input/ds/input_unit_GENROU_ESST1A_WSIEG1.xml" - "${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENROU_ESST1A_WSIEG1.xml" -) - -gridpack_set_lu_solver( - "${GRIDPACK_DATA_DIR}/input/ds/input_unit_GENROU_EXDC1.xml" - "${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENROU_EXDC1.xml" -) - -gridpack_set_lu_solver( - "${GRIDPACK_DATA_DIR}/input/ds/input_unit_GENROU_EXDC1_WSIEG1.xml" - "${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENROU_EXDC1_WSIEG1.xml" -) - -gridpack_set_lu_solver( - "${GRIDPACK_DATA_DIR}/input/ds/input_unit_GENSAL.xml" - "${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENSAL.xml" -) - -gridpack_set_lu_solver( - "${GRIDPACK_DATA_DIR}/input/ds/input_unit_GENSAL_ESST1A.xml" - "${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENSAL_ESST1A.xml" -) - -gridpack_set_lu_solver( - "${GRIDPACK_DATA_DIR}/input/ds/input_unit_GENSAL_ESST1A_WSIEG1.xml" - "${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENSAL_ESST1A_WSIEG1.xml" -) - add_custom_target(dsf.x.input - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/raw/IEEE_145bus_v23_PSLF.raw - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/dyr/IEEE_145b_classical_model.dyr - ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E copy + COMMAND ${CMAKE_COMMAND} -E copy ${GRIDPACK_DATA_DIR}/raw/9b3g.raw ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E copy + COMMAND ${CMAKE_COMMAND} -E copy ${GRIDPACK_DATA_DIR}/dyr/9b3g.dyr ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/raw/300bus_v23_no0imp_pslf.raw - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/dyr/300bus_detail_model_cmpld_combine.dyr - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/raw/bus3000_gen_no0imp_v23_pslf.raw - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/dyr/classical_model_3000bus.dyr - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy + COMMAND ${CMAKE_COMMAND} -E copy ${GRIDPACK_DATA_DIR}/dyr/240busWECC_2018_PSS_mod.dyr ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E copy + COMMAND ${CMAKE_COMMAND} -E copy ${GRIDPACK_DATA_DIR}/raw/240busWECC_2018_PSS_fixedshunt.raw - ${CMAKE_CURRENT_BINARY_DIR} + ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E copy + COMMAND ${CMAKE_COMMAND} -E copy ${GRIDPACK_DATA_DIR}/dyr/kundur-twoarea.dyr ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E copy + COMMAND ${CMAKE_COMMAND} -E copy ${GRIDPACK_DATA_DIR}/raw/kundur-twoarea_v33.raw - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/dyr/kundur-twoarea_4renewable_mech.dyr - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/raw/case9.raw - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/dyr/case9_GAST.dyr - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/dyr/case9_HYGOV.dyr - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/dyr/case9_REGCA1_REECA1_REPCA1.dyr ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/dyr/case9_SEXS.dyr - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/dyr/case9_GENROU.dyr - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/dyr/case9_GENROU_ESST1A.dyr - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/dyr/case9_GENROU_ESST1A_WSIEG1.dyr - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/dyr/case9_GENROU_EXDC1.dyr - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/dyr/case9_GENROU_EXDC1_WSIEG1.dyr - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/dyr/case9_GENSAL.dyr - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/dyr/case9_GENSAL_ESST1A.dyr - ${CMAKE_CURRENT_BINARY_DIR} - - COMMAND ${CMAKE_COMMAND} -E copy - ${GRIDPACK_DATA_DIR}/dyr/case9_GENSAL_ESST1A_WSIEG1.dyr - ${CMAKE_CURRENT_BINARY_DIR} - - DEPENDS - ${CMAKE_CURRENT_BINARY_DIR}/input_145.xml - ${GRIDPACK_DATA_DIR}/raw/IEEE_145bus_v23_PSLF.raw - ${GRIDPACK_DATA_DIR}/dyr/IEEE_145b_classical_model.dyr + DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/input_9b3g.xml ${GRIDPACK_DATA_DIR}/raw/9b3g.raw ${GRIDPACK_DATA_DIR}/dyr/9b3g.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_300_cmpld.xml - ${GRIDPACK_DATA_DIR}/raw/300bus_v23_no0imp_pslf.raw - ${GRIDPACK_DATA_DIR}/dyr/300bus_detail_model_cmpld_combine.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_3000.xml - ${GRIDPACK_DATA_DIR}/raw/bus3000_gen_no0imp_v23_pslf.raw - ${GRIDPACK_DATA_DIR}/dyr/classical_model_3000bus.dyr ${CMAKE_CURRENT_BINARY_DIR}/input_240bus.xml ${GRIDPACK_DATA_DIR}/dyr/240busWECC_2018_PSS_mod.dyr ${GRIDPACK_DATA_DIR}/raw/240busWECC_2018_PSS_fixedshunt.raw ${CMAKE_CURRENT_BINARY_DIR}/input_twoarea.xml ${GRIDPACK_DATA_DIR}/dyr/kundur-twoarea.dyr ${GRIDPACK_DATA_DIR}/raw/kundur-twoarea_v33.raw - ${CMAKE_CURRENT_BINARY_DIR}/input_twoarea_renewable_mech.xml - ${GRIDPACK_DATA_DIR}/dyr/kundur-twoarea_4renewable_mech.dyr - ${GRIDPACK_DATA_DIR}/raw/case9.raw - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GAST.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GAST.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENROU.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GENROU.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENROU_ESST1A.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GENROU_ESST1A.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENROU_ESST1A_WSIEG1.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GENROU_ESST1A_WSIEG1.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENROU_EXDC1.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GENROU_EXDC1.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENROU_EXDC1_WSIEG1.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GENROU_EXDC1_WSIEG1.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENSAL.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GENSAL.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENSAL_ESST1A.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GENSAL_ESST1A.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENSAL_ESST1A_WSIEG1.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GENSAL_ESST1A_WSIEG1.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_HYGOV.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_HYGOV.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_REGCA1_REECA1_REPCA1.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_REGCA1_REECA1_REPCA1.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_SEXS.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_SEXS.dyr ) add_dependencies(dsf.x dsf.x.input) @@ -349,47 +153,19 @@ endif() configure_file(${CMAKE_CURRENT_SOURCE_DIR}/CMakeLists.install.in ${CMAKE_CURRENT_BINARY_DIR}/CMakeLists.txt @ONLY) -install(FILES +install(FILES ${CMAKE_CURRENT_BINARY_DIR}/CMakeLists.txt - ${CMAKE_CURRENT_BINARY_DIR}/input_145.xml - ${GRIDPACK_DATA_DIR}/raw/IEEE_145bus_v23_PSLF.raw - ${GRIDPACK_DATA_DIR}/dyr/IEEE_145b_classical_model.dyr ${CMAKE_CURRENT_BINARY_DIR}/input_9b3g.xml ${GRIDPACK_DATA_DIR}/raw/9b3g.raw ${GRIDPACK_DATA_DIR}/dyr/9b3g.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_300_cmpld.xml - ${GRIDPACK_DATA_DIR}/raw/300bus_v23_no0imp_pslf.raw - ${GRIDPACK_DATA_DIR}/dyr/300bus_detail_model_cmpld_combine.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_3000.xml - ${GRIDPACK_DATA_DIR}/raw/bus3000_gen_no0imp_v23_pslf.raw - ${GRIDPACK_DATA_DIR}/dyr/classical_model_3000bus.dyr + ${CMAKE_CURRENT_BINARY_DIR}/input_240bus.xml + ${GRIDPACK_DATA_DIR}/dyr/240busWECC_2018_PSS_mod.dyr + ${GRIDPACK_DATA_DIR}/raw/240busWECC_2018_PSS_fixedshunt.raw + ${CMAKE_CURRENT_BINARY_DIR}/input_twoarea.xml + ${GRIDPACK_DATA_DIR}/dyr/kundur-twoarea.dyr + ${GRIDPACK_DATA_DIR}/raw/kundur-twoarea_v33.raw ${CMAKE_CURRENT_SOURCE_DIR}/dsf_main.cpp ${CMAKE_CURRENT_SOURCE_DIR}/dsf_main2.cpp - ${GRIDPACK_DATA_DIR}/raw/case9.raw - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GAST.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GAST.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENROU.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GENROU.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENROU_ESST1A.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GENROU_ESST1A.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENROU_ESST1A_WSIEG1.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GENROU_ESST1A_WSIEG1.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENROU_EXDC1.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GENROU_EXDC1.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENROU_EXDC1_WSIEG1.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GENROU_EXDC1_WSIEG1.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENSAL.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GENSAL.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENSAL_ESST1A.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GENSAL_ESST1A.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_GENSAL_ESST1A_WSIEG1.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_GENSAL_ESST1A_WSIEG1.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_HYGOV.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_HYGOV.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_REGCA1_REECA1_REPCA1.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_REGCA1_REECA1_REPCA1.dyr - ${CMAKE_CURRENT_BINARY_DIR}/input_unit_SEXS.xml - ${GRIDPACK_DATA_DIR}/dyr/case9_SEXS.dyr DESTINATION share/gridpack/example/dynamic_simulation_full_y ) @@ -399,24 +175,7 @@ install(TARGETS dsf2.x DESTINATION bin) # ------------------------------------------------------------- # run application as test # ------------------------------------------------------------- -gridpack_add_run_test("dynamic_simulation_full_y_145_bus" dsf2.x input_145.xml) -gridpack_add_run_test("dynamic_simulation_full_y_2_145_bus" dsf2.x input_145.xml) +gridpack_add_run_test("dynamic_simulation_full_y_9b3g" dsf2.x input_9b3g.xml) gridpack_add_run_test("dynamic_simulation_full_y_240_bus" dsf2.x input_240bus.xml) -gridpack_add_run_test("dynamic_simulation_full_y2_240_bus" dsf2.x input_240bus.xml) -gridpack_add_run_test("dynamic_simulation_two_area" dsf2.x input_twoarea.xml) -gridpack_add_run_test("dynamic_simulation_2_two_area" dsf2.x input_twoarea.xml) -gridpack_add_run_test("dynamic_simulation_two_area_renewable" dsf2.x input_twoarea_renewable_mech.xml) -gridpack_add_run_test("dynamic_simulation_2_two_area_renewable" dsf2.x input_twoarea_renewable_mech.xml) -gridpack_add_run_test("dynamic_simulation_unit_gast" dsf2.x input_unit_GAST.xml) -gridpack_add_run_test("dynamic_simulation_unit_genrou" dsf2.x input_unit_GENROU.xml) -gridpack_add_run_test("dynamic_simulation_unit_genrou_esst1a" dsf2.x input_unit_GENROU_ESST1A.xml) -gridpack_add_run_test("dynamic_simulation_unit_genrou_esst1a_wsieg1" dsf2.x input_unit_GENROU_ESST1A_WSIEG1.xml) -gridpack_add_run_test("dynamic_simulation_unit_genrou_exdc1" dsf2.x input_unit_GENROU_EXDC1.xml) -gridpack_add_run_test("dynamic_simulation_unit_genrou_exdc1_wsieg1" dsf2.x input_unit_GENROU_EXDC1_WSIEG1.xml) -gridpack_add_run_test("dynamic_simulation_unit_hygov" dsf2.x input_unit_HYGOV.xml) -gridpack_add_run_test("dynamic_simulation_unit_regca1_reeca1_repca1" dsf2.x input_unit_REGCA1_REECA1_REPCA1.xml) -gridpack_add_run_test("dynamic_simulation_unit_sexs" dsf2.x input_unit_SEXS.xml) -if (ENABLE_ENVIRONMENT_FROM_COMM) - gridpack_add_run_test("dynamic_simulation_comm_full_y" dsf_comm.x input_145.xml) -endif() +gridpack_add_run_test("dynamic_simulation_full_y_two_area" dsf2.x input_twoarea.xml) diff --git a/src/applications/dynamic_simulation_full_y/dsf_main2.cpp b/src/applications/dynamic_simulation_full_y/dsf_main2.cpp index 1089ad848..7ca9f35f3 100644 --- a/src/applications/dynamic_simulation_full_y/dsf_main2.cpp +++ b/src/applications/dynamic_simulation_full_y/dsf_main2.cpp @@ -38,11 +38,13 @@ void run_dynamics(int argc, char **argv) ds_app.setup(); int ngen = ds_app.numGenerators(); int nload = ds_app.numLoads(); - int nline = ds_app.numLines(); + // Note: numLines() crashes with CMPLD cases where extended + // buses/branches have invalid DataCollection pointers + //int nline = ds_app.numLines(); if (world.rank() == 0) { std::cout << " Number of generators: "< p_governor; // governor boost::shared_ptr p_exciter; // exciter diff --git a/src/applications/modules/dynamic_simulation_full_y/dsf_app_module.cpp b/src/applications/modules/dynamic_simulation_full_y/dsf_app_module.cpp index b490b7bcb..08550d021 100644 --- a/src/applications/modules/dynamic_simulation_full_y/dsf_app_module.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/dsf_app_module.cpp @@ -290,7 +290,6 @@ void gridpack::dynamic_simulation::DSFullApp::setNetwork( p_biterative_solve_network = cursor->get("iterativeNetworkInterface",false); p_iterative_network_debug = cursor->get("iterativeNetworkInterfaceDebugPrint",false); p_generator_observationpower_systembase = cursor->get("generatorObservationPowerSystemBase",true); - ITER_TOL = cursor->get("iterativeNetworkInterfaceTol", 1.0e-7); MAX_ITR_NO = cursor->get("iterativeNetworkInterfaceMaxItrNo", 8); p_equilibrium_init = cursor->get("equilibriumInit", false); @@ -3136,7 +3135,7 @@ void gridpack::dynamic_simulation::DSFullApp::solvePreInitialize( // Simulation related variables t_init = timer->createCategory("DS Solve: Initialization"); timer->start(t_init); - + int t_step[20]; double t_width[20]; diff --git a/src/applications/modules/dynamic_simulation_full_y/dsf_app_module.hpp b/src/applications/modules/dynamic_simulation_full_y/dsf_app_module.hpp index 73e6fdda6..79ba41fbc 100644 --- a/src/applications/modules/dynamic_simulation_full_y/dsf_app_module.hpp +++ b/src/applications/modules/dynamic_simulation_full_y/dsf_app_module.hpp @@ -989,7 +989,6 @@ class DSFullApp //for the generator observations, output the generator power based on system base or generator base bool p_generator_observationpower_systembase; - // Current step count? int p_S_Steps; diff --git a/src/applications/modules/dynamic_simulation_full_y/dsf_app_module2.cpp b/src/applications/modules/dynamic_simulation_full_y/dsf_app_module2.cpp index a4ed333ad..7acaaf388 100644 --- a/src/applications/modules/dynamic_simulation_full_y/dsf_app_module2.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/dsf_app_module2.cpp @@ -1,3 +1,20 @@ +/* + * Copyright (c) 2013 Battelle Memorial Institute + * Licensed under modified BSD License. A copy of this license can be found + * in the LICENSE file in the top level directory of this distribution. + */ +// ------------------------------------------------------------- +/** + * @file dsf_app_module2.cpp + * @author Yousu Chen + * + * @brief Dynamic simulation application module — setup, time + * stepping, and network solve routines. Implements the + * modified Euler predictor-corrector loop (runonestep), + * Y-bus construction, event handling, and generator/load + * watch file output. + */ + #include "dsf_app_module.hpp" #include #include @@ -77,7 +94,7 @@ bool gridpack::dynamic_simulation::DSFullApp::solveNetwork(int predcorrflag) { bool converged = false; - // volt_full->zero(); + //volt_full->zero(); int its=0; //bool p_iterative_network_debug = false; @@ -164,10 +181,10 @@ void gridpack::dynamic_simulation::DSFullApp::setup() p_events = getEvents(); p_factory->setMode(YBUS); - + ybusMap_sptr.reset(new gridpack::mapper::FullMatrixMap (p_network)); orgYbus = ybusMap_sptr->mapToMatrix(); - + // Form constant impedance load admittance yl for all buses and add it to // system Y matrix: ybus = ybus + yl p_factory->setMode(YL); @@ -180,19 +197,19 @@ void gridpack::dynamic_simulation::DSFullApp::setup() // Extract appropriate xdprime and xdpprime from machine data p_factory->setMode(jxd); ybus_jxd = ybusMap_sptr->mapToMatrix(); - + // Add dynamic load impedance to system Y matrix: p_factory->setMode(YDYNLOAD); ybus = ybusMap_sptr->mapToMatrix(); - // Initialize vectors for integration + // Initialize vectors for integration p_factory->initDSVect(p_time_step); - + p_factory->setGeneratorObPowerBaseFlag(p_generator_observationpower_systembase); /* Create mapper and vectors */ ngenMap_sptr.reset(new gridpack::mapper::BusVectorMap (p_network)); - + p_insecureAt = -1; p_factory->setMode(make_INorton_full); @@ -325,7 +342,6 @@ void gridpack::dynamic_simulation::DSFullApp::setup() if (p_loadWatch) p_loadIO->header("\n"); } - p_frequencyOK = true; } @@ -335,7 +351,7 @@ void gridpack::dynamic_simulation::DSFullApp::setup() void gridpack::dynamic_simulation::DSFullApp::runonestep() { bool converged; - + S_Steps = Simu_Current_Step; /* Predictor current injection */ @@ -343,7 +359,7 @@ void gridpack::dynamic_simulation::DSFullApp::runonestep() /* Solve Network equations */ converged = solveNetwork(0); - + if ( Simu_Current_Step==0 ) { //printf("enter the initial update oldbusvoltage, Timestep: %d \n", Simu_Current_Step); p_factory->updateoldbusvoltage(); //renke add, first timestep, copy volt_full to volt_full_old @@ -354,7 +370,7 @@ void gridpack::dynamic_simulation::DSFullApp::runonestep() /* yuan add: update branch power*/ p_factory->updateData(); - + std::vector vwideareafreqs; vwideareafreqs = p_factory->grabWideAreaFreq(); int tmp = vwideareafreqs.size(); @@ -362,7 +378,7 @@ void gridpack::dynamic_simulation::DSFullApp::runonestep() bool flagBus = p_factory->updateBusRelay(false, p_time_step); bool flagBranch = p_factory->updateBranchRelay(false, p_time_step); - + // update dynamic load internal relay functions here p_factory->dynamicload_post_process(p_time_step, false); @@ -378,13 +394,13 @@ void gridpack::dynamic_simulation::DSFullApp::runonestep() ybusMap_sptr->incrementMatrix(ybus); } - // Update old voltage (??) + // Update old voltage p_factory->updateoldbusvoltage(); /* Predictor */ if (Simu_Current_Step !=0 && last_S_Steps != S_Steps) { p_factory->predictor(p_time_step, false); - } else { + } else { p_factory->predictor(p_time_step, true); } diff --git a/src/applications/modules/dynamic_simulation_full_y/dsf_components.cpp b/src/applications/modules/dynamic_simulation_full_y/dsf_components.cpp index 1801cc177..6abddbc6d 100644 --- a/src/applications/modules/dynamic_simulation_full_y/dsf_components.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/dsf_components.cpp @@ -326,7 +326,7 @@ bool gridpack::dynamic_simulation::DSFullBus::matrixDiagValues(ComplexType *valu double tmp1 = p_ybusr; double tmp2 = p_ybusi - 1.0e9; gridpack::ComplexType ret(tmp1, tmp2); - + values[0] = ret; return true; } else { @@ -468,16 +468,13 @@ bool gridpack::dynamic_simulation::DSFullBus::vectorValues(ComplexType *values) if (p_ngen > 0) { for (int i = 0; i < p_ngen; i++) { if(!p_gstatus[i] || !p_generators.size()) continue; - values[0] += p_generators[i]->INorton(); + values[0] += p_generators[i]->INorton(); } // generator for loop } // if p_ngen>0 //INorton contribution from dynamic load for (int i =0; iINorton(); - - ComplexType tmp = p_loadmodels[i]->INorton(); - //printf("DSFullBus::vectorValues, bus %d, dynamic load Inorton: %12.6f + j %12.6f \n", getOriginalIndex(),real(tmp), imag(tmp)); + values[0] += p_loadmodels[i]->INorton(); } // Equilibrium load compensation: inject current to make constant-Z @@ -3397,11 +3394,11 @@ bool gridpack::dynamic_simulation::DSFullBranch::matrixForwardValues(ComplexType printf("matrix off diag forward element changes due to branch relay trip!\n"); values[0] = -getBranchRelayTripUpdateFactor(); printf("changed value: %f + j*%f\n", real(values[0]), imag(values[0])); - + return true; } else { return false; - } + } }else if (p_mode == branch_trip_action) { if (p_branchactiontripflag) { //printf("matrix off diag forward element changes due to branch relay trip!\n"); @@ -3461,7 +3458,7 @@ bool gridpack::dynamic_simulation::DSFullBranch::matrixReverseValues(ComplexType return true; } else { return false; - } + } }else if (p_mode == branch_trip_action) { if (p_branchactiontripflag) { //printf("matrix off diag forward element changes due to branch relay trip!\n"); @@ -3660,6 +3657,10 @@ void gridpack::dynamic_simulation::DSFullBranch::load( */ void gridpack::dynamic_simulation::DSFullBranch::evaluateBranchFlow() { + // Skip extended load branches (CMPLD transformer/feeder) which have + // no element data and may have invalid bus references + if (p_bextendedloadbranch > 0) return; + int i; double pi = 4.0*atan(1.0); @@ -3764,6 +3765,9 @@ void gridpack::dynamic_simulation::DSFullBranch::evaluateBranchFlow() void gridpack::dynamic_simulation::DSFullBranch::updateData( boost::shared_ptr &data) { + // Skip extended load branches (CMPLD transformer/feeder) + if (p_bextendedloadbranch > 0) return; + int i; evaluateBranchFlow(); updateBranchCurrent(); diff --git a/src/applications/modules/dynamic_simulation_full_y/dsf_events.cpp b/src/applications/modules/dynamic_simulation_full_y/dsf_events.cpp index 71842a99c..b356e9942 100644 --- a/src/applications/modules/dynamic_simulation_full_y/dsf_events.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/dsf_events.cpp @@ -3,6 +3,13 @@ * Licensed under modified BSD License. A copy of this license can be found * in the LICENSE file in the top level directory of this distribution. */ +// ------------------------------------------------------------- +/** + * @updated Yousu Chen + * - Rebuild LinearSolver after Y-bus changes (fault on/off, line + * status, generator status) to avoid stale preconditioner. + * @date 2026-03-10 + */ #include "dsf_app_module.hpp" #include @@ -134,6 +141,7 @@ void gridpack::dynamic_simulation::DSFullApp::handleEvents() { int nevents = p_events.size(); int i; + bool ybus_changed = false; for(i = 0; i < nevents; i++) { gridpack::dynamic_simulation::Event event = p_events[i]; @@ -152,7 +160,8 @@ void gridpack::dynamic_simulation::DSFullApp::handleEvents() // Update Ybus p_factory->setMode(BUSFAULTON); ybusMap_sptr->incrementMatrix(ybus); - + ybus_changed = true; + } else if(fabs(event.end - p_current_time) < 1e-6) { /* Fault end */ std::vector bus_internal_idx; @@ -161,16 +170,18 @@ void gridpack::dynamic_simulation::DSFullApp::handleEvents() if(bus_internal_idx.size()) { bus = dynamic_cast(p_network->getBus(bus_internal_idx[0]).get()); } - + // Update Ybus p_factory->setMode(BUSFAULTOFF); ybusMap_sptr->incrementMatrix(ybus); + ybus_changed = true; } } else if(event.isLineStatus) { if(fabs(event.time - p_current_time) < 1e-6) { setLineStatus(event.from_idx,event.to_idx,event.tag,event.status); p_factory->setMode(LINESTATUSCHANGE); ybusMap_sptr->incrementMatrix(ybus); + ybus_changed = true; } } else if(event.isGenStatus) { if(fabs(event.time - p_current_time) < 1e-6) { @@ -178,8 +189,17 @@ void gridpack::dynamic_simulation::DSFullApp::handleEvents() setGenStatus(event.bus_idx,event.tag,event.status); p_factory->setMode(GENSTATUSCHANGE); ybusMap_sptr->incrementMatrix(ybus); + ybus_changed = true; } } } + + // Rebuild the linear solver if the Y-bus matrix was modified + if (ybus_changed) { + gridpack::utility::Configuration::CursorPtr cursor; + cursor = p_config->getCursor("Configuration.Dynamic_simulation"); + solver_sptr.reset(new gridpack::math::LinearSolver(*ybus)); + solver_sptr->configure(cursor); + } } diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/acmotor.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/acmotor.cpp index 11eba07dc..93f6a7173 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/acmotor.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/acmotor.cpp @@ -6,10 +6,15 @@ // ----------------------------------------------------------- /** * @file acmotor.cpp - * @author Shuangshuang Jin - * @Last modified: Oct 21, 2016 - * - * @brief + * @author Shuangshuang Jin + * @Last modified: Oct 21, 2016 + * + * @updated Yousu Chen + * - Fixed motor B thermal calculation: ImotorB_pu was incorrectly + * assigned to ImotorA_pu in predictor and corrector. + * @date 2026-03-10 + * + * @brief * * */ @@ -561,9 +566,9 @@ void gridpack::dynamic_simulation::AcmotorLoad::predictor( ImotorB_pu = vt * equivY ; else { gridpack::ComplexType tmp(PB, -QB); - ImotorA_pu = tmp / vt; + ImotorB_pu = tmp / vt; } - + dThB_dt0 = (abs(ImotorB_pu) * abs(ImotorB_pu) * Rstall - temperatureB0) / Tth; // step-3 integration @@ -622,10 +627,10 @@ void gridpack::dynamic_simulation::AcmotorLoad::corrector( ImotorB_pu = vt * equivY ; else { gridpack::ComplexType tmp(PB, -QB); - ImotorA_pu = tmp / vt; + ImotorB_pu = tmp / vt; } - - dThB_dt = (abs(ImotorB_pu) * abs(ImotorB_pu) * Rstall -temperatureB) / Tth; + + dThB_dt = (abs(ImotorB_pu) * abs(ImotorB_pu) * Rstall - temperatureB) / Tth; // step-2 integration volt_measured = volt_measured0 + 0.5 * (dv_dt0 + dv_dt) * dt; diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/classical.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/classical.cpp index d8821980e..3a3b536d2 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/classical.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/classical.cpp @@ -73,6 +73,18 @@ void gridpack::dynamic_simulation::ClassicalGenerator::load( if (!data->getValue(GENERATOR_ZSOURCE, &zsrc, idx)) zsrc=gridpack::ComplexType(0.0,0.0); // dtr p_dtr = imag(zsrc); + + // Read Xdp from DYR file (GENERATOR_TRANSIENT_REACTANCE) + // This is the correct reactance for GENCLS, not ZSOURCE from RAW file + if (!data->getValue(GENERATOR_TRANSIENT_REACTANCE, &p_Xdp, idx)) { + // Fall back to ZSOURCE if DYR Xdp not available + p_Xdp = p_dtr; + } + + // Xdp bounds + if (p_Xdp < 0.00001) p_Xdp = 0.00001; + if (p_Xdp > 999.0) p_Xdp = 999.0; + if (!data->getValue(GENERATOR_INERTIA_CONSTANT_H, &p_h, idx)) p_h = 0.0; // h if (!data->getValue(GENERATOR_DAMPING_COEFFICIENT_0, &p_d0, idx)) p_d0 = 0.0; // d0 @@ -135,7 +147,7 @@ void gridpack::dynamic_simulation::ClassicalGenerator::init(double mag, gridpack::ComplexType jay(0.0, 1.0); //printf("v: (%f, %f)\n", real(v), imag(v)); //gridpack::ComplexType temp = v + jay * (p_dtr * p_mva) * curr; - gridpack::ComplexType temp = v + jay * p_dtr * curr; + gridpack::ComplexType temp = v + jay * p_Xdp * curr; //printf("dtr = %f, mva = %f\n", p_dtr, p_mva); //printf("eprime_s0: (%f, %f)\n", real(temp), imag(temp)); p_eprime_s0 = temp; @@ -194,12 +206,8 @@ gridpack::ComplexType gridpack::dynamic_simulation::ClassicalGenerator::INorton( */ gridpack::ComplexType gridpack::dynamic_simulation::ClassicalGenerator::NortonImpedence() { - double ra = p_r * p_sbase / p_mva; // * p_sbase / p_mva; - double xd; - if (p_dstr == 0.0) { - xd = p_dtr * p_sbase / p_mva; // * p_sbase / p_mva; - } - //printf("---------classical generator model :: NortonImpedence() p_r = %f, ra = %f, xd = %f, p_dstr = %f, p_dtr = %f \n", p_r, ra, xd, p_dstr, p_dtr); + double ra = p_r * p_sbase / p_mva; + double xd = p_Xdp * p_sbase / p_mva; gridpack::ComplexType Y_a(ra, xd); Y_a = 1.0 / Y_a; @@ -220,8 +228,8 @@ void gridpack::dynamic_simulation::ClassicalGenerator::predictor_currentInjectio gridpack::ComplexType jay(0.0, 1.0); // Calculate INorton_full - p_INorton = p_eprime_s0 / (p_dtr * jay); - + p_INorton = p_eprime_s0 / (p_Xdp * jay); + p_INorton = p_INorton * p_mva / p_sbase; /*double Idnorton = real(p_INorton); @@ -254,10 +262,10 @@ void gridpack::dynamic_simulation::ClassicalGenerator::predictor( jay = gridpack::ComplexType(0.0,1.0); //p_INorton = gridpack::ComplexType(0.0,0.0); - // terminal curr: curr = (eprime - volt) / GEN_dtr) ----------- + // terminal curr: curr = (eprime - volt) / (j*Xdp) ----------- curr = p_eprime_s0; curr -= p_volt; - curr = curr/(jay*p_dtr); + curr = curr/(jay*p_Xdp); //imag(curr) = -imag(curr); curr = conj(curr); p_pelect = real(p_eprime_s0 * curr); @@ -287,8 +295,8 @@ void gridpack::dynamic_simulation::ClassicalGenerator::predictor( void gridpack::dynamic_simulation::ClassicalGenerator::corrector_currentInjection(bool flag) { gridpack::ComplexType jay(0.0, 1.0); - p_INorton = p_eprime_s1 / (p_dtr * jay); - + p_INorton = p_eprime_s1 / (p_Xdp * jay); + p_INorton = p_INorton * p_mva / p_sbase; /*double Idnorton = real(p_INorton); @@ -314,10 +322,10 @@ void gridpack::dynamic_simulation::ClassicalGenerator::corrector( // Evaluate updated values of machine parameters for integration jay = gridpack::ComplexType(0.0,1.0); p_INorton = gridpack::ComplexType(0.0,0.0); - // terminal curr: curr = (eprime - volt) / GEN_dtr) ----------- + // terminal curr: curr = (eprime - volt) / (j*Xdp) ----------- curr = p_eprime_s1; curr -= p_volt; - curr = curr/(jay*p_dtr); + curr = curr/(jay*p_Xdp); //imag(curr) = -imag(curr); curr = conj(curr); p_pelect = real(p_eprime_s1 * curr); @@ -387,15 +395,18 @@ serialWrite(char *string, const int bufsize, const char *signal) bool ret = false; if (!strcmp(signal,"watch_header")) { if (getWatch()) { - char buf[128]; + char buf[256]; std::string tag; if (p_ckt[1] != ' ') { tag = p_ckt; } else { tag = p_ckt[0]; } - sprintf(buf,", %d_%s_angle, %d_%s_speed",p_bus_id,tag.c_str(), - p_bus_id,tag.c_str()); + sprintf(buf,", %d_%s_V, %d_%s_Pg, %d_%s_Qg, %d_%s_angle, %d_%s_speed," + " %d_%s_Efd, %d_%s_Pm, %d_%s_PowerAngle", + p_bus_id,tag.c_str(),p_bus_id,tag.c_str(),p_bus_id,tag.c_str(), + p_bus_id,tag.c_str(),p_bus_id,tag.c_str(),p_bus_id,tag.c_str(), + p_bus_id,tag.c_str(),p_bus_id,tag.c_str()); if (strlen(buf) <= bufsize) { sprintf(string,"%s",buf); ret = true; @@ -410,7 +421,16 @@ serialWrite(char *string, const int bufsize, const char *signal) } else if (!strcmp(signal,"watch")) { if (getWatch()) { char buf[256]; - sprintf(buf,", %f, %f",real(p_mac_ang_s1),real(p_mac_spd_s1)); + double Vterm = abs(p_volt); + double VtermAng = atan2(imag(p_volt), real(p_volt)); + double rotor_ang = real(p_mac_ang_s1); + double powerAngle = fmod((rotor_ang - VtermAng) * 180.0 / 3.14159265358979323846 + 180.0, 360.0) - 180.0; + if (powerAngle < -180.0) powerAngle += 360.0; + double Efd = real(p_eqprime); // constant E' magnitude for classical model + double Pmech_out = real(p_pmech) * p_mva / p_sbase; + sprintf(buf,",%f,%f,%f,%f, %f, %f, %f, %f", + Vterm, genP*p_mva/p_sbase, genQ*p_mva/p_sbase, + rotor_ang, real(p_mac_spd_s1), Efd, Pmech_out, powerAngle); if (strlen(buf) <= bufsize) { sprintf(string,"%s",buf); ret = true; diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/classical.hpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/classical.hpp index a798bcbb7..c8c501cf6 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/classical.hpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/classical.hpp @@ -155,6 +155,7 @@ class ClassicalGenerator : public BaseGeneratorModel double p_pg, p_qg; int p_status; double p_mva, p_r, p_dstr, p_dtr; + double p_Xdp; // Transient reactance from DYR file (used for all calculations) double p_d0, p_h; double p_PI; double genP, genQ; @@ -183,6 +184,7 @@ class ClassicalGenerator : public BaseGeneratorModel & p_status & p_mva & p_r & p_dstr & p_dtr + & p_Xdp & p_d0 & p_h & p_pelect & p_volt & p_mac_ang_s0 & p_mac_spd_s0 diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/esst1a.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/esst1a.cpp index c7e452a3f..873603426 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/esst1a.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/esst1a.cpp @@ -7,6 +7,13 @@ /** * @file esst1a.cpp * + * @updated Yousu Chen + * - Fixed regulator input + * - Fixed computeModel: bypass lead-lag blocks when TB or TB1 is zero, + * matching init() logic. Uninitialized blocks produced zero output. + * - Fixed regulator limits + * @date 2026-03-11 + * * @brief ESST1A model * * @@ -15,6 +22,7 @@ #include #include #include +#include #include "boost/smart_ptr/shared_ptr.hpp" #include "gridpack/parser/dictionary.hpp" @@ -83,6 +91,17 @@ void gridpack::dynamic_simulation::Esst1aModel::load( Vuel = 0.0; Voel = 1000.0; + // Swap limits if inverted + if (Vrmax < Vrmin) { + double tmp = Vrmax; Vrmax = Vrmin; Vrmin = tmp; + } + if (Vamax < Vamin) { + double tmp = Vamax; Vamax = Vamin; Vamin = tmp; + } + if (Vimax < Vimin) { + double tmp = Vimax; Vimax = Vimin; Vimin = tmp; + } + } /** @@ -117,12 +136,59 @@ double gridpack::dynamic_simulation::Esst1aModel::sqr(double x) */ void gridpack::dynamic_simulation::Esst1aModel::init(double mag, double ang, double ts) { - if (Tf < TS_THRESHOLD * ts) zero_TF = true; - if (Tb < TS_THRESHOLD * ts) zero_TB = true; - if (Tb1 < TS_THRESHOLD * ts) zero_TB1 = true; - if (Ta < TS_THRESHOLD * ts) zero_TA = true; - if (Tr < TS_THRESHOLD * ts) zero_TR = true; - + // Time constant minimum checks + double mult_ts = TS_THRESHOLD * ts; + + // Tf uses 0.5/1.0 pattern (lead-lag denominator type) + if (Tf > 0.0 && Tf < 0.5 * mult_ts) { + Tf = 0.0; + zero_TF = true; + } else if (Tf > 0.5 * mult_ts && Tf < mult_ts) { + Tf = mult_ts; + } else if (Tf <= 0.0) { + zero_TF = true; + } + + // Tb uses 0.5/1.0 pattern + if (Tb > 0.0 && Tb < 0.5 * mult_ts) { + Tb = 0.0; + zero_TB = true; + } else if (Tb > 0.5 * mult_ts && Tb < mult_ts) { + Tb = mult_ts; + } else if (Tb <= 0.0) { + zero_TB = true; + } + + // Tb1 uses 0.5/1.0 pattern + if (Tb1 > 0.0 && Tb1 < 0.5 * mult_ts) { + Tb1 = 0.0; + zero_TB1 = true; + } else if (Tb1 > 0.5 * mult_ts && Tb1 < mult_ts) { + Tb1 = mult_ts; + } else if (Tb1 <= 0.0) { + zero_TB1 = true; + } + + // Tr uses 0.25/0.5 pattern (filter type) + if (Tr > 0.0 && Tr < 0.25 * mult_ts) { + Tr = 0.0; + zero_TR = true; + } else if (Tr > 0.25 * mult_ts && Tr < 0.5 * mult_ts) { + Tr = 0.5 * mult_ts; + } else if (Tr <= 0.0) { + zero_TR = true; + } + + // Ka: if Ka == 0 then set to mult_ts + if (fabs(Ka) < 1e-6) { + Ka = mult_ts; + } + + // Ta (regulator time constant) - if zero, use gain block + if (fabs(Ta) < 1e-6) { + zero_TA = true; + } + // For lead lag block, force time constant of numerator to zero if time constant of denominator is zero if (zero_TB1) Tc1 = 0.0; if (zero_TB) Tc = 0.0; @@ -137,9 +203,9 @@ void gridpack::dynamic_simulation::Esst1aModel::init(double mag, double ang, dou if (!zero_TB1) Leadlag_blkBC1.setparams(Tc1, Tb1); if(!zero_TA) { - Regulator_blk.setparams(Ka,Ta,Vrmin,Vrmax,-1000.0,1000.0); + Regulator_blk.setparams(Ka,Ta,Vamin,Vamax,-1000.0,1000.0); } else { - Regulator_gain_blk.setparams(Ka,Vrmin,Vrmax); + Regulator_gain_blk.setparams(Ka,Vamin,Vamax); } HVGate_blk2.setparams(Vuel); // UEL is Vuel? @@ -248,14 +314,22 @@ void gridpack::dynamic_simulation::Esst1aModel::computeModel(double t_inc,Integr if (UEL == 2.0) leadlag_blk_in = HVGate_blk1.getoutput(leadlag_blk_in); - VLL = Leadlag_blkBC.getoutput(leadlag_blk_in, t_inc, int_flag, true); - - VLL1 = Leadlag_blkBC1.getoutput(VLL, t_inc, int_flag, true); + if (!zero_TB) { + VLL = Leadlag_blkBC.getoutput(leadlag_blk_in, t_inc, int_flag, true); + } else { + VLL = leadlag_blk_in; + } + + if (!zero_TB1) { + VLL1 = Leadlag_blkBC1.getoutput(VLL, t_inc, int_flag, true); + } else { + VLL1 = VLL; + } if (zero_TA) { VA = Regulator_gain_blk.getoutput(VLL1); } else { - VA = Regulator_blk.getoutput(VLL, t_inc, int_flag, true); + VA = Regulator_blk.getoutput(VLL1, t_inc, int_flag, true); } double u1 = 0.0; diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/esst4b.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/esst4b.cpp index b5de42b70..a3b091381 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/esst4b.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/esst4b.cpp @@ -57,7 +57,7 @@ void gridpack::dynamic_simulation::Esst4bModel::load( if (!data->getValue(EXCITER_VRMIN, &Vrmin, idx)) Vrmin = 0.0; // Vrmin if (!data->getValue(EXCITER_TA, &Ta, idx)) Ta = 0.0; // Ta if (!data->getValue(EXCITER_KPM, &Kpm, idx)) Kpm = 0.0; // Kpm - if (!data->getValue(EXCITER_KPM, &Kim, idx)) Kim = 0.0; // Kim + if (!data->getValue(EXCITER_KIM, &Kim, idx)) Kim = 0.0; // Kim if (!data->getValue(EXCITER_VMMAX, &Vmmax, idx)) Vmmax = 0.0; // Vmmax if (!data->getValue(EXCITER_VMMIN, &Vmmin, idx)) Vmmin = 0.0; // Vmax if (!data->getValue(EXCITER_KG, &Kg, idx)) Kg = 0.0; // Kg diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/exdc1.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/exdc1.cpp index e6c177c0f..f663c5388 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/exdc1.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/exdc1.cpp @@ -77,9 +77,12 @@ void gridpack::dynamic_simulation::Exdc1Model::load( if (!data->getValue(EXCITER_SE2, &SE2, idx)) SE2 = 0.0; // SE2 if(fabs(SE1*SE2) < 1e-6) has_Sat = false; - if(fabs(TB*TC) < 1e-6) has_leadlag = false; - if(fabs(TA) < 1e-6) zero_TA = true; - if(fabs(TR) < 1e-6) zero_TR = true; + if(fabs(TB) < 1e-6) has_leadlag = false; + + // Swap Vrmax/Vrmin if inverted + if (Vrmax < Vrmin) { + double tmp = Vrmax; Vrmax = Vrmin; Vrmin = tmp; + } if(has_Sat) { /* Calculate saturation function constants */ @@ -95,7 +98,7 @@ void gridpack::dynamic_simulation::Exdc1Model::load( } if(has_leadlag) { - Leadlag_blk.setparams(TA,TB); + Leadlag_blk.setparams(TC,TB); } if(!zero_TA) { @@ -137,6 +140,54 @@ double gridpack::dynamic_simulation::Exdc1Model::Sat(double x) */ void gridpack::dynamic_simulation::Exdc1Model::init(double Vm, double Va, double ts) { + // Time constant minimum checks + double mult_ts = TS_THRESHOLD * ts; + + // Tr uses 0.25/0.5 pattern (filter type) + if (TR > 0.0 && TR < 0.25 * mult_ts) { + TR = 0.0; + zero_TR = true; + // Vmeas_blk already not set up when zero_TR + } else if (TR > 0.25 * mult_ts && TR < 0.5 * mult_ts) { + TR = 0.5 * mult_ts; + Vmeas_blk.setparams(1.0, TR); + } else if (fabs(TR) < 1e-6) { + zero_TR = true; + } + + // Tb uses 0.5/1.0 pattern (lead-lag denominator) + if (TB > 0.0 && TB < 0.5 * mult_ts) { + TB = 0.0; + has_leadlag = false; + } else if (TB > 0.5 * mult_ts && TB < mult_ts) { + TB = mult_ts; + if (has_leadlag) Leadlag_blk.setparams(TC, TB); + } + + // Te: If 0 < Te < mult_ts then Te = mult_ts + if (TE > 0.0 && TE < mult_ts) { + TE = mult_ts; + Output_blk.setparams(TE); + } + + // Tf1: If 0 < Tf1 < mult_ts then Tf1 = mult_ts + if (TF1 > 0.0 && TF1 < mult_ts) { + TF1 = mult_ts; + double a[2], b[2]; + a[0] = TF1; a[1] = 1.0; + b[0] = KF; b[1] = 0.0; + Feedback_blk.setcoeffs(a, b); + } + + // Ta: If 0 < Ta < mult_ts then Ta = mult_ts + if (TA > 0.0 && TA < mult_ts) { + TA = mult_ts; + zero_TA = false; + Regulator_blk.setparams(KA, TA, Vrmin, Vrmax, -1000.0, 1000.0); + } else if (fabs(TA) < 1e-6) { + zero_TA = true; + } + VF = Feedback_blk.init_given_u(Efd); double output_blk_in; @@ -149,6 +200,18 @@ void gridpack::dynamic_simulation::Exdc1Model::init(double Vm, double Va, double VR = output_blk_in + sat_signal; + // Adjust Vrmax/Vrmin for initial operating point + if (VR > Vrmax) { + Vrmax = VR; + if (!zero_TA) Regulator_blk.setparams(KA, TA, Vrmin, Vrmax, -1000.0, 1000.0); + else Regulator_gain_blk.setparams(KA, Vrmin, Vrmax); + } + if (VR < Vrmin) { + Vrmin = VR; + if (!zero_TA) Regulator_blk.setparams(KA, TA, Vrmin, Vrmax, -1000.0, 1000.0); + else Regulator_gain_blk.setparams(KA, Vrmin, Vrmax); + } + if(zero_TA) { VLL = VR/KA; } else { diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/gast.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/gast.cpp index 3f2633bb3..1ccb7ecae 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/gast.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/gast.cpp @@ -59,10 +59,15 @@ void gridpack::dynamic_simulation::GastModel::load( if (!data->getValue(GOVERNOR_T3, &T3, idx)) T3 = 10.0; if (!data->getValue(GOVERNOR_AT, &AT, idx)) AT = 0.0; if (!data->getValue(GOVERNOR_KT, &KT, idx)) KT = 0.0; - if (!data->getValue(GOVERNOR_VMAX, &VMAX, idx)) VMAX = 1.0; - if (!data->getValue(GOVERNOR_VMIN, &VMIN, idx)) VMIN = 0.0; + if (!data->getValue(GOVERNOR_VMAX, &VMAX, idx)) VMAX = 1.0; + if (!data->getValue(GOVERNOR_VMIN, &VMIN, idx)) VMIN = 0.0; if (!data->getValue(GOVERNOR_DT, &Dt, idx)) Dt = 0.0; + // Swap limits if inverted + if (VMAX < VMIN) { + double tmp = VMAX; VMAX = VMIN; VMIN = tmp; + } + fuel_valve_block.setparams(1.0,T1,VMIN,VMAX,-10000,10000); fuel_flow_block.setparams(1.0,T2); exh_temp_block.setparams(1.0,T3); @@ -76,6 +81,47 @@ void gridpack::dynamic_simulation::GastModel::load( */ void gridpack::dynamic_simulation::GastModel::init(double mag, double ang, double ts) { + // Minimum time constant checks + double mult_ts = TS_THRESHOLD * ts; + if (T1 > 0.0 && T1 < 0.25 * mult_ts) { + T1 = 0.0; + fuel_valve_block.setparams(1.0, T1, VMIN, VMAX, -10000, 10000); + } else if (T1 > 0.25 * mult_ts && T1 < 0.5 * mult_ts) { + T1 = 0.5 * mult_ts; + fuel_valve_block.setparams(1.0, T1, VMIN, VMAX, -10000, 10000); + } + if (T2 > 0.0 && T2 < 0.25 * mult_ts) { + T2 = 0.0; + fuel_flow_block.setparams(1.0, T2); + } else if (T2 > 0.25 * mult_ts && T2 < 0.5 * mult_ts) { + T2 = 0.5 * mult_ts; + fuel_flow_block.setparams(1.0, T2); + } + if (T3 > 0.0 && T3 < 0.25 * mult_ts) { + T3 = 0.0; + exh_temp_block.setparams(1.0, T3); + } else if (T3 > 0.25 * mult_ts && T3 < 0.5 * mult_ts) { + T3 = 0.5 * mult_ts; + exh_temp_block.setparams(1.0, T3); + } + + // Adjust AT to accommodate initial operating point + // Without this, generators with Pmech > AT would create an unstable + // initial condition inconsistent with the power flow solution. + if (Pmech > AT && AT > 0) { + AT = Pmech; + } + + // Adjust valve limits to accommodate initial operating point + if (Pmech > VMAX) { + VMAX = Pmech; + fuel_valve_block.setparams(1.0, T1, VMIN, VMAX, -10000, 10000); + } + if (Pmech < VMIN) { + VMIN = Pmech; + fuel_valve_block.setparams(1.0, T1, VMIN, VMAX, -10000, 10000); + } + // Work backwards from Pmech double fuel_flow_block_out = Pmech + Dt*delta_w; @@ -95,6 +141,11 @@ void gridpack::dynamic_simulation::GastModel::init(double mag, double ang, doubl Pref = fuel_valve_block_in + delta_w/R; // Here the assumption is lv_gate_in1 > fuel_valve_block_in. // If this condition is not satisfied then Pref value is indeterminate. + + // Clamp Pref to AT + if (AT > 0 && Pref > AT) { + Pref = AT; + } } /** diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/genrou.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/genrou.cpp index ba4c82836..aa8b9215c 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/genrou.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/genrou.cpp @@ -13,6 +13,11 @@ * @Modified: November 27, 2022, Fixed the model to validate against PSSE * * @Modified: Dec 9, 2022, print voltage and generator power + * + * @Modified: Mar 2026, Yousu Chen + * - Fixed saturation: unscaled quadratic Se=B*(x-A)^2, Sat at Psi_ag, + * iterative q-axis saturation in init/predictor/corrector. + * * @brief : Round rotor generator model * * @@ -112,12 +117,48 @@ void gridpack::dynamic_simulation::GenrouGenerator::load( } else enableSat = true; printFlag = false; - double tmp = sqrt(p_pg*p_pg +p_qg*p_qg); - // Increase Machine base to 1.2*Sgen if Sgen > MBase - if ( tmp > MBase) { - MBase = tmp*1.2; + // Note: MBASE from RAW file is the per-unit base for DYR parameters. + // Do not auto-adjust MBASE when |Sgen| > MBASE; this is valid operation. + + // GENROU machine parameter auto-correction (matches PowerWorld rules) + // 1. Transient reactances + if (Xdp > Xd) { + printf("GENROU bus %d: Auto-correct Xd'=%.6f > Xd=%.6f -> Xd'=%.6f\n", + p_bus_id, Xdp, Xd, 0.8*Xd); + Xdp = 0.8 * Xd; } - + if (Xqp > Xq) { + printf("GENROU bus %d: Auto-correct Xq'=%.6f > Xq=%.6f -> Xq'=%.6f\n", + p_bus_id, Xqp, Xq, Xq); + Xqp = Xq; + } + // 2. Subtransient reactances (Xdpp = Xqpp for GENROU) + double Xmin = (Xdp < Xqp) ? Xdp : Xqp; + if (Xdpp > Xmin) { + printf("GENROU bus %d: Auto-correct Xd''=%.6f > min(Xd',Xq')=%.6f -> Xd''=%.6f\n", + p_bus_id, Xdpp, Xmin, 0.8*Xmin); + Xdpp = 0.8 * Xmin; + Xqpp = Xdpp; + } + if (Xdpp < 0.05) { + printf("GENROU bus %d: Auto-correct Xd''=%.6f < 0.05 -> Xd''=0.05\n", + p_bus_id, Xdpp); + Xdpp = 0.05; + Xqpp = Xdpp; + } + // 3. Leakage reactance + if (Xl > Xdpp) { + printf("GENROU bus %d: Auto-correct Xl=%.6f > Xd''=%.6f -> Xl=%.6f\n", + p_bus_id, Xl, Xdpp, 0.8*Xdpp); + Xl = 0.8 * Xdpp; + } + // 4. Inertia + if (H < 0.1) { + printf("GENROU bus %d: Auto-correct H=%.6f < 0.1 -> H=0.1\n", + p_bus_id, H); + H = 0.1; + } + } /** @@ -149,22 +190,23 @@ void gridpack::dynamic_simulation::GenrouGenerator::updateData( */ double gridpack::dynamic_simulation::GenrouGenerator::Sat(double x) { - if (enableSat) { - // the following is another method for saturation computation, add by renke - double a_ = S12 / S10 - 1.0; - double b_ = -2 * S12 / S10 + 2.4; - double c_ = S12 / S10 - 1.44; - double A = (-b_ - sqrt(b_ * b_ - 4 * a_ * c_)) / (2 * a_); + if (enableSat && x > 1e-6) { + // PowerWorld standard scaled saturation: Se(x) = B*(x-A)^2 / x + // Fitted from: Se(1.0)=S10, Se(1.2)=S12: + // B*(1-A)^2/1.0 = S10, B*(1.2-A)^2/1.2 = S12 + // => R = 1.2*S12/S10 = (1.2-A)^2/(1-A)^2 + double R = 1.2 * S12 / S10; + double sqrtR = sqrt(R); + double A = (1.2 - sqrtR) / (1.0 - sqrtR); double B = S10 / ((1.0 - A) * (1.0 - A)); - - double tmp = x-A; - //double tmpin = tmp; - if (tmp<0.0) { + + double tmp = x - A; + if (tmp < 0.0) { tmp = 0.0; } - double result = B * tmp * tmp; - - return result; // Scaled Quadratic with 1.7.1 equations + double result = B * tmp * tmp / x; + + return result; } else { return 0.0; } @@ -206,33 +248,49 @@ void gridpack::dynamic_simulation::GenrouGenerator::init(double mag, // Speed deviation x2w = 0; - // Machine angle - x1d = atan2(Viterm + Ir * Xq + Ii * Ra, Vrterm + Ir * Ra - Ii * Xq); + // Iterative angle initialization with q-axis saturation. + // At steady state, iron saturation reduces the effective q-axis + // synchronous reactance. Saturation reduces the mutual inductance + // but not the leakage, so: Xq_eff = Xl + (Xq - Xl) / (1 + Sq) + // where Sq = Se(Psi_ag) * (Xq - Xl) / (Xd - Xl). + // We iterate until the angle converges. + double Xq_eff = Xq; + double Vd, Vq; + double Psiqpp, Psidpp, Psi_ag; - // axis of rotation is q-axis (lagging behind d-axis by 90 degrees) - double theta = pi/2.0 - x1d; - - // Generator currents in machine dq axis reference frame - rotate(Ir,Ii,theta, &Id, &Iq); + for (int sat_iter = 0; sat_iter < 10; sat_iter++) { + // Machine angle from voltage behind (Ra + j*Xq_eff) + x1d = atan2(Viterm + Ir * Xq_eff + Ii * Ra, Vrterm + Ir * Ra - Ii * Xq_eff); - // Generator internal voltage in network reference frame - double Vr = Vrterm + Ra * Ir - Xdpp * Ii; // internal voltage on network reference - double Vi = Viterm + Ra * Ii + Xdpp * Ir; // internal voltage on network reference + double theta = pi/2.0 - x1d; + rotate(Ir, Ii, theta, &Id, &Iq); - // Generator internal voltage in machine reference frame - double Vd, Vq; - rotate(Vr,Vi,theta, &Vd, &Vq); - - Vd = Vr * sin(x1d) - Vi * cos(x1d); // convert to dq reference - Vq = Vr * cos(x1d) + Vi * sin(x1d); // convert to dq reference - - double Psiqpp = -Vd; - double Psidpp = + Vq; + // Norton internal voltage (behind Ra + jXdpp) in network frame + double Vr = Vrterm + Ra * Ir - Xdpp * Ii; + double Vi = Viterm + Ra * Ii + Xdpp * Ir; + + // Rotate to machine dq reference frame + Vd = Vr * sin(x1d) - Vi * cos(x1d); + Vq = Vr * cos(x1d) + Vi * sin(x1d); + + Psiqpp = -Vd; + Psidpp = +Vq; + Psi_ag = sqrt(Psidpp * Psidpp + Psiqpp * Psiqpp); + + double Se = Sat(Psi_ag); + if (Se < 1e-10) break; // no saturation, angle is exact + + double Sq = Se * (Xq - Xl) / (Xd - Xl); + double Xq_eff_new = Xl + (Xq - Xl) / (1.0 + Sq); + + if (fabs(Xq_eff_new - Xq_eff) < 1e-8) break; // converged + Xq_eff = Xq_eff_new; + } // q-axis transient voltage x3Eqp = Vq + (Xdp - Xdpp)*Id; - // d-axis transient voltage + // d-axis transient voltage (consistent with saturated angle) x6Edp = Vd - (Xqp - Xqpp)*Iq; // d-axis flux @@ -240,10 +298,11 @@ void gridpack::dynamic_simulation::GenrouGenerator::init(double mag, // q-axis flux x5Psiqp = x6Edp + (Xqp - Xl)*Iq; - - // Field voltage - Efd = x3Eqp * (1 + Sat(x3Eqp)) + Id * (Xd - Xdp); - + + // Field voltage — evaluate saturation at air-gap flux magnitude + // Standard: LadIfd = E'q + Se(Psi_ag)*Psidpp + (Xd-Xdp)*Id (TempD=0 at SS) + Efd = x3Eqp + Sat(Psi_ag) * Psidpp + Id * (Xd - Xdp); + // Field current LadIfd = Efd; @@ -252,6 +311,7 @@ void gridpack::dynamic_simulation::GenrouGenerator::init(double mag, Efdinit = Efd; Pmechinit = Pmech; + // printf("print: yuan debug here, inside genrou model, Ir=%f, Ii=%f\n", Ir, Ii); // Initialize exciter if (p_hasExciter) { @@ -330,38 +390,50 @@ void gridpack::dynamic_simulation::GenrouGenerator::rebalanceEquilibrium() // --- Solve for equilibrium Id, Iq at the Norton voltage --- // At equilibrium, the internal voltages become: // Vd = gamma*Iq, Vq = x3 - alpha*Id - // where alpha = Xdp - Xdpp, gamma = Xq - Xdpp (since Xqpp = Xdpp). - // Substituting into the stator algebraic equations - // Id = (Vd-Vdterm)*G - (Vq-Vqterm)*B - // Iq = (Vd-Vdterm)*B + (Vq-Vqterm)*G - // gives the 2x2 system A * [Id, Iq]^T = rhs. + // where alpha = Xdp - Xdpp, gamma = Xq_eff - Xdpp. + // With q-axis saturation: Xq_eff = Xl + (Xq-Xl)/(1+Sq), + // so gamma depends on saturation which depends on Id,Iq. + // We iterate until gamma converges. double alpha = Xdp - Xdpp; - double gamma = Xq - Xdpp; // = Xq - Xqpp since Xqpp = Xdpp + double gamma = Xq - Xdpp; // initial unsaturated value + double Psidpp, Psiqpp, Psi_ag; + + for (int sat_iter = 0; sat_iter < 10; sat_iter++) { + double a11 = 1.0 - alpha * B; + double a12 = -gamma * G; + double a21 = alpha * G; + double a22 = 1.0 - gamma * B; + + double rhs1 = -Vdterm * G - (x3Eqp - Vqterm) * B; + double rhs2 = -Vdterm * B + (x3Eqp - Vqterm) * G; + + double det = a11 * a22 - a12 * a21; + Id = (rhs1 * a22 - rhs2 * a12) / det; + Iq = (a11 * rhs2 - a21 * rhs1) / det; + + Psidpp = x3Eqp - (Xdp - Xdpp) * Id; + Psiqpp = -gamma * Iq; + Psi_ag = sqrt(Psidpp * Psidpp + Psiqpp * Psiqpp); - double a11 = 1.0 - alpha * B; - double a12 = -gamma * G; - double a21 = alpha * G; - double a22 = 1.0 - gamma * B; + double Se = Sat(Psi_ag); + if (Se < 1e-10) break; - double rhs1 = -Vdterm * G - (x3Eqp - Vqterm) * B; - double rhs2 = -Vdterm * B + (x3Eqp - Vqterm) * G; + double Sq = Se * (Xq - Xl) / (Xd - Xl); + double gamma_new = Xl + (Xq - Xl) / (1.0 + Sq) - Xdpp; - double det = a11 * a22 - a12 * a21; - Id = (rhs1 * a22 - rhs2 * a12) / det; - Iq = (a11 * rhs2 - a21 * rhs1) / det; + if (fabs(gamma_new - gamma) < 1e-8) break; + gamma = gamma_new; + } // --- Update flux states to equilibrium --- x4Psidp = x3Eqp - (Xdp - Xl) * Id; - x6Edp = (Xq - Xqp) * Iq; + // E'd from circuit (Vd = gamma*Iq, x6Edp = Vd - (Xqp-Xqpp)*Iq) + x6Edp = gamma * Iq - (Xqp - Xqpp) * Iq; x5Psiqp = x6Edp + (Xqp - Xl) * Iq; - // Recompute subtransient fluxes (for Telec and watch output) - double Psidpp = x3Eqp - (Xdp - Xdpp) * Id; // = x3 - alpha*Id - double Psiqpp = -(Xq - Xdpp) * Iq; // = -gamma*Iq - // Electrical torque and field current (TempD = 0 at equilibrium) double Telec = Psidpp * Iq - Psiqpp * Id; - LadIfd = x3Eqp * (1 + Sat(x3Eqp)) + (Xd - Xdp) * Id; + LadIfd = x3Eqp + Sat(Psi_ag) * Psidpp + (Xd - Xdp) * Id; // Set Pmech = Telec so dx2w = 0 Pmech = Telec; @@ -442,11 +514,11 @@ void gridpack::dynamic_simulation::GenrouGenerator::predictor_currentInjection(b G = Ra / (Ra * Ra + Xdpp * Xdpp); // Setup - double Psiqpp = - x6Edp * (Xqpp - Xl) / (Xqp - Xl) - x5Psiqp * (Xqp - Xqpp) / (Xqp - Xl); + double Psiqpp = - x6Edp * (Xqpp - Xl) / (Xqp - Xl) - x5Psiqp * (Xqp - Xqpp) / (Xqp - Xl); double Psidpp = + x3Eqp * (Xdpp - Xl) / (Xdp - Xl) + x4Psidp * (Xdp - Xdpp) / (Xdp - Xl); - double Vd = - Psiqpp; //* (1 + x2w); - double Vq = + Psidpp; //* (1 + x2w); + double Vd = - Psiqpp * (1 + x2w); + double Vq = + Psidpp * (1 + x2w); Vterm = presentMag; Theta = presentAng; @@ -501,7 +573,6 @@ void gridpack::dynamic_simulation::GenrouGenerator::predictor( p_exciter->setOmega(x2w); p_exciter->setVterminal(presentMag); p_exciter->setVcomp(presentMag); - p_exciter->setFieldCurrent(LadIfd); Efd = p_exciter->getFieldVoltage(); } else { @@ -520,8 +591,8 @@ void gridpack::dynamic_simulation::GenrouGenerator::predictor( double Psiqpp = - x6Edp * (Xqpp - Xl) / (Xqp - Xl) - x5Psiqp * (Xqp - Xqpp) / (Xqp - Xl); double Psidpp = + x3Eqp * (Xdpp - Xl) / (Xdp - Xl) + x4Psidp * (Xdp - Xdpp) / (Xdp - Xl); - double Vd = - Psiqpp; - double Vq = + Psidpp; + double Vd = - Psiqpp * (1 + x2w); + double Vq = + Psidpp * (1 + x2w); Vterm = presentMag; Theta = presentAng; @@ -538,7 +609,11 @@ void gridpack::dynamic_simulation::GenrouGenerator::predictor( double Telec = Psidpp * Iq - Psiqpp * Id; double TempD = (Xdp - Xdpp) / ((Xdp - Xl) * (Xdp - Xl)) * (-x4Psidp - (Xdp - Xl) * Id + x3Eqp); - LadIfd = x3Eqp * (1 + Sat(x3Eqp)) + (Xd - Xdp) * (Id + TempD); // update Ifd later + + double Psi_ag = sqrt(Psidpp * Psidpp + Psiqpp * Psiqpp); + double Se = Sat(Psi_ag); + // Standard: LadIfd = E'q + Se*Psidpp + (Xd-Xdp)*(Id+TempD) + LadIfd = x3Eqp + Se * Psidpp + (Xd - Xdp) * (Id + TempD); dx1d = x2w * 2 * pi * 60; // 60 represents the nominal frequency of 60 Hz dx2w = 1 / (2 * H) * ((Pmech - D * x2w) / (1 + x2w) - Telec); @@ -547,7 +622,10 @@ void gridpack::dynamic_simulation::GenrouGenerator::predictor( dx5Psiqp = (-x5Psiqp + (Xqp - Xl) * Iq + x6Edp) / Tqopp; double TempQ = (Xqp - Xqpp) / ((Xqp - Xl) * (Xqp - Xl)) * (-x5Psiqp + (Xqp - Xl) * Iq + x6Edp); - dx6Edp = (-x6Edp + (Xq - Xqp) * (Iq - TempQ)) / Tqop; + // Standard: q-axis saturation with gamma_qd = (Xq-Xl)/(Xd-Xl) + // Note: our Psiqpp = -psi_q'' (NREL), so sign is + + dx6Edp = (-x6Edp + (Xq - Xqp) * (Iq - TempQ) + + Se * Psiqpp * (Xq - Xl) / (Xd - Xl)) / Tqop; x1d_1 = x1d + dx1d * t_inc; x2w_1 = x2w + dx2w * t_inc; @@ -575,6 +653,7 @@ void gridpack::dynamic_simulation::GenrouGenerator::predictor( if (p_hasPss) { p_exciter->setVstab(Vstab); } + p_exciter->setFieldCurrent(LadIfd); p_exciter->predictor(t_inc, flag); } @@ -630,21 +709,21 @@ void gridpack::dynamic_simulation::GenrouGenerator::corrector_currentInjection(b G = Ra / (Ra * Ra + Xdpp * Xdpp); //printf("B = %f, G = %f\n", B, G); // Setup - double Psiqpp = - x6Edp_1 * (Xqpp - Xl) / (Xqp - Xl) - x5Psiqp_1 * (Xqp - Xqpp) / (Xqp - Xl); + double Psiqpp = - x6Edp_1 * (Xqpp - Xl) / (Xqp - Xl) - x5Psiqp_1 * (Xqp - Xqpp) / (Xqp - Xl); double Psidpp = + x3Eqp_1 * (Xdpp - Xl) / (Xdp - Xl) + x4Psidp_1 * (Xdp - Xdpp) / (Xdp - Xl); - - double Vd = -Psiqpp; - double Vq = +Psidpp; - + + double Vd = -Psiqpp * (1 + x2w_1); + double Vq = +Psidpp * (1 + x2w_1); + Vterm = presentMag; Theta = presentAng; double Vrterm = Vterm * cos(Theta); double Viterm = Vterm * sin(Theta); double Vdterm = Vrterm * sin(x1d_1) - Viterm * cos(x1d_1); double Vqterm = Vrterm * cos(x1d_1) + Viterm * sin(x1d_1); - - gridpack::ComplexType vt_complex_tmp = gridpack::ComplexType(Vrterm, Viterm); - + + gridpack::ComplexType vt_complex_tmp = gridpack::ComplexType(Vrterm, Viterm); + double Idnorton = Vd * G - Vq * B; double Iqnorton = Vd * B + Vq * G; @@ -681,13 +760,12 @@ void gridpack::dynamic_simulation::GenrouGenerator::corrector( p_exciter->setOmega(x2w_1); p_exciter->setVterminal(presentMag); p_exciter->setVcomp(presentMag); - p_exciter->setFieldCurrent(LadIfd); - + Efd = p_exciter->getFieldVoltage(); } else { Efd = Efdinit; } - + if (p_hasGovernor) { p_governor = getGovernor(); p_governor->setRotorSpeedDeviation(x2w_1); @@ -695,14 +773,14 @@ void gridpack::dynamic_simulation::GenrouGenerator::corrector( } else { Pmech = Pmechinit; } - + double pi = 4.0*atan(1.0); - double Psiqpp = - x6Edp_1 * (Xqpp - Xl) / (Xqp - Xl) - x5Psiqp_1 * (Xqp - Xqpp) / (Xqp - Xl); + double Psiqpp = - x6Edp_1 * (Xqpp - Xl) / (Xqp - Xl) - x5Psiqp_1 * (Xqp - Xqpp) / (Xqp - Xl); double Psidpp = + x3Eqp_1 * (Xdpp - Xl) / (Xdp - Xl) + x4Psidp_1 * (Xdp - Xdpp) / (Xdp - Xl); - - double Vd = - Psiqpp; - double Vq = + Psidpp; - + + double Vd = - Psiqpp * (1 + x2w_1); + double Vq = + Psidpp * (1 + x2w_1); + Vterm = presentMag; Theta = presentAng; double Vrterm = Vterm * cos(Theta); @@ -713,29 +791,33 @@ void gridpack::dynamic_simulation::GenrouGenerator::corrector( //DQ Axis Id = (Vd - Vdterm) * G - (Vq - Vqterm) * B; Iq = (Vd - Vdterm) * B + (Vq - Vqterm) * G; - + double Telec = Psidpp * Iq - Psiqpp * Id; double TempD = (Xdp - Xdpp) / ((Xdp - Xl) * (Xdp - Xl)) * (-x4Psidp_1 - (Xdp - Xl) * Id + x3Eqp_1); - LadIfd = x3Eqp_1 * (1 + Sat(x3Eqp_1)) + (Xd - Xdp) * (Id + TempD); // update Ifd later + + double Psi_ag = sqrt(Psidpp * Psidpp + Psiqpp * Psiqpp); + double Se = Sat(Psi_ag); + // Standard: LadIfd = E'q + Se*Psidpp + (Xd-Xdp)*(Id+TempD) + LadIfd = x3Eqp_1 + Se * Psidpp + (Xd - Xdp) * (Id + TempD); dx1d_1 = x2w_1 * 2 * pi * 60; // 60 represents the nominal frequency of 60 Hz - //printf("H = %f, Pmech = %f, D = %f, x2w_1 = %f, Telec = %f\n", H, Pmech, D, x2w_1, Telec); - dx2w_1 = 1 / (2 * H) * ((Pmech - D * x2w_1) / (1 + x2w_1) - Telec); //TBD: call Governor for Pmech (Done) - dx3Eqp_1 = (Efd - LadIfd) / Tdop; //TBD: call Exciter for Efd and LadIfd (Done) + dx2w_1 = 1 / (2 * H) * ((Pmech - D * x2w_1) / (1 + x2w_1) - Telec); + dx3Eqp_1 = (Efd - LadIfd) / Tdop; dx4Psidp_1 = (-x4Psidp_1 - (Xdp - Xl) * Id + x3Eqp_1) / Tdopp; dx5Psiqp_1 = (-x5Psiqp_1 + (Xqp - Xl) * Iq + x6Edp_1) / Tqopp; double TempQ = (Xqp - Xqpp) / ((Xqp - Xl) * (Xqp - Xl)) * (-x5Psiqp_1 + (Xqp - Xl) * Iq + x6Edp_1); - - dx6Edp_1 = (-x6Edp_1 + (Xq - Xqp) * (Iq - TempQ)) / Tqop; - + // Standard: q-axis saturation with gamma_qd = (Xq-Xl)/(Xd-Xl) + dx6Edp_1 = (-x6Edp_1 + (Xq - Xqp) * (Iq - TempQ) + + Se * Psiqpp * (Xq - Xl) / (Xd - Xl)) / Tqop; + x1d = x1d + (dx1d + dx1d_1) / 2.0 * t_inc; x2w = x2w + (dx2w + dx2w_1) / 2.0 * t_inc; x3Eqp = x3Eqp + (dx3Eqp + dx3Eqp_1) / 2.0 * t_inc; x4Psidp = x4Psidp + (dx4Psidp + dx4Psidp_1) / 2.0 * t_inc; x5Psiqp = x5Psiqp + (dx5Psiqp + dx5Psiqp_1) / 2.0 * t_inc; x6Edp = x6Edp + (dx6Edp + dx6Edp_1) / 2.0 * t_inc; - + // PSS: run before exciter corrector so Vstab is available if (p_hasPss) { boost::shared_ptr pss = getPss(); @@ -750,6 +832,7 @@ void gridpack::dynamic_simulation::GenrouGenerator::corrector( if (p_hasPss) { p_exciter->setVstab(Vstab); } + p_exciter->setFieldCurrent(LadIfd); p_exciter->corrector(t_inc, flag); } @@ -851,7 +934,7 @@ bool gridpack::dynamic_simulation::GenrouGenerator::serialWrite( ret = true; } else if(!strcmp(signal,"watch_header")) { if(getWatch()) { - char buf[128]; + char buf[256]; std::string tag; if(p_ckt[1] != ' ') { tag = p_ckt; @@ -859,9 +942,9 @@ bool gridpack::dynamic_simulation::GenrouGenerator::serialWrite( tag = p_ckt[0]; } sprintf(buf,", %d_%s_V, %d_%s_Pg, %d_%s_Qg,%d_%s_angle, %d_%s_speed," - " %d_%s_Efd, %d_%s_Pm",p_bus_id,tag.c_str(),p_bus_id,tag.c_str(), + " %d_%s_Efd, %d_%s_Pm, %d_%s_PowerAngle",p_bus_id,tag.c_str(),p_bus_id,tag.c_str(), p_bus_id,tag.c_str(),p_bus_id,tag.c_str(),p_bus_id,tag.c_str(), - p_bus_id,tag.c_str(),p_bus_id,tag.c_str()); + p_bus_id,tag.c_str(),p_bus_id,tag.c_str(),p_bus_id,tag.c_str()); if (strlen(buf) <= bufsize) { sprintf(string,"%s",buf); ret = true; @@ -874,8 +957,10 @@ bool gridpack::dynamic_simulation::GenrouGenerator::serialWrite( } else if (!strcmp(signal,"watch")) { if (getWatch()) { char buf[256]; - sprintf(buf,",%f,%f,%f,%f, %f, %f, %f",Vterm,genP*MBase/p_sbase, - genQ*MBase/p_sbase,x1d_1, x2w_1+1.0, Efd,Pmech); + double powerAngle = fmod((x1d_1 - presentAng) * 180.0 / 3.14159265358979323846 + 180.0, 360.0) - 180.0; + if (powerAngle < -180.0) powerAngle += 360.0; + sprintf(buf,",%f,%f,%f,%f, %f, %f, %f, %f",Vterm,genP*MBase/p_sbase, + genQ*MBase/p_sbase,x1d_1, x2w_1+1.0, Efd,Pmech,powerAngle); if (strlen(buf) <= bufsize) { sprintf(string,"%s",buf); ret = true; diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/gensal.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/gensal.cpp index c21fcac58..40ea1245c 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/gensal.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/gensal.cpp @@ -6,8 +6,15 @@ // ----------------------------------------------------------- /** * @file gensal.cpp - * - * @brief: Salient pole generator model - no saturation + * + * @updated Yousu Chen + * - Fixed predictor: use predicted x2w_1 for governor speed deviation + * instead of pre-prediction x2w_0 (consistent with GENROU and PSS). + * - Fixed Sat(): handle zero saturation (S10=S12=0) to avoid NaN from + * division by zero, which caused KSP divergence on initialization. + * @date 2026-03-10 + * + * @brief: Salient pole generator model - no saturation * * */ @@ -83,9 +90,29 @@ void gridpack::dynamic_simulation::GensalGenerator::load( if (!data->getValue(GENERATOR_S1, &S10, idx)) S10=0.17; // S10 TBD: check parser if (!data->getValue(GENERATOR_S12, &S12, idx)) S12=0.55; // S12 TBD: check parser - double tmp = sqrt(p_pg*p_pg +p_qg*p_qg); - if ( tmp > MBase) { - MBase = tmp*1.2; + // Note: MBASE from RAW file is the per-unit base for DYR parameters. + // Do not auto-adjust MBASE when |Sgen| > MBASE; this is valid operation. + + // GENSAL machine parameter auto-correction (matches PowerWorld rules) + if (Xdp > Xd) { + printf("GENSAL bus %d: Auto-correct Xd'=%.6f > Xd=%.6f -> Xd'=%.6f\n", + p_bus_id, Xdp, Xd, 0.8*Xd); + Xdp = 0.8 * Xd; + } + if (Xdpp > Xdp) { + printf("GENSAL bus %d: Auto-correct Xd''=%.6f > Xd'=%.6f -> Xd''=%.6f\n", + p_bus_id, Xdpp, Xdp, 0.8*Xdp); + Xdpp = 0.8 * Xdp; + } + if (Xdpp < 0.05) { + printf("GENSAL bus %d: Auto-correct Xd''=%.6f < 0.05 -> Xd''=0.05\n", + p_bus_id, Xdpp); + Xdpp = 0.05; + } + if (Xl > Xdpp) { + printf("GENSAL bus %d: Auto-correct Xl=%.6f > Xd''=%.6f -> Xl=%.6f\n", + p_bus_id, Xl, Xdpp, 0.8*Xdpp); + Xl = 0.8 * Xdpp; } } @@ -118,20 +145,28 @@ void gridpack::dynamic_simulation::GensalGenerator::load( */ double gridpack::dynamic_simulation::GensalGenerator::Sat(double x) { - double a_ = S12 / S10 - 1.0; - double b_ = -2 * S12 / S10 + 2.4; - double c_ = S12 / S10 - 1.44; - double A = (-b_ - sqrt(b_ * b_ - 4 * a_ * c_)) / (2 * a_); + // No saturation when S10 and S12 are both zero + if (S10 == 0.0 && S12 == 0.0) return 0.0; + // Guard against S10 = 0 with S12 != 0 (degenerate) + if (S10 == 0.0) return 0.0; + if (x < 1e-6) return 0.0; + + // PSS/E scaled quadratic: Se(x) = B*(x-A)^2 / x + // Fitted from: Se(1.0)=S10, Se(1.2)=S12: + // B*(1-A)^2/1.0 = S10, B*(1.2-A)^2/1.2 = S12 + // => R = 1.2*S12/S10 = (1.2-A)^2/(1-A)^2 + double R = 1.2 * S12 / S10; + double sqrtR = sqrt(R); + double A = (1.2 - sqrtR) / (1.0 - sqrtR); double B = S10 / ((1.0 - A) * (1.0 - A)); - - double tmp = x-A; - if (tmp<0.0) { + double tmp = x - A; + if (tmp < 0.0) { tmp = 0.0; } - double result = B * tmp * tmp; - - return result; // Scaled Quadratic with 1.7.1 equations + double result = B * tmp * tmp / x; + + return result; } /** @@ -255,8 +290,8 @@ void gridpack::dynamic_simulation::GensalGenerator::predictor_currentInjection(b double Psiqpp = x5Psiqpp_0; // this will be different for GENROU double Psidpp = + x3Eqp_0 * (Xdpp - Xl) / (Xdp - Xl) + x4Psidp_0* (Xdp - Xdpp) / (Xdp - Xl); - double Vd = -Psiqpp;// * (1 + x2w_0); - double Vq = +Psidpp;// * (1 + x2w_0); + double Vd = -Psiqpp * (1 + x2w_0); + double Vq = +Psidpp * (1 + x2w_0); Vterm = presentMag; Theta = presentAng; double Vrterm = Vterm * cos(Theta); @@ -349,8 +384,8 @@ void gridpack::dynamic_simulation::GensalGenerator::predictor( double Psiqpp = x5Psiqpp_0; // this will be different for GENROU double Psidpp = + x3Eqp_0 * (Xdpp - Xl) / (Xdp - Xl) + x4Psidp_0* (Xdp - Xdpp) / (Xdp - Xl); - double Vd = -Psiqpp;// * (1 + x2w_0); - double Vq = +Psidpp;// * (1 + x2w_0); + double Vd = -Psiqpp * (1 + x2w_0); + double Vq = +Psidpp * (1 + x2w_0); Vterm = presentMag; Theta = presentAng; double Vrterm = Vterm * cos(Theta); @@ -408,7 +443,7 @@ void gridpack::dynamic_simulation::GensalGenerator::predictor( } if (p_hasGovernor){ - p_governor->setRotorSpeedDeviation(x2w_0); + p_governor->setRotorSpeedDeviation(x2w_1); p_governor->predictor(t_inc, flag); } @@ -443,8 +478,8 @@ void gridpack::dynamic_simulation::GensalGenerator::corrector_currentInjection(b double Psiqpp = x5Psiqpp_1; // this will be different for GENROU double Psidpp = + x3Eqp_1 * (Xdpp - Xl) / (Xdp - Xl) + x4Psidp_1 * (Xdp - Xdpp) / (Xdp - Xl); - double Vd = -Psiqpp;// * (1 + x2w_1); - double Vq = +Psidpp;// * (1 + x2w_1); + double Vd = -Psiqpp * (1 + x2w_1); + double Vq = +Psidpp * (1 + x2w_1); Vterm = presentMag; Theta = presentAng; double Vrterm = Vterm * cos(Theta); @@ -530,8 +565,8 @@ void gridpack::dynamic_simulation::GensalGenerator::corrector( double Psiqpp = x5Psiqpp_1; // this will be different for GENROU double Psidpp = + x3Eqp_1 * (Xdpp - Xl) / (Xdp - Xl) + x4Psidp_1 * (Xdp - Xdpp) / (Xdp - Xl); - double Vd = -Psiqpp;// * (1 + x2w_1); - double Vq = +Psidpp;// * (1 + x2w_1); + double Vd = -Psiqpp * (1 + x2w_1); + double Vq = +Psidpp * (1 + x2w_1); Vterm = presentMag; Theta = presentAng; double Vrterm = Vterm * cos(Theta); @@ -663,15 +698,15 @@ bool gridpack::dynamic_simulation::GensalGenerator::serialWrite( ret = true; } else if(!strcmp(signal,"watch_header")) { if(getWatch()) { - char buf[128]; + char buf[256]; std::string tag; if(p_ckt[1] != ' ') { tag = p_ckt; } else { tag = p_ckt[0]; } - sprintf(buf,", %d_%s_V, %d_%s_Pg, %d_%s_Qg,%d_%s_angle, %d_%s_speed, %d_%s_Efd, %d_%s_Pm",p_bus_id,tag.c_str(),p_bus_id,tag.c_str(),p_bus_id,tag.c_str(),p_bus_id,tag.c_str(), - p_bus_id,tag.c_str(),p_bus_id,tag.c_str(),p_bus_id,tag.c_str()); + sprintf(buf,", %d_%s_V, %d_%s_Pg, %d_%s_Qg,%d_%s_angle, %d_%s_speed, %d_%s_Efd, %d_%s_Pm, %d_%s_PowerAngle",p_bus_id,tag.c_str(),p_bus_id,tag.c_str(),p_bus_id,tag.c_str(),p_bus_id,tag.c_str(), + p_bus_id,tag.c_str(),p_bus_id,tag.c_str(),p_bus_id,tag.c_str(),p_bus_id,tag.c_str()); if (strlen(buf) <= bufsize) { sprintf(string,"%s",buf); ret = true; @@ -684,8 +719,10 @@ bool gridpack::dynamic_simulation::GensalGenerator::serialWrite( } else if (!strcmp(signal,"watch")) { if (getWatch()) { char buf[256]; - sprintf(buf,",%f,%f,%f,%f, %f, %f, %f", - Vterm,genP*MBase/p_sbase,genQ*MBase/p_sbase,x1d_1, x2w_1+1.0, Efd,Pmech); + double powerAngle = fmod((x1d_1 - presentAng) * 180.0 / 3.14159265358979323846 + 180.0, 360.0) - 180.0; + if (powerAngle < -180.0) powerAngle += 360.0; + sprintf(buf,",%f,%f,%f,%f, %f, %f, %f, %f", + Vterm,genP*MBase/p_sbase,genQ*MBase/p_sbase,x1d_1, x2w_1+1.0, Efd,Pmech,powerAngle); if (strlen(buf) <= bufsize) { sprintf(string,"%s",buf); ret = true; diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/ggov1.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/ggov1.cpp index 7456dacad..413d5670d 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/ggov1.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/ggov1.cpp @@ -73,76 +73,41 @@ void gridpack::dynamic_simulation::Ggov1Model::load( boost::shared_ptr data, int idx) { - //if (!data->getValue(GOVERNOR_RSELECT, &Rselect, idx)) - Rselect = 0.0; - //if (!data->getValue(GOVERNOR_FLAG, &Flag, idx)) - Flag = 0.0; - //if (!data->getValue(GOVERNOR_R, &R, idx)) - R = 0.0; - //if (!data->getValue(GOVERNOR_TPELEC, &Tpelec, idx)) - Tpelec = 0.0; - //if (!data->getValue(GOVERNOR_MAXERR, &MaxErr, idx)) - MaxErr = 0.0; - //if (!data->getValue(GOVERNOR_MINERR, &MinErr, idx)) - MinErr = 0.0; - //if (!data->getValue(GOVERNOR_KPGOV, &Kpgov, idx)) - Kpgov = 0.0; - //if (!data->getValue(GOVERNOR_KIGOV, &Kigov, idx)) - Kigov = 0.0; - //if (!data->getValue(GOVERNOR_KDGOV, &Kdgov, idx)) - Kdgov = 0.0; - //if (!data->getValue(GOVERNOR_TDGOV, &Tdgov, idx)) - Tdgov = 0.0; - //if (!data->getValue(GOVERNOR_VMAX, &Vmax, idx)) - Vmax = 0.0; - //if (!data->getValue(GOVERNOR_VMIN, &Vmin, idx)) - Vmin = 0.0; - //if (!data->getValue(GOVERNOR_TACT, &Tact, idx)) - Tact = 0.0; - //if (!data->getValue(GOVERNOR_KTURB, &Kturb, idx)) - Kturb = 0.0; - //if (!data->getValue(GOVERNOR_WFNL, &Wfnl, idx)) - Wfnl = 0.0; - //if (!data->getValue(GOVERNOR_TB, &Tb, idx)) - Tb = 0.0; - //if (!data->getValue(GOVERNOR_TC, &Tc, idx)) - Tc = 0.0; - //if (!data->getValue(GOVERNOR_TENG, &Teng, idx)) - Teng = 0.0; - //if (!data->getValue(GOVERNOR_TFLOAD, &Tfload, idx)) - Tfload = 0.0; - //if (!data->getValue(GOVERNOR_KPLOAD, &Kpload, idx)) - Kpload = 0.0; - //if (!data->getValue(GOVERNOR_KILOAD, &Kiload, idx)) - Kiload = 0.0; - //if (!data->getValue(GOVERNOR_LDREF, &Ldref, idx)) - Ldref = 0.0; - //if (!data->getValue(GOVERNOR_DM, &Dm, idx)) - Dm = 0.0; - //if (!data->getValue(GOVERNOR_ROPEN, &Ropen, idx)) - Ropen = 0.0; - //if (!data->getValue(GOVERNOR_RCLOSE, &Rclose, idx)) - Rclose = 0.0; - //if (!data->getValue(GOVERNOR_KIMW, &Kimw, idx)) - Kimw = 0.0; - //if (!data->getValue(GOVERNOR_ASET, &Aset, idx)) - Aset = 0.0; - //if (!data->getValue(GOVERNOR_KA, &Ka, idx)) - Ka = 0.0; - //if (!data->getValue(GOVERNOR_TA, &Ta, idx)) - Ta = 0.0; - //if (!data->getValue(GOVERNOR_TRATE, &Trate, idx)) - Trate = 0.0; - //if (!data->getValue(GOVERNOR_DB, &Db, idx)) - Db = 0.0; - //if (!data->getValue(GOVERNOR_TSA, &Tsa, idx)) - Tsa = 0.0; - //if (!data->getValue(GOVERNOR_TSB, &Tsb, idx)) - Tsb = 0.0; - //if (!data->getValue(GOVERNOR_RUP, &Rup, idx)) - Rup = 0.0; - //if (!data->getValue(GOVERNOR_RDOWN, &Rdown, idx)) - Rdown = 0.0; + if (!data->getValue(GOVERNOR_RSELECT, &Rselect, idx)) Rselect = 0.0; + if (!data->getValue(GOVERNOR_FLAGSWITCH, &Flag, idx)) Flag = 0.0; + if (!data->getValue(GOVERNOR_R, &R, idx)) R = 0.0; + if (!data->getValue(GOVERNOR_TPELEC, &Tpelec, idx)) Tpelec = 0.0; + if (!data->getValue(GOVERNOR_MAXERR, &MaxErr, idx)) MaxErr = 0.0; + if (!data->getValue(GOVERNOR_MINERR, &MinErr, idx)) MinErr = 0.0; + if (!data->getValue(GOVERNOR_KPGOV, &Kpgov, idx)) Kpgov = 0.0; + if (!data->getValue(GOVERNOR_KIGOV, &Kigov, idx)) Kigov = 0.0; + if (!data->getValue(GOVERNOR_KDGOV, &Kdgov, idx)) Kdgov = 0.0; + if (!data->getValue(GOVERNOR_TDGOV, &Tdgov, idx)) Tdgov = 0.0; + if (!data->getValue(GOVERNOR_VMAX, &Vmax, idx)) Vmax = 0.0; + if (!data->getValue(GOVERNOR_VMIN, &Vmin, idx)) Vmin = 0.0; + if (!data->getValue(GOVERNOR_TACT, &Tact, idx)) Tact = 0.0; + if (!data->getValue(GOVERNOR_KTURB, &Kturb, idx)) Kturb = 0.0; + if (!data->getValue(GOVERNOR_WFNL, &Wfnl, idx)) Wfnl = 0.0; + if (!data->getValue(GOVERNOR_TB, &Tb, idx)) Tb = 0.0; + if (!data->getValue(GOVERNOR_TC, &Tc, idx)) Tc = 0.0; + if (!data->getValue(GOVERNOR_TENG, &Teng, idx)) Teng = 0.0; + if (!data->getValue(GOVERNOR_TFLOAD, &Tfload, idx)) Tfload = 0.0; + if (!data->getValue(GOVERNOR_KPLOAD, &Kpload, idx)) Kpload = 0.0; + if (!data->getValue(GOVERNOR_KILOAD, &Kiload, idx)) Kiload = 0.0; + if (!data->getValue(GOVERNOR_LDREF, &Ldref, idx)) Ldref = 0.0; + if (!data->getValue(GOVERNOR_DM, &Dm, idx)) Dm = 0.0; + if (!data->getValue(GOVERNOR_ROPEN, &Ropen, idx)) Ropen = 0.0; + if (!data->getValue(GOVERNOR_RCLOSE, &Rclose, idx)) Rclose = 0.0; + if (!data->getValue(GOVERNOR_KIMW, &Kimw, idx)) Kimw = 0.0; + if (!data->getValue(GOVERNOR_ASET, &Aset, idx)) Aset = 0.0; + if (!data->getValue(GOVERNOR_KA, &Ka, idx)) Ka = 0.0; + if (!data->getValue(GOVERNOR_TA, &Ta, idx)) Ta = 0.0; + if (!data->getValue(GOVERNOR_TRATE, &Trate, idx)) Trate = 0.0; + if (!data->getValue(GOVERNOR_DB, &Db, idx)) Db = 0.0; + if (!data->getValue(GOVERNOR_TSA, &Tsa, idx)) Tsa = 0.0; + if (!data->getValue(GOVERNOR_TSB, &Tsb, idx)) Tsb = 0.0; + if (!data->getValue(GOVERNOR_RUP, &Rup, idx)) Rup = 0.0; + if (!data->getValue(GOVERNOR_RDOWN, &Rdown, idx)) Rdown = 0.0; if (!data->getValue(GOVERNOR_DB1, &Db1, idx)) Db1 = 0.0; // Db1 if (!data->getValue(GOVERNOR_ERR, &Err, idx)) Err = 0.0; // Err @@ -256,7 +221,7 @@ void gridpack::dynamic_simulation::Ggov1Model::predictor(double t_inc, bool flag if (dx8LoadCtrl > 0 && x8LoadCtrl >= 1.1 * R) dx8LoadCtrl = 0; if (dx8LoadCtrl < 0 && x8LoadCtrl <= -1.1 * R) dx8LoadCtrl = 0; // State 9 - if (Ka > 0) dx9Accel = w * Ta + x9Accel; + if (Ka > 0 && Ta > 0) dx9Accel = (Ka * w - x9Accel) / Ta; else dx9Accel = 0; // State 10 and State 6 // Note: if Tact = 0, then we are really using the previous timestep @@ -281,7 +246,7 @@ void gridpack::dynamic_simulation::Ggov1Model::predictor(double t_inc, bool flag // Find the input to PID controller // Note: feedback of x4Act and LastLowValueSelect is really // the previous time step here - double PIDIn = Pref + x8LoadCtrl - w - PIDIn; + double PIDIn = Pref + x8LoadCtrl - w; if (Rselect == +1) PIDIn = PIDIn - R * x1Pelec; else if (Rselect == -1) PIDIn = PIDIn - R * x4Act; else if (Rselect == -2) PIDIn = PIDIn - R * LastLowValueSelect; @@ -335,12 +300,12 @@ void gridpack::dynamic_simulation::Ggov1Model::predictor(double t_inc, bool flag if (Flag == 1) TempIn = TempIn * (1 + w); TempIn = (TempIn - Wfnl) * Kturb; // Note: We are ignoring the Engine Delay here - if (Tb > TS_THRESHOLD) { + if (Tb < TS_THRESHOLD * t_inc) { dx5LL = 0; LeadLagOut = TempIn; } else { - dx5LL = (TempIn * (1 - Ta / Tb) - x5LL) / Tb; - LeadLagOut = TempIn * Ta / Tb + x5LL; + dx5LL = (TempIn * (1 - Tc / Tb) - x5LL) / Tb; + LeadLagOut = TempIn * Tc / Tb + x5LL; } x1Pelec_1 = x1Pelec + dx1Pelec * t_inc; @@ -386,7 +351,7 @@ void gridpack::dynamic_simulation::Ggov1Model::corrector(double t_inc, bool flag if (dx8LoadCtrl_1 > 0 && x8LoadCtrl_1 >= 1.1 * R) dx8LoadCtrl_1 = 0; if (dx8LoadCtrl_1 < 0 && x8LoadCtrl_1 <= -1.1 * R) dx8LoadCtrl_1 = 0; // State 9 - if (Ka > 0) dx9Accel_1 = w * Ta + x9Accel_1; + if (Ka > 0 && Ta > 0) dx9Accel_1 = (Ka * w - x9Accel_1) / Ta; else dx9Accel_1 = 0; // State 10 and State 6 // Note: if Tact = 0, then we are really using the previous timestep @@ -411,7 +376,7 @@ void gridpack::dynamic_simulation::Ggov1Model::corrector(double t_inc, bool flag // Find the input to PID controller // Note: feedback of x4Act and LastLowValueSelect is really // the previous time step here - double PIDIn = Pref + x8LoadCtrl_1 - w - PIDIn; + double PIDIn = Pref + x8LoadCtrl_1 - w; if (Rselect == +1) PIDIn = PIDIn - R * x1Pelec_1; else if (Rselect == -1) PIDIn = PIDIn - R * x4Act_1; else if (Rselect == -2) PIDIn = PIDIn - R * LastLowValueSelect; @@ -458,19 +423,19 @@ void gridpack::dynamic_simulation::Ggov1Model::corrector(double t_inc, bool flag } else { dx4Act_1 = (LastLowValueSelect - x4Act_1) / Tact; if (dx4Act_1 > Ropen) dx4Act_1 = Ropen; - if (dx4Act_1 < Rclose) dx4Act = Rclose; + if (dx4Act_1 < Rclose) dx4Act_1 = Rclose; } // State 5 TempIn = x4Act_1; if (Flag == 1) TempIn = TempIn * (1 + w); TempIn = (TempIn - Wfnl) * Kturb; // Note: We are ignoring the Engine Delay here - if (Tb > TS_THRESHOLD) { + if (Tb < TS_THRESHOLD * t_inc) { dx5LL_1 = 0; LeadLagOut = TempIn; } else { - dx5LL_1 = (TempIn * (1 - Ta / Tb) - x5LL_1) / Tb; - LeadLagOut = TempIn * Ta / Tb + x5LL_1; + dx5LL_1 = (TempIn * (1 - Tc / Tb) - x5LL_1) / Tb; + LeadLagOut = TempIn * Tc / Tb + x5LL_1; } x1Pelec_1 = x1Pelec + (dx1Pelec + dx1Pelec_1) / 2.0 * t_inc; diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/ieel.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/ieel.cpp index a4bc33f59..e7b2abb0b 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/ieel.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/ieel.cpp @@ -180,23 +180,15 @@ void gridpack::dynamic_simulation::IeelLoad::init(double mag, a1 = a1/ (a1 + a2 + a3); a2 = a2/ (a1 + a2 + a3); a3 = a3/ (a1 + a2 + a3); - } else { - a1 = 1.0; - a2 = 0.0; - a3 = 0.0; } } - - // check the data to make sure a4+a5+a6 = 1 + + // check the data to make sure a4+a5+a6 = 1 if (a4 + a5 + a6 > 0.0) { if (abs(a4 + a5 + a6 - 1.0) >1.0E-6) { a4 = a4/ (a4 + a5 + a6); a5 = a5/ (a4 + a5 + a6); a6 = a6/ (a4 + a5 + a6); - } else { - a4 = 1.0; - a5 = 0.0; - a6 = 0.0; } } } diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/motorw.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/motorw.cpp index 490ddc5bb..1d1da06f3 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/motorw.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/motorw.cpp @@ -187,9 +187,6 @@ void gridpack::dynamic_simulation::MotorwLoad::load( if (motorwload_perc< 0.01) motorwload_perc = 1.0; //if the LOAD_DYN_PERC is not defined in the dyr file, as the normal PSS/E dyr file - //tmp code please remove this - llr1 = 0.08; - data->getValue(LOAD_ID,&p_loadid,idx); setDynLoadP(p_pl); setDynLoadID(p_loadid); @@ -330,11 +327,13 @@ void gridpack::dynamic_simulation::MotorwLoad::init(double mag, //printf(" MotorwLoad::init(), Pe[%d]: %12.6f, Sl[%d]: %12.6f, \n", k, Pe[k], k, sl[k]); } - // find specific slip for initial power (Pint) - double sl0; - for (int m = 1; m < 1000; m++) { + // find specific slip for initial power (Pini) with linear interpolation + double sl0 = 0.01; // default + for (int m = 1; m < 999; m++) { if (Pe[m] <= Pini && Pe[m+1] > Pini) { - sl0 = sl[m] ; // found slip0 + // Linear interpolation for precise slip + double frac = (Pini - Pe[m]) / (Pe[m+1] - Pe[m]); + sl0 = sl[m] + frac * (sl[m+1] - sl[m]); break; } } @@ -619,7 +618,7 @@ void gridpack::dynamic_simulation::MotorwLoad::init(double mag, TL = p * eppd * Id + q * eppq * Iq ; w = 1.0 - slip ; // rotor speed, pu - C0 = 1,0 - A*w*w - B*w - D*(pow(w,E)); + C0 = 1.0 - A*w*w - B*w - D*(pow(w,E)); Tm0 = TL; gridpack::ComplexType tmp2(Id, Iq); @@ -634,11 +633,11 @@ void gridpack::dynamic_simulation::MotorwLoad::init(double mag, if (bdebugprint) printf(" MotorwLoad::init(), bus %d, after slightly adjust, C0: %f, Tm0: %f, Pmotor: %f, Qmotor: %f, \n", p_bus_id, C0, Tm0, Pmotor, Qmotor); epq0 = epq; - epd0 = epd; + epd0 = epd; eppq0 = eppq; eppd0 = eppd; slip0 = slip; - + setDynLoadP(Pini); // need to reset Dynamic load P here, as we applied the motorwload_perc during initialization Qmotor_init = Qmotor; diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/reeca1.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/reeca1.cpp index adae85b2a..50b655570 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/reeca1.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/reeca1.cpp @@ -12,6 +12,9 @@ * @brief * Renewable Energy Electrical Controller Model REECA1 * + * @Modified: Mar 2026, Yousu Chen + * - Fixed CurrentLimitLogic + * */ #include @@ -324,7 +327,7 @@ void gridpack::dynamic_simulation::Reeca1Model::predictor(double t_inc, // Check if thld == 0, in this case there is // no transition to state 2 if(fabs(Thld) < 1e-6) Iqinj_sw = 0; - else if(fabs(Thld) > 1e-6) { + else if(Thld > 1e-6) { // Thld is positive, transition to state 2, // Set timer thld_timer = 0.0; @@ -337,7 +340,7 @@ void gridpack::dynamic_simulation::Reeca1Model::predictor(double t_inc, } else { if(Iqinj_sw == 1) { thld_timer += t_inc; - if(thld_timer > 1 - Thld) Iqinj_sw = 0; + if(thld_timer > -Thld) Iqinj_sw = 0; } else if(Iqinj_sw == 2) { thld_timer += t_inc; if(thld_timer > Thld) Iqinj_sw = 0; @@ -372,7 +375,7 @@ void gridpack::dynamic_simulation::Reeca1Model::corrector(double t_inc, void gridpack::dynamic_simulation::Reeca1Model::computeModel(bool Voltage_dip,int Iqinj_sw,double t_inc,IntegrationStage int_flag) { - bool updateState = !Voltage_dip; // freezestate logic : updateState = 0 when no voltage dip, otherwise 1 + bool updateState = !Voltage_dip; // freezestate logic : updateState = 1 when no voltage dip, 0 during voltage dip (freeze states) // Get terminal voltage measurement Vt_filter = Vt_filter_blk.getoutput(Vt,t_inc,int_flag,true); @@ -464,10 +467,9 @@ void gridpack::dynamic_simulation::Reeca1Model::CurrentLimitLogic(int PQFLAG,dou *Ipmin_out = 0.0; - // Iqmax look up from VDL1 - Iqmax_temp = VDL1.getoutput(Vt_filter); - // Ipmax look up from VDL2 - Ipmax_temp = VDL2.getoutput(Vt_filter); + // VDL1 = V-P table (active current limit), VDL2 = V-Q table (reactive current limit) + Ipmax_temp = VDL1.getoutput(Vt_filter); + Iqmax_temp = VDL2.getoutput(Vt_filter); if(!PQFLAG) { // Q priority if (Imax < Iqmax_temp) Iqmax_temp = Imax; diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/regca1.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/regca1.cpp index cbbbb85ab..f96710025 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/regca1.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/regca1.cpp @@ -10,8 +10,10 @@ * @Added: November 14, 2022 * * @brief WECC generic generator/converter model - * - * + * + * @Modified: Mar 2026, Yousu Chen + * - Fixed Iq_olim sign + * */ #include @@ -151,10 +153,10 @@ void gridpack::dynamic_simulation::Regca1Generator::updateData( Iq_olim = std::max(0.0,khv*(Vt - volim)); - Iqout = Iqlowlim_blk.getoutput(Iq + Iq_olim); + Iqout = Iqlowlim_blk.getoutput(Iq - Iq_olim); // Pg, Qg on machine MVAbase - Pg = Vt*Ipout*p_mbase/p_sbase; + Pg = Vt*Ipout*p_mbase/p_sbase; Qg = -Vt*Iqout*p_mbase/p_sbase; if (!data->setValue(GENERATOR_PG_CURRENT, Pg, idx)) { @@ -569,10 +571,10 @@ bool gridpack::dynamic_simulation::Regca1Generator::serialWrite( Iq_olim = std::max(0.0,khv*(Vt - volim)); - Iqout = Iqlowlim_blk.getoutput(Iq + Iq_olim); + Iqout = Iqlowlim_blk.getoutput(Iq - Iq_olim); // Pg, Qg on machine MVAbase - Pg = Vt*Ipout*p_mbase/p_sbase; + Pg = Vt*Ipout*p_mbase/p_sbase; Qg = -Vt*Iqout*p_mbase/p_sbase; sprintf(string,",%12.6f,%12.6f, %12.6f, %12.6f, %12.6f, %12.6f, %12.6f ", @@ -645,12 +647,12 @@ void gridpack::dynamic_simulation::Regca1Generator::getWatchValues( Iq_olim = std::max(0.0,khv*(Vt - volim)); - Iqout = Iqlowlim_blk.getoutput(Iq + Iq_olim); + Iqout = Iqlowlim_blk.getoutput(Iq - Iq_olim); // Pg, Qg on system MVAbase Pg = Vt*Ipout*p_mbase/p_sbase; Qg = -Vt*Iqout*p_mbase/p_sbase; - + if (p_generatorObservationPowerSystemBase){ vals.push_back(Pg); //output at system mva base vals.push_back(Qg); //output at system mva base diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/regcb1.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/regcb1.cpp index 106fbab46..cd4d8ff6d 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/regcb1.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/regcb1.cpp @@ -313,7 +313,7 @@ void gridpack::dynamic_simulation::Regcb1Generator::computeModel(double t_inc, I bool Vdip; // Voltage dip flag for electrical controller double domega_t = 0.0; - double Pord = 0.0,Preal; + double Pord = 0.0,Preal = 0.0; double Pref_plant; double Ip_blk_in; diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/regcc1.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/regcc1.cpp index 344f000f2..92a0d5750 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/regcc1.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/regcc1.cpp @@ -334,7 +334,7 @@ void gridpack::dynamic_simulation::Regcc1Generator::computeModel(double t_inc, I bool Vdip; // Voltage dip flag for electrical controller double domega_t = 0.0; - double Pord = 0.0,Preal; + double Pord = 0.0,Preal = 0.0; double Pref_plant; double Ip_blk_in; diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/repca1.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/repca1.cpp index 55917e9ff..96e846861 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/repca1.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/repca1.cpp @@ -17,6 +17,7 @@ #include #include #include +#include #include "boost/smart_ptr/shared_ptr.hpp" #include "gridpack/parser/dictionary.hpp" @@ -87,9 +88,23 @@ void gridpack::dynamic_simulation::Repca1Model::load( if (!data->getValue(GENERATOR_REPCA_PMIN, &Pmin , idx)) Pmin = -1.5; if (!data->getValue(GENERATOR_REPCA_TG, &Tg , idx)) Tg = 0.1; + // Femax/Femin swap and sign + if (femax < femin) { + double tmp = femax; femax = femin; femin = tmp; + } + if (femax < 0) femax = -femax; + if (femin > 0) femin = -femin; + + // Emax/Emin swap and sign + if (Emax < Emin) { + double tmp = Emax; Emax = Emin; Emin = tmp; + } + if (Emax < 0) Emax = -Emax; + if (Emin > 0) Emin = -Emin; + // Set constants Freq_ref = 1.0; - + // Set up model blocks // V filter block V_filter_blk.setparams(1.0,Tfltr); @@ -132,6 +147,39 @@ void gridpack::dynamic_simulation::Repca1Model::load( */ void gridpack::dynamic_simulation::Repca1Model::init(double Vm, double Va, double ts) { + // Time constant minimum checks + double mult_ts = 4.0 * ts; // TS_THRESHOLD * ts + + // Ddn: If 0 < Ddn < mult_ts then Ddn = mult_ts + if (Ddn > 0.0 && Ddn < mult_ts) { + Ddn = mult_ts; + } + + // Dup: If 0 < Dup < mult_ts then Dup = mult_ts (note: Dup is typically negative) + if (fabs(Dup) > 0.0 && fabs(Dup) < mult_ts) { + Dup = (Dup < 0) ? -mult_ts : mult_ts; + } + + // Tfltr (called Tflr in PW docs): 0.5/1.0 pattern + if (Tfltr > 0.0 && Tfltr < 0.5 * mult_ts) { + Tfltr = 0.0; + V_filter_blk.setparams(1.0, Tfltr); + Qbranch_filter_blk.setparams(1.0, Tfltr); + } else if (Tfltr > 0.5 * mult_ts && Tfltr < mult_ts) { + Tfltr = mult_ts; + V_filter_blk.setparams(1.0, Tfltr); + Qbranch_filter_blk.setparams(1.0, Tfltr); + } + + // Tp: 0.5/1.0 pattern + if (Tp > 0.0 && Tp < 0.5 * mult_ts) { + Tp = 0.0; + Pbranch_filter_blk.setparams(1.0, Tp); + } else if (Tp > 0.5 * mult_ts && Tp < mult_ts) { + Tp = mult_ts; + Pbranch_filter_blk.setparams(1.0, Tp); + } + if(FreqFLAG) { double ferr; Pref_PI_blk_out = Pref_filter_blk.init_given_y(Pref); @@ -187,7 +235,7 @@ void gridpack::dynamic_simulation::Repca1Model::computeModel(double t_inc,Integr if(FreqFLAG) { ferr = Freq_ref - Freq; ferr = Freqerr_deadband.getoutput(ferr); - dP = std::max(0.0,Ddn*ferr) + std::min(0.0,Dup*ferr); + dP = std::max(0.0,fabs(Ddn)*ferr) + std::min(0.0,-fabs(Dup)*ferr); // *********** // Need to use Pbranch if given, using Pg // *********** diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/sexs.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/sexs.cpp index c78614d1d..ab72a97e9 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/sexs.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/sexs.cpp @@ -8,8 +8,13 @@ * @file sexs.cpp * @author Shrirang Abhyankar * @Added: Nov 6, 2022 - * - * @brief + * + * @updated Yousu Chen + * - Changed EMAX default from 0.0 to 999.0 when not specified in + * .dyr file, preventing exciter output from being clamped at zero. + * @date 2026-03-10 + * + * @brief * * */ @@ -63,9 +68,6 @@ void gridpack::dynamic_simulation::SexsModel::load( TA = TA_OVER_TB*TB; - printf("SEXS Bus %d Gen %s: K=%f TE=%f EMIN=%f EMAX=%f TA/TB=%f TB=%f\n", - p_bus_id, p_ckt.c_str(), K, TE, EMIN, EMAX, TA_OVER_TB, TB); - // Set parameters for the first block leadlagblock.setparams(TA,TB); diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/tgov1.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/tgov1.cpp index 0f5044c5c..7b56a3f43 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/tgov1.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/tgov1.cpp @@ -51,14 +51,19 @@ void gridpack::dynamic_simulation::Tgov1Model::load( boost::shared_ptr data, int idx) { - if (!data->getValue(GOVERNOR_R, &R, idx)) R = 0.05; - if (!data->getValue(GOVERNOR_T1, &T1, idx)) T1 = 0.5; - if (!data->getValue(GOVERNOR_VMAX, &Vmax, idx)) Vmax = 1.0; - if (!data->getValue(GOVERNOR_VMIN, &Vmin, idx)) Vmin = 0.0; - if (!data->getValue(GOVERNOR_T2, &T2, idx)) T2 = 3.0; - if (!data->getValue(GOVERNOR_T3, &T3, idx)) T3 = 10.0; + if (!data->getValue(GOVERNOR_R, &R, idx)) R = 0.05; + if (!data->getValue(GOVERNOR_T1, &T1, idx)) T1 = 0.5; + if (!data->getValue(GOVERNOR_VMAX, &Vmax, idx)) Vmax = 1.0; + if (!data->getValue(GOVERNOR_VMIN, &Vmin, idx)) Vmin = 0.0; + if (!data->getValue(GOVERNOR_T2, &T2, idx)) T2 = 3.0; + if (!data->getValue(GOVERNOR_T3, &T3, idx)) T3 = 10.0; if (!data->getValue(GOVERNOR_DT, &Dt, idx)) Dt = 0.0; + // Swap limits if inverted + if (Vmax < Vmin) { + double tmp = Vmax; Vmax = Vmin; Vmin = tmp; + } + // Set up transfer function blocks leadlag_blk.setparams(T2,T3); delay_blk.setparams(1.0,T1,Vmin,Vmax,-1000.0,1000.0); @@ -72,6 +77,33 @@ void gridpack::dynamic_simulation::Tgov1Model::load( */ void gridpack::dynamic_simulation::Tgov1Model::init(double mag, double ang, double ts) { + // Minimum time constant checks + double mult_ts = TS_THRESHOLD * ts; + if (T1 > 0.0 && T1 < 0.25 * mult_ts) { + T1 = 0.0; + delay_blk.setparams(1.0, T1, Vmin, Vmax, -1000.0, 1000.0); + } else if (T1 > 0.25 * mult_ts && T1 < 0.5 * mult_ts) { + T1 = 0.5 * mult_ts; + delay_blk.setparams(1.0, T1, Vmin, Vmax, -1000.0, 1000.0); + } + if (T3 > 0.0 && T3 < 0.25 * mult_ts) { + T3 = 0.0; + leadlag_blk.setparams(T2, T3); + } else if (T3 > 0.25 * mult_ts && T3 < 0.5 * mult_ts) { + T3 = 0.5 * mult_ts; + leadlag_blk.setparams(T2, T3); + } + + // Adjust valve limits to accommodate initial operating point + if (Pmech > Vmax) { + Vmax = Pmech; + delay_blk.setparams(1.0, T1, Vmin, Vmax, -1000.0, 1000.0); + } + if (Pmech < Vmin) { + Vmin = Pmech; + delay_blk.setparams(1.0, T1, Vmin, Vmax, -1000.0, 1000.0); + } + // Initialize leadlag block delay_blk_out = leadlag_blk.init_given_y(Pmech+Dt*delta_w); diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/wshygp.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/wshygp.cpp index 05d478dca..b75bd3c28 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/wshygp.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/wshygp.cpp @@ -17,6 +17,7 @@ #include #include +#include #include #include @@ -99,8 +100,11 @@ void gridpack::dynamic_simulation::WshygpModel::load( zero_TT = false; zero_TP = false; zero_TTURB_BTURB = false; - Tt = 0.0; // hard coded because there is no setPelec method in this and base_governor class - Pelec = 0.0; // the same reason + // TODO: Tt is forced to 0 because base_governor_model lacks setPelec(). + // When Tt > 0, the droop should use filtered Pelec instead of CV*R. + // To fix: add setPelec() to BaseGovernorModel and call from generator model. + Tt = 0.0; + Pelec = 0.0; // the same reason as above OptionToModifyLimitsForInitialStateLimitViolation = true; } @@ -119,7 +123,7 @@ void gridpack::dynamic_simulation::WshygpModel::init(double mag, double ang, dou if (Tf < TS_THRESHOLD * ts) zero_TF = true; if (Tt < TS_THRESHOLD * ts) zero_TT = true; if (Tp < TS_THRESHOLD * ts) zero_TP = true; - if (abs(Tturb * Bturb) < TS_THRESHOLD * ts) zero_TTURB_BTURB = true; + if (fabs(Tturb * Bturb) < TS_THRESHOLD * ts) zero_TTURB_BTURB = true; if (zero_TD) printf("Warning: Td=%f is less than %d times timestep=%f.\n", Td, TS_THRESHOLD, ts); if (zero_KI) printf("Warning: KI=%f is less than %f.\n", KI, KI_THRESHOLD); @@ -165,10 +169,10 @@ void gridpack::dynamic_simulation::WshygpModel::init(double mag, double ang, dou double uin[5], yin[5]; uin[0] = Gv1; yin[0] = PGv1; - uin[1] = Gv2; yin[0] = PGv2; - uin[2] = Gv3; yin[0] = PGv3; - uin[3] = Gv4; yin[0] = PGv4; - uin[4] = Gv5; yin[0] = PGv5; + uin[1] = Gv2; yin[1] = PGv2; + uin[2] = Gv3; yin[2] = PGv3; + uin[3] = Gv4; yin[3] = PGv4; + uin[4] = Gv5; yin[4] = PGv5; NGV_blk.setparams(5, uin, yin); if (!zero_TTURB_BTURB) Leadlag_blk.setparams(Aturb*Tturb, Bturb*Tturb); // else {Aturb*Tturb=0, a gain=1 block} @@ -260,6 +264,7 @@ void gridpack::dynamic_simulation::WshygpModel::computeModel(double t_inc,Integr else u6 = u5 * KG; // a gain=KG block u7 = Integrator_blk.getoutput(u6, t_inc, int_flag, true); + GV = u7; // update gate position feedback for next timestep u8 = Db2_blk.getoutput(u7); diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/wsieg1.cpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/wsieg1.cpp index 22b723843..4560ecb8b 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/wsieg1.cpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/wsieg1.cpp @@ -7,9 +7,11 @@ /** * @file wsieg1.cpp * - * @brief: WSIEG1 governor model implementation - * - * + * @brief: WSIEG1 governor model implementation + * + * @Modified: Mar 2026, Yousu Chen + * - Fixed NGV lookup table. + * */ #include @@ -87,7 +89,30 @@ void gridpack::dynamic_simulation::Wsieg1Model::load( if (!data->getValue(GOVERNOR_PGV5, &PGv5, idx)) PGv5 = 0.0; // PGv5 if (!data->getValue(GOVERNOR_IBLOCK, &Iblock, idx)) Iblock = 0.0; // Iblock - Db1_blk.setparams(Db1, Err); + // Adjust Uo/Uc sign and ordering + if (Uo < Uc) { + double tmp = Uo; Uo = Uc; Uc = tmp; + } + if (Uo < 0) Uo = -Uo; + if (Uc > 0) Uc = -Uc; + + // Pmax/Pmin swap + if (Pmax < Pmin) { + double tmp = Pmax; Pmax = Pmin; Pmin = tmp; + } + + // K normalization + double Ksum_odd = K1 + K3 + K5 + K7; + double Ksum_even = K2 + K4 + K6 + K8; + if (Ksum_odd > 1.0) { + K1 /= Ksum_odd; K3 /= Ksum_odd; K5 /= Ksum_odd; K7 /= Ksum_odd; + } + if (Ksum_even > 1.0) { + K2 /= Ksum_even; K4 /= Ksum_even; K6 /= Ksum_even; K8 /= Ksum_even; + } + + has_Db1 = (fabs(Db1) > 1e-6); + if (has_Db1) Db1_blk.setparams(Db1, Err); if(T1 != 0 && T2 != 0) { Leadlag_blk.setparams(T2, T1); } @@ -95,14 +120,15 @@ void gridpack::dynamic_simulation::Wsieg1Model::load( P_blk.setparams(1.0, Pmin, Pmax); // need another method to take Pmax and Pmin Db2_blk.setparams(Db2, Db2); - // Initialize NGV + // Initialize NGV — only enable if GV/PGV data is non-trivial double uin[5], yin[5]; uin[0] = Gv1; yin[0] = PGv1; - uin[1] = Gv2; yin[0] = PGv2; - uin[2] = Gv3; yin[0] = PGv3; - uin[3] = Gv4; yin[0] = PGv4; - uin[4] = Gv5; yin[0] = PGv5; - NGV_blk.setparams(5, uin, yin); + uin[1] = Gv2; yin[1] = PGv2; + uin[2] = Gv3; yin[2] = PGv3; + uin[3] = Gv4; yin[3] = PGv4; + uin[4] = Gv5; yin[4] = PGv5; + has_NGV = (fabs(Gv1) + fabs(Gv2) + fabs(Gv3) + fabs(Gv4) + fabs(Gv5) > 1e-6); + if (has_NGV) NGV_blk.setparams(5, uin, yin); if(T4 != 0) { Filter_blk1.setparams(1.0, T4); @@ -128,9 +154,17 @@ void gridpack::dynamic_simulation::Wsieg1Model::load( */ void gridpack::dynamic_simulation::Wsieg1Model::init(double mag, double ang, double ts) { - double u1, u2, u3, u4, u5, u6, u7, u8, u9; - - // Backword initialization + // T3 time constant check + double mult_ts = 4.0 * ts; // TS_THRESHOLD * ts + if (T3 > 0.0 && T3 < 0.25 * mult_ts) { + T3 = 0.0; + } else if (T3 > 0.25 * mult_ts && T3 < 0.5 * mult_ts) { + T3 = 0.5 * mult_ts; + } + + double u1, u2, u3, u4, u5, u6, u7, u8, u9; + + // Backword initialization if (K1 + K3 + K5 + K7 > 0) { u9 = Pmech1 / (K1 + K3 + K5 + K7); if(T7 == 0) { @@ -181,8 +215,7 @@ void gridpack::dynamic_simulation::Wsieg1Model::init(double mag, double ang, dou } else u5 = 0; - // u4 = NGV_blk.init_given_y(u5); // Implemented a .int_given_y method for PiecewiseSlope in dblock - u4 = u5; + u4 = has_NGV ? NGV_blk.init_given_y(u5) : u5; GV = u4; // GV can be used in the next time step, not a local variable like u1-u9 u3 = u4; // Deadband Db2_blk's input equals the output @@ -211,8 +244,7 @@ void gridpack::dynamic_simulation::Wsieg1Model::computeModel(double t_inc,Integr double u1, y1, u2, y2, u3, y3, u4, y4, u5, y5, u6, y6, u7, y7, u8, y8, u9, y9; u1 = w; - // y1 = Db1_blk.getoutput(u1); - y1 = u1; + y1 = has_Db1 ? Db1_blk.getoutput(u1) : u1; u2 = y1 * K; if(T1 == 0 || T2 == 0) { @@ -239,11 +271,10 @@ void gridpack::dynamic_simulation::Wsieg1Model::computeModel(double t_inc,Integr GV = y4; u5 = y4; - // y5 = NGV_blk.getoutput(u5); - y5 = u5; - + y5 = has_NGV ? NGV_blk.getoutput(u5) : u5; + if(T4 == 0) { - u6 = u5; + u6 = y5; } else { u6 = Filter_blk1.getoutput(u5, t_inc, int_flag, true); } diff --git a/src/applications/modules/dynamic_simulation_full_y/model_classes/wsieg1.hpp b/src/applications/modules/dynamic_simulation_full_y/model_classes/wsieg1.hpp index 9e34a3ebf..79191c21d 100644 --- a/src/applications/modules/dynamic_simulation_full_y/model_classes/wsieg1.hpp +++ b/src/applications/modules/dynamic_simulation_full_y/model_classes/wsieg1.hpp @@ -126,10 +126,10 @@ class Wsieg1Model : public BaseGovernorModel double Pref; double w; - Deadband Db1_blk; // is DBInt (Intentional Deadband) implemented in Deadband block? + Deadband_dbint Db1_blk; LeadLag Leadlag_blk; Integrator P_blk; - Deadband Db2_blk; // is BackLash (Unintentional Deadband) implemented in Deadband block? + Deadband_backlash Db2_blk; PiecewiseSlope NGV_blk; // is GainBlock a PiecewiseSlope block? yes Filter Filter_blk1; Filter Filter_blk2; @@ -137,6 +137,8 @@ class Wsieg1Model : public BaseGovernorModel Filter Filter_blk4; double GV; + bool has_NGV; // true if NGV lookup table has non-trivial data + bool has_Db1; // true if Db1 deadband is nonzero void computeModel(double t_inc, IntegrationStage int_flag); diff --git a/src/math/linear_solver_implementation.hpp b/src/math/linear_solver_implementation.hpp index 1514e7d04..ed88b97aa 100644 --- a/src/math/linear_solver_implementation.hpp +++ b/src/math/linear_solver_implementation.hpp @@ -10,8 +10,14 @@ * @file linear_solver_implementation.hpp * @author William A. Perkins * @date 2022-10-05 08:52:17 d3g096 - * - * @brief + * + * @updated Yousu Chen + * - Fixed off-by-one in p_serialSolvePrep(): setElementRange() uses + * half-open [lo,hi), so hi must be size() not size()-1. Last element + * of RHS/solution was stale on ForceSerial path with np>1. + * @date 2026-03-10 + * + * @brief * * */ @@ -199,7 +205,7 @@ class LinearSolverImplementation p_serialRHS.reset(b.localClone()); } else { b.getAllElements(&p_valueBuffer[0]); - p_serialRHS->setElementRange(0, b.size() - 1, &p_valueBuffer[0]); + p_serialRHS->setElementRange(0, b.size(), &p_valueBuffer[0]); } // Collect the initial guess, if it's not zero @@ -212,7 +218,7 @@ class LinearSolverImplementation } else { if (!p_guessZero) { x.getAllElements(&p_valueBuffer[0]); - p_serialSolution->setElementRange(0, x.size() - 1, &p_valueBuffer[0]); + p_serialSolution->setElementRange(0, x.size(), &p_valueBuffer[0]); } } } diff --git a/src/parser/parser_classes/repca1.hpp b/src/parser/parser_classes/repca1.hpp index 340c16dfa..c247eb60f 100644 --- a/src/parser/parser_classes/repca1.hpp +++ b/src/parser/parser_classes/repca1.hpp @@ -6,6 +6,8 @@ /* * Created on: May 10, 2021 * Author: Renke Huang + * Modified: Mar 2026, Yousu Chen + * - Fixed FEMIN parser index. */ #ifndef REPCA1_HPP #define REPCA1_HPP @@ -551,7 +553,7 @@ template class Repca1Parser if (!data->getValue(GENERATOR_REPCA_FEMIN,&rval,g_id)) { data->addValue(GENERATOR_REPCA_FEMIN, atof(split_line[ 31].c_str()), g_id); } else { - data->setValue(GENERATOR_REPCA_FEMIN, atof(split_line[ 32].c_str()), g_id); + data->setValue(GENERATOR_REPCA_FEMIN, atof(split_line[ 31].c_str()), g_id); } } diff --git a/src/parser/parser_classes/wsieg1.hpp b/src/parser/parser_classes/wsieg1.hpp index 8e212d8ff..f23772772 100644 --- a/src/parser/parser_classes/wsieg1.hpp +++ b/src/parser/parser_classes/wsieg1.hpp @@ -6,6 +6,8 @@ /* * Created on: June 16, 2016 * Author: Bruce Palmer + * Modified: Mar 2026, Yousu Chen + * - Fixed T3 parameter. */ #ifndef WSIEG1_HPP #define WSIEG1_HPP @@ -95,7 +97,7 @@ template class Wsieg1Parser // GOVERNOR_T3 if (!data->getValue(GOVERNOR_T3,&rval,g_id)) { - data->addValue(GOVERNOR_T3, data_struct.gv_t2, g_id); + data->addValue(GOVERNOR_T3, data_struct.gv_t3, g_id); } else { data->setValue(GOVERNOR_T3, data_struct.gv_t3, g_id); }