diff --git a/CHANGELOG.md b/CHANGELOG.md index ed853112cd..81ed51aea8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,7 @@ - Removed `BusToSignalAdapter`. Prefer `BusSignalVoltageOut` or `BusSignalVoltageIn`. - Added `Bus::setFault` to apply or clear a fault to ground directly at a bus. - Removed the `BusFault` class in favor of `Bus::setFault`. +- Added `enable` and `disable` solver events that put a branch in or out of service. ## v0.2 diff --git a/GridKit/Model/PhasorDynamics/Branch/Branch.hpp b/GridKit/Model/PhasorDynamics/Branch/Branch.hpp index 82a5144a21..8588f7869f 100644 --- a/GridKit/Model/PhasorDynamics/Branch/Branch.hpp +++ b/GridKit/Model/PhasorDynamics/Branch/Branch.hpp @@ -74,6 +74,8 @@ namespace GridKit virtual int evaluateJacobian() override final; virtual int verify() const override final; + int setInService(bool in_service); + void setR(RealT R) { R_ = R; @@ -195,6 +197,7 @@ namespace GridKit RealT Bmag_{0.0}; RealT tap_{1.0}; RealT phase_{0.0}; + RealT in_service_{1.0}; ///< 1 in service, 0 out of service IdxT bus1_id_{0}; IdxT bus2_id_{0}; diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp b/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp index f7b1391f17..90d9ea5e28 100644 --- a/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp @@ -184,6 +184,19 @@ namespace GridKit return ret; } + /// Put the branch in or out of service; out of service, it injects no current. + template + int Branch::setInService(bool in_service) + { + in_service_ = ZERO; + if (in_service) + { + in_service_ = ONE; + } + setDerivedParams(); + return 0; + } + template __attribute__((always_inline)) inline void Branch::addAdmittanceContribution( const RealT G, @@ -296,6 +309,9 @@ namespace GridKit /** * @brief Residual contribution of the branch is computed and pushed to the terminal buses. * + * Out-of-service branches retain bus voltages but contribute no current. + * Empty isolated buses and floating series-only islands can make the + * system Jacobian singular; generator-supported islands need not do so. */ template int Branch::evaluateResidual() @@ -514,6 +530,16 @@ namespace GridKit g22_ = g_diag - RealT{0.5} * G_; b22_ = b_diag - RealT{0.5} * B_; + + // An out-of-service branch has zero admittance + g11_ *= in_service_; + b11_ *= in_service_; + g12_ *= in_service_; + b12_ *= in_service_; + g21_ *= in_service_; + b21_ *= in_service_; + g22_ *= in_service_; + b22_ *= in_service_; } } // namespace PhasorDynamics diff --git a/GridKit/Model/PhasorDynamics/Branch/README.md b/GridKit/Model/PhasorDynamics/Branch/README.md index adebd94fdf..c746520bbf 100644 --- a/GridKit/Model/PhasorDynamics/Branch/README.md +++ b/GridKit/Model/PhasorDynamics/Branch/README.md @@ -13,6 +13,16 @@ contributions are oriented entering the adjacent buses. bus 1; both shunts are added outside the $\mathbf{M}$ transformation. - The branch has no solver-owned variables; it contributes current residuals directly to the connected buses. +- Taking a branch out of service can make the system Jacobian singular, which + the model does not detect. See + [Appendix A](#appendix-a-out-of-service-branches) and + [Appendix B](#appendix-b-series-resonance). + +## Block Diagram + +![Branch model block diagram](../../../../docs/Figures/PhasorDynamics/Branch/diagram.png) + +Figure 1: Branch model. ## Model Parameters @@ -94,12 +104,15 @@ The off-nominal transformer transformation uses bus 1 as the tap side: \end{aligned} ``` -The magnetizing and line shunts are added outside the transformation: +The magnetizing and line shunts are added outside the transformation, and the +service status $u$ scales the whole branch: ```math \begin{aligned} \mathbf{Y} &= + u + \left( \mathbf{M}^{\dagger} \mathbf{Y}_0 \mathbf{M} @@ -107,6 +120,7 @@ The magnetizing and line shunts are added outside the transformation: \mathbf{Y}_\mathrm{mag} + \mathbf{Y}_\mathrm{sh} + \right) \end{aligned} ``` @@ -201,22 +215,161 @@ Monitor | Units | Description | Note `p2` | [p.u.] | Active power at bus 2 terminal | Positive entering bus 2 `q2` | [p.u.] | Reactive power at bus 2 terminal | Positive entering bus 2 -Current magnitudes are: +## Appendix A: Out-of-Service Branches + +Taking a branch out of service changes the system Jacobian and can make it +singular. Let the system be described by the DAE +$\mathbf{F}(\mathbf{x}, \dot{\mathbf{x}}, t) = \mathbf{0}$. For a purely +algebraic network, $\mathbf{x} = \mathbf{V}$ are the bus voltages and + +```math +\mathbf{F}(\mathbf{V}) = -\mathbf{I} + \mathbf{Y}\mathbf{V} +``` + +where $\mathbf{I}$ is the current injected by devices and $\mathbf{Y}$ is the +bus admittance matrix. In the three-bus example below, bus 1 has a Norton +generator, a constant $I_G$ in parallel with $Y_G$, and the branch statuses are +$u_{12}$ and $u_{23}$. + +![Three-bus example](../../../../docs/Figures/PhasorDynamics/Branch/oneline.png) + +Since $\mathbf{I} = (I_G, 0, 0)$ is constant, the Jacobian is the admittance +matrix: ```math \begin{aligned} - I_{m1} &= \sqrt{I_{r1}^2 + I_{i1}^2} \\ - I_{m2} &= \sqrt{I_{r2}^2 + I_{i2}^2} + \mathbf{J} &= \dfrac{\partial \mathbf{F}}{\partial \mathbf{V}} = \mathbf{Y} \\ + \mathbf{Y} &= + u_{12} + \begin{bmatrix} + Y_{12} & -Y_{12} & 0 \\ + -Y_{12} & Y_{12} & 0 \\ + 0 & 0 & 0 + \end{bmatrix} + + u_{23} + \begin{bmatrix} + 0 & 0 & 0 \\ + 0 & Y_{23} & -Y_{23} \\ + 0 & -Y_{23} & Y_{23} + \end{bmatrix} + + + \begin{bmatrix} + Y_G & 0 & 0 \\ + 0 & 0 & 0 \\ + 0 & 0 & Y_L + \end{bmatrix} \\ + \det \mathbf{J} &= u_{12} u_{23} Y_{12} Y_{23} (Y_G + Y_L) + (u_{12} Y_{12} + u_{23} Y_{23}) Y_G Y_L \end{aligned} ``` -Complex power at each end is defined as $S=VI^{\ast}$: +### Scenario A.1: Isolated Bus + +![Scenario A.1: isolated bus](../../../../docs/Figures/PhasorDynamics/Branch/scenario1.png) + +With $u_{23} = 0$ and $Y_L = 0$: ```math \begin{aligned} - P_1 &= V_{r1} I_{r1} + V_{i1} I_{i1} \\ - Q_1 &= V_{i1} I_{r1} - V_{r1} I_{i1} \\ - P_2 &= V_{r2} I_{r2} + V_{i2} I_{i2} \\ - Q_2 &= V_{i2} I_{r2} - V_{r2} I_{i2} + \mathbf{J} &= + \begin{bmatrix} + Y_G + Y_{12} & -Y_{12} & 0 \\ + -Y_{12} & Y_{12} & 0 \\ + 0 & 0 & 0 + \end{bmatrix} \\ + \det \mathbf{J} &= 0 \end{aligned} ``` + +The bus 3 row is zero, so $V_3$ is undetermined. + +### Scenario A.2: Floating Island + +![Scenario A.2: floating island](../../../../docs/Figures/PhasorDynamics/Branch/scenario2.png) + +With $u_{12} = 0$ and $Y_L = 0$: + +```math +\begin{aligned} + \mathbf{J} &= + \begin{bmatrix} + Y_G & 0 & 0 \\ + 0 & Y_{23} & -Y_{23} \\ + 0 & -Y_{23} & Y_{23} + \end{bmatrix} \\ + \det \mathbf{J} &= 0 +\end{aligned} +``` + +The island rows sum to zero, so $V_2 = V_3$ is undetermined. (Edge case: Parallel +branches whose taps or phase shifts disagree make the island nonsingular). + +### Scenario A.3: Grounded Island + +![Scenario A.3: grounded island](../../../../docs/Figures/PhasorDynamics/Branch/scenario3.png) + +With $u_{12} = 0$ and $Y_L \neq 0$: + +```math +\begin{aligned} + \mathbf{J} &= + \begin{bmatrix} + Y_G & 0 & 0 \\ + 0 & Y_{23} & -Y_{23} \\ + 0 & -Y_{23} & Y_{23} + Y_L + \end{bmatrix} \\ + \det \mathbf{J} &= Y_G Y_{23} Y_L \neq 0 +\end{aligned} +``` + +The load grounds the island. A fault or generator at bus 3 does the same. + +## Appendix B: Series Resonance + +An outage can also make $\mathbf{J}$ singular without creating an island. In +the two-bus example below, bus 1 has a Norton generator with a capacitive +$Y_G$, as when a capacitor bank outweighs the generator's own admittance, and +the second circuit has status $u$. + +![Two-bus example](../../../../docs/Figures/PhasorDynamics/Branch/resonance.png) + +With $Y_G = j$, $Y_{12} = -j2$ per circuit, and $Y_L = -j2$: + +```math +\begin{aligned} + \mathbf{J} &= + (1 + u) + \begin{bmatrix} + Y_{12} & -Y_{12} \\ + -Y_{12} & Y_{12} + \end{bmatrix} + + + \begin{bmatrix} + Y_G & 0 \\ + 0 & Y_L + \end{bmatrix} \\ + \det \mathbf{J} &= (1 + u) Y_{12} (Y_G + Y_L) + Y_G Y_L = -2u +\end{aligned} +``` + +### Scenario B.1: Circuit Trip + +![Scenario B.1: circuit trip](../../../../docs/Figures/PhasorDynamics/Branch/resonance-trip.png) + +With $u = 0$: + +```math +\begin{aligned} + \mathbf{J} &= + j + \begin{bmatrix} + -1 & 2 \\ + 2 & -4 + \end{bmatrix} \\ + \det \mathbf{J} &= 0 +\end{aligned} +``` + +The network stays connected, but the loop impedance +$1/Y_G + 1/Y_{12} + 1/Y_L = -j + j0.5 + j0.5$ is zero, so the loop through +ground is series resonant. Any resistance makes $\det \mathbf{J} \neq 0$, but a +nearly lossless loop still leaves $\mathbf{J}$ ill-conditioned. diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index f55d8aa68e..9d8c938ce7 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -2,6 +2,7 @@ #include #include +#include #include #include @@ -22,6 +23,9 @@ namespace GridKit template class BusBase; + template + class Branch; + template class SignalNode; @@ -98,9 +102,10 @@ namespace GridKit void setSystemBase(RealT freq_system_base, RealT va_system_base); - BusT* getBus(IdxT bus_id); - SignalNodeT* getSignalNode(IdxT signal_id); - ComponentT* getComponent(IdxT gridkit_component_id); + BusT* getBus(IdxT bus_id); + SignalNodeT* getSignalNode(IdxT signal_id); + ComponentT* getComponent(IdxT gridkit_component_id); + Branch* getBranch(const std::string& id); private: std::vector buses_; @@ -109,6 +114,10 @@ namespace GridKit std::map gridkit_bus_indices_; ///< Map between gridkit_bus_id and bus_id + /// Map between branch id and gridkit_component_id + /// @todo remove this map and use branch IDs to access individual branches. + std::map gridkit_branch_indices_; + bool owns_components_{false}; /// Variable monitor diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 53acd0d42f..0d7e1859b6 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -98,6 +98,8 @@ namespace GridKit auto* branch = new Branch(getBus(bus1_index), getBus(bus2_index), branchdata); + + gridkit_branch_indices_[branchdata.disambiguation_string] = static_cast(components_.size()); addComponent(branch); } @@ -798,5 +800,16 @@ namespace GridKit return components_[gridkit_component_id]; } + /** + * @brief Return pointer to a branch by case-file `id` + */ + template + Branch* + SystemModel::getBranch(const std::string& id) + { + // Should fail if user-provided id is incorrect + return dynamic_cast*>(components_[gridkit_branch_indices_.at(id)]); + } + } // namespace PhasorDynamics } // namespace GridKit diff --git a/application/PhasorDynamics/AnalysisUtilities.hpp b/application/PhasorDynamics/AnalysisUtilities.hpp index d8a3008b39..85ed86c188 100644 --- a/application/PhasorDynamics/AnalysisUtilities.hpp +++ b/application/PhasorDynamics/AnalysisUtilities.hpp @@ -29,7 +29,7 @@ namespace GridKit * @brief Describes an event that is used to modify the simulation at the * given time point * - * ContingencyAnalysis repeats this event at each bus with the same R and X. + * ContingencyAnalysis repeats fault events at each bus with the same R and X. */ struct SystemEvent { @@ -37,7 +37,9 @@ namespace GridKit enum class Type { FAULT_ON, - FAULT_OFF + FAULT_OFF, + IN_SERVICE, + OUT_OF_SERVICE }; /// Time event takes place @@ -50,6 +52,8 @@ namespace GridKit double R{0.0}; /// Fault reactance [p.u.] double X{0.0}; + /// Case-file `id` of the branch to put in or out of service + std::string device; }; /** @@ -137,7 +141,8 @@ namespace GridKit { auto& event = c.events.emplace_back(); raw_event.at("time").get_to(event.time); - event.bus = raw_event.value("bus", INVALID_INDEX); + event.bus = raw_event.value("bus", INVALID_INDEX); + event.device = raw_event.value("device", std::string{}); auto type_str = raw_event.at("type").get(); using EventType = SystemEvent::Type; diff --git a/application/PhasorDynamics/ContingencyAnalysis.cpp b/application/PhasorDynamics/ContingencyAnalysis.cpp index b880f600cf..c1a5b61ee2 100644 --- a/application/PhasorDynamics/ContingencyAnalysis.cpp +++ b/application/PhasorDynamics/ContingencyAnalysis.cpp @@ -7,6 +7,7 @@ #include #endif +#include #include #include #include @@ -56,16 +57,21 @@ TestStatus runStudy(StudyData study_data) ida.runSimulation(event.time, dt_monitor); // Set up run for event (to start at event time) - auto* bus = sys.getBus(event.bus); - int ret = 0; + int ret = 0; switch (event.type) { case EventType::FAULT_ON: - ret = bus->setFault(true, event.R, event.X); + ret = sys.getBus(event.bus)->setFault(true, event.R, event.X); break; case EventType::FAULT_OFF: - ret = bus->clearFault(); + ret = sys.getBus(event.bus)->clearFault(); + break; + case EventType::IN_SERVICE: + sys.getBranch(event.device)->setInService(true); + break; + case EventType::OUT_OF_SERVICE: + sys.getBranch(event.device)->setInService(false); break; } @@ -91,10 +97,13 @@ TestStatus runStudy(StudyData study_data) TestStatus singleFaultStudy(std::size_t fault_id, StudyData study_data) { - // Change bus in schedule to current fault bus + // Retarget fault events while preserving branch-switching events for (auto& event : study_data.events) { - event.bus = study_data.model_data.bus[fault_id].bus_id; + if (event.type == SystemEvent::Type::FAULT_ON || event.type == SystemEvent::Type::FAULT_OFF) + { + event.bus = study_data.model_data.bus[fault_id].bus_id; + } } // Make distinct output files diff --git a/application/PhasorDynamics/DynamicSimulation.cpp b/application/PhasorDynamics/DynamicSimulation.cpp index 02dd547251..351a3a1ed1 100644 --- a/application/PhasorDynamics/DynamicSimulation.cpp +++ b/application/PhasorDynamics/DynamicSimulation.cpp @@ -2,6 +2,7 @@ #include #include +#include #include #include #include @@ -27,9 +28,12 @@ int runApplication(int argc, const char* argv[]) checkCommandLine(argc, "DynamicSimulation"); auto study = parseStudyData(argv[1]); + using EventType = SystemEvent::Type; + for (const auto& event : study.events) { - if (event.bus == GridKit::INVALID_INDEX) + if ((event.type == EventType::FAULT_ON || event.type == EventType::FAULT_OFF) + && event.bus == GridKit::INVALID_INDEX) { Log::error() << "Fault event at time " << event.time << " requires a bus number\n"; @@ -53,8 +57,6 @@ int runApplication(int argc, const char* argv[]) // Start timer real_type start = static_cast(clock()); - using EventType = SystemEvent::Type; - // Initialize simulation for first run auto dt_monitor = study.dt_monitor; real_type final_time = study.tmax; @@ -65,16 +67,21 @@ int runApplication(int argc, const char* argv[]) ida.runSimulation(event.time, dt_monitor); // Set up run for event (to start at event time) - auto* bus = sys.getBus(event.bus); - int ret = 0; + int ret = 0; switch (event.type) { case EventType::FAULT_ON: - ret = bus->setFault(true, event.R, event.X); + ret = sys.getBus(event.bus)->setFault(true, event.R, event.X); break; case EventType::FAULT_OFF: - ret = bus->clearFault(); + ret = sys.getBus(event.bus)->clearFault(); + break; + case EventType::IN_SERVICE: + sys.getBranch(event.device)->setInService(true); + break; + case EventType::OUT_OF_SERVICE: + sys.getBranch(event.device)->setInService(false); break; } diff --git a/application/PhasorDynamics/README.md b/application/PhasorDynamics/README.md index 512ddbc776..a14c5fe7ed 100644 --- a/application/PhasorDynamics/README.md +++ b/application/PhasorDynamics/README.md @@ -29,7 +29,8 @@ Each event group describes a system event that occurs at a given time point Name | Value --------------------|------------------------------------------------------- `time` | A floating point value for time event occurs - `type` | Event type (one of { "fault_on", "fault_off" }) + `type` | Event type (one of { "fault_on", "fault_off", "in_service", "out_of_service" }) `bus` | Number of the faulted bus (omitted for `ContingencyAnalysis`, which faults every bus in turn using the same `R` and `X`) `R` | Fault resistance in per unit (required for "fault_on") `X` | Fault reactance in per unit (required for "fault_on") + `device` | `id` of the branch to put in service ("in_service") or take out of service ("out_of_service") diff --git a/docs/Figures/PhasorDynamics/Branch/diagram.png b/docs/Figures/PhasorDynamics/Branch/diagram.png new file mode 100644 index 0000000000..e81eadabd0 Binary files /dev/null and b/docs/Figures/PhasorDynamics/Branch/diagram.png differ diff --git a/docs/Figures/PhasorDynamics/Branch/diagram.tex b/docs/Figures/PhasorDynamics/Branch/diagram.tex new file mode 100644 index 0000000000..0589bdd1df --- /dev/null +++ b/docs/Figures/PhasorDynamics/Branch/diagram.tex @@ -0,0 +1,42 @@ +\documentclass[tikz,border=12pt]{standalone} +\input{../../EMT/diagram-style.tex} + +\begin{document} +\begin{tikzpicture}[emt diagram] + +% Branch network: shunts remain outside the bus-1 tap transformation. +\node[hblock] (series) at (1.8, 0) {$Y_\mathrm{br}$}; +\node[vblock] (mag) at (-4.8, -1.6) {$Y_\mathrm{mag}$}; +\node[vblock] (shunt1) at (-3.0, -1.6) {$\dfrac{Y_\mathrm{sh}}{2}$}; +\node[vblock] (shunt2) at (4.8, -1.6) {$\dfrac{Y_\mathrm{sh}}{2}$}; + +% Single-line transformer symbol, with the tap on the bus-1 side. +\draw[semithick] (-1.35, 0) circle (0.45); +\draw[semithick] (-0.75, 0) circle (0.45); +\node[sig, anchor=south] at (-1.05, 0.6) {$\tau:1$}; +\node[sig, anchor=north] at (-1.05, -0.6) {$\theta$}; + +\draw[line] (-0.3, 0) -- (series.west); + +\foreach \shunt/\x in {mag/-4.8,shunt1/-3.0,shunt2/4.8} { + \node[signal] at (\x, 0) {}; + \draw[line] (\x, 0) -- (\shunt.north); + \draw[line] (\shunt.south) -- (\x, -3.2); + \draw[semithick] (\x-0.35, -3.20) -- (\x+0.35, -3.20); + \draw[semithick] (\x-0.22, -3.35) -- (\x+0.22, -3.35); + \draw[semithick] (\x-0.09, -3.50) -- (\x+0.09, -3.50); +} + +% Terminal voltages are bus-owned; current arrows point into the buses. +\draw[->, line] (-1.8, 0) -- (-7.5, 0); +\draw[->, line] (series.east) -- (7.5, 0); +\node[signal, label={[sig]above:$V_1$}] at (-6.7, 0) {}; +\node[signal, label={[sig]above:$V_2$}] at (6.7, 0) {}; +\node[sig, anchor=east] at (-7.7, 0) {$I_1$}; +\node[sig, anchor=west] at (7.7, 0) {$I_2$}; + +% Model enclosure. +\draw[enclosure] (-6.0, -4.0) rectangle (6.0, 1.5); + +\end{tikzpicture} +\end{document} diff --git a/docs/Figures/PhasorDynamics/Branch/oneline.png b/docs/Figures/PhasorDynamics/Branch/oneline.png new file mode 100644 index 0000000000..2ad8acb181 Binary files /dev/null and b/docs/Figures/PhasorDynamics/Branch/oneline.png differ diff --git a/docs/Figures/PhasorDynamics/Branch/outage.tex b/docs/Figures/PhasorDynamics/Branch/outage.tex new file mode 100644 index 0000000000..33809e4846 --- /dev/null +++ b/docs/Figures/PhasorDynamics/Branch/outage.tex @@ -0,0 +1,116 @@ +\documentclass[tikz,border=12pt]{standalone} +\input{../../EMT/diagram-style.tex} + +% Out-of-service elements are drawn dashed in red. +\definecolor{outage}{HTML}{D62828} +\tikzset{ + in service/.style = {}, + out of service/.style = {outage, densely dashed}, + bus/.style = {very thick}, + series/.style = {block, minimum width=1.6cm, minimum height=0.9cm}, + shunt/.style = {block, minimum width=1.0cm, minimum height=1.0cm} +} + +\newcommand{\ground}[1]{ + \begin{scope}[shift={#1}, solid] + \draw[semithick] (-0.35, 0.00) -- (0.35, 0.00); + \draw[semithick] (-0.22, -0.15) -- (0.22, -0.15); + \draw[semithick] (-0.09, -0.30) -- (0.09, -0.30); + \end{scope} +} + +% Constant current source at bus 1, with a tap for YG at x = -1.3. +\newcommand{\source}{ + \node[draw, semithick, fill=white, circle, minimum size=1.0cm] + (source) at (-2.7, -1.5) {}; + \draw[->, line] ($(source.south)+(0, 0.18)$) + -- ($(source.north)+(0, -0.18)$); + \node[sig, anchor=east] at ($(source.west)+(-0.15, 0)$) {$I_G$}; + \draw[line] (0, 0) -| (source.north); + \draw[line] (source.south) -- (-2.7, -2.5); + \ground{(-2.7, -2.5)} + \node[signal] at (-1.3, 0) {}; +} + +% Load at the bus at x = #2; #1 is its style. +\newcommand{\load}[2]{ + \begin{scope}[#1] + \node[shunt] (YL) at ({#2 + 1.3}, -1.5) {$Y_L$}; + \draw[line] (#2, 0) -| (YL.north); + \draw[line] (YL.south) -- ({#2 + 1.3}, -2.5); + \ground{({#2 + 1.3}, -2.5)} + \end{scope} +} + +% Bus bars as number/x pairs, drawn last to stay on top of the wires. +\newcommand{\buses}[1]{ + \foreach \k/\x in {#1} { + \draw[bus] (\x, -0.7) -- (\x, 0.7); + \node[sig, anchor=south] at (\x, 0.75) {$V_\k$}; + } +} + +% Four-turn line inductor from (#1, #2) to the next bus. +\newcommand{\inductor}[2]{ + \draw[semithick] (#1, #2) -- ++(1.4, 0) + arc[start angle=180, end angle=0, radius=0.15] + arc[start angle=180, end angle=0, radius=0.15] + arc[start angle=180, end angle=0, radius=0.15] + arc[start angle=180, end angle=0, radius=0.15] + -- ++(1.4, 0); +} + +% Three-bus chain; #1, #2, #3 are the styles of Y12, Y23, and YL. +\newcommand{\chain}[3]{ +\begin{tikzpicture}[emt diagram] +\source +\node[shunt] (YG) at (-1.3, -1.5) {$Y_G$}; +\draw[line] (-1.3, 0) -- (YG.north); +\draw[line] (YG.south) -- (-1.3, -2.5); +\ground{(-1.3, -2.5)} +\begin{scope}[#1] + \node[series] (Y12) at (2, 0) {$Y_{12}$}; + \node[sig, anchor=south] at (2, 0.5) {$u_{12}$}; + \draw[line] (0, 0) -- (Y12.west); + \draw[line] (Y12.east) -- (4, 0); +\end{scope} +\begin{scope}[#2] + \node[series] (Y23) at (6, 0) {$Y_{23}$}; + \node[sig, anchor=south] at (6, 0.5) {$u_{23}$}; + \draw[line] (4, 0) -- (Y23.west); + \draw[line] (Y23.east) -- (8, 0); +\end{scope} +\load{#3}{8} +\buses{1/0,2/4,3/8} +\end{tikzpicture} +} + +% Two buses: capacitive YG and a double circuit; #1 is the lower circuit style. +\newcommand{\resonance}[1]{ +\begin{tikzpicture}[emt diagram] +\source +\draw[line] (-1.3, 0) -- (-1.3, -1.38); +\draw[semithick] (-1.65, -1.38) -- (-0.95, -1.38); +\draw[semithick] (-1.65, -1.62) -- (-0.95, -1.62); +\draw[line] (-1.3, -1.62) -- (-1.3, -2.5); +\ground{(-1.3, -2.5)} +\node[sig, anchor=west] at (-0.85, -1.5) {$Y_G$}; +\inductor{0}{0.35} +\node[sig, anchor=south] at (2, 0.55) {$Y_{12}$}; +\begin{scope}[#1] + \inductor{0}{-0.35} + \node[sig, anchor=north] at (2, -0.45) {$u Y_{12}$}; +\end{scope} +\load{in service}{4} +\buses{1/0,2/4} +\end{tikzpicture} +} + +\begin{document} +\chain{in service}{in service}{in service} +\chain{in service}{out of service}{out of service} +\chain{out of service}{in service}{out of service} +\chain{out of service}{in service}{in service} +\resonance{in service} +\resonance{out of service} +\end{document} diff --git a/docs/Figures/PhasorDynamics/Branch/resonance-trip.png b/docs/Figures/PhasorDynamics/Branch/resonance-trip.png new file mode 100644 index 0000000000..3cfdca2087 Binary files /dev/null and b/docs/Figures/PhasorDynamics/Branch/resonance-trip.png differ diff --git a/docs/Figures/PhasorDynamics/Branch/resonance.png b/docs/Figures/PhasorDynamics/Branch/resonance.png new file mode 100644 index 0000000000..86afe95e29 Binary files /dev/null and b/docs/Figures/PhasorDynamics/Branch/resonance.png differ diff --git a/docs/Figures/PhasorDynamics/Branch/scenario1.png b/docs/Figures/PhasorDynamics/Branch/scenario1.png new file mode 100644 index 0000000000..ce1453f9db Binary files /dev/null and b/docs/Figures/PhasorDynamics/Branch/scenario1.png differ diff --git a/docs/Figures/PhasorDynamics/Branch/scenario2.png b/docs/Figures/PhasorDynamics/Branch/scenario2.png new file mode 100644 index 0000000000..f1b30a2e06 Binary files /dev/null and b/docs/Figures/PhasorDynamics/Branch/scenario2.png differ diff --git a/docs/Figures/PhasorDynamics/Branch/scenario3.png b/docs/Figures/PhasorDynamics/Branch/scenario3.png new file mode 100644 index 0000000000..776efbb420 Binary files /dev/null and b/docs/Figures/PhasorDynamics/Branch/scenario3.png differ diff --git a/examples/PhasorDynamics/DynamicSimulation/README.md b/examples/PhasorDynamics/DynamicSimulation/README.md index ef309ebf15..1fc1f2d986 100644 --- a/examples/PhasorDynamics/DynamicSimulation/README.md +++ b/examples/PhasorDynamics/DynamicSimulation/README.md @@ -5,5 +5,6 @@ inputs. | Example | Description | | --- | --- | +| [ThreeBusBasic](Toy/ThreeBusBasic/README.md) | A three-bus line opened and reclosed. | | [ThreeBusConstantSource](Toy/ThreeBusConstantSource/README.md) | A three-bus constant signal source example. | | [ACTIVSg10k](ACTIVSg10k/README.md) | A short simulation without disturbances using the reusable ACTIVSg10k case. | diff --git a/examples/PhasorDynamics/DynamicSimulation/Toy/CMakeLists.txt b/examples/PhasorDynamics/DynamicSimulation/Toy/CMakeLists.txt index b42f4e33a8..781579c495 100644 --- a/examples/PhasorDynamics/DynamicSimulation/Toy/CMakeLists.txt +++ b/examples/PhasorDynamics/DynamicSimulation/Toy/CMakeLists.txt @@ -1 +1,2 @@ +add_subdirectory(ThreeBusBasic) add_subdirectory(ThreeBusConstantSource) diff --git a/examples/PhasorDynamics/DynamicSimulation/Toy/ThreeBusBasic/CMakeLists.txt b/examples/PhasorDynamics/DynamicSimulation/Toy/ThreeBusBasic/CMakeLists.txt new file mode 100644 index 0000000000..f92d9f3b47 --- /dev/null +++ b/examples/PhasorDynamics/DynamicSimulation/Toy/ThreeBusBasic/CMakeLists.txt @@ -0,0 +1,7 @@ +gridkit_example_add_file(${CMAKE_SOURCE_DIR}/cases/PhasorDynamics/Toy/ThreeBusBasic.case.json) +gridkit_example_add_file(ThreeBusBasic.solver.json) + +add_test( + NAME ThreeBusBasic_line_switching + COMMAND DynamicSimulation ThreeBusBasic.solver.json + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/examples/PhasorDynamics/DynamicSimulation/Toy/ThreeBusBasic/README.md b/examples/PhasorDynamics/DynamicSimulation/Toy/ThreeBusBasic/README.md new file mode 100644 index 0000000000..9c6a3ed708 --- /dev/null +++ b/examples/PhasorDynamics/DynamicSimulation/Toy/ThreeBusBasic/README.md @@ -0,0 +1,8 @@ +# ThreeBusBasic Line Switching + +This study uses the reusable +[ThreeBusBasic case](../../../../../cases/PhasorDynamics/Toy/README.md). +Branch `BR_1_2` is taken out of service at 1 s and put back in service at 5 s. + +![Bus voltage magnitude](figures/ThreeBusBasic.Vm.png) +![Generator speed](figures/ThreeBusBasic.speed.png) diff --git a/examples/PhasorDynamics/DynamicSimulation/Toy/ThreeBusBasic/ThreeBusBasic.solver.json b/examples/PhasorDynamics/DynamicSimulation/Toy/ThreeBusBasic/ThreeBusBasic.solver.json new file mode 100644 index 0000000000..9bfa9704c5 --- /dev/null +++ b/examples/PhasorDynamics/DynamicSimulation/Toy/ThreeBusBasic/ThreeBusBasic.solver.json @@ -0,0 +1,10 @@ +{ + "system_model_file": "ThreeBusBasic.case.json", + "output_file": "ThreeBusBasic.csv", + "dt_monitor": 0.004166666666666667, + "tmax": 10.0, + "events": [ + { "time": 1.0, "type": "out_of_service", "device": "BR_1_2" }, + { "time": 5.0, "type": "in_service", "device": "BR_1_2" } + ] +} diff --git a/examples/PhasorDynamics/DynamicSimulation/Toy/ThreeBusBasic/figures/ThreeBusBasic.Vm.png b/examples/PhasorDynamics/DynamicSimulation/Toy/ThreeBusBasic/figures/ThreeBusBasic.Vm.png new file mode 100644 index 0000000000..554697d13f Binary files /dev/null and b/examples/PhasorDynamics/DynamicSimulation/Toy/ThreeBusBasic/figures/ThreeBusBasic.Vm.png differ diff --git a/examples/PhasorDynamics/DynamicSimulation/Toy/ThreeBusBasic/figures/ThreeBusBasic.speed.png b/examples/PhasorDynamics/DynamicSimulation/Toy/ThreeBusBasic/figures/ThreeBusBasic.speed.png new file mode 100644 index 0000000000..8a4946aa09 Binary files /dev/null and b/examples/PhasorDynamics/DynamicSimulation/Toy/ThreeBusBasic/figures/ThreeBusBasic.speed.png differ diff --git a/tests/IntegrationTests/PhasorDynamics/TenGen/Classical/TenGenClassical.cpp b/tests/IntegrationTests/PhasorDynamics/TenGen/Classical/TenGenClassical.cpp index 2058e3250c..4fcb118460 100644 --- a/tests/IntegrationTests/PhasorDynamics/TenGen/Classical/TenGenClassical.cpp +++ b/tests/IntegrationTests/PhasorDynamics/TenGen/Classical/TenGenClassical.cpp @@ -231,5 +231,20 @@ int main() error_set->display(); success *= error_set->total_error.max_value < 2.6e-4; + // Check that the generator-supported island remains solvable. + success *= branch56.setInService(false) == 0; + success *= ida.initializeSimulation(10.0) == 0; + success *= ida.runSimulation(10.1) == 0; + success *= sys.evaluateResidual() == 0; + // Check algebraic consistency at the output time; IDA controls differential solution error. + for (index_type i = 0; i < sys.size(); ++i) + { + success *= std::isfinite(sys.getResidual().getData()[i]); + if (!sys.tag()[i]) + { + success *= isEqual(sys.getResidual().getData()[i], real_type{0.0}, 1e-5); + } + } + return success.report("TenGenClassical"); } diff --git a/tests/UnitTests/PhasorDynamics/BranchTests.hpp b/tests/UnitTests/PhasorDynamics/BranchTests.hpp index ba09fd4e53..5dea2270af 100644 --- a/tests/UnitTests/PhasorDynamics/BranchTests.hpp +++ b/tests/UnitTests/PhasorDynamics/BranchTests.hpp @@ -1,8 +1,10 @@ #include #include #include +#include #include +#include #include #include #include @@ -243,6 +245,70 @@ namespace GridKit return success.report(__func__); } +#ifdef GRIDKIT_ENABLE_ENZYME + TestOutcome outOfServiceJacobian() + { + // Verify zero currents, fixed sparsity, and restoration after switching. + TestStatus success = true; + + PhasorDynamics::Bus bus1(10.0, 20.0); + PhasorDynamics::Bus bus2(30.0, 40.0); + bus1.allocate(); + bus2.allocate(); + bus1.initialize(); + bus2.initialize(); + for (IdxT i = 0; i < 2; ++i) + { + bus1.setVariableIndex(i, i); + bus1.setResidualIndex(i, i); + bus2.setVariableIndex(i, i + 2); + bus2.setResidualIndex(i, i + 2); + } + + PhasorDynamics::Branch branch(&bus1, &bus2, 2.0, 4.0, 0.2, 1.2); + branch.allocate(); + branch.evaluateJacobian(); + const IdxT nnz = branch.nnz(); + success *= nnz > 0; + auto* jacobian = branch.getCooJacobian(); + const std::vector rows(jacobian->getRowData(), jacobian->getRowData() + nnz); + const std::vector cols(jacobian->getColData(), jacobian->getColData() + nnz); + const std::vector values(jacobian->getValues(), jacobian->getValues() + nnz); + + success *= branch.setInService(false) == 0; + bus1.evaluateResidual(); + bus2.evaluateResidual(); + branch.evaluateResidual(); + success *= isEqual(bus1.Ir(), ScalarT{0.0}); + success *= isEqual(bus1.Ii(), ScalarT{0.0}); + success *= isEqual(bus2.Ir(), ScalarT{0.0}); + success *= isEqual(bus2.Ii(), ScalarT{0.0}); + branch.evaluateJacobian(); + success *= branch.nnz() == nnz; + + jacobian = branch.getCooJacobian(); + for (IdxT i = 0; i < nnz; ++i) + { + success *= jacobian->getRowData()[i] == rows[i]; + success *= jacobian->getColData()[i] == cols[i]; + success *= isEqual(jacobian->getValues()[i], RealT{0.0}); + } + + success *= branch.setInService(true) == 0; + branch.evaluateJacobian(); + success *= branch.nnz() == nnz; + jacobian = branch.getCooJacobian(); + for (IdxT i = 0; i < nnz; ++i) + { + success *= jacobian->getRowData()[i] == rows[i]; + success *= jacobian->getColData()[i] == cols[i]; + success *= isEqual(jacobian->getValues()[i], values[i]); + } + + return success.report(__func__); + } +#endif + TestOutcome parameterSetters() { // Verifies parameter setters refresh derived admittance values. diff --git a/tests/UnitTests/PhasorDynamics/SystemTests.hpp b/tests/UnitTests/PhasorDynamics/SystemTests.hpp index 0d1e8bdcca..e0fd40ae85 100644 --- a/tests/UnitTests/PhasorDynamics/SystemTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemTests.hpp @@ -221,6 +221,44 @@ namespace GridKit return success.report(__func__); } +#ifdef GRIDKIT_ENABLE_ENZYME + TestOutcome isolatedBusJacobian() + { + TestStatus success = true; + + PhasorDynamics::BusInfinite source(1.0, 0.0); + PhasorDynamics::Bus bus(1.0, 0.0); + PhasorDynamics::Branch branch(&source, &bus, 0.1, 0.2, 0.0, 0.0); + PhasorDynamics::SystemModel system; + system.addBus(&source); + system.addBus(&bus); + system.addComponent(&branch); + success *= system.allocate() == 0; + success *= system.initialize() == 0; + success *= system.evaluateJacobian() == 0; + + const IdxT nnz = system.getCsrJacobian()->getNnz(); + success *= system.size() == 2; + success *= nnz > 0; + + success *= branch.setInService(false) == 0; + success *= system.evaluateResidual() == 0; + success *= system.evaluateJacobian() == 0; + success *= isEqual(bus.Ir(), ScalarT{0.0}); + success *= isEqual(bus.Ii(), ScalarT{0.0}); + + // Both voltage variables remain, but their KCL rows are identically zero. + auto* jacobian = system.getCsrJacobian(); + success *= jacobian->getNnz() == nnz; + for (IdxT i = 0; i < nnz; ++i) + { + success *= isEqual(jacobian->getValues()[i], RealT{0.0}); + } + + return success.report(__func__); + } +#endif + TestOutcome reallocateAfterTopologyChange() { TestStatus success = true; diff --git a/tests/UnitTests/PhasorDynamics/runBranchTests.cpp b/tests/UnitTests/PhasorDynamics/runBranchTests.cpp index c0fe130c87..81d5aaf916 100644 --- a/tests/UnitTests/PhasorDynamics/runBranchTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runBranchTests.cpp @@ -16,6 +16,9 @@ int main() result += test.offNominalResidual(); result += test.jacobian(); result += test.offNominalJacobian(); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.outOfServiceJacobian(); +#endif return result.summary(); } diff --git a/tests/UnitTests/PhasorDynamics/runSystemTests.cpp b/tests/UnitTests/PhasorDynamics/runSystemTests.cpp index 71cbe85a36..6a52b1b279 100644 --- a/tests/UnitTests/PhasorDynamics/runSystemTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runSystemTests.cpp @@ -14,6 +14,7 @@ int main() result += test.modelVectorsAliasSystemStorage(); #ifdef GRIDKIT_ENABLE_ENZYME result += test.jacobian(); + result += test.isolatedBusJacobian(); #endif result += test.allocationError();