Experimental OPF Implementation using MATPOWER Input by lukelowry · Pull Request #594 · ORNL/GridKit

lukelowry · 2026-09-28T17:25:11Z

Description

Experimental @superwhiskers @nkoukpaizan which uses MATPOWER input which makes interfacing with PhasorDynamics harder but most OPF data comes in MATPOWER, so this is probably most practical choice

Application OptimalDispatch

Checklist

All tests pass.
Code compiles cleanly with flags -Wall -Wpedantic -Wconversion -Wextra.
The new code follows GridKit™ style guidelines.
There are unit tests for the new code.
The new code is documented.
The feature branch is rebased with respect to the target branch.
The CHANGELOG.md has been updated to reflect the changes. If this is a minor PR that is part of a larger fix already included in the file, state so.

Further comments

Needs some cleaning first

pelesh · 2026-09-28T18:22:04Z

A few comments:

If you have JSON input parser already available, I suggest you use that instead of MATPOWER format. I was thinking you could use some of the existing power flow implementation, but since you took a different route, I suggest you ditch MATPOWER and go straight to JSON.
I would have optimal power flow app in EconomicDispatch subdirectory in examples directory, at the same level as PhasorDynamics.
I would rename Optimization/OptimizationProblem.*pp files (and the class within) into Optimization/SteadyState.*pp. We will likely have Multiperiod and MixedInteger as new optimization modules in the future.

lukelowry · 2026-09-28T20:04:42Z

@pelesh I changed it so it uses the PowerFlow parser with very few changes. I do believe its valuable to have that as the input format only because it makes validation easier and there is practically unlimited data in .m, and OPF data does not always come with dynamic case files.

I believed I addressed your other issues too. I am able to use the example app to change the objective/model, solve for new operating point, and start simulation in steady state at that new point!

pelesh

This is very nice work. I suggest limiting the scope of this PR to providing ACOPF module only. Imho, the integration with dynamic simulations modules should be done after ACOPF is tested and input file format specified.

superwhiskers

tests pass on my end. i think it works quite well, but there are a few implementation questions we should sort out before this moves forward.

additionally, if we are to continue with matpower, are we intending to handle capability curves for generators or the convention of a dispatchable load as a generator with negative real-power injection? if they are outside the intended scope, could we reject them and document that limitation?

i think these conventions, particularly the latter, strengthen the argument for going with our own format, as using an existing format such as matpower requires us to adhere to its conventions, which may not be desirable.

nkoukpaizan

Overall nice work! The component implementations are clear and I like the consolidated ComponentModelImpl.hpp.

A few comments:

For Enzyme, can we re-use some of the existing boilerplate for LowerSparseStorage and use buffers instead of CooEntries wrapping std::vector)? Feel free to amend the wrappers.
Can we bind GridKit::LinearAlgebra::Vector at the component level instead of std::vector? This could be a follow-up PR.
Several methods, e.g., in OptimizationEvaluator* and SteadyState* don't have a Doxygen documentation.

nkoukpaizan · 2026-10-16T22:16:25Z

+      static constexpr RealT POWER_TOL    = 2 * std::numeric_limits<RealT>::epsilon();
+      static constexpr RealT JACOBIAN_TOL = 2 * std::numeric_limits<RealT>::epsilon();
+      static constexpr RealT HESSIAN_TOL  = 9 * std::numeric_limits<RealT>::epsilon();


I would round to 10 * and use the same tolerance.

nkoukpaizan · 2026-10-16T22:18:26Z

+if(GRIDKIT_ENABLE_ENZYME)
+  add_subdirectory(OptimalPowerFlow)
+endif()


I would move the guard into the OptimalPowerFlow directory.

nkoukpaizan · 2026-10-16T22:22:21Z

+      static constexpr RealT DERIVATIVE_TOL = 8 * std::numeric_limits<RealT>::epsilon();
+      static constexpr RealT BALANCE_TOL    = 1.0e-12;


Same here, I would round 8 to 10. And replace 1e-12 with a multiple of std::numeric_limits<RealT>::epsilon();

nkoukpaizan · 2026-10-16T22:25:01Z

