expose values for gens and buses for later interop

Author: lukelowry

GridKit/Model/OptimalPowerFlow/Bus/Bus.hpp

@@ -75,6 +75,18 @@ namespace GridKit
 
       [[gnu::always_inline]] void internalConstraints(const ScalarT* x, const ScalarT* external, ScalarT* g) const override;
 
+      /// Real voltage [p.u.]
+      const ScalarT& Vr() const
+      {
+        return x_.getData()[static_cast<size_t>(BusInternalVariables::VR)];
+      }
+
+      /// Imaginary voltage [p.u.]
+      const ScalarT& Vi() const
+      {
+        return x_.getData()[static_cast<size_t>(BusInternalVariables::VI)];
+      }
+
     private:
       void initializeParameters(const ModelDataT& data);
 

GridKit/Model/OptimalPowerFlow/Generator/Generator.hpp

@@ -44,6 +44,7 @@ namespace GridKit
       using ModelDataT = GeneratorData<RealT, IdxT>;
       using BusT       = Bus<scalar_type, index_type>;
 
+      using ComponentT::x_;
       using ComponentT::x_lower_;
       using ComponentT::x_upper_;
 
@@ -69,6 +70,18 @@ namespace GridKit
       [[gnu::always_inline]] ScalarT objective(const ScalarT* x, const ScalarT* external) const override;
       [[gnu::always_inline]] void    externalConstraints(const ScalarT* x, const ScalarT* external, ScalarT* g) const override;
 
+      /// Active power into the bus [p.u.]
+      const ScalarT& P() const
+      {
+        return x_.getData()[static_cast<size_t>(GeneratorInternalVariables::P)];
+      }
+
+      /// Reactive power into the bus [p.u.]
+      const ScalarT& Q() const
+      {
+        return x_.getData()[static_cast<size_t>(GeneratorInternalVariables::Q)];
+      }
+
     private:
       void initializeParameters(const ModelDataT& data);
       void setDerivedParameters();

GridKit/Model/OptimalPowerFlow/SystemModel.hpp

@@ -7,20 +7,20 @@
 
 #include <map>
 #include <memory>
+#include <string>
 #include <vector>
 
 #include <GridKit/Model/ConfigurationChecks.hpp>
+#include <GridKit/Model/OptimalPowerFlow/Bus/Bus.hpp>
 #include <GridKit/Model/OptimalPowerFlow/Component.hpp>
