Branch switching events for applications

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
 

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};
 

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 <typename scalar_type, typename index_type>
+    int Branch<scalar_type, index_type>::setInService(bool in_service)
+    {
+      in_service_ = ZERO<RealT>;
+      if (in_service)
+      {
+        in_service_ = ONE<RealT>;
+      }
+      setDerivedParams();
+      return 0;
+    }
+
     template <typename scalar_type, typename index_type>
     __attribute__((always_inline)) inline void Branch<scalar_type, index_type>::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 <typename scalar_type, typename index_type>
     int Branch<scalar_type, index_type>::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

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,19 +104,23 @@ 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}
     +
     \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.

GridKit/Model/PhasorDynamics/SystemModel.hpp

@@ -2,6 +2,7 @@
 
 #include <map>
 #include <memory>
+#include <string>
 
 #include <GridKit/Model/PhasorDynamics/Component.hpp>
 #include <GridKit/Model/PhasorDynamics/SignalNode/SignalNodeSet.hpp>
@@ -22,6 +23,9 @@ namespace GridKit
     template <typename scalar_type, typename index_type>
     class BusBase;
 
+    template <typename scalar_type, typename index_type>
+    class Branch;
+
     template <typename scalar_type, typename index_type>
     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<ScalarT, IdxT>* getBranch(const std::string& id);
 
     private:
       std::vector<BusT*>       buses_;
@@ -109,6 +114,10 @@ namespace GridKit
 
       std::map<IdxT, IdxT> 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<std::string, IdxT> gridkit_branch_indices_;
+
       bool owns_components_{false};
 
       /// Variable monitor

GridKit/Model/PhasorDynamics/SystemModelImpl.hpp

@@ -98,6 +98,8 @@ namespace GridKit
         auto* branch = new Branch<ScalarT, IdxT>(getBus(bus1_index),
                                                  getBus(bus2_index),
                                                  branchdata);
+
+        gridkit_branch_indices_[branchdata.disambiguation_string] = static_cast<IdxT>(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 <typename scalar_type, typename index_type>
+    Branch<scalar_type, index_type>*
+    SystemModel<scalar_type, index_type>::getBranch(const std::string& id)
+    {
+      // Should fail if user-provided id is incorrect
+      return dynamic_cast<Branch<ScalarT, IdxT>*>(components_[gridkit_branch_indices_.at(id)]);
+    }
+
   } // namespace PhasorDynamics
 } // namespace GridKit

application/PhasorDynamics/AnalysisUtilities.hpp

@@ -29,15 +29,17 @@ 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
     {
       /// Type of event determines action performed
       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<std::size_t>);
+        event.bus    = raw_event.value("bus", INVALID_INDEX<std::size_t>);
+        event.device = raw_event.value("device", std::string{});
 
         auto type_str   = raw_event.at("type").get<std::string>();
         using EventType = SystemEvent::Type;