+      bool isZero(double value)
+      {
+        return !(value < ZERO<double>) && !(value > ZERO<double>);
+      }


We've been trying to avoid these. If something needs == 0/!= 0, please highlight and we can evaluate.

nkoukpaizan · 2026-10-16T22:27:41Z

+          return true;
+        }
+        std::cerr << std::setprecision(std::numeric_limits<RealT>::max_digits10) << value << " differs from "
+                  << reference << " by " << std::abs(value - reference) / (ONE<RealT> + std::abs(reference)) / std::numeric_limits<RealT>::epsilon()


I wouldn't re-implement Testing::isEqual.

nkoukpaizan · 2026-10-16T22:30:48Z

+        const RealT* x = system.x().getData();
+
+        success *= isEqual(x[1], ZERO<RealT>, BALANCE_TOL);
+        success *= std::abs(x[3]) > BALANCE_TOL && std::abs(x[5]) > BALANCE_TOL;


Is checking the absolute values really needed instead of isEqual?

nkoukpaizan · 2026-10-16T22:48:43Z

    is >> br.fbus >> br.tbus >> br.r >> br.x >> br.b >> br.rateA >> br.rateB
        >> br.rateC >> br.ratio >> br.angle >> br.status >> br.angmin
        >> br.angmax;
-    checkEndOfMatrixRow(is);


Removing checkEndOfMatrixRow partially addresses #587, but we still need to ensure the read was successful (e.g., row not too short).

pelesh

This looks very good and has clean component model structure. Main issues that I would raise are:

Separating equality and inequality constraint evaluation methods.
Selection of coordinates polar vs. Cartesian.
Inheritance patterns should mirror PhasorDynamics module. If we want to change it, then the change should made be across the GridKit modules.

pelesh · 2026-10-16T23:55:04Z

+       * Filled by `mapped_store` and `lower_store`. Assembly sums duplicate
+       * coordinates.
+       */
+      struct CooEntries


What is the purpose of this struct? This needs to be documented.

pelesh · 2026-10-17T00:05:05Z

+      using ScalarT    = scalar_type;
+      using IdxT       = index_type;
+      using RealT      = typename ComponentModelT::RealT;
+      using ModelDataT = BranchData<RealT, IdxT>;


These should not be public. ScalarT, IdxT, etc. are local type aliases.

pelesh · 2026-10-17T00:07:51Z

+      g[0] = p1 * p1 + q1 * q1;
+      g[1] = p2 * p2 + q2 * q2;
+      g[2] = p1;
+      g[3] = q1;
+      g[4] = p2;
+      g[5] = q2;


Need to specify here what are equality and what are inequality constraints.

pelesh · 2026-10-17T00:17:31Z