+#include <GridKit/Model/OptimalPowerFlow/Generator/Generator.hpp>
 #include <GridKit/Model/OptimalPowerFlow/SystemModelData.hpp>
 #include <GridKit/Model/OptimizationEvaluator.hpp>
 
 namespace GridKit
 {
   namespace OptimalPowerFlow
   {
-    template <typename scalar_type, typename index_type>
-    class Bus;
-
     /**
      * @brief AC optimal power flow in Cartesian coordinates
      *
@@ -38,6 +38,7 @@ namespace GridKit
       using IndexVectorT = LinearAlgebra::Vector<IdxT, IdxT>;
       using CsrMatrixT   = typename EvaluatorT::CsrMatrixT;
       using BusT         = Bus<ScalarT, IdxT>;
+      using GeneratorT   = Generator<ScalarT, IdxT>;
       using Block        = typename ComponentT::Block;
       using CooBufferT   = typename ComponentT::CooBuffer;
 
@@ -72,6 +73,12 @@ namespace GridKit
       int evaluateJacobian() override;
       int evaluateHessian(RealT sigma) override;
 
+      /// Bus by number, which views its solution after `allocate()`
+      BusT* getBus(IdxT number);
+
+      /// Generator by case-file `id`, which views its solution after `allocate()`
+      GeneratorT* getGenerator(const std::string& id);
+
     private:
       /// Component entries of one derivative matrix and the CSR entry of each
       struct Assembly
@@ -83,7 +90,7 @@ namespace GridKit
       void addBus(const BusData<RealT, IdxT>& data);
 
       template <typename ModelT, typename DataT>
-      void addDevice(const DataT& data);
+      ModelT* addDevice(const DataT& data);
 
       int compose();
       int bind(ComponentT& component, const Block& offset);
@@ -98,6 +105,9 @@ namespace GridKit
       /// Buses by number
       std::map<IdxT, BusT*> buses_;
 
+      /// Generators by id
+      std::map<std::string, GeneratorT*> generators_;
+
       Assembly jacobian_assembly_;
       Assembly hessian_assembly_;
     };

GridKit/Model/OptimalPowerFlow/SystemModelImpl.hpp

@@ -45,7 +45,7 @@ namespace GridKit
 
       for (const auto& generator : data.generator)
       {
-        addDevice<Generator<ScalarT, IdxT>>(generator);
+        generators_.emplace(generator.id, addDevice<GeneratorT>(generator));
       }
 
       for (const auto& load : data.load)
@@ -234,6 +234,21 @@ namespace GridKit
       return 0;
     }
 
+    /// @throw std::out_of_range for an unknown bus
+    template <typename scalar_type, typename index_type>
+    typename SystemModel<scalar_type, index_type>::BusT* SystemModel<scalar_type, index_type>::getBus(IdxT number)
+    {
+      return buses_.at(number);
+    }
+
+    /// @throw std::out_of_range for an unknown generator
+    template <typename scalar_type, typename index_type>
+    typename SystemModel<scalar_type, index_type>::GeneratorT*
+    SystemModel<scalar_type, index_type>::getGenerator(const std::string& id)
+    {
+      return generators_.at(id);
+    }
+
     template <typename scalar_type, typename index_type>
     void SystemModel<scalar_type, index_type>::addBus(const BusData<RealT, IdxT>& data)
     {
@@ -248,7 +263,7 @@ namespace GridKit
     /// Device from `data`, connected at each port to the bus it names
     template <typename scalar_type, typename index_type>
     template <typename ModelT, typename DataT>
-    void SystemModel<scalar_type, index_type>::addDevice(const DataT& data)
+    ModelT* SystemModel<scalar_type, index_type>::addDevice(const DataT& data)
     {
       using Ports = typename DataT::Ports;
       static_assert(ModelT::NUM_PORT_VARIABLES.size() == Utilities::enum_size<Ports>());
@@ -276,7 +291,9 @@ namespace GridKit
         }
       }
 
+      ModelT* const added = device.get();
       components_.push_back(std::move(device));
+      return added;
     }
 
     /// Size the system vectors and bind each component to its blocks

GridKit/Solver/Optimization/IpoptInterface.cpp

@@ -207,6 +207,7 @@ namespace AnalysisManager
       copyInto(model_->zLower(), z_L);
       copyInto(model_->zUpper(), z_U);
       model_->evaluateObjective();
+      model_->evaluateConstraints();
     }
 
     template <typename scalar_type, typename index_type>

GridKit/Solver/Optimization/IpoptInterface.hpp

@@ -91,7 +91,7 @@ namespace AnalysisManager
                   ipindex*        jCol,
                   ipnumber*       values) override;
 
-      /// Leave the final point, its multipliers, and its objective in the model
+      /// Leave the final point, its multipliers, its objective, and its constraints in the model
       void finalize_solution(SolverReturn               status,
                              ipindex                    n,
                              const ipnumber*            x,

tests/UnitTests/OptimalPowerFlow/IpoptTests.hpp

@@ -22,7 +22,8 @@ namespace GridKit
       using SystemT = OPF::SystemModel<RealT, IdxT>;
 
     public:
-      // Ipopt stops with the constraints and the stationarity within its tolerance. They reach 5e-13 and 2e-13.
+      // Ipopt stops with the constraints and the stationarity within its tolerance. They reach 5e-13 and 2e-13,
+      // and the marginal prices 8e-13.
       static constexpr RealT kTol = static_cast<RealT>(1.0e-10);
 
       TestOutcome reference()
@@ -41,10 +42,9 @@ namespace GridKit
           return success.report(__func__);
         }
 
-        const RealT* x = system.x().getData();
-
-        success *= isEqual(x[1], ZERO<RealT>, kTol);
-        success *= !isEqual(x[3], ZERO<RealT>, kTol) && !isEqual(x[5], ZERO<RealT>, kTol);
+        success *= isEqual(system.getBus(1)->Vi(), ZERO<RealT>, kTol);
+        success *= !isEqual(system.getBus(2)->Vi(), ZERO<RealT>, kTol);
+        success *= !isEqual(system.getBus(3)->Vi(), ZERO<RealT>, kTol);
 
         return success.report(__func__);
       }
@@ -65,22 +65,50 @@ namespace GridKit
           return success.report(__func__);
         }
 
-        success *= system.evaluateConstraints() == 0;
-        if (!success)
-        {
-          return success.report(__func__);
-        }
         const RealT* g = system.g().getData();
         for (const IdxT i : SystemTestData::BALANCE_ROWS)
         {
           success *= isEqual(g[i], ZERO<RealT>, kTol);
         }
         success *= isEqual(g[SystemTestData::FLOW_ROWS[0]], SystemTestData::SMAX * SystemTestData::SMAX, kTol);
 
-        const RealT* x = system.x().getData();
+        const RealT p1 = system.getGenerator("gen_1")->P();
+        const RealT p3 = system.getGenerator("gen_3")->P();
+
+        success *= p1 > p3;
+        success *= p3 > ZERO<RealT>;
 
-        success *= x[6] > x[8];
-        success *= x[8] > ZERO<RealT>;
+        return success.report(__func__);
+      }
+
+      /// The negative of the active power multiplier of a bus is the marginal
+      /// cost of a generator there within its limits
+      TestOutcome marginalPrices()
+      {
+        TestStatus success = true;
+
+        SystemT system(SystemTestData::data());
+        success *= system.allocate() == 0 && system.initialize() == 0;
+        if (!success)
+        {
+          return success.report(__func__);
+        }
+        success *= optimize(system);
+        if (!success)
+        {
+          return success.report(__func__);
+        }
+
+        const auto price = [&system](IdxT number)
+        {
+          const auto P = static_cast<size_t>(OPF::BusInternalConstraints::P);
+          return -system.getBus(number)->lambda().getData()[P];
+        };
+        const auto marginalCost = [&system](const char* id, RealT c1)
+        { return c1 + TWO<RealT> * SystemTestData::C2 * system.getGenerator(id)->P(); };
+
+        success *= isEqual(price(1), marginalCost("gen_1", SystemTestData::C1_GEN_1), kTol);
+        success *= isEqual(price(3), marginalCost("gen_3", SystemTestData::C1_GEN_3), kTol);
 
         return success.report(__func__);
       }

tests/UnitTests/OptimalPowerFlow/SystemTestData.hpp

@@ -20,6 +20,11 @@ namespace GridKit
       static constexpr RealT LOAD_P = 1.2;
       static constexpr RealT LOAD_Q = 0.3;
 
+      /// Linear costs of gen_1 and gen_3, and the quadratic cost of both
+      static constexpr RealT C1_GEN_1 = 10.0;
+      static constexpr RealT C1_GEN_3 = 20.0;
+      static constexpr RealT C2       = 0.01;
+
       /// Constraints of the buses in order, P, Q, and VM2 each, then S1 and S2
       /// of the rated branch, then of the unrated branches
       static constexpr std::array<IdxT, 6> BALANCE_ROWS{0, 1, 3, 4, 6, 7};
@@ -49,8 +54,8 @@ namespace GridKit
         addBranch(data, "branch_1_3", 1, 3, 0.02, 0.2, 0.0);
         data.branch[0].parameters[OPF::BranchParameters::Smax] = SMAX;
 
-        addGenerator(data, "gen_1", 1, 10.0);
-        addGenerator(data, "gen_3", 3, 20.0);
+        addGenerator(data, "gen_1", 1, C1_GEN_1);
+        addGenerator(data, "gen_3", 3, C1_GEN_3);
 
         auto& load                              = data.load.emplace_back();
         load.id                                 = "load_2";
@@ -88,7 +93,7 @@ namespace GridKit
         generator.parameters[OPF::GeneratorParameters::Qmin] = -1.0;
         generator.parameters[OPF::GeneratorParameters::Qmax] = 1.0;
         generator.parameters[OPF::GeneratorParameters::c1]   = c1;
-        generator.parameters[OPF::GeneratorParameters::c2]   = 0.01;
+        generator.parameters[OPF::GeneratorParameters::c2]   = C2;
       }
     };
   } // namespace Testing

tests/UnitTests/OptimalPowerFlow/runIpoptTests.cpp

@@ -8,6 +8,7 @@ int main()
   result += test.reference();
   result += test.solve();
   result += test.multipliers();
+  result += test.marginalPrices();
 
   return result.summary();
 }