grid-opt-alg-lab · samtalki · Jun 15, 2026 · Mar 31, 2026 · Apr 24, 2026 · Apr 24, 2026
diff --git a/.github/workflows/Benchmark.yml b/.github/workflows/Benchmark.yml
@@ -21,3 +21,4 @@ jobs:
           julia-version: '1'
           mode: 'time,memory'
           script: benchmark/benchmarks.jl
+          extra-pkgs: PowerModels
diff --git a/Artifacts.toml b/Artifacts.toml
@@ -0,0 +1,7 @@
+[PGLib_opf]
+git-tree-sha1 = "0e8968a89b6ad43910a8eda4ec30656add35cf91"
+lazy = true
+
+    [[PGLib_opf.download]]
+    sha256 = "f1421ce22f0a7b9de8a8b2111776b496348220192ad24aace392c3bf608706c2"
+    url = "https://github.com/power-grid-lib/pglib-opf/archive/refs/tags/v23.07.tar.gz"
diff --git a/Project.toml b/Project.toml
@@ -4,10 +4,13 @@ version = "0.1.0"
 authors = ["Samuel Talkington", "Michael Klamkin", "Cameron Khanpour"]
 
 [deps]
+ExaModels = "1037b233-b668-4ce9-9b63-f9f681f55dd2"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655"
+NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
+PowerIO = "05ed8b54-f668-4096-9d0d-e8c3dd9dc169"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [weakdeps]
@@ -18,15 +21,19 @@ PowerDiffAPFExt = "AcceleratedDCPowerFlows"
 
 [compat]
 AcceleratedDCPowerFlows = "0.1"
+ExaModels = "0.9, 0.10"
 Ipopt = "1"
 JuMP = "1"
-PowerModels = "0.21"
+LazyArtifacts = "1"
+NLPModelsIpopt = "0.11.2"
+PowerIO = "0.1.3"
 julia = "1.9"
 
 [extras]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ForwardDiff", "Statistics", "Test"]