+     * @tparam model_type - Derived component model
+     */
+    template <typename model_type, typename scalar_type, typename index_type>
+    class ComponentModel : public Component<scalar_type, index_type>
+    {
+    public:


I suggest removing Component class and using the same inheritance structure as in PhasorDynamics module.

pelesh · 2026-10-17T00:23:34Z

+      const auto& terminals = component.terminals();
+      for (IdxT t = 0; t < terminals.size(); ++t)
+      {


I believe that we decided to use term "ports" rather than "terminals". It would be good to use consistent nomenclature across all GridKit modules.

pelesh · 2026-10-17T00:32:11Z

+                                                               Number* x_u,
+                                                               Index   m,


These are deprecated Ipopt types. Use ipnumber and ipindex instead.

pelesh · 2026-10-17T00:36:34Z

+    template <typename scalar_type, typename index_type>
+    class SteadyState : public Ipopt::TNLP
+    {


This is not SteadyState but rather Ipopt specific interface. Perhaps use SteadyState for the namespace and IpoptInterface for the class name.

pelesh · 2026-10-17T00:48:55Z

+    public:
+      using ScalarT     = scalar_type;
+      using IdxT        = index_type;
+      using RealT       = typename ScalarTraits<ScalarT>::RealT;
+      using CooEntriesT = Enzyme::Sparse::CooEntries;


These type aliases should not be public.

pelesh · 2026-10-17T00:53:00Z

+      std::vector<std::unique_ptr<ComponentT>> components_;
+
+      /// Buses by number
+      std::map<IdxT, BusT*> buses_;
+
+      /// Map from component entries to CSR entries
+      std::vector<IdxT> jacobian_map_;
+      std::vector<IdxT> hessian_map_;


I would avoid using STL containers for system level objects. If it is a temporary solution, it should be marked as such.

pelesh · 2026-10-17T01:04:21Z

+      virtual int evaluateConstraints() = 0;
+      virtual int evaluateJacobian()    = 0;


It would be good to separate equality and inequality constraints, as well as their respective Jacobians.

lukelowry · 2026-10-18T13:40:07Z

Thank you all for detailed review, I believe I have addressed all except a distinct few. I will go response by response soon, but the overall design is overall more consistent with PhasorDynamics except for the inheritance structure. I am considering removing the larger ACTIVSg2000 case until after this is merged, please let me know what you think.

superwhiskers

some further changes since i looked a little closer at the ipopt interface code

superwhiskers · 2026-10-19T14:37:15Z

       */
      template <typename ScalarT, typename IdxT>
-      __attribute__((always_inline)) static ScalarT sparse_load(IdxT, IdxT, IdxT*, IdxT*, ScalarT*)
+      __attribute__((always_inline)) static ScalarT sparse_load(IdxT, IdxT, ScalarT, const IdxT*, const IdxT*, IdxT*, IdxT*, ScalarT*, IdxT&)


should still be marked with [[gnu::always_inline]]

superwhiskers · 2026-10-19T14:54:48Z

+      /// Bus `id`, from its number
+      std::string id() const
+      {
+        return std::format("bus {}", number);
+      }


this is a bit nitpicky but if we're just using this to pretty-print the bus, it would be better to either 1) do it at the site of pretty-printing or 2) overload << against an iostream / specialize std::formatter

superwhiskers · 2026-10-19T18:14:56Z

+    /// Every problem with the data
+    template <typename scalar_type, typename index_type>
+    Model::ConfigurationChecks SystemModel<scalar_type, index_type>::verify() const
+    {
+      Model::ConfigurationChecks checks;
+
+      std::unordered_set<const ComponentT*> components;
+      for (const auto& component : components_)
+      {
+        components.insert(component.get());
+      }
+
+      for (const auto& [number, bus] : buses_)
+      {
+        checks.check(number != INVALID_INDEX<IdxT>,
+                     std::format("{}: invalid bus number", bus->id()));
+      }
+
+      std::set<std::string> ids;
+      for (const auto& component : components_)
+      {
+        checks.check(ids.insert(component->id()).second, std::format("{}: duplicate id", component->id()));
+
+        const auto component_checks = component->verify();
+        for (const auto& error : component_checks.errors())
+        {
+          checks.fail(std::format("{}: {}", component->id(), error));
+        }
+
+        const auto& ports = component->ports();
+        for (IdxT port = 0; port < ports.size(); ++port)
+        {
+          checks.check(components.contains(ports[port].component),
+                       std::format("{}: port {} must connect to a component in this system", component->id(), port));
+        }
+      }
+
+      return checks;
+    }


the way this is laid out makes me think we could be checking some of these things at construction time and throwing an exception instead of looping over the entire system after constructing it, but it may not be right to deviate from PhasorDynamics that way in this pull request

this would have the benefit of guaranteeing that any SystemModel is valid, making validity an invariant of the SystemModel structure as there is no way that i see to add more buses, components onto an existing SystemModel

superwhiskers · 2026-10-19T18:30:27Z

+/**
+ * @file IpoptInterface.cpp
+ * @brief Ipopt problem over a `Model::OptimizationEvaluator`.
+ */


given we already have an ipopt interface under GridKit/Solver/Optimization this generic name is confusing as it isn't immediately obvious how this is different (it uses the OptimizationEvaluator interface)

Yes, we need to sort out the naming so that we distinguish between dynamic and steady state optimization interfaces. Note that Model::OptimizationEvaluator is the model interface to optimization solvers and can in principle work with different solvers. It is different from IpoptInterface, which contains code specific to Ipopt.

superwhiskers · 2026-10-19T18:35:28Z

+    }
+
+    // Available template instantiations
+    template class IpoptInterface<double, size_t>;


Suggested change

template class IpoptInterface<double, size_t>;

template class IpoptInterface<Ipopt::Number, Ipopt::Index>;

i'm going to guess that since we're interfacing with ipopt here, it would be more proper to instantiate based upon the types it is compatible with, as is done in the existing interface code

It would in general but that may cause type mismatch elsewhere in the code. This is one of poorly addressed issues in GridKit. To make everything consistent, we would need to build all dependencies with the same types and to type check at compile time verifying types are matching.

In this case I would set

template class IpoptInterface<RealT, IdxT>;

and check that types match Ipopt::Number and Ipopt::Index.

lukelowry requested review from nkoukpaizan and superwhiskers September 28, 2026 17:25

pelesh self-requested a review September 28, 2026 18:13

pelesh requested a review from abirchfield September 28, 2026 18:22

lukelowry force-pushed the lukel/OPF-NEW-dev branch 4 times, most recently from 95913b8 to 5c4fb92 Compare October 10, 2026 21:02

lukelowry force-pushed the lukel/OPF-NEW-dev branch from 5c4fb92 to 49ad14b Compare October 13, 2026 22:04

pelesh reviewed Oct 14, 2026

View reviewed changes

lukelowry added 9 commits October 14, 2026 23:54

Initial OPF implementation

860af60

Minimize Churn on PhasorDynamics

81ba74e

clenaup matpower for both models

dd18dcd

rename to EconomicDispatch and downgrade to an example app

51996ad

ipopt cmake settings

2b40e6b

remove integration with PhasorDynamics Apps

583d1af

cleanup uncessary tests

d428593

BUG: Application symlink issue

b74ff99

Enhance and cleanup

a6e3ae1

lukelowry force-pushed the lukel/OPF-NEW-dev branch from 979e671 to a6e3ae1 Compare October 15, 2026 07:09

platofrm dependent correction

d5c7d68

superwhiskers requested changes Oct 16, 2026

View reviewed changes

nkoukpaizan reviewed Oct 16, 2026

View reviewed changes

nkoukpaizan assigned lukelowry Oct 16, 2026

nkoukpaizan added enhancement New feature or request new model new solver labels Oct 16, 2026

pelesh reviewed Oct 17, 2026

View reviewed changes

lukelowry added 5 commits October 18, 2026 02:42

Read OPF cases from JSON and rebuild components on ports

680fdde

Adopt the standard model interface and start from flat

8e91d10

Tag inequality constraints, store multipliers, and simplify parsing

f679ffd

Remove island checks

fe85267

distill documentaiton

5bc7404

lukelowry force-pushed the lukel/OPF-NEW-dev branch from 38dc4b5 to 5bc7404 Compare October 18, 2026 07:54

remove edge case

657ac9f

lukelowry requested a review from pelesh October 18, 2026 13:40

superwhiskers requested changes Oct 19, 2026

View reviewed changes

adjust inheritance and minor type issues

9ebec1b

		static constexpr RealT DERIVATIVE_TOL = 8 * std::numeric_limits<RealT>::epsilon();
		static constexpr RealT BALANCE_TOL = 1.0e-12;

		virtual int evaluateConstraints() = 0;
		virtual int evaluateJacobian() = 0;

	template class IpoptInterface<double, size_t>;
	template class IpoptInterface<Ipopt::Number, Ipopt::Index>;

Conversation

lukelowry commented Sep 28, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Description

Checklist

Further comments

Uh oh!

pelesh commented Sep 28, 2026

Uh oh!

lukelowry commented Sep 28, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

pelesh left a comment

Choose a reason for hiding this comment

Uh oh!

superwhiskers left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

nkoukpaizan left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pelesh left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

lukelowry commented Oct 18, 2026

Uh oh!

superwhiskers left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

lukelowry commented Sep 28, 2026 •

edited

Loading

lukelowry commented Sep 28, 2026 •

edited

Loading