diff --git a/CHANGELOG.md b/CHANGELOG.md index 4b07b80..a78dd66 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,85 @@ # Changelog +## Unreleased + +Tracks **PowerIO 0.11** (powerio C ABI 7). The `[compat]` bound is now +`PowerIO = "0.11"`; the binding gates its ABI handshake on equality, so it and its +binaries move together and the bound cannot be relaxed. + +PowerIO 0.11 replaces the accessor-and-JSON-payload layer with typed element +tables, so this is a rewrite of the seam rather than a version bump. PowerDiff is a +consumer of PowerIO: every electrical quantity PowerIO states is now read from it, +and nothing that PowerIO computes is computed again here. + +### Changed + +- **Breaking:** `parse_file` and `parse_matpower` return a + `PowerIO.PioModule{PowerIO.BalancedNetwork}` rather than a bare network. The + module carries the reader's diagnostics, the source record and the history + alongside the case, and it is what `PowerIO.emit` writes back out. Every network + constructor accepts the module or the network inside it, so `DCNetwork(net)`, + `ACNetwork(net)`, `DCOPFProblem(net)`, `ACOPFProblem(net)`, + `DCPowerFlowState(net)` and `calc_demand_vector(net)` are unchanged at the call + site. +- **Breaking:** a failure inside PowerIO reaches the caller as a + `PowerIO.PowerIOError` instead of being flattened into an `ArgumentError`. It + carries the diagnostic `code` and the records behind it, both of which wrapping + discarded. PowerDiff's own refusals stay `ArgumentError`. +- **Breaking:** `parse_matpower_struct` is removed. It was a compatibility alias + for `parse_matpower` with no callers. +- **Breaking:** `network_findings` is removed, one release after it was added and + before it ever shipped. It existed because PowerIO 0.9 reached its findings only + through a normalize pass PowerDiff had to own, deduplicate and label. In 0.11 they + are `m.diagnostics`, a property of the module `parse_file` already returns, and + a `Diagnostic` is a record with a `code` and a `severity` rather than a line of + text. Wrapping a property in an exported function is surface, not integration. +- The series conductance and susceptance are read back from PowerIO's terminal + admittance coefficients rather than derived from `r` and `x`, and each branch + terminal's charging admittance is carried on its own side. PowerDiff previously + summed the two sides, split the total evenly and discarded the charging + conductance, which lost fidelity on every source that states the terminals + separately. Nothing in `src/` inverts a branch impedance any more. +- The branch-by-bus incidence matrix is stated once and shared by both network + types, which previously assembled one each. +- Out-of-service and isolated elements are selected out here rather than by + PowerIO. `to_powerdata` is unfiltered in 0.11: every row carries a `status` and + its source row number, which is exactly the `IDMapping` index, so the two-pass + reconciliation that used to recover those numbers is gone. Isolated (`type == 4`) + buses and everything standing on them are dropped, as before. +- A generator's cost model is read off the element rather than inferred from the + shape of a converted row, so a piecewise linear cost is refused for what it is. +- The format vocabulary is PowerIO's rather than a copy of it: an unrecognized + token goes to PowerIO, which answers with what the linked library actually reads, + and a reader PowerIO gains is reachable without a PowerDiff release. PowerDiff's + historical short spellings still resolve. A bare `json` is still refused as + ambiguous. The hand-maintained distribution-format blocklist is replaced by a + check on what the source actually parsed to, which also covers time series, + scenario sets and calculation instances. +- `docs/powerio-integration.md` moves to `docs/src/` and is published; it stated the + ingest contract but sat outside the docs build, so it was never rendered. + +### Fixed + +- A branch whose source states no thermal limit takes the synthesized limit whether + the source spells that `0` or `Inf`. PowerIO 0.11 carries MATPOWER's `rate_a == 0` + out as `Inf`, which the previous `rate_a > 0` test accepted, so every unrated + branch would have reached the solver with an unbounded flow. Stock IEEE 300 leaves + all 411 branches unrated. +- Generator reactive limits may be absent. PowerIO carries a bound the case does not + state as `±Inf` rather than refusing the case, and stock case9241pegase leaves them + off seven generators. PowerDiff leaves the bound off the solver model and its KKT + complementarity row reads `ρ = 0`, the multiplier of a constraint that is not + there and the value a solver reports for a bound it was never given. Without this + the upgrade would have turned a clean refusal into `0 * Inf` — a `NaN` in the + residual and an `Inf` in the Jacobian. The KKT sparsity pattern is unchanged. + +### Added + +- A test that the ingest still reproduces MATPOWER's network, comparing + `ptdf_matrix` against `PowerModels.calc_basic_ptdf_matrix` on three PGLib cases. + It shares no code with the ingest and pins the topology, the branch susceptances + and the reference bus together. + ## 0.1.0 First release on the PowerIO 0.6.x data layer, and the first cut through the diff --git a/Project.toml b/Project.toml index 7abd798..07c1fa7 100644 --- a/Project.toml +++ b/Project.toml @@ -21,7 +21,7 @@ JuMP = "1" LazyArtifacts = "1" NLPModelsIpopt = "0.11.2" PowerModels = "0.21" -PowerIO = "0.8.2" +PowerIO = "0.11" julia = "1.10" [extras] diff --git a/README.md b/README.md index fe46892..5c8ecd1 100644 --- a/README.md +++ b/README.md @@ -33,7 +33,7 @@ Pkg.add(url="https://github.com/grid-opt-alg-lab/PowerDiff.jl.git") ```julia using PowerDiff -# Parse a supported PowerIO case into a PowerIO.BalancedNetwork +# Parse a case into a PowerIO module net = parse_file("case14.m") dc_net = DCNetwork(net) d = calc_demand_vector(net) @@ -61,13 +61,24 @@ See the [Getting Started guide](https://samueltalkington.com/research/powerdiff/ ## 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`. +PowerDiff reads files through PowerIO, and `parse_file` reads every transmission +format the linked PowerIO library does — MATPOWER `.m`, PSS/E `.raw`, PowerWorld, +PowerModels JSON, Egret JSON, pandapower, PyPSA, PSLF, gridfm, GO Challenge 3 and +the rest. The format tokens are PowerIO's, so a reader PowerIO gains works here at +once. + +A path's format is inferred from its extension; a stream has no extension, so pass +`from` (MATPOWER is assumed otherwise). A bare `json` names a container rather than +a reader, so name the one you mean: `from=:powermodels`, `:egret`, `:pandapower`. + +`parse_file` returns a `PowerIO.PioModule`. Beyond the case itself it carries +`m.diagnostics`, the reader's findings as records you can branch on by `code` and +`severity`; `m.sources[1].format`, the reader that ran; and enough for +`PowerIO.emit(m, "psse", path)` to write the case out again. ## Dependencies -- [PowerIO.jl](https://github.com/eigenergy/PowerIO.jl) — Parser and data layer (see `docs/powerio-integration.md`) +- [PowerIO.jl](https://github.com/eigenergy/PowerIO.jl) — Parser and data layer (see [PowerIO Integration](https://samueltalkington.com/research/powerdiff/powerio-integration/)) - [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 diff --git a/docs/make.jl b/docs/make.jl index bcd303b..14ce4f0 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -46,6 +46,7 @@ makedocs( "AC Power Flow" => "math/ac-power-flow.md", "AC Optimal Power Flow" => "math/ac-opf.md", ], + "PowerIO Integration" => "powerio-integration.md", "Advanced Topics" => "advanced.md", "API Reference" => "api.md", ], diff --git a/docs/powerio-integration.md b/docs/powerio-integration.md deleted file mode 100644 index 3f307ab..0000000 --- a/docs/powerio-integration.md +++ /dev/null @@ -1,29 +0,0 @@ -# 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.BalancedNetwork` 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/api.md b/docs/src/api.md index 0d24952..ca66b56 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -5,7 +5,6 @@ ```@docs parse_file parse_matpower -parse_matpower_struct get_path ``` diff --git a/docs/src/getting-started.md b/docs/src/getting-started.md index 4a1f187..43a5e16 100644 --- a/docs/src/getting-started.md +++ b/docs/src/getting-started.md @@ -7,14 +7,23 @@ This guide walks through the main workflows: DC power flow, DC OPF with LMP anal ```julia using PowerDiff -# Parse a supported PowerIO case into a PowerIO.BalancedNetwork +# Parse a case into a PowerIO module 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. +PowerDiff reads files through PowerIO, and `parse_file` reads every transmission +format the linked PowerIO library does — MATPOWER `.m`, PSS/E `.raw`, PowerWorld, +PowerModels JSON, Egret JSON, pandapower, PyPSA, PSLF, gridfm, GO Challenge 3 and the +rest. The format tokens are PowerIO's, so a reader PowerIO gains works here at once. + +A path's format is inferred from its extension; a stream has no extension, so pass +`from` (MATPOWER is assumed otherwise). A bare `json` names a container rather than a +reader, so name the one you mean: `from=:powermodels`, `:egret`, `:pandapower`. +The former PowerModels dictionary constructors were removed. + +`parse_file` returns a `PowerIO.PioModule`, which carries the reader's findings on +`m.diagnostics` alongside the case. See +[PowerIO Integration](powerio-integration.md) for the whole seam. ## Interactive Exploration diff --git a/docs/src/index.md b/docs/src/index.md index ce4b94e..35cf500 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -25,7 +25,7 @@ Pkg.add(url="https://github.com/grid-opt-alg-lab/PowerDiff.jl.git") ```julia using PowerDiff -# Parse a supported PowerIO case into a PowerIO.BalancedNetwork +# Parse a case into a PowerIO module net = parse_file("case14.m") dc_net = DCNetwork(net) d = calc_demand_vector(net) diff --git a/docs/src/powerio-integration.md b/docs/src/powerio-integration.md new file mode 100644 index 0000000..36ebbf2 --- /dev/null +++ b/docs/src/powerio-integration.md @@ -0,0 +1,90 @@ +# PowerIO Integration + +PowerIO is PowerDiff's parser and data layer, and the only one. PowerDiff has no +parser backend switch and no second representation of a case, and it derives no +electrical quantity that PowerIO already states. + +## Parsing + +`PowerDiff.parse_file(path)` returns a `PowerIO.PioModule{PowerIO.BalancedNetwork}`. +The module carries more than the network: `m.value` is the case, `m.diagnostics` +holds the reader's findings, `m.sources[1].format` names the reader that ran, and +`m.history` records what has been done to it. Pass the module, or the network +inside it, to [`DCNetwork`](@ref) or [`ACNetwork`](@ref). + +For paths the extension selects the reader unless `from` is given. For streams +pass `from`, because a stream has no extension; MATPOWER is assumed when neither +is given. + +The format vocabulary is PowerIO's rather than a copy of it. PowerDiff normalizes +its own historical short spellings (`:m`, `:raw`, `:aux`, `:pm`, `:powermodels`, +`:egret`) and hands every other token straight over, so a reader PowerIO gains is +reachable without a PowerDiff release, and an unknown token is answered by PowerIO +with what the linked library actually reads. Two refusals stay local because +PowerIO cannot make them: a bare `json` names a container rather than a reader, and +a source that parses to anything other than a balanced transmission network is +refused by naming what it holds. + +A failure inside PowerIO reaches the caller as a `PowerIO.PowerIOError`, which +carries the diagnostic `code` and the records behind it. PowerDiff does not wrap +it; wrapping would discard both. + +## What PowerIO states and what PowerDiff adds + +The network constructors read `PowerIO.to_powerdata`. It supplies per-unit powers +and ratings, radian branch angles and angle bounds, per-bus aggregated load and +shunt, right-aligned per-unit polynomial cost coefficients, each branch terminal's +own charging admittance, and the four terminal admittance coefficients from which +the series conductance and susceptance are read back exactly. + +It is **unfiltered**: every generator and branch row carries a `status`, isolated +buses are present, and each row's `i` is its source row number. PowerDiff selects +what it models -- rows whose status is set, on buses that are not isolated -- and +that source row number becomes the [`IDMapping`](@ref) index, so an out-of-service +row leaves a gap instead of renumbering the rows that remain. + +On top of that selection sit the four modeling decisions PowerIO leaves to its +consumer: + +- polynomial cost interpretation, reading the cost model off the generator and + refusing piecewise linear costs. Higher-than-quadratic polynomials are refused by + `to_powerdata` itself. +- a finite thermal limit when the source states none. MATPOWER spells that `0` and + PowerIO carries it out as `Inf`; both take the same synthesized limit, the largest + flow the endpoint voltage limits and the angle window physically admit. +- default angle-difference bounds, mapping MATPOWER's `±360` and `0, 0` spellings to + a `±60°` window. +- refusal of storage and HVDC records, which PowerDiff does not model. + +## What PowerDiff computes and why it is not a duplicate + +Two objects look like something PowerIO offers and are not. + +`admittance_matrix(net, sw)` assembles `Y` from `g`, `b`, `tap`, `shift` and the +four terminal terms on every call. `PowerIO.calc_admittance_matrix` tabulates `Y` +for one network at its parsed values with a boolean in-service flag; PowerDiff +needs `Y` as a function of parameters a caller perturbs, because `sw` is continuous +in `[0, 1]` and `calc_sensitivity(state, :vm, :b)` and its nine siblings +differentiate through `g` and `b`. The same argument applies to the susceptance +Laplacian `B = AᵀWA` with `W = -b .* sw`, and to the branch flow matrix. + +The branch-by-bus incidence matrix is stated once, in `_incidence_matrix`, and +shared by both network types. `PowerIO.calc_incidence_matrix` covers in-service +branches in table order against all buses and returns a bare sparse matrix with no +branch or bus index map, so relabeling it into PowerDiff's sorted-source-id space +costs more than stating it. + +## Reading what PowerIO found + +`m.diagnostics` is a `Vector{PowerIO.Diagnostic}`, available as soon as a case is +parsed and without emitting it again. A `Diagnostic` is a record, not a line of +text: branch on `d.code` and `d.severity`, and treat `d.message` as prose that +carries no stability promise. + +```julia +m = parse_file("case14.m") +any(d -> d.severity === :warning, m.diagnostics) +``` + +PowerDiff logs none of it. Writing a case back out in another format goes through +PowerIO directly, `PowerIO.emit(m, "psse", path)`. diff --git a/examples/sens_examples/post_event_restoration.jl b/examples/sens_examples/post_event_restoration.jl index 76de1ec..8c133d5 100644 --- a/examples/sens_examples/post_event_restoration.jl +++ b/examples/sens_examples/post_event_restoration.jl @@ -24,9 +24,8 @@ const BASE_MVA = 100.0 """Solve the IEEE300 DC OPF at supplied continuous branch switching states.""" function solve_at_switching(network_data, demand, switching) - network = DCNetwork(network_data) - network.sw .= switching - problem = DCOPFProblem(network, demand) + problem = DCOPFProblem(DCNetwork(network_data), demand) + update_switching!(problem, switching) solution = solve!(problem) return problem, solution end diff --git a/examples/sens_examples/post_event_restoration.md b/examples/sens_examples/post_event_restoration.md index d69ee4c..23e8d63 100644 --- a/examples/sens_examples/post_event_restoration.md +++ b/examples/sens_examples/post_event_restoration.md @@ -10,9 +10,8 @@ julia --project=. -e 'using Pkg; Pkg.instantiate()' julia --project=. examples/sens_examples/post_event_restoration.jl ``` -PowerDiff currently requires PowerIO 0.8.2 or later within the compatible 0.8 -release series. If an older untracked `Manifest.toml` pins PowerIO 0.1.x, it -must be regenerated from the current `Project.toml` before running the example. +The example resolves whatever `Project.toml` pins. If an untracked `Manifest.toml` +holds an older PowerIO, regenerate it from `Project.toml` before running. ## Example Problem Setup diff --git a/src/PowerDiff.jl b/src/PowerDiff.jl index 797577c..37157b7 100644 --- a/src/PowerDiff.jl +++ b/src/PowerDiff.jl @@ -103,7 +103,7 @@ export Sensitivity, silence export operand_symbols, parameter_symbols export jvp, vjp, jvp!, vjp!, dict_to_vec, vec_to_dict export kkt_layout, kkt_dims, kkt_indices -export parse_file, parse_matpower, parse_matpower_struct, get_path +export parse_file, parse_matpower, get_path # DC Power Flow Types export DCNetwork, DCPowerFlowState diff --git a/src/prob/ac_opf_solve.jl b/src/prob/ac_opf_solve.jl index 6d1f96e..8b29b51 100644 --- a/src/prob/ac_opf_solve.jl +++ b/src/prob/ac_opf_solve.jl @@ -137,8 +137,11 @@ function _extract_ac_opf_solution(prob::ACOPFProblem{JuMPBackend}) # Extract dual variables - generation bounds (inequality) ρ_pg_lb = [dual(LowerBoundRef(prob.pg[i])) for i in 1:k] ρ_pg_ub = [dual(UpperBoundRef(prob.pg[i])) for i in 1:k] - ρ_qg_lb = [dual(LowerBoundRef(prob.qg[i])) for i in 1:k] - ρ_qg_ub = [dual(UpperBoundRef(prob.qg[i])) for i in 1:k] + # An absent reactive limit has no bound reference and no multiplier. Zero is not a + # stand-in there: it is the value the KKT system pins the multiplier of an absent + # constraint to, so the solved duals and the residual agree by construction. + ρ_qg_lb = [has_lower_bound(prob.qg[i]) ? dual(LowerBoundRef(prob.qg[i])) : 0.0 for i in 1:k] + ρ_qg_ub = [has_upper_bound(prob.qg[i]) ? dual(UpperBoundRef(prob.qg[i])) : 0.0 for i in 1:k] # Extract dual variables - flow variable bounds (inequality) σ_p_fr_lb = zeros(m) diff --git a/src/prob/kkt_ac_opf.jl b/src/prob/kkt_ac_opf.jl index 5e46586..0b8fbb9 100644 --- a/src/prob/kkt_ac_opf.jl +++ b/src/prob/kkt_ac_opf.jl @@ -28,6 +28,34 @@ # Helpers # ============================================================================= +# --- Bound complementarity, with or without a stated bound --------------------- +# +# A case need not state every bound: PowerIO 0.9 carries an absent one as `±Inf`, +# which is how MATPOWER, PowerModels, pandapower and PyPSA all spell "no limit". +# PowerDiff's KKT layout is fixed, with one complementarity row per bound in a +# structure that must not move with the data, so an absent bound cannot be dropped +# from the system — and it cannot be carried either: `ρ * (x - (-Inf))` is `0 * Inf`, +# a `NaN` in the residual and an `Inf` in the Jacobian. +# +# So the row stays and states the right thing. A constraint that is not there has no +# multiplier, so its row reads `ρ = 0` — exactly the value a solver reports for a +# bound it was never given, which is what keeps the solved duals and this residual in +# agreement. Its derivative is `∂/∂ρ = 1` and `∂/∂x = 0`, both constants, so the +# fixed-regime derivative of an absent bound is zero, as it should be. +# +# Only the generator reactive limits reach here unbounded; `_network_data` refuses a +# non-finite value anywhere else, and `ACNetwork` refuses one on the bounds a JuMP +# variable must carry. + +_lb_complementarity(ρ, x, lb) = _absent_bound(lb) ? ρ : ρ * (x - lb) +_ub_complementarity(ρ, x, ub) = _absent_bound(ub) ? ρ : ρ * (ub - x) + +_lb_dcomp_dx(ρ, lb) = _absent_bound(lb) ? zero(ρ) : ρ +_ub_dcomp_dx(ρ, ub) = _absent_bound(ub) ? zero(ρ) : -ρ + +_lb_dcomp_ddual(x, lb) = _absent_bound(lb) ? one(x) : x - lb +_ub_dcomp_ddual(x, ub) = _absent_bound(ub) ? one(x) : ub - x + """Return sorted reference bus indices for the problem.""" _ref_bus_indices(prob::ACOPFProblem) = prob.data.ref_bus_keys @@ -827,14 +855,17 @@ function calc_kkt_jacobian(prob::ACOPFProblem; sol::Union{ACOPFSolution,Nothing} -vars.rho_pg_ub[i]) _add_sparse_entry!(row_idxs, col_idxs, vals, idx.rho_pg_ub[i], idx.rho_pg_ub[i], constants.pmax[i] - vars.pg[i]) + # An absent bound contributes `∂/∂qg = 0` and `∂/∂ρ = 1`. The zero is dropped + # by `_add_sparse_entry!`, the same treatment a stated bound already gets when + # its multiplier is zero because it is not binding. _add_sparse_entry!(row_idxs, col_idxs, vals, idx.rho_qg_lb[i], idx.qg[i], - vars.rho_qg_lb[i]) + _lb_dcomp_dx(vars.rho_qg_lb[i], constants.qmin[i])) _add_sparse_entry!(row_idxs, col_idxs, vals, idx.rho_qg_lb[i], idx.rho_qg_lb[i], - vars.qg[i] - constants.qmin[i]) + _lb_dcomp_ddual(vars.qg[i], constants.qmin[i])) _add_sparse_entry!(row_idxs, col_idxs, vals, idx.rho_qg_ub[i], idx.qg[i], - -vars.rho_qg_ub[i]) + _ub_dcomp_dx(vars.rho_qg_ub[i], constants.qmax[i])) _add_sparse_entry!(row_idxs, col_idxs, vals, idx.rho_qg_ub[i], idx.rho_qg_ub[i], - constants.qmax[i] - vars.qg[i]) + _ub_dcomp_ddual(vars.qg[i], constants.qmax[i])) end @inbounds for l in 1:m @@ -962,8 +993,8 @@ function kkt(z::AbstractVector, prob::ACOPFProblem, sw::AbstractVector; # Generation bounds K[idx.rho_pg_lb] .= vars.rho_pg_lb .* (pg .- pmin) K[idx.rho_pg_ub] .= vars.rho_pg_ub .* (pmax .- pg) - K[idx.rho_qg_lb] .= vars.rho_qg_lb .* (qg .- qmin) - K[idx.rho_qg_ub] .= vars.rho_qg_ub .* (qmax .- qg) + K[idx.rho_qg_lb] .= _lb_complementarity.(vars.rho_qg_lb, qg, qmin) + K[idx.rho_qg_ub] .= _ub_complementarity.(vars.rho_qg_ub, qg, qmax) # Flow variable bounds (reduced-space) K[idx.sig_p_fr_lb] .= vars.sig_p_fr_lb .* (p_fr .+ rate_a) diff --git a/src/types/ac_network.jl b/src/types/ac_network.jl index 3740a42..b30ecf7 100644 --- a/src/types/ac_network.jl +++ b/src/types/ac_network.jl @@ -207,7 +207,19 @@ function ACNetwork(net::Dict{String,<:Any}; idx_slack::Union{Nothing,Int}=nothin throw(ArgumentError("dictionary constructors were removed; parse a network file with PowerDiff.parse_file")) end -ACNetwork(net::PowerIO.BalancedNetwork; idx_slack::Union{Nothing,Int}=nothing) = +# Only the reactive generator limits model an absent bound end to end. A non-finite +# value on any other bound would reach JuMP as an infinite variable bound and the KKT +# system as `0 * Inf`; name the row and the field instead. +function _require_finite_bounds(values, field::Symbol) + for (i, v) in enumerate(values) + isfinite(v) || throw(ArgumentError( + "ACNetwork: row $i has non-finite `$field` ($v); only generator reactive " * + "limits (`qmin`/`qmax`) may be left unbounded")) + end + return nothing +end + +ACNetwork(net::PowerIOSource; idx_slack::Union{Nothing,Int}=nothing) = ACNetwork(_network_data(net); idx_slack=idx_slack) # Build from PowerDiff network tables (see `_network_data`). The `PowerIO.BalancedNetwork` @@ -222,7 +234,6 @@ function ACNetwork(data::NamedTuple; idx_slack::Union{Nothing,Int}=nothing) branch_tbl = Dict(branch.index => branch for branch in data.branch) gen_tbl = Dict(gen.index => gen for gen in data.gen) - A = spzeros(n_branch, n_bus) incidences = Vector{Tuple{Int,Int}}(undef, n_branch) f_bus = Vector{Int}(undef, n_branch) t_bus = Vector{Int}(undef, n_branch) @@ -242,40 +253,34 @@ function ACNetwork(data::NamedTuple; idx_slack::Union{Nothing,Int}=nothing) g = zeros(n_branch) b = zeros(n_branch) + A = _incidence_matrix(branch_tbl, id_map) + for orig_id in id_map.branch_ids branch = branch_tbl[orig_id] l = id_map.branch_to_idx[orig_id] fb = id_map.bus_to_idx[branch.f_bus] tb = id_map.bus_to_idx[branch.t_bus] - A[l, fb] = 1.0 - A[l, tb] = -1.0 incidences[l] = (fb, tb) f_bus[l] = fb t_bus[l] = tb br_r[l] = branch.br_r br_x[l] = branch.br_x - br_b[l] = branch.br_b - # MATPOWER models line charging as a single symmetric susceptance with no - # charging conductance, and `_network_data` already folded the two PowerIO - # sides into `br_b`. Split it evenly and leave g_fr/g_to at zero (initialized - # above). Asymmetric b_fr != b_to or nonzero g_fr/g_to from a non MATPOWER - # source would be averaged/dropped here; thread the per-side values through - # if that fidelity is ever needed. - b_fr[l] = branch.br_b / 2 - b_to[l] = branch.br_b / 2 - tap[l] = iszero(branch.tap) ? 1.0 : branch.tap + # Each terminal states its own charging admittance. A MATPOWER source splits + # one symmetric susceptance evenly and states no charging conductance; a + # source that states the terminals separately keeps that split. + g_fr[l] = branch.g_fr + b_fr[l] = branch.b_fr + g_to[l] = branch.g_to + b_to[l] = branch.b_to + br_b[l] = branch.b_fr + branch.b_to + tap[l] = branch.tap shift[l] = branch.shift tm[l] = tap[l]^2 angmin[l] = branch.angmin angmax[l] = branch.angmax rate_a[l] = branch.rate_a - z2 = branch.br_r^2 + branch.br_x^2 - if z2 > 1e-10 - g[l] = branch.br_r / z2 - b[l] = -branch.br_x / z2 - else - _SILENCE_WARNINGS[] || @warn "Branch $orig_id has near-zero impedance; treating as open." - end + g[l] = branch.g + b[l] = branch.b end g_shunt = zeros(n_bus) @@ -334,6 +339,15 @@ function ACNetwork(data::NamedTuple; idx_slack::Union{Nothing,Int}=nothing) vm_min = [bus_tbl[id].vmin for id in id_map.bus_ids] vm_max = [bus_tbl[id].vmax for id in id_map.bus_ids] + # Reactive generator limits may be absent (`±Inf`): the solver model leaves the + # bound off and the KKT system pins its multiplier to zero. Every other bound here + # becomes a JuMP variable bound and a complementarity row that needs a number, so + # say which one is missing rather than letting it reach the solver. + _require_finite_bounds(pmin, :pmin) + _require_finite_bounds(pmax, :pmax) + _require_finite_bounds(vm_min, :vmin) + _require_finite_bounds(vm_max, :vmax) + return ACNetwork( n_bus, n_branch, sparse(A), incidences, g, b, g_shunt, b_shunt, ones(n_branch), trues(n_branch), idx_slack, vm_min, vm_max, diff --git a/src/types/ac_opf_problem.jl b/src/types/ac_opf_problem.jl index d4ca223..e4b2ca4 100644 --- a/src/types/ac_opf_problem.jl +++ b/src/types/ac_opf_problem.jl @@ -500,9 +500,15 @@ function _build_jump_model(network::ACNetwork, data::ACOPFData, optimizer, silen @variable(model, va[1:n]) @variable(model, constants.vmin[i] <= vm[i in 1:n] <= constants.vmax[i], start=1.0) - # Generation variables + # Generation variables. Reactive limits may be absent (`±Inf`, PowerIO's spelling + # for a case that states none); JuMP takes only finite bounds, so an absent one is + # left off the variable rather than passed through as an infinite bound. @variable(model, constants.pmin[i] <= pg[i in 1:n_gen] <= constants.pmax[i]) - @variable(model, constants.qmin[i] <= qg[i in 1:n_gen] <= constants.qmax[i]) + @variable(model, qg[1:n_gen]) + for i in 1:n_gen + _absent_bound(constants.qmin[i]) || set_lower_bound(qg[i], constants.qmin[i]) + _absent_bound(constants.qmax[i]) || set_upper_bound(qg[i], constants.qmax[i]) + end # Branch flow variables @variable(model, -arc_fmax[i] <= p[i in 1:length(data.arcs)] <= arc_fmax[i]) @@ -741,5 +747,5 @@ function ACOPFProblem(pm_data::Dict; kwargs...) throw(ArgumentError("dictionary constructors were removed; parse a network file with PowerDiff.parse_file")) end -ACOPFProblem(net::PowerIO.BalancedNetwork; kwargs...) = ACOPFProblem(ACNetwork(net); kwargs...) +ACOPFProblem(net::PowerIOSource; kwargs...) = ACOPFProblem(ACNetwork(net); kwargs...) ACOPFProblem(data::NamedTuple; kwargs...) = ACOPFProblem(ACNetwork(data); kwargs...) diff --git a/src/types/dc_network.jl b/src/types/dc_network.jl index 9e9dc63..30c96c3 100644 --- a/src/types/dc_network.jl +++ b/src/types/dc_network.jl @@ -185,85 +185,84 @@ const DEFAULT_SHED_COST_MULTIPLIER = 10 # PowerIO input and network-table construction # ============================================================================= # -# PowerIO is the parser and data layer. `PowerIO.parse_*` reads MATPOWER/PSSE/etc. -# and `PowerIO.to_powerdata` returns normalized, per-unit, status/isolated-filtered -# data with the reference bus inferred (`type == 3`), source bus ids on `bus_i`, -# loads/shunts aggregated per bus, and polynomial costs collapsed and rescaled. -# These thin wrappers return a `PowerIO.BalancedNetwork`, and `_network_data` turns one into -# the network tables the DCNetwork and ACNetwork constructors consume. The only -# logic beyond re-keying to source bus ids is the OPF solver modeling PowerIO leaves -# to the consumer: polynomial cost interpretation, finite flow limits, default -# angle difference bounds, and rejection of records PowerDiff does not model. - -""" - parse_file(path::String; library=nothing, from=nothing, filetype=nothing) -> PowerIO.BalancedNetwork - parse_file(io::IO; from="matpower", filetype=nothing) -> PowerIO.BalancedNetwork - -Parse a supported PowerIO network file into a `PowerIO.BalancedNetwork`. - -For paths, PowerIO infers the format from the extension unless `from` is given. -For streams, pass `from` (or `filetype`) because there is no extension. JSON formats -are ambiguous; use `from=:egret` or `from=:powermodels`. - -Supported format tokens are PowerIO's tokens: `:matpower` / `:m`, `:psse` / `:raw`, -`:powerworld` / `:aux`, `:powermodels`, and `:egret`. -Pass the result to [`DCNetwork`](@ref) / [`ACNetwork`](@ref). -""" -function parse_file(io::Union{IO,String}; library=nothing, filetype=nothing, from=nothing, kwargs...) +# PowerIO is the parser and data layer, and the only one: PowerDiff has no parser +# backend switch and no second representation of a case. `PowerIO.parse` reads every +# format the linked library ships and returns a `PowerIO.PioModule{T}` carrying the +# parsed value, the reader's diagnostics, and the source record. +# +# `PowerIO.to_powerdata` states a balanced network in per unit, with radian branch +# angles, per-bus load and shunt aggregation, rescaled polynomial costs, a `status` +# flag and a source row number `i` on every row, and the four terminal admittance +# coefficients `c1..c8` per branch. It is unfiltered: out-of-service rows and +# isolated buses reach the caller. +# +# `_network_data` selects what PowerDiff models -- in-service elements on buses that +# are not isolated -- and adds the OPF modeling PowerIO leaves to its consumers: +# polynomial cost interpretation, a finite thermal limit when the source states +# none, default angle-difference bounds, and refusal of records PowerDiff has no +# model for. It derives no electrical quantity: the series admittance and the +# terminal charging terms are read off PowerIO's coefficients. + +""" + parse_file(path::String; library=nothing, from=nothing, filetype=nothing) -> PowerIO.PioModule + parse_file(io::IO; from="matpower", filetype=nothing) -> PowerIO.PioModule + +Parse a network into a `PowerIO.PioModule{PowerIO.BalancedNetwork}`, in any +transmission format the linked PowerIO library reads. + +For paths, PowerIO infers the format from the extension unless `from` is given. For +streams, pass `from` (or `filetype`), because a stream has no extension; MATPOWER is +assumed when neither is given. A bare `json` names a container rather than a reader, +so name the reader (`from=:powermodels`, `:egret`, `:pandapower`, `:goc3`, `:surge`, +`:opfdata`). + +`from` takes PowerIO's own format tokens, and the vocabulary is PowerIO's rather +than a copy of it: an unrecognized token goes to PowerIO, which answers with what +the linked library actually reads, so a reader PowerIO gains is usable here at once. +PowerDiff's historical short spellings (`:m`, `:raw`, `:aux`, `:pm`, `:powermodels`, +`:egret`) resolve to the tokens they always meant. + +A source that parses to anything other than a balanced transmission network -- a +distribution case, a time series, a scenario set, a calculation instance -- is +refused by naming what it holds. + +The returned module carries more than the network: `m.value` is the case, +`m.diagnostics` is what the reader retained or had to assume, `m.sources[1].format` +is the reader that ran, and `PowerIO.emit(m, format, path)` writes the case out +again. Pass the module to [`DCNetwork`](@ref) or [`ACNetwork`](@ref). +""" +function parse_file(source::Union{IO,AbstractString}; library=nothing, filetype=nothing, + from=nothing, kwargs...) isempty(kwargs) || throw(ArgumentError( "unsupported parse_file keyword(s): $(join(string.(keys(kwargs)), ", "))")) fmt = _powerio_format_hint(from, filetype) - if io isa String - resolved = _resolve_case_path(io, library) - try - return isnothing(fmt) ? PowerIO.parse_file(resolved) : PowerIO.parse_file(resolved; from=fmt) - catch e - e isa ArgumentError && rethrow() - throw(ArgumentError("PowerDiff.parse_file: " * sprint(showerror, e))) - end + m = if source isa AbstractString + PowerIO.parse(_resolve_case_path(source, library); format=fmt) else - fmt = isnothing(fmt) ? "matpower" : fmt - try - return PowerIO.parse_file(io, fmt) - catch e - e isa ArgumentError && rethrow() - throw(ArgumentError("PowerDiff.parse_file: " * sprint(showerror, e))) - end + PowerIO.parse(source; format=something(fmt, "matpower")) end + return _require_balanced(m) end """ - parse_matpower(io::IO) -> PowerIO.BalancedNetwork - parse_matpower(file::String; library=nothing) -> PowerIO.BalancedNetwork + parse_matpower(io::IO) -> PowerIO.PioModule + parse_matpower(file::String; library=nothing) -> PowerIO.PioModule -Parse MATPOWER v2 data into a `PowerIO.BalancedNetwork`. +Parse MATPOWER v2 data into a `PowerIO.PioModule{PowerIO.BalancedNetwork}`. """ -function parse_matpower(io::IO) - try - return PowerIO.parse_file(io, "matpower") - catch e - e isa ArgumentError && rethrow() - throw(ArgumentError("PowerDiff.parse_matpower: " * sprint(showerror, e))) - end -end - -function parse_matpower(file::String; library=nothing) - resolved = _resolve_case_path(file, library) - isfile(resolved) || throw(ArgumentError("invalid MATPOWER file $resolved")) - try - return PowerIO.parse_file(String(resolved); from="matpower") - catch e - e isa ArgumentError && rethrow() - throw(ArgumentError("PowerDiff.parse_matpower: " * sprint(showerror, e))) - end -end +parse_matpower(io::IO) = _require_balanced(PowerIO.parse(io; format="matpower")) -""" - parse_matpower_struct(file::String; library=nothing) +parse_matpower(file::String; library=nothing) = + _require_balanced(PowerIO.parse(_resolve_case_path(file, library); format="matpower")) -Compatibility alias for [`parse_matpower`](@ref). -""" -parse_matpower_struct(file::String; library=nothing) = parse_matpower(file; library=library) +# PowerIO returns one of twenty value kinds. Only a balanced transmission network +# has the buses, branches and generators PowerDiff differentiates, so the type check +# is the refusal: it covers distribution cases, series carriers and calculation +# instances alike, without a list here that trails what PowerIO parses. +_require_balanced(m::PowerIO.PioModule{PowerIO.BalancedNetwork}) = m +_require_balanced(m::PowerIO.PioModule) = throw(ArgumentError( + "PowerDiff models balanced transmission networks; this source parsed to a " * + "$(nameof(typeof(m.value)))")) _resolve_case_path(path::AbstractString, ::Nothing) = String(path) _resolve_case_path(path::AbstractString, library) = joinpath(get_path(library), path) @@ -278,152 +277,195 @@ function _powerio_format_hint(from, filetype) return f1 end +# The short spellings PowerDiff has always accepted, mapped to PowerIO's own token. +# This is a courtesy layer, not a gate: a token absent from it passes through +# untouched, and PowerIO answers for it. +const _FORMAT_ALIASES = Dict( + "m" => "matpower", + "raw" => "psse", + "aux" => "powerworld", + "pm" => "powermodels-json", + "powermodels" => "powermodels-json", + "egret" => "egret-json", +) + function _format_token(x) s = lowercase(String(x)) startswith(s, ".") && (s = s[2:end]) + isempty(s) && throw(ArgumentError("network format hint is empty")) + # A bare `json` names a container, not a reader, and PowerIO has several. Say so + # here: routing it would pick one by guess. s == "json" && throw(ArgumentError( - "JSON input is ambiguous; pass from=:egret or from=:powermodels")) - s in ("m", "matpower") && return "matpower" - s in ("raw", "psse") && return "psse" - s in ("aux", "powerworld") && return "powerworld" - s in ("pm", "powermodels", "powermodels-json") && return "powermodels-json" - s in ("egret", "egret-json") && return "egret-json" - throw(ArgumentError( - "unsupported network format $x (expected matpower, psse/raw, powerworld/aux, powermodels-json, or egret-json)")) + "JSON input is ambiguous; name the reader (from=:powermodels, from=:egret, " * + "from=:pandapower, from=:goc3, from=:surge, from=:opfdata)")) + return get(_FORMAT_ALIASES, s, s) end """ + _network_data(m::PowerIO.PioModule) -> NamedTuple _network_data(net::PowerIO.BalancedNetwork) -> NamedTuple -Build PowerDiff network tables from `PowerIO.to_powerdata(net)`. - -`to_powerdata` does per-unit scaling, status/isolated filtering, per-bus -load/shunt aggregation, reference-bus inference (`type == 3`), source bus ids on -`bus_i`, and polynomial cost collapse/rescaling, returning dense file-order rows. -This adapter keys bus references back to source bus ids (so [`IDMapping`](@ref)'s -sorted ordering is preserved) and applies the OPF modeling PowerIO leaves to the -consumer: polynomial cost interpretation (rejecting PWL and higher-than-quadratic), -a finite flow limit fallback when `rate_a == 0`, default angle difference bounds, -and rejection of storage / HVDC records that PowerDiff does not model. - -The returned `bus`/`gen`/`branch` rows mirror the field names the network -constructors expect, with loads/shunts already folded into per-bus `pd/qd/gs/bs`. -`shunt` re-exposes those bus shunts as a table (one `(; index, shunt_bus, gs, bs)` -record per bus with a nonzero shunt admittance) for callers that want shunt records. - -Bus rows carry the source bus id on `bus_i`, so [`IDMapping`](@ref)`.bus_ids` -(and any bus-indexed sensitivity `row_to_id`) map back to the input network. -Generator and branch `index` values are source row numbers among the unfiltered -PowerIO rows, so out-of-service rows leave gaps instead of renumbering active rows. -""" -function _network_data(net) - # Reject records PowerDiff does not model. Both guards read the raw network so - # they stay consistent: to_powerdata's filtered output drops out-of-service - # records, which would silently accept a file that declares them. - isempty(PowerIO.hvdc(net)) || throw(ArgumentError( +Build PowerDiff network tables from a parsed PowerIO network. + +`PowerIO.to_powerdata` supplies every value: per-unit powers and ratings, radian +branch angles and angle bounds, per-bus aggregated load and shunt, right-aligned +per-unit polynomial cost coefficients, and the four terminal admittance +coefficients from which the series conductance and susceptance are read. It is +unfiltered, so this selects what PowerDiff models: rows whose `status` is set, on +buses that are not isolated (`type == 4`). Its row field `i` is the source row +number, which becomes the [`IDMapping`](@ref) index, so out-of-service rows leave +gaps rather than renumbering the rows that remain. + +On top of that selection sit the four modeling decisions PowerIO leaves to its +consumer: polynomial cost interpretation (rejecting piecewise linear), a +synthesized thermal limit when the source states none, default angle-difference +bounds, and refusal of storage and HVDC records. + +Bus rows carry the source bus id on `bus_i`, so [`IDMapping`](@ref)`.bus_ids` and +any bus-indexed sensitivity `row_to_id` map back to the input network. +""" +# A parsed case reaches PowerDiff either as the module `parse_file` returns or as the +# network inside it. Both name the same case, so both construct the same networks. +const PowerIOSource = Union{PowerIO.PioModule{PowerIO.BalancedNetwork},PowerIO.BalancedNetwork} + +_network_data(m::PowerIO.PioModule{PowerIO.BalancedNetwork}) = _network_data(m.value) + +function _network_data(net::PowerIO.BalancedNetwork) + isempty(net.hvdc) || throw(ArgumentError( "PowerDiff does not support HVDC/dcline records; remove or convert dcline before parsing")) - isempty(PowerIO.storage(net)) || throw(ArgumentError( + isempty(net.storage) || throw(ArgumentError( "PowerDiff does not support storage records; remove or convert storage before parsing")) + pd = PowerIO.to_powerdata(net) - isempty(pd.bus) && throw(ArgumentError("network has no active buses")) - isempty(pd.gen) && throw(ArgumentError("network has no active generators")) - isempty(pd.branch) && throw(ArgumentError("network has no active branches")) - - orig = [Int(b.bus_i) for b in pd.bus] # dense file-order index -> source bus id - gen_source_rows, branch_source_rows = _active_source_rows(net, pd) - - buses = [(; bus_i=orig[i], bus_type=Int(b.type), - pd=Float64(b.pd), qd=Float64(b.qd), gs=Float64(b.gs), bs=Float64(b.bs), - vm=Float64(b.vm), va=Float64(b.va), vmin=Float64(b.vmin), vmax=Float64(b.vmax)) - for (i, b) in enumerate(pd.bus)] - - # Costs come straight from to_powerdata's gen rows (already per-unit and - # right-aligned). Map dense `gen.bus` to the source bus id via `orig`. - gens = [(; index=gen_source_rows[j], gen_bus=orig[g.bus], - pg=Float64(g.pg), qg=Float64(g.qg), qmin=Float64(g.qmin), qmax=Float64(g.qmax), - vg=Float64(g.vg), pmin=Float64(g.pmin), pmax=Float64(g.pmax), cost=_poly_cost(g)) - for (j, g) in enumerate(pd.gen)] - - branches = [_branch_row(branch_source_rows[l], br, orig, buses) for (l, br) in enumerate(pd.branch)] - all(br.rate_a > 0 for br in branches) || throw(ArgumentError( - "branches must have positive thermal limits after normalization")) - - # to_powerdata folds shunts into per-bus gs/bs (which the constructors consume). - # Re-expose them as a table, one record per bus with a nonzero shunt admittance, - # for callers that want shunt records back. - shunt_buses = [b for b in buses if b.gs != 0.0 || b.bs != 0.0] - shunts = [(; index=i, shunt_bus=b.bus_i, gs=b.gs, bs=b.bs) for (i, b) in enumerate(shunt_buses)] - - return (; name=PowerIO.network_name(net), baseMVA=Float64(pd.baseMVA), - bus=buses, gen=gens, branch=branches, shunt=shunts) -end + isempty(pd.bus) && throw(ArgumentError("network has no buses")) -function _active_source_rows(net, pd) - raw = PowerIO.to_powerdata(net; filtered=false) - kept_bus_ids = Set(Int(b.bus_i) for b in pd.bus) - raw_bus_id = Dict(Int(b.i) => Int(b.bus_i) for b in raw.bus) + # `to_powerdata` states the cost coefficients but not which cost model produced + # them, so the model comes off the element table. One pass, read once. + gens = collect(net.generators) - gen_rows = Int[] - for (row, gen) in enumerate(raw.gen) - status = hasproperty(gen, :status) ? Int(gen.status) != 0 : true - bus_id = get(raw_bus_id, Int(gen.bus), nothing) - status && bus_id in kept_bus_ids && push!(gen_rows, row) - end + # An isolated bus is out of service, and so is every element on it. + live = [Int(b.type) != 4 for b in pd.bus] + keep_bus = findall(live) + isempty(keep_bus) && throw(ArgumentError("network has no in-service buses")) - branch_rows = Int[] - for (row, br) in enumerate(raw.branch) - status = hasproperty(br, :status) ? Int(br.status) != 0 : true - f_id = get(raw_bus_id, Int(br.f_bus), nothing) - t_id = get(raw_bus_id, Int(br.t_bus), nothing) - status && f_id in kept_bus_ids && t_id in kept_bus_ids && push!(branch_rows, row) - end + orig = [Int(b.bus_i) for b in pd.bus] # table position -> source bus id + buses = [_bus_row(orig[i], pd.bus[i]) for i in keep_bus] + vmax_by_id = Dict(b.bus_i => b.vmax for b in buses) + + kept_gen = [Int(g.i) for g in pd.gen if Int(g.status) != 0 && live[g.bus]] + isempty(kept_gen) && throw(ArgumentError("network has no in-service generators")) + gen_rows = [_gen_row(pd.gen[j], orig[pd.gen[j].bus], gens[j]) for j in kept_gen] - length(gen_rows) == length(pd.gen) || throw(ArgumentError( - "PowerDiff could not map active generators back to source rows")) - length(branch_rows) == length(pd.branch) || throw(ArgumentError( - "PowerDiff could not map active branches back to source rows")) + kept_branch = [Int(br.i) for br in pd.branch + if Int(br.status) != 0 && live[br.f_bus] && live[br.t_bus]] + isempty(kept_branch) && throw(ArgumentError("network has no in-service branches")) + branches = [_branch_row(pd.branch[l], orig, vmax_by_id) for l in kept_branch] - return gen_rows, branch_rows + return (; name = net.name, + baseMVA = Float64(pd.baseMVA), + bus = buses, gen = gen_rows, branch = branches) end -# Build one PowerDiff branch row from a to_powerdata branch: map dense f_bus/t_bus to -# source ids, default the angle window, and synthesize a finite rate_a when the -# source leaves it at 0 (unlimited), using the endpoint buses' vmax limits. -function _branch_row(l, br, orig, buses) +# ============================================================================= +# Absent numeric bounds +# ============================================================================= +# +# PowerIO passes an absent bound through as `±Inf` instead of refusing the case: +# `Inf` is how MATPOWER, PowerModels, pandapower and PyPSA all spell "no limit", and +# stock pglib cases carry it (case9241pegase leaves the reactive limits off seven +# generators). PowerDiff's KKT layout is fixed, with one complementarity row per +# bound, so `Inf` cannot simply flow in: `ρ * (qg - qmin)` with `qmin == -Inf` is +# `0 * Inf`, a `NaN` in the residual and an `Inf` in the Jacobian. +# +# Reactive generator limits model absence properly -- the bound is left off the +# solver model and its complementarity row reads `ρ = 0`, the multiplier of a +# constraint that is not there (see `_lb_complementarity` and its derivatives in +# `prob/kkt_ac_opf.jl`). A branch that states no rating takes a synthesized one. +# Everywhere else PowerIO rejects a non-finite value itself, naming the element and +# the field, so nothing here repeats that check; `_require_finite_bounds` covers the +# caller-built tables PowerIO never sees. + +""" + _absent_bound(v) -> Bool + +Whether a variable bound is absent. PowerIO spells an absent bound `±Inf`. +""" +_absent_bound(v::Real) = !isfinite(v) + +_bus_row(bus_id, b) = (; + bus_i = bus_id, bus_type = Int(b.type), + pd = Float64(b.pd), qd = Float64(b.qd), + gs = Float64(b.gs), bs = Float64(b.bs), + vm = Float64(b.vm), va = Float64(b.va), + vmin = Float64(b.vmin), vmax = Float64(b.vmax), +) + +_gen_row(g, bus_id, element) = (; + index = Int(g.i), gen_bus = bus_id, + pg = Float64(g.pg), qg = Float64(g.qg), + qmin = Float64(g.qmin), qmax = Float64(g.qmax), + vg = Float64(g.vg), + pmin = Float64(g.pmin), pmax = Float64(g.pmax), + cost = _poly_cost(g, element), +) + +# One PowerDiff branch row from a `to_powerdata` branch row: default the angle +# window, read the series admittance off PowerIO's terminal coefficients, and +# synthesize a finite rate_a when the source states no thermal limit. +# +# "States no limit" is `rate_a == 0`, MATPOWER's spelling, or a non-finite `rate_a`, +# which is how PowerIO carries an unbounded rating out of the formats that write +# one. Both mean the same thing and both take the same synthesized limit, the +# largest flow the endpoint voltage limits and the angle window physically admit -- +# a bound by construction, not an invented rating that could bind. +function _branch_row(br, orig, vmax_by_id) angmin, angmax = _normalize_angle_bounds(Float64(br.angmin), Float64(br.angmax)) - rate_a = br.rate_a > 0 ? Float64(br.rate_a) : + f_bus = orig[br.f_bus] + t_bus = orig[br.t_bus] + g_to = Float64(br.g_to) + b_to = Float64(br.b_to) + # `to_powerdata`'s eighth and seventh coefficients are the to-side terminal + # admittance added to the series admittance, so the series part reads back + # exactly. Nothing here inverts an impedance. + g = Float64(br.c7) - g_to + b = Float64(br.c8) - b_to + raw_rate_a = Float64(br.rate_a) + rate_a = (isfinite(raw_rate_a) && raw_rate_a > 0) ? raw_rate_a : _fallback_rate_a(Float64(br.br_r), Float64(br.br_x), angmin, angmax, - buses[br.f_bus].vmax, buses[br.t_bus].vmax) - return (; index=l, f_bus=orig[br.f_bus], t_bus=orig[br.t_bus], - br_r=Float64(br.br_r), br_x=Float64(br.br_x), br_b=Float64(br.b_fr + br.b_to), - rate_a=rate_a, rate_b=Float64(br.rate_b), rate_c=Float64(br.rate_c), - tap=Float64(br.tap), shift=Float64(br.shift), angmin=angmin, angmax=angmax) + vmax_by_id[f_bus], vmax_by_id[t_bus]) + return (; index = Int(br.i), f_bus = f_bus, t_bus = t_bus, + br_r = Float64(br.br_r), br_x = Float64(br.br_x), + g = g, b = b, + g_fr = Float64(br.g_fr), b_fr = Float64(br.b_fr), + g_to = g_to, b_to = b_to, + rate_a = rate_a, + rate_b = _unlimited_as_zero(br.rate_b), rate_c = _unlimited_as_zero(br.rate_c), + tap = Float64(br.tap), shift = Float64(br.shift), + angmin = angmin, angmax = angmax) end -# Interpret a PowerIO gen row's polynomial cost as PowerDiff's (quadratic, linear, -# constant) tuple. to_powerdata returns polynomial (model 2) costs as a right-aligned, -# per-unit (cq, cl, cc) triple and rejects higher-than-quadratic itself. A generator -# with no gencost row comes back as `model_poly == false` with `n == 0` (cost-free); -# piecewise-linear (model 1) is `model_poly == false` with `n > 0` and is unsupported. -function _poly_cost(g) - if !g.model_poly - Int(g.n) == 0 && return (0.0, 0.0, 0.0) - throw(ArgumentError( - "piecewise linear generator costs are not supported; convert model 1 costs to polynomial model 2 before parsing")) - end - # to_powerdata right-aligns the (quadratic, linear, constant) triple, but guard the - # indexing so a model-2 cost shorter than 3 terms (purely linear/constant) zero-pads - # the missing leading coefficients instead of throwing a BoundsError. +# `rate_b` / `rate_c` are carried through untouched and unused by either formulation. +# Keep them numeric so a caller reading the tables never meets an `Inf`: an unbounded +# rating reads as `0`, which is the same "no limit" spelling `rate_a` arrives in. +_unlimited_as_zero(x) = (v = Float64(x); isfinite(v) ? v : 0.0) + +# Interpret a generator's polynomial cost as PowerDiff's (quadratic, linear, +# constant) tuple. `to_powerdata` returns a model 2 cost as a right-aligned per-unit +# triple and rejects higher-than-quadratic itself; the cost model comes off the +# element, which states it verbatim. A generator with no cost record is cost-free. +function _poly_cost(g, element) + cost = element.cost + cost === nothing && return (0.0, 0.0, 0.0) + Int(cost.model) == 2 || throw(ArgumentError( + "generator $(Int(g.i)) states cost model $(Int(cost.model)); PowerDiff models " * + "polynomial (model 2) costs. Convert piecewise linear costs before parsing")) c = g.c - cq = length(c) >= 3 ? Float64(c[end-2]) : 0.0 - cl = length(c) >= 2 ? Float64(c[end-1]) : 0.0 - cc = length(c) >= 1 ? Float64(c[end]) : 0.0 - return (cq, cl, cc) + return (Float64(c[1]), Float64(c[2]), Float64(c[3])) end # PowerDiff's OPF needs a finite thermal limit on every branch. When the source -# leaves rate_a == 0 (unlimited), synthesize one from the bus voltage limits and -# the branch impedance / angle window, matching the previous native parser. +# states none, synthesize one from the bus voltage limits and the branch impedance +# and angle window. function _fallback_rate_a(r::Float64, x::Float64, angmin::Float64, angmax::Float64, fr_vmax::Float64, to_vmax::Float64) theta_max = max(abs(angmin), abs(angmax)) @@ -444,6 +486,26 @@ function _normalize_angle_bounds(angmin::Float64, angmax::Float64) return angmin, angmax end +# The branch-by-bus incidence matrix over PowerDiff's own index space: one row per +# branch in `id_map.branch_ids` order, one column per bus in sorted source-id order, +# `+1` at the from bus and `-1` at the to bus. Both network types read it from here, +# so the two describe the same graph by construction. +# +# PowerIO's `calc_incidence_matrix` covers in-service branches in table order against +# all buses and returns a bare sparse matrix with no branch or bus index map +# (eigenergy/PowerIO.jl#2 in the issues this port filed), so relabeling it into this +# space costs more than stating it. +function _incidence_matrix(branch_tbl, id_map::IDMapping) + m = length(id_map.branch_ids) + A = spzeros(Float64, m, length(id_map.bus_ids)) + for orig_id in id_map.branch_ids + br = branch_tbl[orig_id] + row = id_map.branch_to_idx[orig_id] + A[row, id_map.bus_to_idx[br.f_bus]] = 1.0 + A[row, id_map.bus_to_idx[br.t_bus]] = -1.0 + end + return A +end # ============================================================================= # DCNetwork Constructors # ============================================================================= @@ -468,7 +530,7 @@ net = parse_file("case14.m") dc_net = DCNetwork(net) ``` """ -DCNetwork(net::PowerIO.BalancedNetwork; tau::Float64=DEFAULT_TAU, ref_bus::Union{Nothing,Int}=nothing) = +DCNetwork(net::PowerIOSource; tau::Float64=DEFAULT_TAU, ref_bus::Union{Nothing,Int}=nothing) = DCNetwork(_network_data(net); tau=tau, ref_bus=ref_bus) # Build from PowerDiff network tables (see `_network_data`). The `PowerIO.BalancedNetwork` @@ -484,16 +546,7 @@ function DCNetwork(data::NamedTuple; tau::Float64=DEFAULT_TAU, ref_bus::Union{No branch_tbl = Dict(branch.index => branch for branch in data.branch) gen_tbl = Dict(gen.index => gen for gen in data.gen) - # Incidence matrix A (m × n) from active branches using id_map translation - A = spzeros(m, n) - for orig_id in id_map.branch_ids - br = branch_tbl[orig_id] - row = id_map.branch_to_idx[orig_id] - f_col = id_map.bus_to_idx[br.f_bus] - t_col = id_map.bus_to_idx[br.t_bus] - A[row, f_col] = 1.0 - A[row, t_col] = -1.0 - end + A = _incidence_matrix(branch_tbl, id_map) # Generator-bus incidence matrix G_inc (n × k) G_inc = spzeros(n, k) @@ -504,20 +557,9 @@ function DCNetwork(data::NamedTuple; tau::Float64=DEFAULT_TAU, ref_bus::Union{No G_inc[row, col] = 1.0 end - # Branch susceptances: b = imag(1/z) - b = zeros(m) - for orig_id in id_map.branch_ids - br = branch_tbl[orig_id] - idx = id_map.branch_to_idx[orig_id] - r = br.br_r - x = br.br_x - z2 = r^2 + x^2 - if z2 > 1e-10 - b[idx] = -x / z2 - else - _SILENCE_WARNINGS[] || @warn "Branch $(orig_id) has near-zero impedance (|z|² = $(z2)); treating as open (zero admittance)." - end - end + # Series susceptance per branch, as PowerIO states it: `imag(1/z)`, negative for + # an inductive branch, so `W = -b .* sw` weights the Laplacian positively. + b = [Float64(branch_tbl[id_map.branch_ids[i]].b) for i in 1:m] # All branches initially active sw = ones(m) @@ -628,7 +670,7 @@ function calc_demand_vector(network::DCNetwork) return copy(network.demand) end -calc_demand_vector(net::PowerIO.BalancedNetwork) = calc_demand_vector(_network_data(net)) +calc_demand_vector(net::PowerIOSource) = calc_demand_vector(_network_data(net)) calc_demand_vector(data::NamedTuple) = calc_demand_vector(data, IDMapping(data)) function calc_demand_vector(data::NamedTuple, id_map::IDMapping) @@ -927,7 +969,7 @@ Construct DCPowerFlowState from a parsed PowerIO network. If `d` is not provided, extracts demand from the network. If `g` is not provided, aggregates generation from gen data to buses. """ -function DCPowerFlowState(net::PowerIO.BalancedNetwork; g::Union{Nothing,AbstractVector}=nothing, d::Union{Nothing,AbstractVector}=nothing) +function DCPowerFlowState(net::PowerIOSource; g::Union{Nothing,AbstractVector}=nothing, d::Union{Nothing,AbstractVector}=nothing) net = DCNetwork(net) if isnothing(d) diff --git a/src/types/dc_opf_problem.jl b/src/types/dc_opf_problem.jl index 3894014..fccde9c 100644 --- a/src/types/dc_opf_problem.jl +++ b/src/types/dc_opf_problem.jl @@ -330,7 +330,7 @@ function DCOPFProblem(pm_data::Dict{String,<:Any}; kwargs...) throw(ArgumentError("dictionary constructors were removed; parse a network file with PowerDiff.parse_file")) end -DCOPFProblem(net::PowerIO.BalancedNetwork; kwargs...) = DCOPFProblem(_network_data(net); kwargs...) +DCOPFProblem(net::PowerIOSource; kwargs...) = DCOPFProblem(_network_data(net); kwargs...) function DCOPFProblem(data::NamedTuple; d::Union{Nothing,AbstractVector}=nothing, tau::Float64=DEFAULT_TAU, kwargs...) network = DCNetwork(data; tau=tau) diff --git a/test/common.jl b/test/common.jl index cf4a3f2..68dcdfc 100644 --- a/test/common.jl +++ b/test/common.jl @@ -58,9 +58,20 @@ pd_bus(bus_i, bus_type; pd=0.0, qd=0.0, gs=0.0, bs=0.0, vm=1.0, va=0.0, vmin=0.9 (; bus_i, bus_type, pd, qd, gs, bs, vm, va, vmin, vmax) pd_gen(index, gen_bus; pg=0.0, qg=0.0, qmin=0.0, qmax=0.0, vg=1.0, pmin=0.0, pmax=0.0, cost=(0.0, 0.0, 0.0)) = (; index, gen_bus, pg, qg, qmin, qmax, vg, pmin, pmax, cost) -pd_branch(index, f_bus, t_bus; br_r, br_x, br_b=0.0, rate_a=Inf, rate_b=0.0, rate_c=0.0, - tap=1.0, shift=0.0, angmin=-pi / 3, angmax=pi / 3) = - (; index, f_bus, t_bus, br_r, br_x, br_b, rate_a, rate_b, rate_c, tap, shift, angmin, angmax) +# `_network_data` reads the series admittance and each terminal's charging +# admittance off PowerIO rather than deriving them, so a hand-built branch row has +# to state them. This is fixture arithmetic: PowerDiff itself never inverts an +# impedance. `br_b` keeps MATPOWER's whole-line charging spelling and splits evenly, +# which is what a MATPOWER source states; pass `b_fr`/`b_to` to state each terminal. +function pd_branch(index, f_bus, t_bus; br_r, br_x, br_b=0.0, + g_fr=0.0, b_fr=br_b / 2, g_to=0.0, b_to=br_b / 2, + rate_a=Inf, rate_b=0.0, rate_c=0.0, + tap=1.0, shift=0.0, angmin=-pi / 3, angmax=pi / 3) + y = (iszero(br_r) && iszero(br_x)) ? zero(ComplexF64) : inv(complex(br_r, br_x)) + return (; index, f_bus, t_bus, br_r, br_x, g=real(y), b=imag(y), + g_fr, b_fr, g_to, b_to, + rate_a, rate_b, rate_c, tap, shift, angmin, angmax) +end pd_case(bus, gen, branch; name="case", baseMVA=100.0) = (; name, baseMVA, bus, gen, branch) """ diff --git a/test/runtests.jl b/test/runtests.jl index 265acb1..3d7ff30 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -899,6 +899,7 @@ include("test_acpf_va_flow.jl") include("test_parameter_transforms.jl") include("test_parser_parity.jl") include("test_non_matpower_parsers.jl") +include("test_powerio_integration.jl") include("test_ac_opf_exa_backend.jl") include("test_ac_opf_all_sens.jl") include("test_ac_topology_sens.jl") diff --git a/test/test_non_matpower_parsers.jl b/test/test_non_matpower_parsers.jl index 417070e..d14621e 100644 --- a/test/test_non_matpower_parsers.jl +++ b/test/test_non_matpower_parsers.jl @@ -8,6 +8,8 @@ function _assert_network_tables_compatible(actual, baseline; label, demand_rtol= @test sum(calc_demand_vector(actual)) ≈ sum(calc_demand_vector(baseline)) rtol=demand_rtol atol=1e-8 @test all(isfinite, [br.rate_a for br in actual.branch]) @test all(>(0), [br.rate_a for br in actual.branch]) + @test sum(br.b_fr + br.b_to for br in actual.branch) ≈ + sum(br.b_fr + br.b_to for br in baseline.branch) rtol=1e-6 atol=1e-8 end end @@ -37,31 +39,39 @@ end baseline = PowerDiff._network_data(matpower_net) @testset "PowerModels JSON" begin - text, warnings = PowerIO.to_format(matpower_net, "powermodels-json") - @test isempty(warnings) + # `emit` returns an `EmitResult`: the artifacts it wrote, the fidelity of the + # conversion, and the writer's findings as `Diagnostic` records. + result = PowerIO.emit(matpower_net, "powermodels-json") + @test result.layout == "file" + @test isempty(result.diagnostics) + text = result.text + @test text isa String parsed = PowerDiff.parse_file(IOBuffer(text); from=:powermodels) - @test PowerIO.source_format(parsed) == "PowerModelsJson" + @test parsed.sources[1].format == "powermodels-json" data = PowerDiff._network_data(parsed) _assert_network_tables_compatible(data, baseline; label="PowerModels JSON") _assert_constructs_and_solves(data; label="PowerModels JSON") mktempdir() do dir path = joinpath(dir, "case14.json") - write(path, text) + written = PowerIO.emit(matpower_net, "powermodels-json", path) + @test written.layout == "file" + @test isfile(written.artifacts[1].path) parsed_path = PowerDiff.parse_file(path; from="powermodels-json") - @test PowerIO.source_format(parsed_path) == "PowerModelsJson" + @test parsed_path.sources[1].format == "powermodels-json" _assert_network_tables_compatible( PowerDiff._network_data(parsed_path), baseline; label="PowerModels JSON path") end end @testset "Egret JSON" begin - text, warnings = PowerIO.to_format(matpower_net, "egret-json") - @test isempty(warnings) + result = PowerIO.emit(matpower_net, "egret-json") + @test isempty(result.diagnostics) + text = result.text parsed = PowerDiff.parse_file(IOBuffer(text); from=:egret) - @test PowerIO.source_format(parsed) == "EgretJson" + @test parsed.sources[1].format == "egret-json" data = PowerDiff._network_data(parsed) _assert_network_tables_compatible(data, baseline; label="Egret JSON") _assert_constructs_and_solves(data; label="Egret JSON") @@ -70,20 +80,26 @@ end end @testset "PSS/E RAW" begin - text, warnings = PowerIO.to_format(matpower_net, "psse") - @test warnings isa AbstractVector + result = PowerIO.emit(matpower_net, "psse") + # A `Diagnostic` is a record, not a line of text: branch on `code` and + # `severity`, never on the message, which carries no stability promise. + @test result.diagnostics isa Vector{PowerIO.Diagnostic} + @test all(d -> d.severity in (:error, :warning, :remark, :note), result.diagnostics) + @test !any(d -> d.severity === :error, result.diagnostics) + text = result.text parsed = PowerDiff.parse_file(IOBuffer(text); from=:psse) - @test PowerIO.source_format(parsed) == "Psse" + @test parsed.sources[1].format == "psse" data = PowerDiff._network_data(parsed) _assert_network_tables_compatible(data, baseline; label="PSS/E RAW", demand_rtol=1e-5) _assert_constructs_and_solves(data; label="PSS/E RAW") mktempdir() do dir path = joinpath(dir, "case14.raw") - write(path, text) + written = PowerIO.emit(matpower_net, "psse", path) + @test isfile(written.artifacts[1].path) parsed_path = PowerDiff.parse_file(path) - @test PowerIO.source_format(parsed_path) == "Psse" + @test parsed_path.sources[1].format == "psse" _assert_network_tables_compatible( PowerDiff._network_data(parsed_path), baseline; label="PSS/E RAW path", demand_rtol=1e-5) end diff --git a/test/test_parser_parity.jl b/test/test_parser_parity.jl index 4c6b3c6..453cc0f 100644 --- a/test/test_parser_parity.jl +++ b/test/test_parser_parity.jl @@ -12,9 +12,11 @@ mpc.areas = [1 1]; mpc.bus_name = ['one'; 'two']; """ -# `parse_file`/`parse_matpower` return a PowerIO.BalancedNetwork; `_network_data` applies -# PowerDiff's normalization (per-unit via PowerIO, cost right-align, rate_a fallback, -# angle defaults, storage/HVDC rejection) into the tables the constructors consume. +# `parse_file`/`parse_matpower` return a `PowerIO.PioModule{PowerIO.BalancedNetwork}`. +# `_network_data` reads `PowerIO.to_powerdata` off it -- per unit, radian angles, +# rescaled costs, a source row number and a status on every row -- selects the +# elements PowerDiff models, and applies the modeling PowerIO leaves to its consumer +# (cost interpretation, the rate_a fallback, angle defaults, storage/HVDC refusal). _inline_data() = PowerDiff._network_data(PowerDiff.parse_matpower(IOBuffer(_INLINE_CASE))) @testset "MATPOWER Parser Semantics" begin @@ -33,10 +35,14 @@ _inline_data() = PowerDiff._network_data(PowerDiff.parse_matpower(IOBuffer(_INLI # Loads and shunts are aggregated into per-bus values. @test data.bus[1].pd == 0.5 @test data.bus[1].gs == 0.01 - # Shunts are also re-exposed as a per-bus table (bus 1: Gs=1, Bs=-2 -> 0.01, -0.02 pu). - @test length(data.shunt) == 1 - @test data.shunt[1].shunt_bus == 1 - @test data.shunt[1].gs ≈ 0.01 && data.shunt[1].bs ≈ -0.02 + @test data.bus[1].bs ≈ -0.02 + @test data.bus[2].gs == 0.0 && data.bus[2].bs == 0.0 + # The series admittance and each terminal's charging admittance come from + # PowerIO; `br_r`/`br_x` are carried alongside for the rate_a fallback. + @test data.branch[1].g ≈ 0.02 / (0.02^2 + 0.2^2) + @test data.branch[1].b ≈ -0.2 / (0.02^2 + 0.2^2) + @test data.branch[1].b_fr ≈ 0.005 && data.branch[1].b_to ≈ 0.005 + @test data.branch[1].g_fr == 0.0 && data.branch[1].g_to == 0.0 @test data.branch[1].tap == 1.0 @test data.branch[1].rate_a > 0 @test data.branch[1].angmin ≈ -π / 3 @@ -46,7 +52,11 @@ _inline_data() = PowerDiff._network_data(PowerDiff.parse_matpower(IOBuffer(_INLI @testset "Multiline arrays and artifact path" begin parsed = PowerDiff.parse_file("pglib_opf_case14_ieee.m"; library=:pglib) - @test parsed isa PowerIO.BalancedNetwork + @test parsed isa PowerIO.PioModule{PowerIO.BalancedNetwork} + @test parsed.sources[1].format == "matpower" + # Reader findings travel with the module as records, not log output. + @test parsed.diagnostics isa Vector{PowerIO.Diagnostic} + @test !any(d -> d.severity === :error, parsed.diagnostics) nd = PowerDiff._network_data(parsed) @test length(nd.bus) == 14 @test length(nd.branch) == 20 @@ -54,28 +64,62 @@ _inline_data() = PowerDiff._network_data(PowerDiff.parse_matpower(IOBuffer(_INLI end @testset "Rejected inputs" begin - @test_throws ArgumentError PowerDiff.parse_file("case.json") + # PowerDiff's own refusals: a container token that names no reader, an + # unknown keyword, an unknown artifact. @test_throws ArgumentError PowerDiff.parse_file(IOBuffer(_INLINE_CASE); filetype="json") @test_throws ArgumentError PowerDiff.parse_file(IOBuffer(_INLINE_CASE); unsupported=true) @test_throws ArgumentError PowerDiff.get_path(:unknown) - # Modeling-level rejections happen when the network tables are built. - unsupported = replace(_INLINE_CASE, "mpc.areas = [1 1];" => "mpc.storage = [1 1];") + # PowerIO's refusals reach the caller as PowerIOError, carrying the code and + # the diagnostics that say why. Wrapping them would discard both. + missing_file = joinpath(@__DIR__, "no_such_case.m") + @test_throws PowerIO.PowerIOError PowerDiff.parse_file(missing_file) + err = try + PowerDiff.parse_file(missing_file) + catch e + e + end + @test err.code == "READ.IO.OPEN" + + # Modeling-level rejections happen when the network tables are built, on a + # case PowerIO parses without complaint. + storage_row = "mpc.storage = [1 0 0 0 100 50 50 0.9 0.9 100 -50 50 0 0 0 0 1];" + unsupported = replace(_INLINE_CASE, "mpc.areas = [1 1];" => storage_row) @test_throws ArgumentError PowerDiff._network_data(PowerDiff.parse_matpower(IOBuffer(unsupported))) invalid = replace(_INLINE_CASE, "0.01 0.1" => "NaN 0.1") @test_throws ArgumentError PowerDiff._network_data(PowerDiff.parse_matpower(IOBuffer(invalid))) - pwl = replace(_INLINE_CASE, "2 0 0 3 0.01 2 3" => "1 0 0 3 0.01 2 3") + # A well-formed piecewise linear cost: PowerDiff models polynomial costs, and + # reads the model off the generator rather than inferring it from the row. + pwl = replace(_INLINE_CASE, + "mpc.gencost = [2 0 0 3 0.01 2 3; 2 0 0 3 0.02 3 4];" => + "mpc.gencost = [1 0 0 3 0 0 50 500 100 1200; 1 0 0 3 0 0 25 300 50 700];") @test_throws ArgumentError PowerDiff._network_data(PowerDiff.parse_matpower(IOBuffer(pwl))) - quartic = replace(_INLINE_CASE, "2 0 0 3 0.01 2 3" => "2 0 0 4 1 0.01 2 3") + quartic = replace(_INLINE_CASE, + "mpc.gencost = [2 0 0 3 0.01 2 3; 2 0 0 3 0.02 3 4];" => + "mpc.gencost = [2 0 0 4 1 0.01 2 3; 2 0 0 4 1 0.02 3 4];") @test_throws ArgumentError PowerDiff._network_data(PowerDiff.parse_matpower(IOBuffer(quartic))) end @testset "Parser contract" begin - @test PowerDiff.parse_file(IOBuffer(_INLINE_CASE)) isa PowerIO.BalancedNetwork + @test PowerDiff.parse_file(IOBuffer(_INLINE_CASE)) isa PowerIO.PioModule{PowerIO.BalancedNetwork} @test_throws ArgumentError PowerDiff.parse_file(IOBuffer(_INLINE_CASE); backend=:native) + + # PowerIO owns the format vocabulary: a token PowerDiff has never named + # passes through for PowerIO to resolve or refuse, so a reader PowerIO gains + # is reachable without a PowerDiff release. The historical short spellings + # still resolve to the tokens they always meant. + @test PowerDiff._format_token(:m) == "matpower" + @test PowerDiff._format_token(:raw) == "psse" + @test PowerDiff._format_token(:pm) == "powermodels-json" + @test PowerDiff._format_token(:egret) == "egret-json" + for token in (:xiidm, :cgmes, :pandapower, :gridfm, Symbol("pypsa-csv")) + @test PowerDiff._format_token(token) == lowercase(String(token)) + end + @test_throws ArgumentError PowerDiff._format_token(:json) + @test_throws ArgumentError PowerDiff._format_token("") end end @@ -100,7 +144,10 @@ function _assert_netdata_equal(a, b, label) end for (x, y) in zip(a.branch, b.branch) @test x.f_bus == y.f_bus && x.t_bus == y.t_bus - @test x.br_r ≈ y.br_r && x.br_x ≈ y.br_x && x.br_b ≈ y.br_b + @test x.br_r ≈ y.br_r && x.br_x ≈ y.br_x + @test x.g ≈ y.g && x.b ≈ y.b + @test x.g_fr ≈ y.g_fr && x.b_fr ≈ y.b_fr + @test x.g_to ≈ y.g_to && x.b_to ≈ y.b_to @test x.rate_a ≈ y.rate_a @test x.tap ≈ y.tap && x.shift ≈ y.shift && x.angmin ≈ y.angmin && x.angmax ≈ y.angmax end @@ -110,19 +157,16 @@ end @testset "PowerIO parser path and IO parity" begin # PowerIO is the only parser/data layer. Path parsing and IO parsing must land on # the same PowerDiff network tables after normalization. - if !PowerIO.library_available() - @info "libpowerio_capi not found (set POWERIO_CAPI to a local build); skipping parser parity" - @test_skip false - else - cases = filter(c -> isfile(joinpath(PD_PGLIB_DIR, c)), - ["pglib_opf_case5_pjm.m", "pglib_opf_case14_ieee.m", "pglib_opf_case30_ieee.m"]) - @test !isempty(cases) - for c in cases - path_data = PowerDiff._network_data(PowerDiff.parse_file(c; library=:pglib)) - io_data = PowerDiff._network_data( - PowerDiff.parse_file(IOBuffer(read(joinpath(PD_PGLIB_DIR, c), String)))) - _assert_netdata_equal(path_data, io_data, c) - end + # PowerIO resolves libpowerio_capi from a bundled artifact, so this path is + # reachable wherever the rest of the suite is. + cases = filter(c -> isfile(joinpath(PD_PGLIB_DIR, c)), + ["pglib_opf_case5_pjm.m", "pglib_opf_case14_ieee.m", "pglib_opf_case30_ieee.m"]) + @test !isempty(cases) + for c in cases + path_data = PowerDiff._network_data(PowerDiff.parse_file(c; library=:pglib)) + io_data = PowerDiff._network_data( + PowerDiff.parse_file(IOBuffer(read(joinpath(PD_PGLIB_DIR, c), String)))) + _assert_netdata_equal(path_data, io_data, c) end end diff --git a/test/test_powerio_integration.jl b/test/test_powerio_integration.jl new file mode 100644 index 0000000..88712e7 --- /dev/null +++ b/test/test_powerio_integration.jl @@ -0,0 +1,292 @@ +# Copyright 2026 Samuel Talkington and contributors +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# The PowerDiff/PowerIO seam: format routing, the one-pass ingest and its findings, +# and the bounds a case may leave unstated. + +import JuMP +import PowerIO + +@testset "PowerIO format routing" begin + # PowerIO owns the format vocabulary. PowerDiff normalizes its own historical + # short spellings and hands everything else straight over, so a reader PowerIO + # ships is reachable without a PowerDiff release. + @test PowerDiff._format_token(:m) == "matpower" + @test PowerDiff._format_token(".m") == "matpower" + @test PowerDiff._format_token(:raw) == "psse" + @test PowerDiff._format_token(:aux) == "powerworld" + @test PowerDiff._format_token(:pm) == "powermodels-json" + @test PowerDiff._format_token(:powermodels) == "powermodels-json" + @test PowerDiff._format_token(:egret) == "egret-json" + + # Tokens PowerDiff has never named: passed through for PowerIO to resolve, not + # refused here. This is the property that keeps the two in step, and it covers + # the distribution readers too -- what refuses those is the value they parse to. + for token in (:psse34, :psse35, :pandapower, :pypsa, :pslf, :epc, :pwb, + :gridfm, :goc3, :surge, :opfdata, :xiidm, :cgmes, :dss, :pmd, + Symbol("pypsa-csv")) + @test PowerDiff._format_token(token) == lowercase(String(token)) + end + + # Two hints must agree, and the alias layer is what lets them. + @test PowerDiff._powerio_format_hint(:matpower, ".m") == "matpower" + @test_throws ArgumentError PowerDiff._powerio_format_hint(:matpower, :psse) + + # A bare `json` names a container, not a reader. + @test_throws ArgumentError PowerDiff._format_token(:json) + @test_throws ArgumentError PowerDiff._format_token(".json") + + @test_throws ArgumentError PowerDiff._format_token("") +end + +@testset "Only a balanced transmission network is modeled" begin + # PowerIO returns one of twenty value kinds, and only a balanced transmission + # network has what PowerDiff differentiates. The refusal is the type of what was + # parsed, so one rule covers distribution cases, series carriers and calculation + # instances alike, with no token list here to trail what PowerIO reads. + net = PowerDiff.parse_file("pglib_opf_case14_ieee.m"; library=:pglib) + @test PowerDiff._require_balanced(net) === net + + dss = """ + New Circuit.mini basekv=12.47 pu=1.0 phases=3 bus1=sourcebus + New Linecode.lc1 nphases=3 r1=0.1 x1=0.2 c1=0 units=km + New Line.l1 bus1=sourcebus bus2=b1 linecode=lc1 length=1 units=km + New Load.ld1 bus1=b1 phases=3 kV=12.47 kW=100 kvar=30 + Set voltagebases=[12.47] + Calcvoltagebases + Solve + """ + distribution = PowerIO.parse(IOBuffer(dss); format="dss", name="mini.dss") + @test distribution isa PowerIO.PioModule{PowerIO.MulticonductorNetwork} + + err = try + PowerDiff._require_balanced(distribution) + nothing + catch caught + caught + end + @test err isa ArgumentError + @test occursin("balanced transmission networks", sprint(showerror, err)) + @test occursin("MulticonductorNetwork", sprint(showerror, err)) + + @test_throws ArgumentError PowerDiff.parse_file(IOBuffer(dss); from=:dss) +end + +@testset "PowerIO states the network, PowerDiff selects from it" begin + net = PowerDiff.parse_file("pglib_opf_case14_ieee.m"; library=:pglib) + tables = PowerDiff._network_data(net) + + # Both network types read the same tables, so they cannot describe different + # cases, and both accept the module or the network inside it. + dc = DCNetwork(net) + ac = ACNetwork(net) + @test dc.n == ac.n + @test dc.m == ac.m + @test dc.id_map.bus_ids == ac.id_map.bus_ids + @test dc.id_map.branch_ids == ac.id_map.branch_ids + @test DCNetwork(net.value).id_map.bus_ids == dc.id_map.bus_ids + + # `to_powerdata` is unfiltered: row `i` is the source row number and `status` + # says whether the row is in service, so an out-of-service row leaves a gap + # rather than renumbering the rows that remain. + raw = PowerIO.to_powerdata(net) + @test [Int(g.i) for g in raw.gen if Int(g.status) != 0] == [g.index for g in tables.gen] + @test [Int(b.i) for b in raw.branch if Int(b.status) != 0] == [br.index for br in tables.branch] + + # The series admittance is read off PowerIO's terminal coefficients rather than + # derived, and the incidence matrix labels the same edges the tables name. + for (l, br) in enumerate(tables.branch) + @test br.g ≈ Float64(raw.branch[br.index].c7) - br.g_to + @test br.b ≈ Float64(raw.branch[br.index].c8) - br.b_to + @test dc.A[l, dc.id_map.bus_to_idx[br.f_bus]] == 1.0 + @test dc.A[l, dc.id_map.bus_to_idx[br.t_bus]] == -1.0 + end + @test dc.b == [br.b for br in tables.branch] + @test ac.g == [br.g for br in tables.branch] + + # Reader findings are records on the module, not log output and not a wrapper. + @test net.diagnostics isa Vector{PowerIO.Diagnostic} + @test all(d -> d.code isa String && !isempty(d.code), net.diagnostics) + @test net.sources[1].format == "matpower" + @test net.producer.name isa String +end + +@testset "A source that states no thermal limit" begin + # MATPOWER spells an absent rating `0`, and PowerIO carries that out of + # `to_powerdata` as `Inf`. Both mean the same thing, both take the synthesized + # limit, and neither may reach the solver as an unbounded flow. + template = """ + function mpc = case_rate + mpc.version = '2'; + mpc.baseMVA = 100; + mpc.bus = [1 3 0 0 0 0 1 1.0 0 230 1 1.1 0.9; 2 1 50 10 0 0 1 1.0 0 230 1 1.1 0.9]; + mpc.gen = [1 60 0 100 -100 1 100 1 150 0]; + mpc.branch = [1 2 0.01 0.1 0.02 %RATE% %RATE% %RATE% 0 0 1 -60 60]; + mpc.gencost = [2 0 0 3 0.01 2 3]; + """ + absent = PowerDiff.parse_matpower(IOBuffer(replace(template, "%RATE%" => "0"))) + huge = PowerDiff.parse_matpower(IOBuffer(replace(template, "%RATE%" => "1e6"))) + + nd = PowerDiff._network_data(absent) + @test all(br -> isfinite(br.rate_a) && br.rate_a > 0, nd.branch) + @test all(isfinite, DCNetwork(absent).fmax) + + a = solve!(DCOPFProblem(DCNetwork(absent))) + h = solve!(DCOPFProblem(DCNetwork(huge))) + @test a.pg ≈ h.pg atol=1e-6 + @test a.f ≈ h.f atol=1e-6 + @test a.psh ≈ h.psh atol=1e-6 +end + +@testset "Isolated buses do not reach the model" begin + # `to_powerdata` no longer drops a `type == 4` bus, so PowerDiff drops it and + # every branch and generator standing on it. Bus 3 sits in the middle of the bus + # table, so a missed re-index would move branch 2 onto the wrong bus. + isolated = """ + function mpc = case_isolated + mpc.version = '2'; + mpc.baseMVA = 100; + mpc.bus = [1 3 0 0 0 0 1 1.0 0 230 1 1.1 0.9; 2 1 30 5 0 0 1 1.0 0 230 1 1.1 0.9; 3 4 25 5 0 0 1 1.0 0 230 1 1.1 0.9; 4 1 20 5 0 0 1 1.0 0 230 1 1.1 0.9]; + mpc.gen = [1 60 0 100 -100 1 100 1 150 0]; + mpc.branch = [1 2 0.01 0.1 0.02 100 100 100 0 0 1 -60 60; 2 4 0.02 0.2 0.01 100 100 100 0 0 1 -60 60; 2 3 0.03 0.3 0.01 100 100 100 0 0 1 -60 60]; + mpc.gencost = [2 0 0 3 0.01 2 3]; + """ + nd = PowerDiff._network_data(PowerDiff.parse_matpower(IOBuffer(isolated))) + @test [b.bus_i for b in nd.bus] == [1, 2, 4] + @test all(b -> b.bus_type != 4, nd.bus) + # The branch to the isolated bus goes with it; the survivors keep their source + # row numbers and land on the right buses. + @test [br.index for br in nd.branch] == [1, 2] + @test (nd.branch[2].f_bus, nd.branch[2].t_bus) == (2, 4) + # Demand at an isolated bus leaves the model with the bus. + @test sum(b.pd for b in nd.bus) ≈ 0.50 +end + +@testset "Terminal charging is carried, not averaged" begin + # Each terminal states its own charging admittance. Summing the two and splitting + # the total evenly, as PowerDiff once did, loses an asymmetric source's fidelity. + net = ACNetwork(pd_case([pd_bus(1, 3), pd_bus(2, 1; pd=0.2)], + [pd_gen(1, 1; pmax=2.0, cost=(1.0, 1.0, 0.0))], + [pd_branch(1, 1, 2; br_r=0.01, br_x=0.1, rate_a=2.0, + g_fr=0.001, b_fr=0.004, g_to=0.002, b_to=0.012)])) + @test net.b_fr == [0.004] && net.b_to == [0.012] + @test net.g_fr == [0.001] && net.g_to == [0.002] + Y = Matrix(admittance_matrix(net)) + y = inv(complex(0.01, 0.1)) + @test Y[1, 1] ≈ y + complex(0.001, 0.004) + @test Y[2, 2] ≈ y + complex(0.002, 0.012) +end + +@testset "Absent reactive limits" begin + # PowerIO 0.9 carries a bound the case does not state as ±Inf rather than refusing + # the case (stock case9241pegase leaves the reactive limits off seven generators). + # PowerDiff leaves the bound off the solver model and pins its multiplier to zero, + # which is what a solver reports for a bound it was never given. + buses = [pd_bus(1, 3; vmax=1.1, vmin=0.9), pd_bus(2, 1; pd=0.3, qd=0.1, vmax=1.1, vmin=0.9)] + gens = [pd_gen(1, 1; pg=0.3, qmin=-Inf, qmax=Inf, pmax=2.0, pmin=0.0, cost=(0.0, 1.0, 0.0))] + branches = [pd_branch(1, 1, 2; br_r=0.01, br_x=0.1, rate_a=2.0)] + data = pd_case(buses, gens, branches; name="absent_q_limits") + + net = ACNetwork(data) + @test net.qmin == [-Inf] + @test net.qmax == [Inf] + @test PowerDiff._absent_bound(net.qmin[1]) + @test PowerDiff._absent_bound(net.qmax[1]) + @test !PowerDiff._absent_bound(net.pmax[1]) + + prob = ACOPFProblem(net; silent=true) + @test !JuMP.has_lower_bound(prob.qg[1]) + @test !JuMP.has_upper_bound(prob.qg[1]) + @test JuMP.has_lower_bound(prob.pg[1]) + + sol = solve!(prob) + @test all(isfinite, sol.qg) + # An absent bound has no multiplier. Zero is the value, not a placeholder. + @test sol.rho_qg_lb == [0.0] + @test sol.rho_qg_ub == [0.0] + + # The KKT residual and Jacobian stay finite: the two rows read `ρ = 0`. + idx = kkt_indices(prob) + K = kkt(flatten_variables(sol, prob), prob) + @test all(isfinite, K) + @test K[idx.rho_qg_lb[1]] == sol.rho_qg_lb[1] + @test K[idx.rho_qg_ub[1]] == sol.rho_qg_ub[1] + + J = calc_kkt_jacobian(prob; sol=sol) + @test all(isfinite, nonzeros(J)) + @test J[idx.rho_qg_lb[1], idx.rho_qg_lb[1]] == 1.0 + @test J[idx.rho_qg_ub[1], idx.rho_qg_ub[1]] == 1.0 + @test J[idx.rho_qg_lb[1], idx.qg[1]] == 0.0 + @test J[idx.rho_qg_ub[1], idx.qg[1]] == 0.0 + + # A stated bound is untouched by any of this. + bounded = ACNetwork(pd_case(buses, + [pd_gen(1, 1; pg=0.3, qmin=-1.0, qmax=1.0, pmax=2.0, pmin=0.0, cost=(0.0, 1.0, 0.0))], + branches; name="stated_q_limits")) + bprob = ACOPFProblem(bounded; silent=true) + @test JuMP.has_lower_bound(bprob.qg[1]) + bsol = solve!(bprob) + bidx = kkt_indices(bprob) + bJ = calc_kkt_jacobian(bprob; sol=bsol) + @test bJ[bidx.rho_qg_lb[1], bidx.rho_qg_lb[1]] ≈ bsol.qg[1] - bounded.qmin[1] + @test all(isfinite, kkt(flatten_variables(bsol, bprob), bprob)) +end + +@testset "Non-finite values PowerDiff cannot model" begin + # Only the reactive limits model an absent bound. Everywhere else a non-finite + # value is named, with the row and the field, rather than reaching a factorization. + buses = [pd_bus(1, 3), pd_bus(2, 1; pd=0.3)] + branches = [pd_branch(1, 1, 2; br_r=0.01, br_x=0.1, rate_a=2.0)] + + for (field, gen) in ( + (:pmax, pd_gen(1, 1; pmax=Inf, pmin=0.0, cost=(0.0, 1.0, 0.0))), + (:pmin, pd_gen(1, 1; pmax=2.0, pmin=-Inf, cost=(0.0, 1.0, 0.0))), + ) + err = try + ACNetwork(pd_case(buses, [gen], branches)) + nothing + catch caught + caught + end + @test err isa ArgumentError + @test occursin(String(field), sprint(showerror, err)) + end + + err = try + ACNetwork(pd_case([pd_bus(1, 3; vmax=Inf), pd_bus(2, 1; pd=0.3)], + [pd_gen(1, 1; pmax=2.0, pmin=0.0, cost=(0.0, 1.0, 0.0))], + branches)) + nothing + catch caught + caught + end + @test err isa ArgumentError + @test occursin("vmax", sprint(showerror, err)) +end + +@testset "The ingest still reproduces MATPOWER's network" begin + # An oracle that shares no code with the ingest: PowerModels parses the same + # files independently, and its basic-network PTDF is a function of the topology, + # the branch susceptances and the reference bus alone. Agreement pins all three. + PowerModels.silence() + for c in ("pglib_opf_case14_ieee.m", "pglib_opf_case30_ieee.m", "pglib_opf_case118_ieee.m") + path = joinpath(PD_PGLIB_DIR, c) + isfile(path) || continue + net = PowerDiff.parse_file(path) + state = DCPowerFlowState(DCNetwork(net), calc_demand_vector(net)) + pm = PowerModels.calc_basic_ptdf_matrix( + PowerModels.make_basic_network(PowerModels.parse_file(path))) + @test maximum(abs, ptdf_matrix(state) .- pm) < 1e-8 + end +end