+test = ["ForwardDiff", "PowerModels", "Statistics", "Test"]
diff --git a/README.md b/README.md
@@ -31,9 +31,9 @@ Pkg.add(url="https://github.com/grid-opt-alg-lab/PowerDiff.jl.git")
 ## Quick Start
 
 ```julia
-using PowerDiff, PowerModels
+using PowerDiff
 
-# Load network (make_basic_network is optional)
+# Parse a supported PowerIO case into a PowerIO.Network
 net = parse_file("case14.m")
 dc_net = DCNetwork(net)
 d = calc_demand_vector(net)
@@ -59,12 +59,18 @@ See the [Getting Started guide](https://samueltalkington.com/research/powerdiff/
 - [Advanced Topics](https://samueltalkington.com/research/powerdiff/advanced/) — Type hierarchy, caching, solver configuration
 - [API Reference](https://samueltalkington.com/research/powerdiff/api/) — Full docstring reference
 
+## Input Format
+
+PowerDiff reads files through PowerIO. `parse_file` supports MATPOWER `.m`,
+PSS/E `.raw`, PowerWorld `.aux`, PowerModels JSON, and Egret JSON. For streams,
+pass `from`; JSON streams need `from=:egret` or `from=:powermodels`.
+
 ## Dependencies
 
-- [PowerModels.jl](https://github.com/lanl-ansi/PowerModels.jl) — Power system modeling
+- [PowerIO.jl](https://github.com/eigenergy/PowerIO.jl) — Parser and data layer (see `docs/powerio-integration.md`)
 - [JuMP.jl](https://github.com/jump-dev/JuMP.jl) — Optimization modeling
+- [ExaModels.jl](https://github.com/exanauts/ExaModels.jl) — Alternative optimization modeling for GPU parallelization
 - [Ipopt.jl](https://github.com/jump-dev/Ipopt.jl) — Default solver for DC and AC OPF
-- [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) — Automatic differentiation
 
 ## License
 

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
@@ -1,29 +1,40 @@
 using BenchmarkTools
 using PowerDiff
+import PowerModels
 
-const PM = PowerDiff.PM
+const PM = PowerModels
 const SUITE = BenchmarkGroup()
 
 PM.silence()
 PowerDiff.silence()
 
+function _parse_benchmark_case(case_path)
+    if isdefined(PowerDiff, :parse_file)
+        return PowerDiff.parse_file(case_path)
+    end
+    return PM.make_basic_network(PM.parse_file(case_path))
+end
+
 function _load_benchmark_case()
     pm_dir = joinpath(dirname(pathof(PM)), "..", "test", "data", "matpower")
     for case_name in ("case30.m", "case24.m", "case14.m", "case9.m", "case5.m")
         case_path = joinpath(pm_dir, case_name)
         isfile(case_path) || continue
-        raw = PM.parse_file(case_path)
-        return case_name, PM.make_basic_network(raw)
+        net_data = _parse_benchmark_case(case_path)
+        return case_name, case_path, net_data
     end
     error("No bundled PowerModels MATPOWER benchmark case found")
 end
 
-case_name, net_data = _load_benchmark_case()
+case_name, case_path, net_data = _load_benchmark_case()
 prob = DCOPFProblem(net_data)
 sol = solve!(prob)
 ac_prob = ACOPFProblem(deepcopy(net_data); silent=true)
 ac_sol = solve!(ac_prob)
 
+SUITE["parser"] = BenchmarkGroup()
+SUITE["parser"][case_name] = @benchmarkable _parse_benchmark_case($case_path)
+
 SUITE["dc_opf"] = BenchmarkGroup()
 SUITE["dc_opf"]["kkt_jacobian"] = BenchmarkGroup()
 kkt_suite = SUITE["dc_opf"]["kkt_jacobian"][case_name] = BenchmarkGroup()

diff --git a/docs/Project.toml b/docs/Project.toml
@@ -3,7 +3,6 @@
 
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655"
 PowerDiff = "4fa8226c-b122-4e48-8217-6f318ba8ef74"
 
 [sources]

diff --git a/docs/powerio-integration.md b/docs/powerio-integration.md
@@ -0,0 +1,29 @@
+# PowerIO Parser Contract
+
+PowerIO is PowerDiff's parser and data layer. PowerDiff does not expose a parser
+backend switch.
+
+`PowerDiff.parse_file(path)` resolves the path and returns a `PowerIO.Network` via
+`PowerIO.parse_file`. PowerIO infers path formats from extensions unless `from` is
+given. `PowerDiff.parse_file(io)` uses MATPOWER by default because streams have no
+extension; pass `from` for PSS/E RAW, PowerWorld AUX, PowerModels JSON, or Egret
+JSON. JSON streams are ambiguous, so use `from=:egret` or `from=:powermodels`.
+Pass the result to [`DCNetwork`](@ref) or [`ACNetwork`](@ref).
+
+The network constructors build directly from `PowerIO.to_powerdata(net)`, which
+already returns normalized data: per-unit scaling by `base_mva`, degree-to-radian
+conversion, out-of-service and isolated-element filtering, bus-type inference,
+per-bus load/shunt aggregation, and polynomial cost rescaling. PowerDiff layers on
+only the OPF modeling it owns:
+
+- polynomial cost interpretation: it reads the constant, linear, and quadratic
+  coefficients straight from `to_powerdata`'s generator rows (already per-unit and
+  right-aligned). PWL costs are rejected; higher-order polynomials are rejected by
+  `to_powerdata` itself. A generator with no cost record is treated as cost-free.
+- a finite `rate_a` fallback when the source leaves the thermal limit at `0`
+- default angle-difference bounds
+
+PowerDiff rejects networks carrying storage or HVDC/dcline records, which it does
+not model.
+
+The parser tests assert path and IO parity through this single PowerIO path.
diff --git a/docs/src/advanced.md b/docs/src/advanced.md
@@ -40,9 +40,9 @@ Stores the DC network topology and parameters.
 | `ref_bus` | `Int` | Reference bus index (sequential) |
 | `tau` | `Float64` | Regularization parameter |
 | `id_map` | `IDMapping` | Bidirectional element ID mapping (original ↔ sequential) |
-| `ref` | `Union{Nothing,Dict}` | Stored `build_ref` result (nothing for programmatic constructors) |
 
-Construct from PowerModels data: `DCNetwork(net)` (accepts both basic and non-basic networks) or with explicit parameters: `DCNetwork(n, m, k, A, G_inc, b; ...)`.
+Construct from a parsed MATPOWER network with `DCNetwork(parse_file("case14.m"))`, or
+with explicit parameters: `DCNetwork(n, m, k, A, G_inc, b; ...)`.
 
 ### ACNetwork
 
@@ -59,9 +59,7 @@ Stores the AC network with vectorized admittance representation.
 | `is_switchable` | `BitVector` | Which branches can be switched |
 | `idx_slack` | `Int` | Slack bus index (sequential) |
 | `vm_min`, `vm_max` | `Vector{Float64}` | Voltage magnitude limits per bus |
-| `i_max` | `Vector{Float64}` | Branch current magnitude limits |
 | `id_map` | `IDMapping` | Bidirectional element ID mapping (original ↔ sequential) |
-| `ref` | `Union{Nothing,Dict}` | Stored `build_ref` result (nothing for Y-matrix constructors) |
 
 ## Sensitivity Caching
 
@@ -100,10 +98,13 @@ prob = DCOPFProblem(dc_net, d; optimizer=HiGHS.Optimizer)
 
 ### AC OPF
 
-Default solver is Ipopt. The `silent` keyword suppresses solver output:
+The default `:jump` backend uses Ipopt. The opt-in CPU `:exa` backend uses
+ExaModels and NLPModelsIpopt. Custom JuMP optimizer objects are accepted only
+by `:jump`.
 
 ```julia
 prob = ACOPFProblem(net; silent=true)
