implemetnation of branch outage minimal

Author: lukelowry

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## v.03
 
+- 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 setEnabled(bool enabled);
+
       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 enabled_{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>::setEnabled(bool enabled)
+    {
+      enabled_ = ZERO<RealT>;
+      if (enabled)
+      {
+        enabled_ = 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,
@@ -514,6 +527,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_ *= enabled_;
+      b11_ *= enabled_;
+      g12_ *= enabled_;
+      b12_ *= enabled_;
+      g21_ *= enabled_;
+      b21_ *= enabled_;
+      g22_ *= enabled_;
+      b22_ *= enabled_;
     }
 
   } // namespace PhasorDynamics

GridKit/Model/PhasorDynamics/Branch/README.md

@@ -94,19 +94,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}
 ```
 

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>
@@ -25,6 +26,9 @@ namespace GridKit
     template <typename scalar_type, typename index_type>
     class BusFault;
 
+    template <typename scalar_type, typename index_type>
+    class Branch;
+
     template <typename scalar_type, typename index_type>
     class SignalNode;
 
@@ -105,6 +109,7 @@ namespace GridKit
       BusT*                    getBus(IdxT bus_id);
       SignalNodeT*             getSignalNode(IdxT signal_id);
       ComponentT*              getComponent(IdxT gridkit_component_id);
+      Branch<ScalarT, IdxT>*   getBranch(const std::string& id);
       BusFault<ScalarT, IdxT>* getBusFault(IdxT fault_id);
 
     private:
@@ -115,6 +120,8 @@ namespace GridKit
       std::map<IdxT, IdxT> gridkit_bus_indices_;   ///< Map between gridkit_bus_id and bus_id
       std::map<IdxT, IdxT> gridkit_fault_indices_; ///< Map between fault_id and component_id
 
+      std::map<std::string, IdxT> gridkit_branch_indices_; ///< Map between branch id and gridkit_component_id
+
       bool owns_components_{false};
 
       /// Variable monitor

GridKit/Model/PhasorDynamics/SystemModelImpl.hpp

@@ -110,6 +110,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);
       }
 
@@ -830,6 +832,17 @@ 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)]);
+    }
+
     /**
      * @brief Return pointer to a bus fault model
      *

application/PhasorDynamics/AnalysisUtilities.hpp

@@ -11,6 +11,7 @@
 #include <magic_enum/magic_enum.hpp>
 #include <nlohmann/json.hpp>
 
+#include <GridKit/Constants.hpp>
 #include <GridKit/Model/PhasorDynamics/SystemModelData.hpp>
 #include <GridKit/Solver/Dynamic/Ida.hpp>
 #include <GridKit/Testing/TestHelpers.hpp>
@@ -34,7 +35,9 @@ namespace GridKit
       enum class Type
       {
         FAULT_ON,
-        FAULT_OFF
+        FAULT_OFF,
+        ENABLE,
+        DISABLE
       };
 
       /// Time event takes place
@@ -43,6 +46,8 @@ namespace GridKit
       Type        type;
       /// ID of element used in event (e.g., bus fault id)
       std::size_t element_id;
+      /// Case-file `id` of the device to enable or disable
+      std::string device;
     };
 
     /**
@@ -130,7 +135,8 @@ namespace GridKit
       {
         auto& event = c.events.emplace_back();
         raw_event.at("time").get_to(event.time);
-        raw_event.at("element_id").get_to(event.element_id);
+        event.element_id = raw_event.value("element_id", 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;

application/PhasorDynamics/ContingencyAnalysis.cpp

@@ -7,6 +7,7 @@
 #include <omp.h>
 #endif
 
+#include <GridKit/Model/PhasorDynamics/Branch/Branch.hpp>
 #include <GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp>
 #include <GridKit/Model/PhasorDynamics/SystemModel.hpp>
 #include <GridKit/Solver/Dynamic/Ida.hpp>
@@ -64,6 +65,12 @@ TestStatus runStudy(StudyData study_data)
     case EventType::FAULT_OFF:
       sys.getBusFault(event.element_id)->setStatus(false);
       break;
+    case EventType::ENABLE:
+      sys.getBranch(event.device)->setEnabled(true);
+      break;
+    case EventType::DISABLE:
+      sys.getBranch(event.device)->setEnabled(false);
+      break;
     }
 
     // Re-initialize simulation at event time

application/PhasorDynamics/DynamicSimulation.cpp

@@ -2,6 +2,7 @@
 #include <filesystem>
 #include <fstream>
 
+#include <GridKit/Model/PhasorDynamics/Branch/Branch.hpp>
 #include <GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp>
 #include <GridKit/Model/PhasorDynamics/SystemModel.hpp>
 #include <GridKit/Solver/Dynamic/Ida.hpp>
@@ -63,6 +64,12 @@ int runApplication(int argc, const char* argv[])
     case EventType::FAULT_OFF:
       sys.getBusFault(event.element_id)->setStatus(false);
       break;
+    case EventType::ENABLE:
+      sys.getBranch(event.device)->setEnabled(true);
+      break;
+    case EventType::DISABLE:
+      sys.getBranch(event.device)->setEnabled(false);
+      break;
     }
 
     // Re-initialize simulation at event time

application/PhasorDynamics/README.md

@@ -29,5 +29,6 @@ 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", "enable", "disable" })
   `element_id`       | An integer value referencing the element associated with the event (e.g., bus fault id)
+  `device`           | `id` of the branch to put in service ("enable") or take out of service ("disable")