+exa_prob = ACOPFProblem(net; backend=:exa, silent=true)
 ```
 
 ## KKT System Access (Qualified)
@@ -124,7 +125,7 @@ idx = PD.kkt_indices(dc_net)           # Named index ranges
 
 # AC OPF — same unified API
 z = PD.flatten_variables(sol, ac_prob)
-J = PD.calc_kkt_jacobian(ac_prob)       # Dense Jacobian via ForwardDiff
+J = PD.calc_kkt_jacobian(ac_prob)       # Sparse analytical Jacobian
 dim = PD.kkt_dims(ac_prob)             # KKT dimension
 idx = PD.kkt_indices(ac_prob)          # Named index ranges
 ```

diff --git a/docs/src/api.md b/docs/src/api.md
@@ -1,5 +1,14 @@
 # API Reference
 
+## Parser
+
+```@docs
+parse_file
+parse_matpower
+parse_matpower_struct
+get_path
+```
+
 ## Sensitivity Interface
 
 ```@docs

diff --git a/docs/src/getting-started.md b/docs/src/getting-started.md
@@ -6,13 +6,16 @@ This guide walks through the main workflows: DC power flow, DC OPF with LMP anal
 
 ```julia
 using PowerDiff
-using PowerModels
 
-# Load a MATPOWER case (make_basic_network is optional)
-net = PowerModels.parse_file("case14.m")
-# OR: net = PowerModels.make_basic_network(raw)  # sequential IDs
+# Parse a supported PowerIO case into a PowerIO.Network
+net = parse_file("case14.m")
 ```
 
+PowerDiff reads files through PowerIO. `parse_file` supports MATPOWER `.m`,
+PSS/E `.raw`, PowerWorld `.aux`, PowerModels JSON, and Egret JSON. For streams,
+pass `from`; JSON streams need `from=:egret` or `from=:powermodels`. The former
+PowerModels dictionary constructors were removed.
+
 ## Interactive Exploration
 
 For a pre-loaded REPL session with case14, run:
@@ -109,11 +112,13 @@ sens * ones(size(sens,2)) # matrix-vector product
 
 ## AC Power Flow
 
-AC power flow sensitivities require a solved AC power flow solution.
+AC power flow sensitivities require complex bus voltages from an external AC
+power flow solve.
 
 ```julia
-PowerModels.compute_ac_pf!(net)
-ac_state = ACPowerFlowState(net)
+ac_net = ACNetwork(net)
+v = external_solver_voltage_vector
+ac_state = ACPowerFlowState(ac_net, v)
 
 # Voltage and current sensitivities
 dvm_dp = calc_sensitivity(ac_state, :vm, :p)   # d|V|/dp (n x n)

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -23,9 +23,9 @@ Pkg.add(url="https://github.com/grid-opt-alg-lab/PowerDiff.jl.git")
 ## Quick Example
 
 ```julia
-using PowerDiff, PowerModels
+using PowerDiff
 
-# Load network (make_basic_network is optional)
+# Parse a supported PowerIO case into a PowerIO.Network
 net = parse_file("case14.m")
 dc_net = DCNetwork(net)
 d = calc_demand_vector(net)

diff --git a/docs/src/math/ac-opf.md b/docs/src/math/ac-opf.md
@@ -72,7 +72,7 @@ q_{g,\min,i} \leq q_{g,i} &\leq q_{g,\max,i} & (\rho_{qg,lb,i}, \rho_{qg,ub,i})
 The implementation uses a **reduced-space** formulation where branch flows ``p_{fr}``, ``q_{fr}``, ``p_{to}``, ``q_{to}`` are treated as functions of the voltage state ``(\theta, |V|)`` rather than separate primal variables. This means:
 
 - The flow definition constraints are eliminated
-- Stationarity conditions automatically include all chain-rule terms via ForwardDiff on the reduced Lagrangian
+- Stationarity conditions include all reduced-space chain-rule terms analytically
 - Flow bound complementary slackness uses the computed flow expressions
 
 ## KKT System
@@ -89,7 +89,7 @@ with total dimension ``6n + 12m + 6k + n_{\text{ref}}``, where ``n`` is the numb
 
 ### KKT Conditions
 
-1. **Stationarity** (``2n + 2k`` conditions): Computed via `ForwardDiff.gradient` on the reduced-space Lagrangian ``\mathcal{L}(\theta, |V|, p_g, q_g)`` which automatically handles chain-rule terms through flow expressions.
+1. **Stationarity** (``2n + 2k`` conditions): Assembled analytically from the reduced-space Lagrangian ``\mathcal{L}(\theta, |V|, p_g, q_g)`` and branch-flow derivatives.
 
 2. **Primal feasibility** (``2n + n_{\text{ref}}`` conditions):
    - Active power balance at each bus
@@ -107,11 +107,15 @@ with total dimension ``6n + 12m + 6k + n_{\text{ref}}``, where ``n`` is the numb
 
 ### KKT Jacobian
 
-The KKT Jacobian ``\partial K / \partial z`` is computed via ForwardDiff (dense matrix). This is applied to the KKT operator which includes stationarity via the reduced-space Lagrangian, making it self-consistent with the flow chain-rule terms.
+The KKT Jacobian ``\partial K / \partial z`` is assembled analytically as a
+sparse matrix. Branch-flow derivatives and Hessians are evaluated only where
+the reduced-space KKT terms require them.
 
 ### Parameter Jacobians
 
-For each parameter ``p``, the parameter Jacobian ``\partial K / \partial p`` is also computed via ForwardDiff. The parameter is injected into the KKT operator through keyword arguments, and ForwardDiff propagates dual numbers through the entire computation.
+For each parameter ``p``, the parameter Jacobian ``\partial K / \partial p`` is
+also assembled analytically. Single-column, JVP, and VJP paths avoid
+materializing the full parameter Jacobian when only one direction is needed.
 
 The full derivative is:
 

diff --git a/examples/apf_integration.jl b/examples/apf_integration.jl
@@ -81,8 +81,10 @@ end
 
 # --- Step 2: PD for economic sensitivity ---
 println("\nStep 2: DC OPF + sensitivity analysis via PD")
-dc_net = DCNetwork(pm_data)
-d = calc_demand_vector(pm_data)
+typed_data = parse_file(case_path)
+dc_net = DCNetwork(typed_data)
+dc_net.fmax .*= 0.2
+d = calc_demand_vector(typed_data)
 prob = DCOPFProblem(dc_net, d)
 PowerDiff.solve!(prob)
 

diff --git a/examples/interactive_repl.jl b/examples/interactive_repl.jl
@@ -12,15 +12,13 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-using PowerDiff, PowerModels
-const PM = PowerModels
+using PowerDiff
 
-case_path = joinpath(dirname(pathof(PM)), "..", "test", "data", "matpower", "case14.m")
-pm_data = PM.parse_file(case_path)
-PM.make_basic_network!(pm_data)
+case_path = joinpath(get_path(:pglib), "pglib_opf_case14_ieee.m")
+typed_data = parse_file(case_path)
 
-net = DCNetwork(pm_data)
-d = PowerDiff.calc_demand_vector(pm_data)
+net = DCNetwork(typed_data)
+d = PowerDiff.calc_demand_vector(typed_data)
 
 pf = DCPowerFlowState(net, d)
 

diff --git a/experiments/ipp_market_planning.jl b/experiments/ipp_market_planning.jl
@@ -571,16 +571,12 @@ function run_tier2(; outdir::String=@__DIR__,
     println("="^65)
 
     case_path = joinpath(dirname(pathof(PM)), "..", "test", "data", "matpower", "case14.m")
-    raw = PM.parse_file(case_path)
-    PM.make_basic_network!(raw)        # populates rate_a defaults
+    net = DCNetwork(parse_file(case_path))
     # case14 ships with very loose flow limits (max loading ~15% on the binding
-    # line at default rate_a). Scale by 0.10 to bring 2 lines to the bound and
+    # line at default rate_a). Scale fmax by 0.10 to bring 2 lines to the bound and
     # produce meaningful LMP variation. This is the "stressed" scenario IPPs
     # actually care about — peak loading, outages, etc.
-    for (_, br) in raw["branch"]
-        br["rate_a"] *= fmax_scale
-    end
-    net = DCNetwork(raw)
+    net.fmax .*= fmax_scale
     # Break generator degeneracy. Without this the KKT system is singular at
     # the optimum (multiple gens at upper bound), Tikhonov regularization kicks
     # in, and the matrix free VJP returns essentially zero gradient.
@@ -728,9 +724,8 @@ const RTS_LOAD_CSV = expanduser("~/Datasets/RTS-GMLC/RTS_Data/timeseries_data_fi
 function load_rts_gmlc()
     isfile(RTS_PATH) || error("RTS_GMLC.m not at $RTS_PATH")
     raw = PM.parse_file(RTS_PATH)
-    if !isempty(raw["dcline"])
-        empty!(raw["dcline"])           # PowerModels DC line workaround
-    end
+    # PowerDiff does not model HVDC/dclines; drop them (RTS-GMLC ships DC lines).
+    isempty(raw["dcline"]) || empty!(raw["dcline"])
     PM.make_basic_network!(raw)         # populates rate_a defaults & sequential IDs
     return raw
 end
@@ -765,11 +760,10 @@ function run_tier4(; outdir::String=@__DIR__,
     println("="^65)
 
     raw = load_rts_gmlc()
+    # Bridge the (dcline-free) PowerModels dict into PowerDiff via PowerIO.
+    net = DCNetwork(PowerDiff.PowerIO.from_powermodels(raw))
     # Tighten flow limits so congestion is meaningful (RTS-GMLC ships generous limits)
-    for (_, br) in raw["branch"]
-        br["rate_a"] *= 0.5
-    end
-    net = DCNetwork(raw)
+    net.fmax .*= 0.5
     # Break gen degeneracy so the KKT system is non singular at the optimum.
     for i in eachindex(net.gmax)
         if net.gmax[i] > 0.01