diff --git a/.github/workflows/Benchmark.yml b/.github/workflows/Benchmark.yml index b27acb3..05525c3 100644 --- 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 new file mode 100644 index 0000000..2f0c6ec --- /dev/null +++ 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 index 36cf030..e9e48ed 100644 --- 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 index a90cc43..45c5e2f 100644 --- 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 index 655f585..b9a1071 100644 --- 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 index f669d0d..dd362f6 100644 --- 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 new file mode 100644 index 0000000..74dfb21 --- /dev/null +++ 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 index 090d02e..3545400 100644 --- 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 index e2c2892..6ad29c5 100644 --- 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 index 97c9bea..8a85bf8 100644 --- 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 index b5009f8..3353838 100644 --- 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 index d50b56c..b4ecedf 100644 --- 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 index 16a261d..aa6d63a 100644 --- 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 index 958c509..cf10b45 100644 --- 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 index d2ea7a1..d439ccb 100644 --- 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 diff --git a/src/PowerDiff.jl b/src/PowerDiff.jl index 508493d..703d1aa 100644 --- a/src/PowerDiff.jl +++ b/src/PowerDiff.jl @@ -18,8 +18,9 @@ using LinearAlgebra using SparseArrays using JuMP using Ipopt -import PowerModels -const PM = PowerModels +using ExaModels +using NLPModelsIpopt +using PowerIO const MOI = JuMP.MOI @@ -28,13 +29,14 @@ const MOI = JuMP.MOI # ============================================================================= const _SILENCE_WARNINGS = Ref(false) +include("artifacts.jl") + """ silence() Suppress all warning messages from PowerDiff for the rest of the session. -Warnings from other packages (PowerModels, JuMP, Ipopt, etc.) are not affected. -To suppress PowerModels output, also call `PowerModels.silence()`. +Warnings from other packages (JuMP, Ipopt, etc.) are not affected. """ function silence() _SILENCE_WARNINGS[] = true @@ -100,6 +102,7 @@ export calc_sensitivity, calc_sensitivity_column export Sensitivity, silence export operand_symbols, parameter_symbols export jvp, vjp, jvp!, vjp!, dict_to_vec, vec_to_dict, kkt_dims +export parse_file, parse_matpower, parse_matpower_struct, get_path # DC Power Flow Types export DCNetwork, DCPowerFlowState diff --git a/src/artifacts.jl b/src/artifacts.jl new file mode 100644 index 0000000..4d32077 --- /dev/null +++ b/src/artifacts.jl @@ -0,0 +1,12 @@ +using LazyArtifacts + +""" + get_path(library::Symbol) + +Resolve an artifact-backed case library bundled with PowerDiff. Currently only +`:pglib` (PGLib-OPF) is available. +""" +function get_path(library::Symbol) + library == :pglib && return joinpath(artifact"PGLib_opf", "pglib-opf-23.07") + throw(ArgumentError("unsupported library $library")) +end diff --git a/src/prob/ac_opf_solve.jl b/src/prob/ac_opf_solve.jl index 279f69a..6d1f96e 100644 --- a/src/prob/ac_opf_solve.jl +++ b/src/prob/ac_opf_solve.jl @@ -18,6 +18,31 @@ # # Functions for solving AC OPF problems and updating parameters. +""" + _check_solve_status(stats, label::String) + +Check the NLPModelsIpopt solve status and throw informative errors for common failure modes. +""" +function _check_solve_status(stats, label::String) + status = getproperty(stats, :status) + status == :first_order && return status + if status == :acceptable + # Solved only to the looser acceptable tolerance, not first-order optimality. + # The duals feed calc_sensitivity, so a nonzero KKT residual here degrades sensitivities. + _SILENCE_WARNINGS[] || @warn "$label converged only to acceptable tolerance, not first-order optimality; duals may carry a nonzero KKT residual and degrade sensitivities" + return status + end + if status == :infeasible + error("$label is infeasible. Check that demand is feasible given generator capacities and network constraints.") + elseif status == :unbounded + error("$label is unbounded. Check cost coefficients and variable bounds.") + elseif status in (:max_iter, :max_time) + error("$label solver reached $status. Try increasing solver limits or simplifying the problem.") + else + error("$label failed with status: $status") + end +end + """ solve!(prob::ACOPFProblem) @@ -31,7 +56,7 @@ ACOPFSolution containing optimal primal and dual variables. # Throws Error if optimization does not converge to optimal/locally optimal solution. """ -function solve!(prob::ACOPFProblem) +function solve!(prob::ACOPFProblem{JuMPBackend}) # Invalidate sensitivity cache since we're re-solving invalidate!(prob.cache) @@ -47,31 +72,49 @@ function solve!(prob::ACOPFProblem) return sol end +function solve!(prob::ACOPFProblem{ExaBackend}) + # Invalidate sensitivity cache since we're re-solving + invalidate!(prob.cache) + + # Match the JuMP backend's Ipopt tolerance so backend correspondence tests + # compare like-for-like solver targets. + result = NLPModelsIpopt.ipopt(prob.model; print_level = prob._silent ? 0 : 5, tol=1e-6) + + _check_solve_status(result, "AC OPF") + + sol = _extract_ac_opf_solution(prob, result) + + # Cache the solution for sensitivity computations + prob.cache.solution = sol + + return sol +end + """ Extract solution from solved AC OPF problem. """ -function _extract_ac_opf_solution(prob::ACOPFProblem) +function _extract_ac_opf_solution(prob::ACOPFProblem{JuMPBackend}) n = prob.network.n m = prob.network.m k = prob.n_gen - ref = prob.ref # Extract primal variables va_val = value.(prob.va) vm_val = value.(prob.vm) pg_val = value.(prob.pg) qg_val = value.(prob.qg) + p_arr = value.(prob.p) + q_arr = value.(prob.q) - p_val = Dict(arc => value(prob.p[arc]) for arc in keys(prob.p)) - q_val = Dict(arc => value(prob.q[arc]) for arc in keys(prob.q)) + p_val = Dict(prob.data.arcs[i] => p_arr[i] for i in eachindex(prob.data.arcs)) + q_val = Dict(prob.data.arcs[i] => q_arr[i] for i in eachindex(prob.data.arcs)) # Extract dual variables - power balance (equality) ν_p_bal = [dual(prob.cons.p_bal[i]) for i in 1:n] ν_q_bal = [dual(prob.cons.q_bal[i]) for i in 1:n] # Extract dual variables - reference bus (equality) - ref_bus_keys = sort(collect(keys(ref[:ref_buses]))) - ν_ref_bus = [dual(prob.cons.ref_bus[i]) for i in ref_bus_keys] + ν_ref_bus = [dual(prob.cons.ref_bus[i]) for i in eachindex(prob.cons.ref_bus)] # Extract dual variables - flow definition equations (equality) ν_p_fr = [dual(prob.cons.p_fr[l]) for l in 1:m] @@ -84,8 +127,8 @@ function _extract_ac_opf_solution(prob::ACOPFProblem) λ_thermal_to = [dual(prob.cons.thermal_to[l]) for l in 1:m] # Extract dual variables - angle difference limits (inequality) - λ_angle_lb = [dual(prob.cons.angle_diff[l][1]) for l in 1:m] - λ_angle_ub = [dual(prob.cons.angle_diff[l][2]) for l in 1:m] + λ_angle_lb = [dual(prob.cons.angle_diff_lb[l]) for l in 1:m] + λ_angle_ub = [dual(prob.cons.angle_diff_ub[l]) for l in 1:m] # Extract dual variables - voltage bounds (inequality) μ_vm_lb = [dual(LowerBoundRef(prob.vm[i])) for i in 1:n] @@ -108,12 +151,8 @@ function _extract_ac_opf_solution(prob::ACOPFProblem) σ_q_to_ub = zeros(m) for l in 1:m - branch = ref[:branch][l] - f_bus = branch["f_bus"] - t_bus = branch["t_bus"] - f_idx = (l, f_bus, t_bus) - t_idx = (l, t_bus, f_bus) - + f_idx = prob.data.arc_from_idx[l] + t_idx = prob.data.arc_to_idx[l] σ_p_fr_lb[l] = dual(LowerBoundRef(prob.p[f_idx])) σ_p_fr_ub[l] = dual(UpperBoundRef(prob.p[f_idx])) σ_q_fr_lb[l] = dual(LowerBoundRef(prob.q[f_idx])) @@ -145,6 +184,71 @@ function _extract_ac_opf_solution(prob::ACOPFProblem) ) end +function _extract_ac_opf_solution(prob::ACOPFProblem{ExaBackend}, result) + m = prob.network.m + + va_val = ExaModels.solution(result, prob.va) + vm_val = ExaModels.solution(result, prob.vm) + pg_val = ExaModels.solution(result, prob.pg) + qg_val = ExaModels.solution(result, prob.qg) + p_arr = ExaModels.solution(result, prob.p) + q_arr = ExaModels.solution(result, prob.q) + + p_val = Dict(prob.data.arcs[i] => p_arr[i] for i in eachindex(prob.data.arcs)) + q_val = Dict(prob.data.arcs[i] => q_arr[i] for i in eachindex(prob.data.arcs)) + + # ExaModels reports minimization-constraint multipliers with the opposite + # sign of JuMP's duals; ACOPFSolution uses the JuMP/KKT convention. + ν_p_bal = -ExaModels.multipliers(result, prob.cons.p_bal) + ν_q_bal = -ExaModels.multipliers(result, prob.cons.q_bal) + ν_ref_bus = -ExaModels.multipliers(result, prob.cons.ref_bus) + ν_p_fr = -ExaModels.multipliers(result, prob.cons.p_fr) + ν_p_to = -ExaModels.multipliers(result, prob.cons.p_to) + ν_q_fr = -ExaModels.multipliers(result, prob.cons.q_fr) + ν_q_to = -ExaModels.multipliers(result, prob.cons.q_to) + λ_thermal_fr = -ExaModels.multipliers(result, prob.cons.thermal_fr) + λ_thermal_to = -ExaModels.multipliers(result, prob.cons.thermal_to) + λ_angle_lb = -ExaModels.multipliers(result, prob.cons.angle_diff_lb) + λ_angle_ub = -ExaModels.multipliers(result, prob.cons.angle_diff_ub) + μ_vm_lb = ExaModels.multipliers_L(result, prob.vm) + μ_vm_ub = -ExaModels.multipliers_U(result, prob.vm) + ρ_pg_lb = ExaModels.multipliers_L(result, prob.pg) + ρ_pg_ub = -ExaModels.multipliers_U(result, prob.pg) + ρ_qg_lb = ExaModels.multipliers_L(result, prob.qg) + ρ_qg_ub = -ExaModels.multipliers_U(result, prob.qg) + + p_lb = ExaModels.multipliers_L(result, prob.p) + p_ub = -ExaModels.multipliers_U(result, prob.p) + q_lb = ExaModels.multipliers_L(result, prob.q) + q_ub = -ExaModels.multipliers_U(result, prob.q) + σ_p_fr_lb = p_lb[prob.data.arc_from_idx] + σ_p_fr_ub = p_ub[prob.data.arc_from_idx] + σ_q_fr_lb = q_lb[prob.data.arc_from_idx] + σ_q_fr_ub = q_ub[prob.data.arc_from_idx] + σ_p_to_lb = p_lb[prob.data.arc_to_idx] + σ_p_to_ub = p_ub[prob.data.arc_to_idx] + σ_q_to_lb = q_lb[prob.data.arc_to_idx] + σ_q_to_ub = q_ub[prob.data.arc_to_idx] + + return ACOPFSolution( + va = va_val, vm = vm_val, + pg = pg_val, qg = qg_val, + p = p_val, q = q_val, + nu_p_bal = ν_p_bal, nu_q_bal = ν_q_bal, + nu_ref_bus = ν_ref_bus, + nu_p_fr = ν_p_fr, nu_p_to = ν_p_to, nu_q_fr = ν_q_fr, nu_q_to = ν_q_to, + lam_thermal_fr = λ_thermal_fr, lam_thermal_to = λ_thermal_to, + lam_angle_lb = λ_angle_lb, lam_angle_ub = λ_angle_ub, + mu_vm_lb = μ_vm_lb, mu_vm_ub = μ_vm_ub, + rho_pg_lb = ρ_pg_lb, rho_pg_ub = ρ_pg_ub, rho_qg_lb = ρ_qg_lb, rho_qg_ub = ρ_qg_ub, + sig_p_fr_lb = σ_p_fr_lb, sig_p_fr_ub = σ_p_fr_ub, + sig_q_fr_lb = σ_q_fr_lb, sig_q_fr_ub = σ_q_fr_ub, + sig_p_to_lb = σ_p_to_lb, sig_p_to_ub = σ_p_to_ub, + sig_q_to_lb = σ_q_to_lb, sig_q_to_ub = σ_q_to_ub, + objective = result.objective + ) +end + """ update_switching!(prob::ACOPFProblem, sw::AbstractVector) @@ -166,8 +270,8 @@ function update_switching!(prob::ACOPFProblem, sw::AbstractVector) # Update network switching state prob.network.sw .= sw - # Rebuild JuMP model with new switching coefficients - _rebuild_jump_model!(prob) + # Rebuild the model with new switching coefficients + _rebuild_model!(prob) return prob end diff --git a/src/prob/dc_opf.jl b/src/prob/dc_opf.jl index 525bb50..50d31c9 100644 --- a/src/prob/dc_opf.jl +++ b/src/prob/dc_opf.jl @@ -28,11 +28,11 @@ const COMPLEMENTARITY_SNAP_TOL = 1e-6 Check JuMP model termination status and throw informative errors for common failure modes. """ -function _check_solve_status(model, label::String) +function _check_solve_status(model::JuMP.Model, label::String) status = termination_status(model) status in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) && return status if status == MOI.ALMOST_LOCALLY_SOLVED - @warn "$label converged at acceptable tolerance (ALMOST_LOCALLY_SOLVED)" + _SILENCE_WARNINGS[] || @warn "$label converged at acceptable tolerance (ALMOST_LOCALLY_SOLVED)" return status end if status == MOI.INFEASIBLE @@ -96,6 +96,7 @@ function solve!(prob::DCOPFProblem) μ_ub = -dual.(prob.cons.shed_ub) γ_lb = dual.(prob.cons.phase_diff_lb) γ_ub = -dual.(prob.cons.phase_diff_ub) + η_ref = dual(prob.cons.ref) # Post-process phase angle difference duals for strict complementarity. # Interior point solvers leave gamma ≈ 1e-8 for non binding constraints. @@ -151,7 +152,7 @@ function solve!(prob::DCOPFProblem) end B_r_factor = prob.cache.b_r_factor - sol = DCOPFSolution(θ_val, g_val, f_val, psh_val, ν_bal, ν_flow, λ_ub, λ_lb, ρ_ub, ρ_lb, μ_lb, μ_ub, γ_lb, γ_ub, obj, B_r_factor) + sol = DCOPFSolution(θ_val, g_val, f_val, psh_val, ν_bal, ν_flow, λ_ub, λ_lb, ρ_ub, ρ_lb, μ_lb, μ_ub, γ_lb, γ_ub, η_ref, obj, B_r_factor) # Cache the solution for sensitivity computations prob.cache.solution = sol diff --git a/src/prob/kkt_ac_opf.jl b/src/prob/kkt_ac_opf.jl index dad2a8f..be8a99e 100644 --- a/src/prob/kkt_ac_opf.jl +++ b/src/prob/kkt_ac_opf.jl @@ -29,7 +29,7 @@ # ============================================================================= """Return sorted reference bus indices for the problem.""" -_ref_bus_indices(prob::ACOPFProblem) = sort(collect(keys(prob.ref[:ref_buses]))) +_ref_bus_indices(prob::ACOPFProblem) = prob.data.ref_bus_keys # ============================================================================= # Dimension Calculations @@ -54,7 +54,7 @@ Total: 6n + 12m + 6k + n_ref """ function kkt_dims(prob::ACOPFProblem) n, m, k = prob.network.n, prob.network.m, prob.n_gen - n_ref = length(prob.ref[:ref_buses]) + n_ref = length(prob.data.ref_bus_keys) return 6n + 12m + 6k + n_ref end @@ -131,7 +131,7 @@ function kkt_indices(n::Int, m::Int, k::Int, n_ref::Int) end function kkt_indices(prob::ACOPFProblem) - n_ref = length(prob.ref[:ref_buses]) + n_ref = length(prob.data.ref_bus_keys) kkt_indices(prob.network.n, prob.network.m, prob.n_gen, n_ref) end @@ -198,13 +198,14 @@ end Compute all branch flows given voltage state and switching state. Returns vectors of p_fr, q_fr, p_to, q_to indexed by branch number. -Flow equations match PowerModels' `constraint_ohms_yt_from`/`constraint_ohms_yt_to` -(polar form with `calc_branch_y`/`calc_branch_t` decomposition, `tm = tap^2`). +Flow equations use the polar pi-model: series admittance `g + jb`, line-charging shunts +`g_fr/b_fr/g_to/b_to`, complex tap (`tap`/`shift`), and `tm = tap^2`. These match the AC OPF +flow constraints assembled by the solver backends. The switching variable sw_l multiplies each flow, so sw_l=0 means the branch contributes zero flow (open), sw_l=1 means full flow (closed). """ -function _compute_branch_flows(va, vm, net::ACNetwork, ref, sw; constants=nothing) +function _compute_branch_flows(va, vm, net::ACNetwork, sw; constants) m = net.m T = promote_type(eltype(va), eltype(vm), eltype(sw)) p_fr = zeros(T, m) @@ -213,30 +214,17 @@ function _compute_branch_flows(va, vm, net::ACNetwork, ref, sw; constants=nothin q_to = zeros(T, m) for l in 1:m - if isnothing(constants) - branch = ref[:branch][l] - f_bus = branch["f_bus"] - t_bus = branch["t_bus"] - g_br, b_br = PM.calc_branch_y(branch) - tr, ti = PM.calc_branch_t(branch) - g_fr_shunt = branch["g_fr"] - b_fr_shunt = branch["b_fr"] - g_to_shunt = branch["g_to"] - b_to_shunt = branch["b_to"] - tm = branch["tap"]^2 - else - f_bus = constants.f_bus[l] - t_bus = constants.t_bus[l] - g_br = constants.g_br[l] - b_br = constants.b_br[l] - tr = constants.tr[l] - ti = constants.ti[l] - g_fr_shunt = constants.g_fr[l] - b_fr_shunt = constants.b_fr[l] - g_to_shunt = constants.g_to[l] - b_to_shunt = constants.b_to[l] - tm = constants.tm[l] - end + f_bus = constants.f_bus[l] + t_bus = constants.t_bus[l] + g_br = constants.g_br[l] + b_br = constants.b_br[l] + tr = constants.tr[l] + ti = constants.ti[l] + g_fr_shunt = constants.g_fr[l] + b_fr_shunt = constants.b_fr[l] + g_to_shunt = constants.g_to[l] + b_to_shunt = constants.b_to[l] + tm = constants.tm[l] sw_l = sw[l] @@ -490,8 +478,8 @@ end Power balance residuals. """ function _power_balance_residuals(va, vm, pg, qg, p_fr, q_fr, p_to, q_to, - net::ACNetwork, ref, prob::ACOPFProblem; - pd=nothing, qd=nothing, constants=nothing) + net::ACNetwork, prob::ACOPFProblem; + pd=nothing, qd=nothing, constants) n = net.n m = net.m _et(x) = isnothing(x) ? Float64 : eltype(x) @@ -504,8 +492,8 @@ function _power_balance_residuals(va, vm, pg, qg, p_fr, q_fr, p_to, q_to, q_flow_sum = zeros(T, n) for l in 1:m - fb = isnothing(constants) ? ref[:branch][l]["f_bus"] : constants.f_bus[l] - tb = isnothing(constants) ? ref[:branch][l]["t_bus"] : constants.t_bus[l] + fb = constants.f_bus[l] + tb = constants.t_bus[l] p_flow_sum[fb] += p_fr[l] p_flow_sum[tb] += p_to[l] q_flow_sum[fb] += q_fr[l] @@ -516,17 +504,16 @@ function _power_balance_residuals(va, vm, pg, qg, p_fr, q_fr, p_to, q_to, pg_sum = zeros(T, n) qg_sum = zeros(T, n) for i in 1:prob.n_gen - bus_idx = isnothing(constants) ? ref[:gen][i]["gen_bus"] : constants.gen_bus[i] + bus_idx = constants.gen_bus[i] pg_sum[bus_idx] += pg[i] qg_sum[bus_idx] += qg[i] end for i in 1:n - gs_i = isnothing(constants) ? sum(ref[:shunt][s]["gs"] for s in ref[:bus_shunts][i]; init=0.0) : constants.gs[i] - bs_i = isnothing(constants) ? sum(ref[:shunt][s]["bs"] for s in ref[:bus_shunts][i]; init=0.0) : constants.bs[i] - - pd_i = isnothing(pd) ? sum(ref[:load][l]["pd"] for l in ref[:bus_loads][i]; init=0.0) : pd[i] - qd_i = isnothing(qd) ? sum(ref[:load][l]["qd"] for l in ref[:bus_loads][i]; init=0.0) : qd[i] + gs_i = constants.gs[i] + bs_i = constants.bs[i] + pd_i = isnothing(pd) ? constants.pd[i] : pd[i] + qd_i = isnothing(qd) ? constants.qd[i] : qd[i] K_p_bal[i] = p_flow_sum[i] + gs_i * vm[i]^2 - pg_sum[i] + pd_i K_q_bal[i] = q_flow_sum[i] - bs_i * vm[i]^2 - qg_sum[i] + qd_i @@ -918,7 +905,6 @@ function kkt(z::AbstractVector, prob::ACOPFProblem, sw::AbstractVector; end vars = unflatten_variables(z, idx) net = prob.network - ref = prob.ref n, m, k = net.n, net.m, prob.n_gen va, vm = vars.va, vars.vm @@ -928,12 +914,11 @@ function kkt(z::AbstractVector, prob::ACOPFProblem, sw::AbstractVector; T = promote_type(eltype(z), eltype(sw), _et(pd), _et(qd), _et(cq), _et(cl), _et(fmax)) # Compute branch flows as functions of voltages - p_fr, q_fr, p_to, q_to = _compute_branch_flows(va, vm, net, ref, sw; constants=constants) + p_fr, q_fr, p_to, q_to = _compute_branch_flows(va, vm, net, sw; constants=constants) - # Materialize rate_a vector once to avoid repeated isnothing checks - rate_a = isnothing(fmax) ? T[ref[:branch][l]["rate_a"] for l in 1:m] : fmax - cq_vec = isnothing(cq) ? T[ref[:gen][i]["cost"][1] for i in 1:k] : cq - cl_vec = isnothing(cl) ? T[ref[:gen][i]["cost"][2] for i in 1:k] : cl + rate_a = isnothing(fmax) ? T.(constants.fmax) : fmax + cq_vec = isnothing(cq) ? T.(constants.cq) : cq + cl_vec = isnothing(cl) ? T.(constants.cl) : cl # Pre-allocate KKT residual vector K = fill(T(NaN), last(idx.sig_q_to_ub)) @@ -949,13 +934,12 @@ function kkt(z::AbstractVector, prob::ACOPFProblem, sw::AbstractVector; # Power balance K_p_bal, K_q_bal = _power_balance_residuals(va, vm, pg, qg, p_fr, q_fr, p_to, q_to, - net, ref, prob; pd=pd, qd=qd, - constants=constants) + net, prob; pd=pd, qd=qd, constants=constants) K[idx.nu_p_bal] = K_p_bal K[idx.nu_q_bal] = K_q_bal # Reference bus: va[ref_bus] == 0 - rbk = isnothing(constants) ? _ref_bus_indices(prob) : constants.ref_bus_keys + rbk = constants.ref_bus_keys for (j, ref_bus_idx) in enumerate(rbk) K[idx.nu_ref_bus[j]] = va[ref_bus_idx] end @@ -967,30 +951,16 @@ function kkt(z::AbstractVector, prob::ACOPFProblem, sw::AbstractVector; # Upper bounds: L -= λ*(x - ub), CS = λ*(ub - x) = 0 (negated residual) # Both are valid; the sign flip cancels in implicit differentiation. - # Use pre-extracted bounds when available to avoid repeated ref lookups - if isnothing(constants) - vmin = T[ref[:bus][i]["vmin"] for i in 1:n] - vmax = T[ref[:bus][i]["vmax"] for i in 1:n] - pmin = T[ref[:gen][i]["pmin"] for i in 1:k] - pmax = T[ref[:gen][i]["pmax"] for i in 1:k] - qmin = T[ref[:gen][i]["qmin"] for i in 1:k] - qmax = T[ref[:gen][i]["qmax"] for i in 1:k] - f_bus_idx = [ref[:branch][l]["f_bus"] for l in 1:m] - t_bus_idx = [ref[:branch][l]["t_bus"] for l in 1:m] - angmin = T[ref[:branch][l]["angmin"] for l in 1:m] - angmax = T[ref[:branch][l]["angmax"] for l in 1:m] - else - vmin = constants.vmin - vmax = constants.vmax - pmin = constants.pmin - pmax = constants.pmax - qmin = constants.qmin - qmax = constants.qmax - f_bus_idx = constants.f_bus - t_bus_idx = constants.t_bus - angmin = constants.angmin - angmax = constants.angmax - end + vmin = constants.vmin + vmax = constants.vmax + pmin = constants.pmin + pmax = constants.pmax + qmin = constants.qmin + qmax = constants.qmax + f_bus_idx = constants.f_bus + t_bus_idx = constants.t_bus + angmin = constants.angmin + angmax = constants.angmax # Thermal limits K[idx.lam_thermal_fr] .= vars.lam_thermal_fr .* (p_fr.^2 .+ q_fr.^2 .- rate_a.^2) @@ -1030,32 +1000,30 @@ kkt(z::AbstractVector, prob::ACOPFProblem) = kkt(z, prob, prob.network.sw) # Parameter Extraction Functions # ============================================================================= +@inline function _require_kkt_constants(prob::ACOPFProblem) + constants = prob.cache.kkt_constants + isnothing(constants) && error( + "ACOPFProblem cache invariant violated: kkt_constants are missing. " * + "ACOPFProblem constructors and rebuilds should populate prob.cache.kkt_constants." + ) + return constants +end + """Extract per-bus aggregated load values for a given key ("pd" or "qd").""" function _extract_bus_load(prob::ACOPFProblem, key::String) - ref = prob.ref - n = prob.network.n - vals = zeros(n) - for i in 1:n - for lid in ref[:bus_loads][i] - vals[i] += ref[:load][lid][key] - end - end - return vals + constants = _require_kkt_constants(prob) + return key == "pd" ? copy(constants.pd) : copy(constants.qd) end """ Extract per-generator cost coefficient at a given index (1=quadratic, 2=linear). -AC OPF construction runs PowerModels cost standardization before this point, so -the analytical path assumes finite numeric coefficients. +MATPOWER parsing standardizes costs before this point, so the analytical path +assumes finite numeric coefficients. """ function _extract_gen_cost(prob::ACOPFProblem, cost_idx::Int) - k = prob.n_gen - vals = zeros(k) - for i in 1:k - vals[i] = prob.ref[:gen][i]["cost"][cost_idx] - end - return vals + constants = _require_kkt_constants(prob) + return cost_idx == 1 ? copy(constants.cq) : copy(constants.cl) end _extract_bus_pd(prob::ACOPFProblem) = _extract_bus_load(prob, "pd") @@ -1064,79 +1032,21 @@ _extract_gen_cq(prob::ACOPFProblem) = _extract_gen_cost(prob, 1) _extract_gen_cl(prob::ACOPFProblem) = _extract_gen_cost(prob, 2) """ -Extract per-branch flow limits (`rate_a`) from the problem's ref. +Extract per-branch flow limits (`rate_a`) from cached constants. -AC OPF construction runs PowerModels thermal-limit preprocessing before this -point, so the analytical path assumes finite numeric limits. +MATPOWER parsing normalizes thermal limits before this point, so the analytical +path assumes finite numeric limits. """ function _extract_branch_fmax(prob::ACOPFProblem) - m = prob.network.m - fmax = zeros(m) - for l in 1:m - fmax[l] = prob.ref[:branch][l]["rate_a"] - end - return fmax + return copy(_require_kkt_constants(prob).fmax) end """ -Pre-extract all constant data from the problem's ref for efficient analytical +Pre-extract all constant data from the problem for efficient analytical KKT assembly and repeated sensitivity evaluation. """ function _extract_kkt_constants(prob::ACOPFProblem) - ref = prob.ref - n, m, k = prob.network.n, prob.network.m, prob.n_gen - - # Branch electrical parameters - bp_g = Vector{Float64}(undef, m) - bp_b = Vector{Float64}(undef, m) - bp_tr = Vector{Float64}(undef, m) - bp_ti = Vector{Float64}(undef, m) - bp_g_fr = Vector{Float64}(undef, m) - bp_b_fr = Vector{Float64}(undef, m) - bp_g_to = Vector{Float64}(undef, m) - bp_b_to = Vector{Float64}(undef, m) - bp_tm = Vector{Float64}(undef, m) - bp_f_bus = Vector{Int}(undef, m) - bp_t_bus = Vector{Int}(undef, m) - bp_angmin = Vector{Float64}(undef, m) - bp_angmax = Vector{Float64}(undef, m) - for l in 1:m - branch = ref[:branch][l] - bp_g[l], bp_b[l] = PM.calc_branch_y(branch) - bp_tr[l], bp_ti[l] = PM.calc_branch_t(branch) - bp_g_fr[l] = branch["g_fr"] - bp_b_fr[l] = branch["b_fr"] - bp_g_to[l] = branch["g_to"] - bp_b_to[l] = branch["b_to"] - bp_tm[l] = branch["tap"]^2 - bp_f_bus[l] = branch["f_bus"] - bp_t_bus[l] = branch["t_bus"] - bp_angmin[l] = branch["angmin"] - bp_angmax[l] = branch["angmax"] - end - - return ( - # Branch parameters - g_br = bp_g, b_br = bp_b, tr = bp_tr, ti = bp_ti, - g_fr = bp_g_fr, b_fr = bp_b_fr, g_to = bp_g_to, b_to = bp_b_to, - tm = bp_tm, f_bus = bp_f_bus, t_bus = bp_t_bus, - angmin = bp_angmin, angmax = bp_angmax, - # Bus bounds - vmin = Float64[ref[:bus][i]["vmin"] for i in 1:n], - vmax = Float64[ref[:bus][i]["vmax"] for i in 1:n], - # Gen bounds and parameters - pmin = Float64[ref[:gen][i]["pmin"] for i in 1:k], - pmax = Float64[ref[:gen][i]["pmax"] for i in 1:k], - qmin = Float64[ref[:gen][i]["qmin"] for i in 1:k], - qmax = Float64[ref[:gen][i]["qmax"] for i in 1:k], - gen_bus = Int[ref[:gen][i]["gen_bus"] for i in 1:k], - cc = Float64[ref[:gen][i]["cost"][3] for i in 1:k], - # Shunt parameters (aggregated per bus) - gs = Float64[sum(ref[:shunt][s]["gs"] for s in ref[:bus_shunts][i]; init=0.0) for i in 1:n], - bs = Float64[sum(ref[:shunt][s]["bs"] for s in ref[:bus_shunts][i]; init=0.0) for i in 1:n], - # Reference bus - ref_bus_keys = sort(collect(keys(ref[:ref_buses]))), - ) + return prob.data.constants end # ============================================================================= @@ -1297,11 +1207,7 @@ function calc_kkt_jacobian_param(prob::ACOPFProblem, sol::ACOPFSolution, param:: p0 = _AC_PARAM_EXTRACT[param](prob) Jp = zeros(Float64, kkt_dims(prob), length(p0)) - constants = prob.cache.kkt_constants - if isnothing(constants) - constants = _extract_kkt_constants(prob) - prob.cache.kkt_constants = constants - end + constants = _require_kkt_constants(prob) if param === :sw return _fill_ac_param_jacobian_sw!(Jp, prob, sol, idx, constants) @@ -1458,11 +1364,7 @@ function _calc_ac_kkt_param_column(prob::ACOPFProblem, sol::ACOPFSolution, param p0 = _AC_PARAM_EXTRACT[param](prob) Kcol = zeros(Float64, kkt_dims(prob)) - constants = prob.cache.kkt_constants - if isnothing(constants) - constants = _extract_kkt_constants(prob) - prob.cache.kkt_constants = constants - end + constants = _require_kkt_constants(prob) if param === :sw prim = _branch_flow_gradients(vars.va, vars.vm, prob.network.sw, constants, col_idx) diff --git a/src/prob/kkt_dc_opf.jl b/src/prob/kkt_dc_opf.jl index aaf07df..53f97b2 100644 --- a/src/prob/kkt_dc_opf.jl +++ b/src/prob/kkt_dc_opf.jl @@ -351,12 +351,10 @@ Flatten solution primal and dual variables into a single vector for KKT evaluati # Variable ordering [va; pg; f; psh; lam_lb; lam_ub; gamma_lb; gamma_ub; rho_lb; rho_ub; mu_lb; mu_ub; nu_bal; nu_flow; eta] -where eta is the dual for the reference bus constraint (set to 0). +where eta is the dual for the reference bus constraint `va[ref_bus] == 0`, packed +from the recovered `sol.eta_ref` (the KKT operator includes this entry). """ function flatten_variables(sol::DCOPFSolution, prob::DCOPFProblem) - # Reference bus dual (typically not needed, set to 0) - η_ref = dual(prob.cons.ref) - return vcat( sol.va, sol.pg, @@ -372,7 +370,7 @@ function flatten_variables(sol::DCOPFSolution, prob::DCOPFProblem) sol.mu_ub, sol.nu_bal, sol.nu_flow, - [η_ref] + [sol.eta_ref] ) end diff --git a/src/sens/current.jl b/src/sens/current.jl index 6655092..a55f640 100644 --- a/src/sens/current.jl +++ b/src/sens/current.jl @@ -39,7 +39,7 @@ where I_ℓ is the current on branch ℓ connecting buses i and j. # Arguments - `v::Vector{ComplexF64}`: Voltage phasors at all buses - `Y::AbstractMatrix{ComplexF64}`: Bus admittance matrix -- `branch_data::Dict`: PowerModels branch dictionary +- `branch_data::Dict`: branch dictionary keyed by sequential branch index, each entry holding `index`/`f_bus`/`t_bus` # Keyword Arguments - `idx_slack::Int=1`: Index of the slack (reference) bus @@ -54,7 +54,7 @@ NamedTuple with fields: # Example ```julia -state = ACPowerFlowState(pm_net) +state = ACPowerFlowState(net, v) sens = calc_sensitivity(state, :im, :p) # How does current on line 2 change when active power at bus 3 increases? dI_dp = sens[2, 3] @@ -116,15 +116,10 @@ end """ calc_current_power_sensitivities(net::Dict; full=true) -Compute current-power sensitivities from a solved PowerModels network. - -Accepts both basic and non-basic networks. For non-basic networks, constructs -an ACPowerFlowState internally which handles ID translation. +Reject the removed dictionary wrapper with a migration hint. """ function calc_current_power_sensitivities(net::Dict; full::Bool=true) - state = ACPowerFlowState(net) - !isnothing(state.branch_data) || throw(ArgumentError("Failed to extract branch data from network")) - return calc_current_power_sensitivities(state; full=full) + throw(ArgumentError("dictionary wrappers were removed; construct ACPowerFlowState(ACNetwork(data), v)")) end """ @@ -133,11 +128,21 @@ end Compute current-power sensitivities from an ACPowerFlowState. This method provides a unified interface consistent with DC OPF sensitivities. -Requires that `state.branch_data` is not nothing. +Requires either `state.net` or `state.branch_data`. """ function calc_current_power_sensitivities(state::ACPowerFlowState; full::Bool=true) - !isnothing(state.branch_data) || throw(ArgumentError("ACPowerFlowState must have branch_data for current sensitivities")) - return calc_current_power_sensitivities(state.v, state.Y, state.branch_data; idx_slack=state.idx_slack, full=full) + ∂v_∂p, _, _ = calc_voltage_active_power_sensitivities( + state.v, state.Y; idx_slack=state.idx_slack, full=full) + ∂v_∂q, _, _ = calc_voltage_reactive_power_sensitivities( + state.v, state.Y; idx_slack=state.idx_slack, full=full) + ∂I_∂p = _branch_current_from_dv(∂v_∂p, state) + ∂I_∂q = _branch_current_from_dv(∂v_∂q, state) + return ( + dI_dp=∂I_∂p, + dI_dq=∂I_∂q, + dIm_dp=_current_magnitude_from_dI(∂I_∂p, state), + dIm_dq=_current_magnitude_from_dI(∂I_∂q, state), + ) end # ============================================================================= @@ -152,7 +157,7 @@ Compute sensitivity of branch active power flows w.r.t. power injections. Uses product rule: P_ℓ = Re(v_f · conj(I_ℓ)), so ∂P_ℓ/∂p_k = Re(∂v_f/∂p_k · conj(I_ℓ) + v_f · conj(∂I_ℓ/∂p_k)) -Requires `state.branch_data` to be set. +Requires either `state.net` or `state.branch_data`. # Returns NamedTuple with: @@ -160,41 +165,13 @@ NamedTuple with: - `df_dq`: ∂P_flow/∂q (m × n) """ function calc_branch_flow_power_sensitivities(state::ACPowerFlowState) - !isnothing(state.branch_data) || throw(ArgumentError("ACPowerFlowState must have branch_data for flow sensitivities")) - - v = state.v - Y = state.Y - n = state.n - m = state.m - - # Get complex voltage sensitivities (full=true for indexing) - ∂v_∂p, _, _ = calc_voltage_active_power_sensitivities(v, Y; idx_slack=state.idx_slack, full=true) - ∂v_∂q, _, _ = calc_voltage_reactive_power_sensitivities(v, Y; idx_slack=state.idx_slack, full=true) - - # Get complex current sensitivities - cur_sens = calc_current_power_sensitivities(state; full=true) - ∂I_∂p = cur_sens.dI_dp - ∂I_∂q = cur_sens.dI_dq - - df_dp = zeros(Float64, m, n) - df_dq = zeros(Float64, m, n) - - for (_, br) in state.branch_data - ℓ = br["index"] - f_bus = br["f_bus"] - t_bus = br["t_bus"] - - Y_ft = Y[f_bus, t_bus] - I_ℓ = Y_ft * (v[f_bus] - v[t_bus]) - v_f = v[f_bus] - - for k in 1:n - # Product rule: ∂P_ℓ/∂p_k = Re(∂v_f/∂p_k · conj(I_ℓ) + v_f · conj(∂I_ℓ/∂p_k)) - df_dp[ℓ, k] = real(∂v_∂p[f_bus, k] * conj(I_ℓ) + v_f * conj(∂I_∂p[ℓ, k])) - df_dq[ℓ, k] = real(∂v_∂q[f_bus, k] * conj(I_ℓ) + v_f * conj(∂I_∂q[ℓ, k])) - end - end - - return (df_dp=df_dp, df_dq=df_dq) + ∂v_∂p, _, _ = calc_voltage_active_power_sensitivities( + state.v, state.Y; idx_slack=state.idx_slack, full=true) + ∂v_∂q, _, _ = calc_voltage_reactive_power_sensitivities( + state.v, state.Y; idx_slack=state.idx_slack, full=true) + return ( + df_dp=_branch_flow_from_dv(∂v_∂p, state), + df_dq=_branch_flow_from_dv(∂v_∂q, state), + ) end diff --git a/src/sens/interface.jl b/src/sens/interface.jl index 49cc442..513ac61 100644 --- a/src/sens/interface.jl +++ b/src/sens/interface.jl @@ -430,19 +430,41 @@ function _voltage_phasor_single_dir(state::ACPowerFlowState, param::Symbol) return ∂v end +function _branch_records(state::ACPowerFlowState) + if !isnothing(state.net) + net = state.net + return [(l, net.f_bus[l], net.t_bus[l]) for l in 1:net.m] + end + isnothing(state.branch_data) && throw(ArgumentError( + "ACPowerFlowState must have an ACNetwork or branch_data for branch sensitivities")) + return sort!([(br["index"], br["f_bus"], br["t_bus"]) for br in values(state.branch_data)]) +end + +@inline function _state_branch_current_coefficients(state::ACPowerFlowState, l::Int, f_bus::Int, t_bus::Int) + if !isnothing(state.net) + return _branch_current_coefficients(state.net, l) + end + yft = state.Y[f_bus, t_bus] + return yft, -yft +end + +function _state_branch_currents(state::ACPowerFlowState) + currents = zeros(ComplexF64, state.m) + for (l, f_bus, t_bus) in _branch_records(state) + yff, yft = _state_branch_current_coefficients(state, l, f_bus, t_bus) + currents[l] = yff * state.v[f_bus] + yft * state.v[t_bus] + end + return currents +end + """Compute branch current phasor sensitivity from pre-computed voltage phasor sensitivity.""" function _branch_current_from_dv(∂v, state::ACPowerFlowState) - !isnothing(state.branch_data) || throw(ArgumentError( - "ACPowerFlowState must have branch_data for current sensitivities")) m = state.m ncols = size(∂v, 2) ∂I = zeros(ComplexF64, m, ncols) - for (_, br) in state.branch_data - ℓ = br["index"] - f_bus = br["f_bus"] - t_bus = br["t_bus"] - Y_ft = state.Y[f_bus, t_bus] - ∂I[ℓ, :] = Y_ft .* (∂v[f_bus, :] .- ∂v[t_bus, :]) + for (l, f_bus, t_bus) in _branch_records(state) + yff, yft = _state_branch_current_coefficients(state, l, f_bus, t_bus) + ∂I[l, :] = yff .* ∂v[f_bus, :] .+ yft .* ∂v[t_bus, :] end return ∂I end @@ -450,40 +472,37 @@ end """ Compute topology direct current term from explicit dependence on branch admittance. -Returns a diagonal m × m matrix: ∂I_ℓ/∂(param)_e = scale · ΔV_ℓ · δ_{ℓe}, where -scale is -1 for :g and -j for :b. Parallel branches (multiple branches between the -same bus pair) are not supported — consistent with the Y[f,t] indexing used elsewhere. +Returns a diagonal m × m matrix. Each diagonal entry is the explicit derivative +of the from-side pi-model branch current with respect to that branch's series +conductance or susceptance, including taps, phase shifts, and switching. """ function _topology_branch_current_direct(state::ACPowerFlowState, param::Symbol) net = state.net isnothing(net) && throw(ArgumentError( "ACPowerFlowState must have an ACNetwork (net field) for topology sensitivities")) - ΔV = net.A * state.v - scale = param === :g ? ComplexF64(-1.0, 0.0) : ComplexF64(0.0, -1.0) - return Diagonal(scale .* ΔV) + direct = zeros(ComplexF64, net.m) + for l in 1:net.m + tap = net.tap[l] * cis(net.shift[l]) + value = net.sw[l] * (state.v[net.f_bus[l]] / abs2(tap) - state.v[net.t_bus[l]] / conj(tap)) + direct[l] = param === :g ? value : im * value + end + return Diagonal(direct) end """Project complex branch current sensitivities to current magnitude sensitivities.""" function _current_magnitude_from_dI(∂I, state::ACPowerFlowState) - !isnothing(state.branch_data) || throw(ArgumentError( - "ACPowerFlowState must have branch_data for current sensitivities")) - v = state.v m = state.m ncols = size(∂I, 2) ∂Im = zeros(Float64, m, ncols) n_suppressed = 0 - for (_, br) in state.branch_data - ℓ = br["index"] - f_bus = br["f_bus"] - t_bus = br["t_bus"] - Y_ft = state.Y[f_bus, t_bus] - I_ℓ = Y_ft * (v[f_bus] - v[t_bus]) - if abs(I_ℓ) <= VOLTAGE_ZERO_TOL + currents = _state_branch_currents(state) + for l in 1:m + if abs(currents[l]) <= VOLTAGE_ZERO_TOL n_suppressed += 1 continue end - ∂Im[ℓ, :] = real.(∂I[ℓ, :] .* conj(I_ℓ)) ./ abs(I_ℓ) + ∂Im[l, :] = real.(∂I[l, :] .* conj(currents[l])) ./ abs(currents[l]) end if n_suppressed > 0 @debug "Current magnitude sensitivity: $n_suppressed branches had |I| < $VOLTAGE_ZERO_TOL; their ∂|I| rows are zero." @@ -504,20 +523,13 @@ end """Project voltage/current sensitivities to branch active power flow sensitivities.""" function _branch_flow_from_dv_dI(∂v, ∂I, state::ACPowerFlowState) - !isnothing(state.branch_data) || throw(ArgumentError( - "ACPowerFlowState must have branch_data for flow sensitivities")) v = state.v m = state.m ncols = size(∂I, 2) df = zeros(Float64, m, ncols) - for (_, br) in state.branch_data - ℓ = br["index"] - f_bus = br["f_bus"] - t_bus = br["t_bus"] - Y_ft = state.Y[f_bus, t_bus] - I_ℓ = Y_ft * (v[f_bus] - v[t_bus]) - v_f = v[f_bus] - df[ℓ, :] = real.(∂v[f_bus, :] .* conj(I_ℓ) .+ v_f .* conj.(∂I[ℓ, :])) + currents = _state_branch_currents(state) + for (l, f_bus, _) in _branch_records(state) + df[l, :] = real.(∂v[f_bus, :] .* conj(currents[l]) .+ v[f_bus] .* conj.(∂I[l, :])) end return df end diff --git a/src/sens/topology_ac.jl b/src/sens/topology_ac.jl index c3a55b1..44fe732 100644 --- a/src/sens/topology_ac.jl +++ b/src/sens/topology_ac.jl @@ -49,42 +49,37 @@ function _build_topology_rhs(state::ACPowerFlowState, param::Symbol) net = state.net isnothing(net) && throw(ArgumentError( "ACPowerFlowState must have an ACNetwork (net field) for topology " * - "sensitivities (:g, :b). Construct via ACPowerFlowState(pm_net::Dict) " * - "or ACPowerFlowState(net::ACNetwork, v).")) + "sensitivities (:g, :b). Construct via ACPowerFlowState(net::ACNetwork, v).")) n = state.n m = state.m v = state.v - A = net.A # m × n sparse incidence matrix - ns = _non_slack_indices(n, state.idx_slack) d = length(ns) - - # Edge voltage drops: ΔV = A * V (m-vector, complex) - ΔV = A * v - - # M = Diag(conj(V_r)) * A_r' * Diag(ΔV) - # where _r denotes reduced (slack rows removed) - # A_r' is n×m with slack row removed → d×m - # M is d × m complex - At_r = transpose(A)[ns, :] # d × m sparse transpose slice - - M = Diagonal(conj.(v[ns])) * At_r * Diagonal(ΔV) - - # Assemble RHS based on parameter - # Sign convention: the power flow equations are F = [P; Q] - [p_spec; q_spec] = 0 - # where P = Re(conj(V)·Y·V) and Q = -Im(conj(V)·Y·V). - # The Jacobian solve is: J_v · dv = -∂F/∂param = -∂[P;Q]/∂param. - # We return the RHS = ∂[P;Q]/∂param, and the caller negates in the solve. RHS = Matrix{Float64}(undef, 2d, m) - if param === :g - RHS[1:d, :] = real.(M) # ∂P/∂g - RHS[d+1:2d, :] = -imag.(M) # ∂Q/∂g - else # :b - RHS[1:d, :] = -imag.(M) # ∂P/∂b - RHS[d+1:2d, :] = -real.(M) # ∂Q/∂b + fill!(RHS, 0.0) + reduced_idx = Dict(bus => i for (i, bus) in enumerate(ns)) + for l in 1:m + fb, tb = net.f_bus[l], net.t_bus[l] + tap = net.tap[l] * cis(net.shift[l]) + scale = net.sw[l] * (param === :g ? 1.0 : im) + dYff = scale / abs2(tap) + dYft = -scale / conj(tap) + dYtf = -scale / tap + dYtt = scale + if haskey(reduced_idx, fb) + i = reduced_idx[fb] + value = conj(v[fb]) * (dYff * v[fb] + dYft * v[tb]) + RHS[i, l] = real(value) + RHS[d + i, l] = -imag(value) + end + if haskey(reduced_idx, tb) + i = reduced_idx[tb] + value = conj(v[tb]) * (dYtf * v[fb] + dYtt * v[tb]) + RHS[i, l] = real(value) + RHS[d + i, l] = -imag(value) + end end - return RHS end diff --git a/src/sens/vjp_jvp.jl b/src/sens/vjp_jvp.jl index aa04fa0..c643e93 100644 --- a/src/sens/vjp_jvp.jl +++ b/src/sens/vjp_jvp.jl @@ -300,11 +300,7 @@ end function _ac_kkt_context(prob::ACOPFProblem) sol = _ensure_ac_solved!(prob) idx = kkt_indices(prob) - constants = prob.cache.kkt_constants - if isnothing(constants) - constants = _extract_kkt_constants(prob) - prob.cache.kkt_constants = constants - end + constants = _require_kkt_constants(prob) fmax = _extract_branch_fmax(prob) return (; sol, idx, constants, fmax) end diff --git a/src/sens/voltage.jl b/src/sens/voltage.jl index 3cab109..435b28a 100644 --- a/src/sens/voltage.jl +++ b/src/sens/voltage.jl @@ -98,14 +98,10 @@ end """ calc_voltage_power_sensitivities(net::Dict; full=true) -Compute voltage-power sensitivities from a solved PowerModels network. - -Accepts both basic and non-basic networks. For non-basic networks, constructs -an ACPowerFlowState internally which handles ID translation. +Reject the removed dictionary wrapper with a migration hint. """ function calc_voltage_power_sensitivities(net::Dict; full::Bool=true) - state = ACPowerFlowState(net) - return calc_voltage_power_sensitivities(state; full=full) + throw(ArgumentError("dictionary wrappers were removed; construct ACPowerFlowState(ACNetwork(data), v)")) end """ diff --git a/src/types/ac_network.jl b/src/types/ac_network.jl index 0f8cde3..8904ba4 100644 --- a/src/types/ac_network.jl +++ b/src/types/ac_network.jl @@ -25,14 +25,9 @@ AC network data with vectorized admittance representation. Provides a unified interface for AC power flow and sensitivity analysis, -analogous to `DCNetwork` for DC formulations. Uses edge-based conductance -and susceptance vectors for differentiable admittance matrix construction. - -The admittance matrix is reconstructed as: - Y = A' * Diag(g + j*b) * A + Diag(g_shunt + j*b_shunt) - -For switching-aware formulation: - Y(sw) = A' * Diag((g + j*b) .* sw) * A + Diag(g_shunt + j*b_shunt) +analogous to `DCNetwork` for DC formulations. Each branch contributes its +from/from, from/to, to/from, and to/to pi-model coefficients, including line +charging, transformer taps, phase shifts, switching, and parallel lines. # Fields - `n`: Number of buses @@ -41,15 +36,14 @@ For switching-aware formulation: - `incidences`: Edge list [(i,j), ...] for each branch (sequential indices) - `g`: Branch conductances - `b`: Branch susceptances (note: typically negative for inductive lines) -- `g_shunt`: Shunt conductances per bus (from shunts + line charging) +- `g_shunt`: Shunt conductances per bus - `b_shunt`: Shunt susceptances per bus - `sw`: Branch switching states ∈ [0,1]^m - `is_switchable`: Which branches can be switched - `idx_slack`: Slack bus index (sequential) - `vm_min`, `vm_max`: Voltage magnitude limits per bus -- `i_max`: Branch current magnitude limits - `id_map`: Bidirectional mapping between original and sequential element IDs -- `ref`: PowerModels build_ref dictionary (nothing for Y-matrix constructors) +- typed branch, bus, and generator arrays used by PF/OPF constructors """ struct ACNetwork <: AbstractPowerNetwork # Dimensions @@ -76,13 +70,45 @@ struct ACNetwork <: AbstractPowerNetwork # Limits vm_min::Vector{Float64} vm_max::Vector{Float64} - i_max::Vector{Float64} # ID mapping id_map::IDMapping - # PowerModels reference (nothing for Y-matrix constructors) - ref::Union{Nothing,Dict} + # Branch parameters + f_bus::Vector{Int} + t_bus::Vector{Int} + br_r::Vector{Float64} + br_x::Vector{Float64} + br_b::Vector{Float64} + g_fr::Vector{Float64} + b_fr::Vector{Float64} + g_to::Vector{Float64} + b_to::Vector{Float64} + tap::Vector{Float64} + shift::Vector{Float64} + tm::Vector{Float64} + angmin::Vector{Float64} + angmax::Vector{Float64} + rate_a::Vector{Float64} + + # Bus injections and shunts + pd::Vector{Float64} + qd::Vector{Float64} + gs::Vector{Float64} + bs::Vector{Float64} + + # Generator data + pg::Vector{Float64} + qg::Vector{Float64} + gen_bus::Vector{Int} + pmin::Vector{Float64} + pmax::Vector{Float64} + qmin::Vector{Float64} + qmax::Vector{Float64} + cq::Vector{Float64} + cl::Vector{Float64} + cc::Vector{Float64} + ref_bus_keys::Vector{Int} end # ============================================================================= @@ -95,8 +121,8 @@ end AC power flow solution with full injection tracking. Provides a common interface for AC sensitivity computations, analogous to -`DCPowerFlowState` for DC power flow. Can be constructed from a PowerModels -network or from raw voltage/admittance data. +`DCPowerFlowState` for DC power flow. Can be constructed from an `ACNetwork` +and externally solved voltages, or from raw voltage/admittance data. # Fields - `net`: ACNetwork reference (optional, provides access to edge-level data) @@ -108,7 +134,7 @@ network or from raw voltage/admittance data. - `pd`: Real power demand per bus - `qg`: Reactive power generation per bus - `qd`: Reactive power demand per bus -- `branch_data`: Branch dictionary with sequential indices (required for :im sensitivity) +- `branch_data`: Optional branch dictionary for raw voltage/admittance states - `idx_slack`: Index of the slack (reference) bus - `n`: Number of buses - `m`: Number of branches @@ -116,7 +142,6 @@ network or from raw voltage/admittance data. # Constructors - `ACPowerFlowState(v, Y; ...)`: From voltage phasors and admittance matrix - `ACPowerFlowState(net::ACNetwork, v; ...)`: From ACNetwork and voltage solution -- `ACPowerFlowState(pm_net::Dict)`: From solved PowerModels network """ struct ACPowerFlowState <: AbstractPowerFlowState net::Union{ACNetwork, Nothing} @@ -176,111 +201,145 @@ end """ ACNetwork(net::Dict; idx_slack=nothing) -Construct ACNetwork from a PowerModels network dictionary. - -Accepts both basic and non-basic networks. Non-basic networks (with arbitrary -bus/branch/gen IDs) are automatically translated to sequential indices internally. -The original IDs are preserved in `id_map` for result interpretation. - -The `build_ref` result is stored on the network for reuse by downstream -constructors (e.g. `ACOPFProblem`, `ACPowerFlowState`), avoiding redundant -`deepcopy` + `build_ref` calls. - -# Arguments -- `net`: PowerModels network dictionary (basic or non-basic) -- `idx_slack`: Slack bus index (if not specified, uses reference bus from data) +Reject the removed dictionary API with a migration hint. """ function ACNetwork(net::Dict{String,<:Any}; idx_slack::Union{Nothing,Int}=nothing) - pm_data, ref, id_map = _prepare_network_data(net) + throw(ArgumentError("dictionary constructors were removed; parse a network file with PowerDiff.parse_file")) +end + +ACNetwork(net::PowerIO.Network; idx_slack::Union{Nothing,Int}=nothing) = + ACNetwork(_network_data(net); idx_slack=idx_slack) +# Build from PowerDiff network tables (see `_network_data`). The `PowerIO.Network` +# method runs PowerDiff's modeling deltas; this assumes the tables are already +# normalized, so programmatic callers can supply ready values directly. +function ACNetwork(data::NamedTuple; idx_slack::Union{Nothing,Int}=nothing) + id_map = IDMapping(data) n_bus = length(id_map.bus_ids) n_branch = length(id_map.branch_ids) + n_gen = length(id_map.gen_ids) + bus_tbl = Dict(bus.bus_i => bus for bus in data.bus) + branch_tbl = Dict(branch.index => branch for branch in data.branch) + gen_tbl = Dict(gen.index => gen for gen in data.gen) - # Build incidence matrix and edge list A = spzeros(n_branch, n_bus) incidences = Vector{Tuple{Int,Int}}(undef, n_branch) - - for (orig_id, br) in ref[:branch] - ix = id_map.branch_to_idx[orig_id] - f_idx = id_map.bus_to_idx[br["f_bus"]] - t_idx = id_map.bus_to_idx[br["t_bus"]] - A[ix, f_idx] = 1.0 - A[ix, t_idx] = -1.0 - incidences[ix] = (f_idx, t_idx) - end - - # Compute individual branch admittances from impedance + f_bus = Vector{Int}(undef, n_branch) + t_bus = Vector{Int}(undef, n_branch) + br_r = Vector{Float64}(undef, n_branch) + br_x = Vector{Float64}(undef, n_branch) + br_b = Vector{Float64}(undef, n_branch) + g_fr = zeros(n_branch) + b_fr = Vector{Float64}(undef, n_branch) + g_to = zeros(n_branch) + b_to = Vector{Float64}(undef, n_branch) + tap = Vector{Float64}(undef, n_branch) + shift = Vector{Float64}(undef, n_branch) + tm = Vector{Float64}(undef, n_branch) + angmin = Vector{Float64}(undef, n_branch) + angmax = Vector{Float64}(undef, n_branch) + rate_a = Vector{Float64}(undef, n_branch) g = zeros(n_branch) b = zeros(n_branch) - for (orig_id, br) in ref[:branch] - ix = id_map.branch_to_idx[orig_id] - r = br["br_r"] - x = br["br_x"] - - # Branch admittance: y = 1/(r + jx) = (r - jx)/(r² + x²) - z2 = r^2 + x^2 + 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 + 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[ix] = r / z2 - b[ix] = -x / z2 + g[l] = branch.br_r / z2 + b[l] = -branch.br_x / z2 else - _SILENCE_WARNINGS[] || @warn "Branch $(orig_id) has near-zero impedance (|z|² = $(z2)); treating as open (zero admittance)." + _SILENCE_WARNINGS[] || @warn "Branch $orig_id has near-zero impedance; treating as open." end end - # Shunt admittances: use PM.calc_admittance_matrix to get the full Y matrix, - # then extract shunts from diagonal minus branch contributions - am = PM.calc_admittance_matrix(pm_data) - g_shunt = zeros(n_bus) b_shunt = zeros(n_bus) - - for i in 1:n_bus - orig_bus_id = id_map.bus_ids[i] - # Find PM's internal index for this bus - pm_idx = am.bus_to_idx[orig_bus_id] - - # Y_ii = sum of all admittances connected to bus i - y_sum = am.matrix[pm_idx, pm_idx] - - # Subtract branch contributions - for (orig_br_id, br) in ref[:branch] - br_idx = id_map.branch_to_idx[orig_br_id] - if br["f_bus"] == orig_bus_id || br["t_bus"] == orig_bus_id - y_sum -= g[br_idx] + im * b[br_idx] - end - end - - g_shunt[i] = real(y_sum) - b_shunt[i] = imag(y_sum) + pd = zeros(n_bus) + qd = zeros(n_bus) + gs = zeros(n_bus) + bs = zeros(n_bus) + # to_powerdata aggregates loads/shunts into per-bus values (per-unit). + for bus in data.bus + i = id_map.bus_to_idx[bus.bus_i] + pd[i] += bus.pd + qd[i] += bus.qd + gs[i] += bus.gs + bs[i] += bus.bs + g_shunt[i] += bus.gs + b_shunt[i] += bus.bs end - # All branches active by default - sw = ones(n_branch) - is_switchable = trues(n_branch) - - # Find slack bus - if isnothing(idx_slack) - orig_ref = first(keys(ref[:ref_buses])) - idx_slack = id_map.bus_to_idx[orig_ref] + pg = zeros(n_bus) + qg = zeros(n_bus) + gen_bus = Vector{Int}(undef, n_gen) + pmin = Vector{Float64}(undef, n_gen) + pmax = Vector{Float64}(undef, n_gen) + qmin = Vector{Float64}(undef, n_gen) + qmax = Vector{Float64}(undef, n_gen) + cq = Vector{Float64}(undef, n_gen) + cl = Vector{Float64}(undef, n_gen) + cc = Vector{Float64}(undef, n_gen) + for orig_id in id_map.gen_ids + gen = gen_tbl[orig_id] + j = id_map.gen_to_idx[orig_id] + i = id_map.bus_to_idx[gen.gen_bus] + gen_bus[j] = i + pg[i] += gen.pg + qg[i] += gen.qg + pmin[j] = gen.pmin + pmax[j] = gen.pmax + qmin[j] = gen.qmin + qmax[j] = gen.qmax + cq[j], cl[j], cc[j] = gen.cost end - # Voltage limits (iterate in sequential order) - vm_min = [get(ref[:bus][id_map.bus_ids[i]], "vmin", 0.9) for i in 1:n_bus] - vm_max = [get(ref[:bus][id_map.bus_ids[i]], "vmax", 1.1) for i in 1:n_bus] - - # Current limits (from rate_a if available) - i_max = [get(ref[:branch][id_map.branch_ids[i]], "rate_a", Inf) for i in 1:n_branch] + if !isnothing(idx_slack) + 1 <= idx_slack <= n_bus || throw(ArgumentError( + "idx_slack=$idx_slack is not a valid bus index (1:$n_bus)")) + ref_bus_keys = [idx_slack] + else + ref_bus_keys = [id_map.bus_to_idx[id] for id in id_map.bus_ids if bus_tbl[id].bus_type == 3] + if isempty(ref_bus_keys) + _SILENCE_WARNINGS[] || @warn "No reference bus (type 3) in the network; defaulting to bus 1 as slack. Pass `idx_slack` to choose explicitly." + push!(ref_bus_keys, 1) + end + idx_slack = first(ref_bus_keys) + end + 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] return ACNetwork( - n_bus, n_branch, - sparse(A), incidences, - g, b, g_shunt, b_shunt, - sw, is_switchable, - idx_slack, - vm_min, vm_max, i_max, - id_map, - ref + 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, + id_map, f_bus, t_bus, br_r, br_x, br_b, g_fr, b_fr, g_to, b_to, + tap, shift, tm, angmin, angmax, rate_a, pd, qd, gs, bs, pg, qg, + gen_bus, pmin, pmax, qmin, qmax, cq, cl, cc, ref_bus_keys ) end @@ -290,7 +349,7 @@ end Construct ACNetwork from a complex admittance matrix. Extracts edge-based representation from the full admittance matrix. -Useful for direct construction without PowerModels. +Useful for direct construction from a raw admittance matrix. """ function ACNetwork(Y::AbstractMatrix{<:Complex}; idx_slack::Int=1) n = size(Y, 1) @@ -334,17 +393,20 @@ function ACNetwork(Y::AbstractMatrix{<:Complex}; idx_slack::Int=1) is_switchable = trues(m) vm_min = fill(0.9, n) vm_max = fill(1.1, n) - i_max = fill(Inf, m) - return ACNetwork( n, m, A, edges, g, b, g_shunt, b_shunt, sw, is_switchable, idx_slack, - vm_min, vm_max, i_max, - IDMapping(n, m, 0, 0), - nothing + vm_min, vm_max, + IDMapping(n, m, 0), + [edge[1] for edge in edges], [edge[2] for edge in edges], + zeros(m), zeros(m), zeros(m), zeros(m), zeros(m), zeros(m), zeros(m), + ones(m), zeros(m), ones(m), fill(-π, m), fill(π, m), fill(Inf, m), + zeros(n), zeros(n), zeros(n), zeros(n), + zeros(n), zeros(n), Int[], Float64[], Float64[], Float64[], + Float64[], Float64[], Float64[], Float64[], [idx_slack] ) end @@ -359,10 +421,7 @@ Reconstruct the bus admittance matrix Y from vectorized representation. Y = A' * Diag(g + j*b) * A + Diag(g_shunt + j*b_shunt) """ -function admittance_matrix(net::ACNetwork) - W = Diagonal(net.g .+ im .* net.b) - return transpose(net.A) * W * net.A + Diagonal(net.g_shunt .+ im .* net.b_shunt) -end +admittance_matrix(net::ACNetwork) = admittance_matrix(net, net.sw) """ admittance_matrix(net::ACNetwork, sw::AbstractVector) → SparseMatrixCSC{ComplexF64} @@ -372,8 +431,36 @@ Reconstruct admittance matrix with switching states. Y(sw) = A' * Diag((g + j*b) .* sw) * A + Diag(g_shunt + j*b_shunt) """ function admittance_matrix(net::ACNetwork, sw::AbstractVector) - W = Diagonal((net.g .+ im .* net.b) .* sw) - return transpose(net.A) * W * net.A + Diagonal(net.g_shunt .+ im .* net.b_shunt) + length(sw) == net.m || throw(DimensionMismatch("switching vector must have length $(net.m)")) + rows = collect(1:net.n) + cols = collect(1:net.n) + vals = ComplexF64.(net.g_shunt .+ im .* net.b_shunt) + sizehint!(rows, net.n + 4net.m) + sizehint!(cols, net.n + 4net.m) + sizehint!(vals, net.n + 4net.m) + for l in 1:net.m + yff, yft, ytf, ytt = _branch_admittance_coefficients(net, l) + fb, tb = net.f_bus[l], net.t_bus[l] + append!(rows, (fb, fb, tb, tb)) + append!(cols, (fb, tb, fb, tb)) + append!(vals, sw[l] .* (yff, yft, ytf, ytt)) + end + return sparse(rows, cols, vals, net.n, net.n) +end + +@inline function _branch_admittance_coefficients(net::ACNetwork, l::Int) + y = net.g[l] + im * net.b[l] + tap = net.tap[l] * cis(net.shift[l]) + yff = (y + net.g_fr[l] + im * net.b_fr[l]) / abs2(tap) + yft = -y / conj(tap) + ytf = -y / tap + ytt = y + net.g_to[l] + im * net.b_to[l] + return yff, yft, ytf, ytt +end + +@inline function _branch_current_coefficients(net::ACNetwork, l::Int) + yff, yft, _, _ = _branch_admittance_coefficients(net, l) + return net.sw[l] * yff, net.sw[l] * yft end # ============================================================================= @@ -435,22 +522,25 @@ q_polar(net::ACNetwork, vm::AbstractVector, δ::AbstractVector) = """ branch_current(net::ACNetwork, v::AbstractVector{<:Complex}) → Vector{ComplexF64} -Complex branch currents: I_branch = Diag(y) * A * v +Complex branch currents injected from the from-side bus. """ function branch_current(net::ACNetwork, v::AbstractVector{<:Complex}) - W = Diagonal(net.g .+ im .* net.b) - return W * net.A * v + return [ + let (yff, yft) = _branch_current_coefficients(net, l) + yff * v[net.f_bus[l]] + yft * v[net.t_bus[l]] + end + for l in 1:net.m + ] end """ branch_power(net::ACNetwork, v::AbstractVector{<:Complex}) → Vector{ComplexF64} -Complex branch power flows: S_branch = diag(A*v) * conj(I_branch) +Complex branch power flows injected from the from-side bus. """ function branch_power(net::ACNetwork, v::AbstractVector{<:Complex}) I = branch_current(net, v) - ΔV = net.A * v # Voltage difference across each branch - return ΔV .* conj.(I) + return v[net.f_bus] .* conj.(I) end # ============================================================================= @@ -467,7 +557,7 @@ Construct ACPowerFlowState from ACNetwork and voltage solution. - `v`: Complex voltage phasors from power flow solution # Keyword Arguments -- `pg`, `pd`, `qg`, `qd`: Generation and demand vectors (default to zeros) +- `pg`, `pd`, `qg`, `qd`: Generation and demand vectors (default to network values) """ function ACPowerFlowState( net::ACNetwork, @@ -483,11 +573,10 @@ function ACPowerFlowState( # Build admittance matrix from network Y = admittance_matrix(net) - # Default to zeros if not provided - pg_vec = isnothing(pg) ? zeros(n) : pg - pd_vec = isnothing(pd) ? zeros(n) : pd - qg_vec = isnothing(qg) ? zeros(n) : qg - qd_vec = isnothing(qd) ? zeros(n) : qd + pg_vec = isnothing(pg) ? copy(net.pg) : pg + pd_vec = isnothing(pd) ? copy(net.pd) : pd + qg_vec = isnothing(qg) ? copy(net.qg) : qg + qd_vec = isnothing(qd) ? copy(net.qd) : qd p_net = pg_vec - pd_vec q_net = qg_vec - qd_vec @@ -503,76 +592,10 @@ end """ ACPowerFlowState(pm_net::Dict) -Construct ACPowerFlowState from a solved PowerModels network. - -Accepts both basic and non-basic networks. Extracts voltage solution and -injection data. Creates an ACNetwork internally for access to edge-level data. -The network must have a solved power flow. +Reject the removed dictionary API with a migration hint. """ function ACPowerFlowState(pm_net::Dict) - # Create ACNetwork from the PowerModels data (handles basic/non-basic, stores ref) - net = ACNetwork(pm_net) - id_map = net.id_map - ref = net.ref # Reuse stored ref — no extra build_ref - - n = net.n - m = net.m - - # Extract bus voltages in sequential order - v = Vector{ComplexF64}(undef, n) - for i in 1:n - orig_id = id_map.bus_ids[i] - bus = ref[:bus][orig_id] - vm_val = get(bus, "vm", 1.0) - va_val = get(bus, "va", 0.0) - v[i] = vm_val * cis(va_val) - end - - Y = admittance_matrix(net) - - # Extract generation and demand in sequential bus order - pg = zeros(n) - qg = zeros(n) - for (orig_id, gen_ids) in ref[:bus_gens] - bus_idx = id_map.bus_to_idx[orig_id] - for gen_id in gen_ids - gen = ref[:gen][gen_id] - pg[bus_idx] += get(gen, "pg", 0.0) - qg[bus_idx] += get(gen, "qg", 0.0) - end - end - - pd = zeros(n) - qd = zeros(n) - for (orig_id, load_ids) in ref[:bus_loads] - bus_idx = id_map.bus_to_idx[orig_id] - for load_id in load_ids - load = ref[:load][load_id] - pd[bus_idx] += get(load, "pd", 0.0) - qd[bus_idx] += get(load, "qd", 0.0) - end - end - - p_net = pg - pd - q_net = qg - qd - - # Build branch_data with sequential indices for current sensitivity - seq_branch = Dict{String,Any}() - for (orig_id, br) in ref[:branch] - seq_idx = id_map.branch_to_idx[orig_id] - seq_br = copy(br) - seq_br["index"] = seq_idx - seq_br["f_bus"] = id_map.bus_to_idx[br["f_bus"]] - seq_br["t_bus"] = id_map.bus_to_idx[br["t_bus"]] - seq_branch[string(seq_idx)] = seq_br - end - - return ACPowerFlowState( - net, v, Y, - p_net, q_net, - pg, pd, qg, qd, - seq_branch, net.idx_slack, n, m - ) + throw(ArgumentError("dictionary constructors were removed; construct ACPowerFlowState(ACNetwork(data), v)")) end """ diff --git a/src/types/ac_opf_problem.jl b/src/types/ac_opf_problem.jl index b9cd9e6..fd6f9a3 100644 --- a/src/types/ac_opf_problem.jl +++ b/src/types/ac_opf_problem.jl @@ -167,7 +167,6 @@ Clear all cached AC sensitivity data. Called when problem parameters change. function invalidate!(cache::ACSensitivityCache) cache.solution = nothing cache.kkt_factor = nothing - cache.kkt_constants = nothing cache.dz_dsw = nothing cache.dz_dd = nothing cache.dz_dqd = nothing @@ -177,44 +176,134 @@ function invalidate!(cache::ACSensitivityCache) return nothing end -# ============================================================================= -# ACOPFProblem -# ============================================================================= +abstract type AbstractACOPFBackend end +struct JuMPBackend <: AbstractACOPFBackend end +struct ExaBackend <: AbstractACOPFBackend end + +struct ACOPFData{C} + arcs::Vector{NTuple{3,Int}} + arc_from_idx::Vector{Int} + arc_to_idx::Vector{Int} + bus_arc_idxs::Vector{Vector{Int}} + bus_gen_idxs::Vector{Vector{Int}} + ref_bus_keys::Vector{Int} + constants::C +end + +struct ACJuMPConstraints + ref_bus::Vector{ConstraintRef} + p_bal::Vector{ConstraintRef} + q_bal::Vector{ConstraintRef} + p_fr::Vector{ConstraintRef} + q_fr::Vector{ConstraintRef} + p_to::Vector{ConstraintRef} + q_to::Vector{ConstraintRef} + thermal_fr::Vector{ConstraintRef} + thermal_to::Vector{ConstraintRef} + angle_diff_lb::Vector{ConstraintRef} + angle_diff_ub::Vector{ConstraintRef} +end + +struct ACExaConstraints{RB,PB,QB,PF,QF,PT,QT,TF,TT,AL,AU} + ref_bus::RB + p_bal::PB + q_bal::QB + p_fr::PF + q_fr::QF + p_to::PT + q_to::QT + thermal_fr::TF + thermal_to::TT + angle_diff_lb::AL + angle_diff_ub::AU +end + +struct ACExaBusRecord + i::Int + pd::Float64 + qd::Float64 + gs::Float64 + bs::Float64 +end + +struct ACExaGenRecord + i::Int + bus::Int + cost1::Float64 + cost2::Float64 + cost3::Float64 +end + +struct ACExaArcRecord + i::Int + bus::Int +end + +struct ACExaBranchRecord + i::Int + f_idx::Int + t_idx::Int + f_bus::Int + t_bus::Int + c1::Float64 + c2::Float64 + c3::Float64 + c4::Float64 + c5::Float64 + c6::Float64 + c7::Float64 + c8::Float64 + sw::Float64 + angmin_scaled::Float64 + angmax_scaled::Float64 + rate_a_sq::Float64 +end + +const _EXA_HAS_FUNCTIONAL_API = isdefined(ExaModels, :add_var) + +_exa_core(; backend=nothing) = _EXA_HAS_FUNCTIONAL_API ? + ExaModels.ExaCore(; backend=backend, concrete=Val(true)) : + ExaModels.ExaCore(; backend=backend) + +_exa_add_var(core, dims...; kwargs...) = _EXA_HAS_FUNCTIONAL_API ? + ExaModels.add_var(core, dims...; kwargs...) : + (core, ExaModels.variable(core, dims...; kwargs...)) + +_exa_add_obj(core, expr; kwargs...) = _EXA_HAS_FUNCTIONAL_API ? + ExaModels.add_obj(core, expr; kwargs...) : + (core, ExaModels.objective(core, expr; kwargs...)) + +_exa_add_con(core, expr; kwargs...) = _EXA_HAS_FUNCTIONAL_API ? + ExaModels.add_con(core, expr; kwargs...) : + (core, ExaModels.constraint(core, expr; kwargs...)) + +_exa_add_con!(core, con, expr) = _EXA_HAS_FUNCTIONAL_API ? + ExaModels.add_con!(core, con, expr) : + (core, ExaModels.constraint!(core, con, expr)) + +_exa_model(core; kwargs...) = ExaModels.ExaModel(core; kwargs...) """ ACOPFProblem <: AbstractOPFProblem -Polar coordinate AC OPF wrapped around a JuMP model. - -# Fields -- `model`: JuMP Model -- `network`: ACNetwork data -- `va`, `vm`: Variable references for voltage angles and magnitudes -- `pg`, `qg`: Variable references for active and reactive generation -- `p`, `q`: Dict of branch flow variable references (keyed by arc tuple) -- `cons`: Named tuple of constraint references -- `ref`: PowerModels-style reference dictionary (sequential keys) -- `gen_buses`: Generator bus indices (maps generator index to bus index) -- `n_gen`: Number of generators -- `cache`: ACSensitivityCache for caching KKT derivatives -- `_optimizer`: Optimizer factory for model rebuilds (internal) -- `_silent`: Whether to suppress solver output (internal) +Polar coordinate AC OPF backed by either a JuMP model or an ExaModels model. """ -mutable struct ACOPFProblem <: AbstractOPFProblem - model::JuMP.Model +mutable struct ACOPFProblem{B<:AbstractACOPFBackend,M,D,VA,VM,PG,QG,P,Q,C,O} <: AbstractOPFProblem + model::M network::ACNetwork - va::Vector{VariableRef} - vm::Vector{VariableRef} - pg::Vector{VariableRef} - qg::Vector{VariableRef} - p::Dict{Tuple{Int,Int,Int}, VariableRef} - q::Dict{Tuple{Int,Int,Int}, VariableRef} - cons::NamedTuple - ref::Dict{Symbol, Any} + data::D + va::VA + vm::VM + pg::PG + qg::QG + p::P + q::Q + cons::C gen_buses::Vector{Int} n_gen::Int cache::ACSensitivityCache - _optimizer::Any + _backend::B + _optimizer::O _silent::Bool end @@ -225,306 +314,432 @@ end """ ACOPFProblem(network::ACNetwork; optimizer=Ipopt.Optimizer, silent=true) -Build a polar AC OPF problem from an ACNetwork that was constructed from a -PowerModels dictionary (i.e., `network.ref` must not be nothing). +Build a polar AC OPF problem from an ACNetwork. Accepts both basic and non-basic networks; internally remaps to sequential indices. # Arguments -- `network`: ACNetwork containing topology, admittances, and stored `ref` -- `optimizer`: JuMP-compatible optimizer (default: Ipopt) +- `network`: ACNetwork containing topology, admittances, and typed branch/gen data +- `backend`: `:jump` (default) or CPU `:exa` +- `optimizer`: JuMP-compatible optimizer for `backend=:jump` (default: Ipopt) - `silent`: Suppress solver output (default: true) # Example ```julia -pm_data = PowerModels.parse_file("case5.m") -prob = ACOPFProblem(pm_data) +data = parse_file("case5.m") +prob = ACOPFProblem(data) solve!(prob) ``` """ function ACOPFProblem( network::ACNetwork; + backend::Symbol=:jump, optimizer=Ipopt.Optimizer, silent::Bool=true ) - isnothing(network.ref) && error( - "ACOPFProblem requires an ACNetwork constructed from a PowerModels dict " * - "(network.ref must not be nothing). Use ACOPFProblem(pm_data::Dict) instead.") - - id_map = network.id_map - - # Remap ref to sequential keys so JuMP/KKT code can iterate 1:n, 1:m, 1:k - seq_ref = _remap_ref_to_sequential(network.ref, id_map) - - n_gen = length(seq_ref[:gen]) - gen_buses = [seq_ref[:gen][i]["gen_bus"] for i in 1:n_gen] + backend_tag = _ac_backend_tag(backend) + backend_tag isa ExaBackend && optimizer !== Ipopt.Optimizer && throw(ArgumentError( + "backend=:exa uses NLPModelsIpopt directly and does not accept a custom optimizer")) + _validate_acopf_network(network) + data = _build_acopf_data(network) + return _acopf_problem(network, data, backend_tag; optimizer=optimizer, silent=silent) +end - prob = ACOPFProblem( - JuMP.Model(), network, - VariableRef[], VariableRef[], VariableRef[], VariableRef[], - Dict{Tuple{Int,Int,Int}, VariableRef}(), Dict{Tuple{Int,Int,Int}, VariableRef}(), - (;), seq_ref, gen_buses, n_gen, ACSensitivityCache(), optimizer, silent - ) - _rebuild_jump_model!(prob) - return prob +function _ac_backend_tag(backend::Symbol) + backend == :jump && return JuMPBackend() + backend == :exa && return ExaBackend() + throw(ArgumentError("unsupported ACOPF backend :$backend (expected :jump or :exa)")) end -""" -Remap element dict entries to sequential indices, updating bus ID fields. -""" -function _remap_element_dict(ref_dict, to_idx::Dict{Int,Int}, bus_to_idx::Dict{Int,Int}, bus_fields::Vector{String}) - seq = Dict{Int, Any}() - for (orig_id, elem) in ref_dict - idx = to_idx[orig_id] - new_elem = deepcopy(elem) - new_elem["index"] = idx - for f in bus_fields - haskey(new_elem, f) && (new_elem[f] = bus_to_idx[new_elem[f]]) - end - seq[idx] = new_elem - end - return seq +function _validate_acopf_network(network::ACNetwork) + isempty(network.gen_bus) && throw(ArgumentError( + "ACOPFProblem requires generator, cost, demand, voltage limit, and finite branch limit data; " * + "ACNetwork values built from a raw admittance matrix are power flow networks only")) + all(isfinite, network.rate_a) || throw(ArgumentError( + "ACOPFProblem requires finite branch rate_a limits")) + all(>(0), network.rate_a) || throw(ArgumentError( + "ACOPFProblem requires positive branch rate_a limits")) + return nothing end -""" -Remap a `build_ref` result to sequential 1-based keys so that the JuMP model -and KKT code can iterate `for i in 1:n` without change. -""" -function _remap_ref_to_sequential(ref::Dict, id_map::IDMapping) - seq_bus = _remap_element_dict(ref[:bus], id_map.bus_to_idx, id_map.bus_to_idx, ["bus_i"]) - seq_branch = _remap_element_dict(ref[:branch], id_map.branch_to_idx, id_map.bus_to_idx, ["f_bus", "t_bus"]) - seq_gen = _remap_element_dict(ref[:gen], id_map.gen_to_idx, id_map.bus_to_idx, ["gen_bus"]) - seq_load = _remap_element_dict(ref[:load], id_map.load_to_idx, id_map.bus_to_idx, ["load_bus"]) - seq_shunt = _remap_element_dict(ref[:shunt], id_map.shunt_to_idx, id_map.bus_to_idx, ["shunt_bus"]) - - # Remap arcs: (l, i, j) → (seq_l, seq_i, seq_j) - function _remap_arc(arc) - (l, i, j) = arc - return (id_map.branch_to_idx[l], id_map.bus_to_idx[i], id_map.bus_to_idx[j]) +function _build_acopf_data(network::ACNetwork) + n, m = network.n, network.m + k = length(network.gen_bus) + g_br = copy(network.g) + b_br = copy(network.b) + tr = network.tap .* cos.(network.shift) + ti = network.tap .* sin.(network.shift) + g_fr = copy(network.g_fr) + b_fr = copy(network.b_fr) + g_to = copy(network.g_to) + b_to = copy(network.b_to) + tm = copy(network.tm) + f_bus = copy(network.f_bus) + t_bus = copy(network.t_bus) + angmin = copy(network.angmin) + angmax = copy(network.angmax) + fmax = copy(network.rate_a) + arcs = Vector{NTuple{3,Int}}(undef, 2m) + arc_from_idx = Vector{Int}(undef, m) + arc_to_idx = Vector{Int}(undef, m) + bus_arc_idxs = [Int[] for _ in 1:n] + + for l in 1:m + from_idx = 2l - 1 + to_idx = 2l + arc_from_idx[l] = from_idx + arc_to_idx[l] = to_idx + arcs[from_idx] = (l, f_bus[l], t_bus[l]) + arcs[to_idx] = (l, t_bus[l], f_bus[l]) + push!(bus_arc_idxs[f_bus[l]], from_idx) + push!(bus_arc_idxs[t_bus[l]], to_idx) end - seq_arcs = [_remap_arc(a) for a in ref[:arcs]] - seq_arcs_from = [_remap_arc(a) for a in ref[:arcs_from]] - seq_arcs_to = [_remap_arc(a) for a in ref[:arcs_to]] - - # Remap bus_arcs, bus_gens, bus_loads, bus_shunts to sequential bus keys - seq_bus_arcs = Dict{Int, Vector{Tuple{Int,Int,Int}}}() - for (orig_bus, arcs) in ref[:bus_arcs] - seq_bus_arcs[id_map.bus_to_idx[orig_bus]] = [_remap_arc(a) for a in arcs] + vmin = copy(network.vm_min) + vmax = copy(network.vm_max) + gs = copy(network.gs) + bs = copy(network.bs) + pd = copy(network.pd) + qd = copy(network.qd) + ref_bus_keys = copy(network.ref_bus_keys) + pmin = copy(network.pmin) + pmax = copy(network.pmax) + qmin = copy(network.qmin) + qmax = copy(network.qmax) + gen_bus = copy(network.gen_bus) + cq = copy(network.cq) + cl = copy(network.cl) + cc = copy(network.cc) + bus_gen_idxs = [Int[] for _ in 1:n] + + for i in 1:k + push!(bus_gen_idxs[gen_bus[i]], i) end - seq_bus_gens = Dict{Int, Vector{Int}}() - for (orig_bus, gen_ids) in ref[:bus_gens] - seq_bus_gens[id_map.bus_to_idx[orig_bus]] = [id_map.gen_to_idx[g] for g in gen_ids] - end + constants = ( + g_br = g_br, b_br = b_br, tr = tr, ti = ti, + g_fr = g_fr, b_fr = b_fr, g_to = g_to, b_to = b_to, + tm = tm, f_bus = f_bus, t_bus = t_bus, + angmin = angmin, angmax = angmax, + vmin = vmin, vmax = vmax, + pmin = pmin, pmax = pmax, qmin = qmin, qmax = qmax, + gen_bus = gen_bus, cq = cq, cl = cl, cc = cc, + fmax = fmax, gs = gs, bs = bs, pd = pd, qd = qd, + ref_bus_keys = ref_bus_keys, + ) - seq_bus_loads = Dict{Int, Vector{Int}}() - for (orig_bus, load_ids) in ref[:bus_loads] - seq_bus_loads[id_map.bus_to_idx[orig_bus]] = [id_map.load_to_idx[l] for l in load_ids] - end + return ACOPFData(arcs, arc_from_idx, arc_to_idx, bus_arc_idxs, bus_gen_idxs, ref_bus_keys, constants) +end - seq_bus_shunts = Dict{Int, Vector{Int}}() - for (orig_bus, shunt_ids) in ref[:bus_shunts] - seq_bus_shunts[id_map.bus_to_idx[orig_bus]] = [id_map.shunt_to_idx[s] for s in shunt_ids] - end +function _acopf_problem(network::ACNetwork, data::ACOPFData, ::JuMPBackend; optimizer, silent::Bool) + model, va, vm, pg, qg, p, q, cons = _build_jump_model(network, data, optimizer, silent) + cache = ACSensitivityCache() + cache.kkt_constants = data.constants + gen_buses = copy(data.constants.gen_bus) + return ACOPFProblem(model, network, data, va, vm, pg, qg, p, q, cons, + gen_buses, length(gen_buses), cache, JuMPBackend(), optimizer, silent) +end - # Remap ref_buses - seq_ref_buses = Dict{Int, Any}() - for (orig_bus, val) in ref[:ref_buses] - seq_ref_buses[id_map.bus_to_idx[orig_bus]] = val - end +function _acopf_problem(network::ACNetwork, data::ACOPFData, ::ExaBackend; optimizer, silent::Bool) + model, va, vm, pg, qg, p, q, cons = _build_examodel(network, data, optimizer, silent) + cache = ACSensitivityCache() + cache.kkt_constants = data.constants + gen_buses = copy(data.constants.gen_bus) + return ACOPFProblem(model, network, data, va, vm, pg, qg, p, q, cons, + gen_buses, length(gen_buses), cache, ExaBackend(), optimizer, silent) +end - return Dict{Symbol, Any}( - :bus => seq_bus, - :gen => seq_gen, - :branch => seq_branch, - :load => seq_load, - :shunt => seq_shunt, - :arcs => seq_arcs, - :arcs_from => seq_arcs_from, - :arcs_to => seq_arcs_to, - :bus_arcs => seq_bus_arcs, - :bus_gens => seq_bus_gens, - :bus_loads => seq_bus_loads, - :bus_shunts => seq_bus_shunts, - :ref_buses => seq_ref_buses, - ) +function _rebuild_model!(prob::ACOPFProblem{JuMPBackend}) + data = _build_acopf_data(prob.network) + model, va, vm, pg, qg, p, q, cons = _build_jump_model(prob.network, data, prob._optimizer, prob._silent) + prob.data = data + prob.model = model + prob.va = va + prob.vm = vm + prob.pg = pg + prob.qg = qg + prob.p = p + prob.q = q + prob.cons = cons + prob.gen_buses = copy(data.constants.gen_bus) + prob.n_gen = length(prob.gen_buses) + prob.cache.kkt_constants = data.constants + return nothing end -""" - _rebuild_jump_model!(prob::ACOPFProblem) +function _rebuild_model!(prob::ACOPFProblem{ExaBackend}) + data = _build_acopf_data(prob.network) + model, va, vm, pg, qg, p, q, cons = _build_examodel(prob.network, data, prob._optimizer, prob._silent) + prob.data = data + prob.model = model + prob.va = va + prob.vm = vm + prob.pg = pg + prob.qg = qg + prob.p = p + prob.q = q + prob.cons = cons + prob.gen_buses = copy(data.constants.gen_bus) + prob.n_gen = length(prob.gen_buses) + prob.cache.kkt_constants = data.constants + return nothing +end -Build (or rebuild) the JuMP model from current network parameters. -Called by the constructor and by `update_switching!` after mutating `network.sw`. -""" -function _rebuild_jump_model!(prob::ACOPFProblem) - network = prob.network - ref = prob.ref +function _build_jump_model(network::ACNetwork, data::ACOPFData, optimizer, silent::Bool) + constants = data.constants n, m = network.n, network.m - n_gen = prob.n_gen + n_gen = length(constants.gen_bus) + arc_fmax = [constants.fmax[arc[1]] for arc in data.arcs] # Create model - model = JuMP.Model(prob._optimizer) - prob._silent && set_silent(model) + model = JuMP.Model(optimizer) + silent && set_silent(model) set_optimizer_attribute(model, "tol", 1e-6) # Voltage variables - @variable(model, va[i in 1:n]) - @variable(model, ref[:bus][i]["vmin"] <= vm[i in 1:n] <= ref[:bus][i]["vmax"], start=1.0) + @variable(model, va[1:n]) + @variable(model, constants.vmin[i] <= vm[i in 1:n] <= constants.vmax[i], start=1.0) # Generation variables - @variable(model, ref[:gen][i]["pmin"] <= pg[i in 1:n_gen] <= ref[:gen][i]["pmax"]) - @variable(model, ref[:gen][i]["qmin"] <= qg[i in 1:n_gen] <= ref[:gen][i]["qmax"]) + @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]) # Branch flow variables - p = Dict{Tuple{Int,Int,Int}, VariableRef}() - q = Dict{Tuple{Int,Int,Int}, VariableRef}() - for (l, i, j) in ref[:arcs] - p[(l,i,j)] = @variable(model, base_name="p[$l,$i,$j]") - q[(l,i,j)] = @variable(model, base_name="q[$l,$i,$j]") - set_lower_bound(p[(l,i,j)], -ref[:branch][l]["rate_a"]) - set_upper_bound(p[(l,i,j)], ref[:branch][l]["rate_a"]) - set_lower_bound(q[(l,i,j)], -ref[:branch][l]["rate_a"]) - set_upper_bound(q[(l,i,j)], ref[:branch][l]["rate_a"]) - end + @variable(model, -arc_fmax[i] <= p[i in 1:length(data.arcs)] <= arc_fmax[i]) + @variable(model, -arc_fmax[i] <= q[i in 1:length(data.arcs)] <= arc_fmax[i]) # Objective: minimize generation cost (quadratic) @objective(model, Min, - sum(gen["cost"][1]*pg[i]^2 + gen["cost"][2]*pg[i] + gen["cost"][3] - for (i, gen) in ref[:gen]) + sum(constants.cq[i] * pg[i]^2 + constants.cl[i] * pg[i] + constants.cc[i] for i in 1:n_gen) ) # Reference bus constraint - ref_bus_con = @constraint(model, [i in keys(ref[:ref_buses])], va[i] == 0) + ref_bus = [@constraint(model, va[i] == 0) for i in data.ref_bus_keys] # Nodal power balance constraints - p_bal_cons = Vector{ConstraintRef}(undef, n) - q_bal_cons = Vector{ConstraintRef}(undef, n) - - for (i, bus) in ref[:bus] - bus_loads = [ref[:load][l] for l in ref[:bus_loads][i]] - bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][i]] - + p_bal = Vector{ConstraintRef}(undef, n) + q_bal = Vector{ConstraintRef}(undef, n) + for i in 1:n # Active power balance - p_bal_cons[i] = @constraint(model, - sum(p[a] for a in ref[:bus_arcs][i]) == - sum(pg[g] for g in ref[:bus_gens][i]) - - sum(load["pd"] for load in bus_loads) - - sum(shunt["gs"] for shunt in bus_shunts) * vm[i]^2 + p_bal[i] = @constraint(model, + sum(p[j] for j in data.bus_arc_idxs[i]) == + sum(pg[g] for g in data.bus_gen_idxs[i]) - constants.pd[i] - constants.gs[i] * vm[i]^2 ) # Reactive power balance - q_bal_cons[i] = @constraint(model, - sum(q[a] for a in ref[:bus_arcs][i]) == - sum(qg[g] for g in ref[:bus_gens][i]) - - sum(load["qd"] for load in bus_loads) + - sum(shunt["bs"] for shunt in bus_shunts) * vm[i]^2 + q_bal[i] = @constraint(model, + sum(q[j] for j in data.bus_arc_idxs[i]) == + sum(qg[g] for g in data.bus_gen_idxs[i]) - constants.qd[i] + constants.bs[i] * vm[i]^2 ) end + p_fr = Vector{ConstraintRef}(undef, m) + q_fr = Vector{ConstraintRef}(undef, m) + p_to = Vector{ConstraintRef}(undef, m) + q_to = Vector{ConstraintRef}(undef, m) + thermal_fr = Vector{ConstraintRef}(undef, m) + thermal_to = Vector{ConstraintRef}(undef, m) + angle_diff_lb = Vector{ConstraintRef}(undef, m) + angle_diff_ub = Vector{ConstraintRef}(undef, m) + # Branch power flow constraints and thermal limits - p_fr_cons = Dict{Int, ConstraintRef}() - q_fr_cons = Dict{Int, ConstraintRef}() - p_to_cons = Dict{Int, ConstraintRef}() - q_to_cons = Dict{Int, ConstraintRef}() - thermal_fr_cons = Dict{Int, ConstraintRef}() - thermal_to_cons = Dict{Int, ConstraintRef}() - angle_diff_cons = Dict{Int, Vector{ConstraintRef}}() - - for (l, branch) in ref[:branch] - f_idx = (l, branch["f_bus"], branch["t_bus"]) - t_idx = (l, branch["t_bus"], branch["f_bus"]) - - p_fr = p[f_idx] - q_fr = q[f_idx] - p_to = p[t_idx] - q_to = q[t_idx] - - vm_fr = vm[branch["f_bus"]] - vm_to = vm[branch["t_bus"]] - va_fr = va[branch["f_bus"]] - va_to = va[branch["t_bus"]] - - # Branch parameters (incorporating switching state sw) - g_br, b_br = PM.calc_branch_y(branch) - tr, ti = PM.calc_branch_t(branch) - g_fr_shunt = branch["g_fr"] - b_fr_shunt = branch["b_fr"] - g_to_shunt = branch["g_to"] - b_to_shunt = branch["b_to"] - tm = branch["tap"]^2 - - # Scale by switching state + for l in 1:m + f_idx = data.arc_from_idx[l] + t_idx = data.arc_to_idx[l] + fb = constants.f_bus[l] + tb = constants.t_bus[l] sw_l = network.sw[l] # AC Power Flow Constraints (from side) - p_fr_cons[l] = @constraint(model, - p_fr == sw_l * ((g_br + g_fr_shunt)/tm * vm_fr^2 + - (-g_br*tr + b_br*ti)/tm * (vm_fr * vm_to * cos(va_fr - va_to)) + - (-b_br*tr - g_br*ti)/tm * (vm_fr * vm_to * sin(va_fr - va_to))) + p_fr[l] = @constraint(model, + p[f_idx] == sw_l * ((constants.g_br[l] + constants.g_fr[l]) / constants.tm[l] * vm[fb]^2 + + (-constants.g_br[l] * constants.tr[l] + constants.b_br[l] * constants.ti[l]) / constants.tm[l] * + (vm[fb] * vm[tb] * cos(va[fb] - va[tb])) + + (-constants.b_br[l] * constants.tr[l] - constants.g_br[l] * constants.ti[l]) / constants.tm[l] * + (vm[fb] * vm[tb] * sin(va[fb] - va[tb]))) ) - q_fr_cons[l] = @constraint(model, - q_fr == sw_l * (-(b_br + b_fr_shunt)/tm * vm_fr^2 - - (-b_br*tr - g_br*ti)/tm * (vm_fr * vm_to * cos(va_fr - va_to)) + - (-g_br*tr + b_br*ti)/tm * (vm_fr * vm_to * sin(va_fr - va_to))) + q_fr[l] = @constraint(model, + q[f_idx] == sw_l * (-(constants.b_br[l] + constants.b_fr[l]) / constants.tm[l] * vm[fb]^2 - + (-constants.b_br[l] * constants.tr[l] - constants.g_br[l] * constants.ti[l]) / constants.tm[l] * + (vm[fb] * vm[tb] * cos(va[fb] - va[tb])) + + (-constants.g_br[l] * constants.tr[l] + constants.b_br[l] * constants.ti[l]) / constants.tm[l] * + (vm[fb] * vm[tb] * sin(va[fb] - va[tb]))) ) # AC Power Flow Constraints (to side) - p_to_cons[l] = @constraint(model, - p_to == sw_l * ((g_br + g_to_shunt) * vm_to^2 + - (-g_br*tr - b_br*ti)/tm * (vm_to * vm_fr * cos(va_to - va_fr)) + - (-b_br*tr + g_br*ti)/tm * (vm_to * vm_fr * sin(va_to - va_fr))) + p_to[l] = @constraint(model, + p[t_idx] == sw_l * ((constants.g_br[l] + constants.g_to[l]) * vm[tb]^2 + + (-constants.g_br[l] * constants.tr[l] - constants.b_br[l] * constants.ti[l]) / constants.tm[l] * + (vm[tb] * vm[fb] * cos(va[tb] - va[fb])) + + (-constants.b_br[l] * constants.tr[l] + constants.g_br[l] * constants.ti[l]) / constants.tm[l] * + (vm[tb] * vm[fb] * sin(va[tb] - va[fb]))) ) - q_to_cons[l] = @constraint(model, - q_to == sw_l * (-(b_br + b_to_shunt) * vm_to^2 - - (-b_br*tr + g_br*ti)/tm * (vm_to * vm_fr * cos(va_fr - va_to)) + - (-g_br*tr - b_br*ti)/tm * (vm_to * vm_fr * sin(va_to - va_fr))) + q_to[l] = @constraint(model, + q[t_idx] == sw_l * (-(constants.b_br[l] + constants.b_to[l]) * vm[tb]^2 - + (-constants.b_br[l] * constants.tr[l] + constants.g_br[l] * constants.ti[l]) / constants.tm[l] * + (vm[tb] * vm[fb] * cos(va[fb] - va[tb])) + + (-constants.g_br[l] * constants.tr[l] - constants.b_br[l] * constants.ti[l]) / constants.tm[l] * + (vm[tb] * vm[fb] * sin(va[tb] - va[fb]))) ) # Angle difference limits - angle_diff_cons[l] = [ - @constraint(model, sw_l * (va_fr - va_to) >= sw_l * branch["angmin"]), - @constraint(model, sw_l * (va_fr - va_to) <= sw_l * branch["angmax"]) - ] + angle_diff_lb[l] = @constraint(model, sw_l * (va[fb] - va[tb]) >= sw_l * constants.angmin[l]) + angle_diff_ub[l] = @constraint(model, sw_l * (va[fb] - va[tb]) <= sw_l * constants.angmax[l]) # Thermal limits (apparent power) - thermal_fr_cons[l] = @constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2) - thermal_to_cons[l] = @constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2) + thermal_fr[l] = @constraint(model, p[f_idx]^2 + q[f_idx]^2 <= constants.fmax[l]^2) + thermal_to[l] = @constraint(model, p[t_idx]^2 + q[t_idx]^2 <= constants.fmax[l]^2) end - prob.model = model - prob.va = collect(va) - prob.vm = collect(vm) - prob.pg = collect(pg) - prob.qg = collect(qg) - prob.p = p - prob.q = q - prob.cons = ( - ref_bus = ref_bus_con, - p_bal = p_bal_cons, - q_bal = q_bal_cons, - p_fr = p_fr_cons, - q_fr = q_fr_cons, - p_to = p_to_cons, - q_to = q_to_cons, - thermal_fr = thermal_fr_cons, - thermal_to = thermal_to_cons, - angle_diff = angle_diff_cons + cons = ACJuMPConstraints(ref_bus, p_bal, q_bal, p_fr, q_fr, p_to, q_to, + thermal_fr, thermal_to, angle_diff_lb, angle_diff_ub) + return model, collect(va), collect(vm), collect(pg), collect(qg), collect(p), collect(q), cons +end + +function _build_ac_exa_records(network::ACNetwork, data::ACOPFData) + constants = data.constants + n, m = network.n, network.m + k = length(constants.gen_bus) + + bus = Vector{ACExaBusRecord}(undef, n) + for i in 1:n + bus[i] = ACExaBusRecord(i, constants.pd[i], constants.qd[i], constants.gs[i], constants.bs[i]) + end + + gen = Vector{ACExaGenRecord}(undef, k) + for i in 1:k + gen[i] = ACExaGenRecord(i, constants.gen_bus[i], constants.cq[i], constants.cl[i], constants.cc[i]) + end + + arc = Vector{ACExaArcRecord}(undef, length(data.arcs)) + for i in eachindex(data.arcs) + arc[i] = ACExaArcRecord(i, data.arcs[i][2]) + end + + branch = Vector{ACExaBranchRecord}(undef, m) + for l in 1:m + sw = network.sw[l] + inv_tm = 1.0 / constants.tm[l] + c1 = sw * ((-constants.g_br[l] * constants.tr[l] - constants.b_br[l] * constants.ti[l]) * inv_tm) + c2 = sw * ((-constants.b_br[l] * constants.tr[l] + constants.g_br[l] * constants.ti[l]) * inv_tm) + c3 = sw * ((-constants.g_br[l] * constants.tr[l] + constants.b_br[l] * constants.ti[l]) * inv_tm) + c4 = sw * ((-constants.b_br[l] * constants.tr[l] - constants.g_br[l] * constants.ti[l]) * inv_tm) + c5 = sw * ((constants.g_br[l] + constants.g_fr[l]) * inv_tm) + c6 = sw * ((constants.b_br[l] + constants.b_fr[l]) * inv_tm) + c7 = sw * (constants.g_br[l] + constants.g_to[l]) + c8 = sw * (constants.b_br[l] + constants.b_to[l]) + branch[l] = ACExaBranchRecord( + l, + data.arc_from_idx[l], + data.arc_to_idx[l], + constants.f_bus[l], + constants.t_bus[l], + c1, + c2, + c3, + c4, + c5, + c6, + c7, + c8, + sw, + sw * constants.angmin[l], + sw * constants.angmax[l], + constants.fmax[l]^2, + ) + end + + return ( + bus = bus, + gen = gen, + arc = arc, + branch = branch, + ref_bus_keys = copy(data.ref_bus_keys), + nonpositive_lcon = fill(-Inf, m), + nonpositive_ucon = zeros(m), ) +end - return nothing +function _build_examodel(network::ACNetwork, data::ACOPFData, optimizer, silent::Bool) + constants = data.constants + n = network.n + n_gen = length(constants.gen_bus) + arc_fmax = [constants.fmax[arc[1]] for arc in data.arcs] + exa = _build_ac_exa_records(network, data) + + core = _exa_core() + core, va = _exa_add_var(core, n) + core, vm = _exa_add_var(core, n; start=ones(n), lvar=constants.vmin, uvar=constants.vmax) + core, pg = _exa_add_var(core, n_gen; lvar=constants.pmin, uvar=constants.pmax) + core, qg = _exa_add_var(core, n_gen; lvar=constants.qmin, uvar=constants.qmax) + core, p = _exa_add_var(core, length(data.arcs); lvar=-arc_fmax, uvar=arc_fmax) + core, q = _exa_add_var(core, length(data.arcs); lvar=-arc_fmax, uvar=arc_fmax) + + core, _ = _exa_add_obj(core, + g.cost1 * pg[g.i]^2 + g.cost2 * pg[g.i] + g.cost3 for g in exa.gen) + + core, ref_bus = _exa_add_con(core, va[i] for i in exa.ref_bus_keys) + core, p_fr = _exa_add_con(core, + p[b.f_idx] - b.c5 * vm[b.f_bus]^2 - + b.c3 * (vm[b.f_bus] * vm[b.t_bus] * cos(va[b.f_bus] - va[b.t_bus])) - + b.c4 * (vm[b.f_bus] * vm[b.t_bus] * sin(va[b.f_bus] - va[b.t_bus])) + for b in exa.branch) + core, q_fr = _exa_add_con(core, + q[b.f_idx] + b.c6 * vm[b.f_bus]^2 + + b.c4 * (vm[b.f_bus] * vm[b.t_bus] * cos(va[b.f_bus] - va[b.t_bus])) - + b.c3 * (vm[b.f_bus] * vm[b.t_bus] * sin(va[b.f_bus] - va[b.t_bus])) + for b in exa.branch) + core, p_to = _exa_add_con(core, + p[b.t_idx] - b.c7 * vm[b.t_bus]^2 - + b.c1 * (vm[b.t_bus] * vm[b.f_bus] * cos(va[b.t_bus] - va[b.f_bus])) - + b.c2 * (vm[b.t_bus] * vm[b.f_bus] * sin(va[b.t_bus] - va[b.f_bus])) + for b in exa.branch) + core, q_to = _exa_add_con(core, + q[b.t_idx] + b.c8 * vm[b.t_bus]^2 + + b.c2 * (vm[b.t_bus] * vm[b.f_bus] * cos(va[b.t_bus] - va[b.f_bus])) - + b.c1 * (vm[b.t_bus] * vm[b.f_bus] * sin(va[b.t_bus] - va[b.f_bus])) + for b in exa.branch) + core, angle_diff_lb = _exa_add_con(core, + b.angmin_scaled - b.sw * va[b.f_bus] + b.sw * va[b.t_bus] for b in exa.branch; + lcon=exa.nonpositive_lcon, ucon=exa.nonpositive_ucon) + core, angle_diff_ub = _exa_add_con(core, + b.sw * va[b.f_bus] - b.sw * va[b.t_bus] - b.angmax_scaled for b in exa.branch; + lcon=exa.nonpositive_lcon, ucon=exa.nonpositive_ucon) + core, thermal_fr = _exa_add_con(core, + p[b.f_idx]^2 + q[b.f_idx]^2 - b.rate_a_sq for b in exa.branch; + lcon=exa.nonpositive_lcon, ucon=exa.nonpositive_ucon) + core, thermal_to = _exa_add_con(core, + p[b.t_idx]^2 + q[b.t_idx]^2 - b.rate_a_sq for b in exa.branch; + lcon=exa.nonpositive_lcon, ucon=exa.nonpositive_ucon) + + core, p_bal = _exa_add_con(core, + b.pd + b.gs * vm[b.i]^2 for b in exa.bus) + core, _ = _exa_add_con!(core, p_bal, a.bus => p[a.i] for a in exa.arc) + core, _ = _exa_add_con!(core, p_bal, g.bus => -pg[g.i] for g in exa.gen) + + core, q_bal = _exa_add_con(core, + b.qd - b.bs * vm[b.i]^2 for b in exa.bus) + core, _ = _exa_add_con!(core, q_bal, a.bus => q[a.i] for a in exa.arc) + core, _ = _exa_add_con!(core, q_bal, g.bus => -qg[g.i] for g in exa.gen) + + model = _exa_model(core) + cons = ACExaConstraints(ref_bus, p_bal, q_bal, p_fr, q_fr, p_to, q_to, + thermal_fr, thermal_to, angle_diff_lb, angle_diff_ub) + return model, va, vm, pg, qg, p, q, cons end """ ACOPFProblem(pm_data::Dict; kwargs...) -Convenience constructor: build ACOPFProblem directly from PowerModels dict. +Reject the removed dictionary API with a migration hint. Accepts both basic and non-basic networks. """ function ACOPFProblem(pm_data::Dict; kwargs...) - network = ACNetwork(pm_data) - return ACOPFProblem(network; kwargs...) + throw(ArgumentError("dictionary constructors were removed; parse a network file with PowerDiff.parse_file")) end + +ACOPFProblem(net::PowerIO.Network; 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 b4ae34e..ed98782 100644 --- a/src/types/dc_network.jl +++ b/src/types/dc_network.jl @@ -39,7 +39,8 @@ topology sensitivity analysis. - `ref_bus`: Reference bus index (phase angle = 0) - `tau`: Regularization parameter for strong convexity - `id_map`: Bidirectional mapping between original and sequential element IDs -- `ref`: PowerModels build_ref dictionary (nothing for programmatic constructors) +- `demand`: Real power demand aggregated per bus +- `pg_init`: Initial real generation aggregated per bus """ struct DCNetwork <: AbstractPowerNetwork n::Int @@ -57,10 +58,11 @@ struct DCNetwork <: AbstractPowerNetwork cq::Vector{Float64} cl::Vector{Float64} c_shed::Vector{Float64} + demand::Vector{Float64} + pg_init::Vector{Float64} ref_bus::Int tau::Float64 id_map::IDMapping - ref::Union{Nothing,Dict} end # ============================================================================= @@ -83,6 +85,7 @@ Solution container for DC OPF problem, storing both primal and dual variables. - `rho_ub`, `rho_lb`: Generator upper/lower bound duals - `mu_lb`, `mu_ub`: Load shedding lower/upper bound duals - `gamma_lb`, `gamma_ub`: Phase angle difference lower/upper bound duals +- `eta_ref`: Reference bus constraint dual (`va[ref_bus] == 0`) - `objective`: Optimal objective value - `B_r_factor`: Cached factorization of reduced susceptance matrix `B[non_ref, non_ref]` """ @@ -101,6 +104,7 @@ struct DCOPFSolution{F<:Factorization{Float64}} <: AbstractOPFSolution mu_ub::Vector{Float64} gamma_lb::Vector{Float64} gamma_ub::Vector{Float64} + eta_ref::Float64 objective::Float64 B_r_factor::F end @@ -148,63 +152,336 @@ const DEFAULT_TAU = 1e-2 # constraints prevent delivery. 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.Network`, 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.Network + parse_file(io::IO; from="matpower", filetype=nothing) -> PowerIO.Network + +Parse a supported PowerIO network file into a `PowerIO.Network`. + +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...) + 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 + 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 + end +end + +""" + parse_matpower(io::IO) -> PowerIO.Network + parse_matpower(file::String; library=nothing) -> PowerIO.Network + +Parse MATPOWER v2 data into a `PowerIO.Network`. +""" +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_struct(file::String; library=nothing) + +Compatibility alias for [`parse_matpower`](@ref). +""" +parse_matpower_struct(file::String; library=nothing) = parse_matpower(file; library=library) + +_resolve_case_path(path::AbstractString, ::Nothing) = String(path) +_resolve_case_path(path::AbstractString, library) = joinpath(get_path(library), path) + +_powerio_format_hint(::Nothing, ::Nothing) = nothing +_powerio_format_hint(from, ::Nothing) = _format_token(from) +_powerio_format_hint(::Nothing, filetype) = _format_token(filetype) +function _powerio_format_hint(from, filetype) + f1 = _format_token(from) + f2 = _format_token(filetype) + f1 == f2 || throw(ArgumentError("conflicting parse format hints: from=$from and filetype=$filetype")) + return f1 +end + +function _format_token(x) + s = lowercase(String(x)) + startswith(s, ".") && (s = s[2:end]) + 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)")) +end + +""" + _network_data(net::PowerIO.Network) -> 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( + "PowerDiff does not support HVDC/dcline records; remove or convert dcline before parsing")) + isempty(PowerIO.storage(net)) || 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 + +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) + + 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 + + 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 + + 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")) + + return gen_rows, branch_rows +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) + angmin, angmax = _normalize_angle_bounds(Float64(br.angmin), Float64(br.angmax)) + rate_a = br.rate_a > 0 ? Float64(br.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) +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. + 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) +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. +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)) + zmag = hypot(r, x) + ymag = iszero(zmag) ? 0.0 : inv(zmag) + cmax = sqrt(fr_vmax^2 + to_vmax^2 - 2fr_vmax * to_vmax * cos(theta_max)) + return ymag * max(fr_vmax, to_vmax) * cmax +end + +# Default angle-difference bounds (radians in, radians out). MATPOWER angmin == angmax +# == 0 means unbounded; treat ±90° or wider and the zero case as a ±60° window, the +# MATPOWER/PowerModels convention. PowerIO's `to_powerdata` already converts to radians. +function _normalize_angle_bounds(angmin::Float64, angmax::Float64) + pad = deg2rad(60.0) + angmin <= -pi / 2 && (angmin = -pad) + angmax >= pi / 2 && (angmax = pad) + iszero(angmin) && iszero(angmax) && return (-pad, pad) + return angmin, angmax +end + # ============================================================================= # DCNetwork Constructors # ============================================================================= """ - DCNetwork(net::Dict; tau=DEFAULT_TAU, ref_bus=nothing) + DCNetwork(net::Dict; kwargs...) -Construct a DCNetwork from a PowerModels network dictionary. +Reject the removed dictionary API with a migration hint. +""" +function DCNetwork(net::Dict{String,<:Any}; kwargs...) + throw(ArgumentError("dictionary constructors were removed; parse a network file with PowerDiff.parse_file")) +end -Accepts both basic and non-basic networks. Non-basic networks (with arbitrary -bus/branch/gen IDs) are automatically translated to sequential indices internally. -The original IDs are preserved in `id_map` for result interpretation. +""" + DCNetwork(net::PowerIO.Network; tau=DEFAULT_TAU, ref_bus=nothing) -The `build_ref` result is stored on the network for reuse by downstream -constructors (e.g. `DCOPFProblem`, `DCPowerFlowState`), avoiding redundant -`deepcopy` + `build_ref` calls. +Construct a DCNetwork from a parsed PowerIO network. # Example ```julia -raw = PowerModels.parse_file("case14.m") -dc_net = DCNetwork(raw) # non-basic OK -# or -net = PowerModels.make_basic_network(raw) -dc_net = DCNetwork(net) # basic also OK +net = parse_file("case14.m") +dc_net = DCNetwork(net) ``` """ -function DCNetwork(net::Dict; tau::Float64=DEFAULT_TAU, ref_bus::Union{Nothing,Int}=nothing) - pm_data, ref, id_map = _prepare_network_data(net) +DCNetwork(net::PowerIO.Network; 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.Network` +# method runs PowerDiff's modeling deltas; this assumes the tables are already +# normalized, so programmatic callers can supply ready values directly. +function DCNetwork(data::NamedTuple; tau::Float64=DEFAULT_TAU, ref_bus::Union{Nothing,Int}=nothing) + id_map = IDMapping(data) n = length(id_map.bus_ids) m = length(id_map.branch_ids) k = length(id_map.gen_ids) + bus_tbl = Dict(bus.bus_i => bus for bus in data.bus) + 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 ref[:branch] using id_map translation + # Incidence matrix A (m × n) from active branches using id_map translation A = spzeros(m, n) - for (orig_id, br) in ref[:branch] + 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"]] + 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 # Generator-bus incidence matrix G_inc (n × k) G_inc = spzeros(n, k) - for (orig_id, gen) in ref[:gen] + for orig_id in id_map.gen_ids + gen = gen_tbl[orig_id] col = id_map.gen_to_idx[orig_id] - row = id_map.bus_to_idx[gen["gen_bus"]] + row = id_map.bus_to_idx[gen.gen_bus] G_inc[row, col] = 1.0 end # Branch susceptances: b = imag(1/z) b = zeros(m) - for (orig_id, br) in ref[:branch] + 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"] + r = br.br_r + x = br.br_x z2 = r^2 + x^2 if z2 > 1e-10 b[idx] = -x / z2 @@ -217,25 +494,36 @@ function DCNetwork(net::Dict; tau::Float64=DEFAULT_TAU, ref_bus::Union{Nothing,I sw = ones(m) # Limits (iterate in sequential order via sorted IDs) - fmax = [ref[:branch][id_map.branch_ids[i]]["rate_a"] for i in 1:m] - gmax = [ref[:gen][id_map.gen_ids[i]]["pmax"] for i in 1:k] - gmin = [ref[:gen][id_map.gen_ids[i]]["pmin"] for i in 1:k] + fmax = [branch_tbl[id_map.branch_ids[i]].rate_a for i in 1:m] + gmax = [gen_tbl[id_map.gen_ids[i]].pmax for i in 1:k] + gmin = [gen_tbl[id_map.gen_ids[i]].pmin for i in 1:k] # Phase angle difference limits - angmax = [ref[:branch][id_map.branch_ids[i]]["angmax"] for i in 1:m] - angmin = [ref[:branch][id_map.branch_ids[i]]["angmin"] for i in 1:m] + angmax = [branch_tbl[id_map.branch_ids[i]].angmax for i in 1:m] + angmin = [branch_tbl[id_map.branch_ids[i]].angmin for i in 1:m] # Cost coefficients (assumes polynomial cost with at least 2 terms) - cq = [ref[:gen][id_map.gen_ids[i]]["cost"][1] for i in 1:k] - cl = [ref[:gen][id_map.gen_ids[i]]["cost"][2] for i in 1:k] - - # Load-shedding cost: high penalty to discourage shedding when feasible - marginal_cost_ub = max(maximum(2cq .* gmax .+ cl), 1.0) + cq = [gen_tbl[id_map.gen_ids[i]].cost[1] for i in 1:k] + cl = [gen_tbl[id_map.gen_ids[i]].cost[2] for i in 1:k] + demand = calc_demand_vector(data, id_map) + pg_init = _calc_generation_vector(data, id_map) + + # Load-shedding cost: high penalty to discourage shedding when feasible. + # Guard the reduction so a generator-free network (valid for pure DC power flow + # built via the NamedTuple constructor) falls back to a unit marginal cost + # instead of `maximum` throwing on an empty collection. + marginal_cost_ub = k == 0 ? 1.0 : max(maximum(2cq .* gmax .+ cl), 1.0) c_shed = fill(DEFAULT_SHED_COST_MULTIPLIER * marginal_cost_ub, n) # Reference bus (translate original ID to sequential index) if isnothing(ref_bus) - orig_ref = first(keys(ref[:ref_buses])) + ref_candidates = [id for id in id_map.bus_ids if bus_tbl[id].bus_type == 3] + if isempty(ref_candidates) + _SILENCE_WARNINGS[] || @warn "No reference bus (type 3) in the network; defaulting to bus $(id_map.bus_ids[1]) as slack. Pass `ref_bus` to choose explicitly." + orig_ref = id_map.bus_ids[1] + else + orig_ref = ref_candidates[1] + end ref_bus = id_map.bus_to_idx[orig_ref] else # If user provided an original bus ID, translate it; validate the result @@ -248,7 +536,7 @@ function DCNetwork(net::Dict; tau::Float64=DEFAULT_TAU, ref_bus::Union{Nothing,I end return DCNetwork(n, m, k, A, G_inc, b, sw, fmax, gmax, gmin, angmax, angmin, - cq, cl, c_shed, ref_bus, tau, id_map, ref) + cq, cl, c_shed, demand, pg_init, ref_bus, tau, id_map) end """ @@ -269,10 +557,14 @@ function DCNetwork( cq::AbstractVector=zeros(k), cl::AbstractVector=zeros(k), c_shed::AbstractVector=fill(1e4, n), + demand::AbstractVector=zeros(n), + pg_init::AbstractVector=zeros(n), ref_bus::Int=1, tau::Float64=DEFAULT_TAU ) length(c_shed) == n || throw(DimensionMismatch("c_shed length $(length(c_shed)) must match number of buses $n")) + length(demand) == n || throw(DimensionMismatch("demand length $(length(demand)) must match number of buses $n")) + length(pg_init) == n || throw(DimensionMismatch("pg_init length $(length(pg_init)) must match number of buses $n")) all(c_shed .> 0) || throw(ArgumentError("c_shed must be strictly positive at all buses")) return DCNetwork( n, m, k, @@ -282,9 +574,9 @@ function DCNetwork( Float64.(angmax), Float64.(angmin), Float64.(cq), Float64.(cl), Float64.(c_shed), + Float64.(demand), Float64.(pg_init), ref_bus, tau, - IDMapping(n, m, k, 0), - nothing + IDMapping(n, m, k) ) end @@ -292,45 +584,24 @@ end # DCNetwork Helper Functions # ============================================================================= -""" - calc_demand_vector(net::Dict) - -Extract demand vector from PowerModels network dictionary. - -Works with both basic and non-basic networks. For non-basic networks, -uses `PM.build_ref` to resolve load-bus mappings and returns a vector -in sequential bus order (matching `IDMapping` from `DCNetwork(net)`). -""" -function calc_demand_vector(net::Dict) - _, ref, id_map = _prepare_network_data(net) - return _calc_demand_vector(ref, id_map) -end - """ calc_demand_vector(network::DCNetwork) -Extract demand vector from a DCNetwork that was constructed from a PowerModels dict. - -Uses the stored `ref` to avoid redundant `build_ref` calls. +Extract demand vector from a DCNetwork. """ function calc_demand_vector(network::DCNetwork) - isnothing(network.ref) && error( - "DCNetwork was not constructed from a PowerModels dict; cannot extract demand. " * - "Provide the demand vector explicitly.") - return _calc_demand_vector(network.ref, network.id_map) + return copy(network.demand) end -""" -Internal demand vector extraction from ref + id_map (avoids redundant build_ref). -""" -function _calc_demand_vector(ref::Dict, id_map::IDMapping) - n = length(id_map.bus_ids) - d = zeros(n) - for (bus_orig_id, load_ids) in ref[:bus_loads] - bus_idx = id_map.bus_to_idx[bus_orig_id] - for load_id in load_ids - d[bus_idx] += ref[:load][load_id]["pd"] - end +calc_demand_vector(net::PowerIO.Network) = 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) + # to_powerdata already aggregates loads into per-bus demand (per-unit). Index by + # the sorted IDMapping so demand aligns even when original bus IDs are unsorted. + d = zeros(length(id_map.bus_ids)) + for bus in data.bus + d[id_map.bus_to_idx[bus.bus_i]] += bus.pd end return d end @@ -377,18 +648,13 @@ function _factorize_B_r(net::DCNetwork) end """ -Aggregate generation to bus-level vector (uses ref + id_map). +Aggregate generation to bus-level vector. """ -function _calc_generation_vector(ref::Dict, id_map::IDMapping) +function _calc_generation_vector(data::NamedTuple, id_map::IDMapping) n = length(id_map.bus_ids) g = zeros(n) - for (bus_orig_id, gen_ids) in ref[:bus_gens] - bus_idx = id_map.bus_to_idx[bus_orig_id] - for gen_id in gen_ids - gen_data = ref[:gen][gen_id] - pg = get(gen_data, "pg", (gen_data["pmin"] + gen_data["pmax"]) / 2) - g[bus_idx] += pg - end + for gen in data.gen + g[id_map.bus_to_idx[gen.gen_bus]] += gen.pg end return g end @@ -474,25 +740,31 @@ function DCPowerFlowState(net::DCNetwork, d::AbstractVector{<:Real}) end """ - DCPowerFlowState(net::Dict; g=nothing, d=nothing) + DCPowerFlowState(net::Dict; kwargs...) -Construct DCPowerFlowState from PowerModels network dictionary. +Reject the removed dictionary API with a migration hint. +""" +function DCPowerFlowState(net::Dict{String,<:Any}; kwargs...) + throw(ArgumentError("dictionary constructors were removed; construct DCPowerFlowState(DCNetwork(data), g, d)")) +end + +""" + DCPowerFlowState(net::PowerIO.Network; g=nothing, d=nothing) -Accepts both basic and non-basic networks. +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(pm_net::Dict; g::Union{Nothing,AbstractVector}=nothing, d::Union{Nothing,AbstractVector}=nothing) - net = DCNetwork(pm_net) +function DCPowerFlowState(net::PowerIO.Network; g::Union{Nothing,AbstractVector}=nothing, d::Union{Nothing,AbstractVector}=nothing) + net = DCNetwork(net) - # Extract demand if not provided (reuses stored ref — no extra build_ref) if isnothing(d) - d = _calc_demand_vector(net.ref, net.id_map) + d = net.demand end # Aggregate generation to buses if not provided if isnothing(g) - g = _calc_generation_vector(net.ref, net.id_map) + g = net.pg_init end return DCPowerFlowState(net, g, d) diff --git a/src/types/dc_opf_problem.jl b/src/types/dc_opf_problem.jl index a03602c..e35ee87 100644 --- a/src/types/dc_opf_problem.jl +++ b/src/types/dc_opf_problem.jl @@ -125,7 +125,7 @@ B-θ formulation of DC OPF wrapped around a JuMP model. - `_optimizer`: Optimizer factory for model rebuilds (internal) - `_silent`: Whether to suppress solver output (internal) """ -mutable struct DCOPFProblem <: AbstractOPFProblem +mutable struct DCOPFProblem{O} <: AbstractOPFProblem model::JuMP.Model network::DCNetwork va::Vector{VariableRef} @@ -135,7 +135,7 @@ mutable struct DCOPFProblem <: AbstractOPFProblem d::Vector{Float64} cons::NamedTuple cache::DCSensitivityCache - _optimizer::Any + _optimizer::O _silent::Bool end @@ -205,14 +205,13 @@ function _rebuild_jump_model!(prob::DCOPFProblem) # Create model model = JuMP.Model(prob._optimizer) prob._silent && set_silent(model) - # Tighten Ipopt tolerances for accurate dual recovery (needed by sensitivity analysis) + # Tighten Ipopt tolerances for accurate dual recovery (needed by sensitivity analysis). if _is_ipopt_optimizer(prob._optimizer) set_optimizer_attribute(model, "tol", 1e-10) set_optimizer_attribute(model, "acceptable_tol", 1e-8) set_optimizer_attribute(model, "max_cpu_time", 30.0) end - # Variables @variable(model, va[1:n]) @variable(model, pg[1:k]) @variable(model, f[1:m]) @@ -279,13 +278,12 @@ end Build a B-θ DC OPF problem from a DCNetwork. -If `d` is not provided and the network was constructed from a PowerModels dict, -demand is extracted from the stored `ref`. +If `d` is not provided, demand is read from the network's typed cache. # Example ```julia -net = DCNetwork(pm_data) -prob = DCOPFProblem(net) # demand extracted from stored ref +net = DCNetwork(data) +prob = DCOPFProblem(net) # demand extracted from network data prob = DCOPFProblem(net; d=d) # explicit demand ``` """ @@ -297,17 +295,20 @@ function DCOPFProblem(network::DCNetwork; d::Union{Nothing,AbstractVector}=nothi end """ - DCOPFProblem(net::Dict; d=nothing, kwargs...) + DCOPFProblem(pm_data::Dict; kwargs...) -Convenience constructor: build DCOPFProblem directly from PowerModels dict. - -Accepts both basic and non-basic networks. -If `d` is not provided, extracts demand from the network data. +Reject the removed dictionary API with a migration hint. """ -function DCOPFProblem(net::Dict; d::Union{Nothing,AbstractVector}=nothing, tau::Float64=DEFAULT_TAU, kwargs...) - network = DCNetwork(net; tau=tau) +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.Network; 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) if isnothing(d) - d = _calc_demand_vector(network.ref, network.id_map) + d = calc_demand_vector(network) end return DCOPFProblem(network, d; kwargs...) end diff --git a/src/types/id_mapping.jl b/src/types/id_mapping.jl index f53f6e4..eb42784 100644 --- a/src/types/id_mapping.jl +++ b/src/types/id_mapping.jl @@ -2,152 +2,68 @@ # # 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. - -# ============================================================================= -# IDMapping: Bidirectional index mapping between original and sequential IDs -# ============================================================================= -# -# Translates between original PowerModels element IDs (which may be arbitrary -# integers like [1,2,3,4,10]) and sequential 1-based indices used internally. """ IDMapping -Bidirectional mapping between original PowerModels element IDs and sequential -1-based indices used for internal computation. - -Constructed from a `PM.build_ref` dictionary or as an identity mapping for -direct/programmatic constructors. - -# Fields -- `bus_ids::Vector{Int}`: Sorted original bus IDs; index position = sequential index -- `branch_ids::Vector{Int}`: Sorted original branch IDs -- `gen_ids::Vector{Int}`: Sorted original generator IDs -- `load_ids::Vector{Int}`: Sorted original load IDs -- `shunt_ids::Vector{Int}`: Sorted original shunt IDs -- `bus_to_idx::Dict{Int,Int}`: Original bus ID → sequential index -- `branch_to_idx::Dict{Int,Int}`: Original branch ID → sequential index -- `gen_to_idx::Dict{Int,Int}`: Original generator ID → sequential index -- `load_to_idx::Dict{Int,Int}`: Original load ID → sequential index -- `shunt_to_idx::Dict{Int,Int}`: Original shunt ID → sequential index +Bidirectional mapping between original network element IDs and sequential +1-based indices used for internal computation. Loads and shunts are aggregated +per bus, so only bus, branch, and generator IDs are tracked. """ struct IDMapping bus_ids::Vector{Int} branch_ids::Vector{Int} gen_ids::Vector{Int} - load_ids::Vector{Int} - shunt_ids::Vector{Int} bus_to_idx::Dict{Int,Int} branch_to_idx::Dict{Int,Int} gen_to_idx::Dict{Int,Int} - load_to_idx::Dict{Int,Int} - shunt_to_idx::Dict{Int,Int} - function IDMapping(bus_ids, branch_ids, gen_ids, load_ids, shunt_ids, - bus_to_idx, branch_to_idx, gen_to_idx, load_to_idx, shunt_to_idx) - issorted(bus_ids) || throw(ArgumentError("bus_ids must be sorted")) - issorted(branch_ids) || throw(ArgumentError("branch_ids must be sorted")) - issorted(gen_ids) || throw(ArgumentError("gen_ids must be sorted")) - issorted(load_ids) || throw(ArgumentError("load_ids must be sorted")) - issorted(shunt_ids) || throw(ArgumentError("shunt_ids must be sorted")) - length(bus_ids) == length(bus_to_idx) || throw(ArgumentError( - "bus_ids length ($(length(bus_ids))) must match bus_to_idx length ($(length(bus_to_idx)))")) - length(branch_ids) == length(branch_to_idx) || throw(ArgumentError( - "branch_ids length ($(length(branch_ids))) must match branch_to_idx length ($(length(branch_to_idx)))")) - length(gen_ids) == length(gen_to_idx) || throw(ArgumentError( - "gen_ids length ($(length(gen_ids))) must match gen_to_idx length ($(length(gen_to_idx)))")) - length(load_ids) == length(load_to_idx) || throw(ArgumentError( - "load_ids length ($(length(load_ids))) must match load_to_idx length ($(length(load_to_idx)))")) - length(shunt_ids) == length(shunt_to_idx) || throw(ArgumentError( - "shunt_ids length ($(length(shunt_ids))) must match shunt_to_idx length ($(length(shunt_to_idx)))")) - new(bus_ids, branch_ids, gen_ids, load_ids, shunt_ids, - bus_to_idx, branch_to_idx, gen_to_idx, load_to_idx, shunt_to_idx) + function IDMapping(bus_ids, branch_ids, gen_ids, + bus_to_idx, branch_to_idx, gen_to_idx) + for (ids, mapping, label) in ( + (bus_ids, bus_to_idx, "bus"), + (branch_ids, branch_to_idx, "branch"), + (gen_ids, gen_to_idx, "generator"), + ) + issorted(ids) || throw(ArgumentError("$label IDs must be sorted")) + length(ids) == length(mapping) || throw(ArgumentError( + "$label ID count must match mapping size")) + end + new(bus_ids, branch_ids, gen_ids, bus_to_idx, branch_to_idx, gen_to_idx) end end """ - IDMapping(ref::Dict) - -Construct IDMapping from a `PM.build_ref` reference dictionary. + IDMapping(data::NamedTuple) -Extracts sorted keys from `ref[:bus]`, `ref[:branch]`, `ref[:gen]`, `ref[:load]`, -and `ref[:shunt]` (if present) and builds inverse dictionaries. +Construct an ID mapping from PowerDiff network tables (see `_network_data`). """ -function IDMapping(ref::Dict) - for key in (:bus, :branch, :gen, :load) - haskey(ref, key) || throw(ArgumentError("ref missing required key :$key")) - end - isempty(keys(ref[:bus])) && throw(ArgumentError("Network has no buses")) - - bus_ids = sort(collect(keys(ref[:bus]))) - branch_ids = sort(collect(keys(ref[:branch]))) - gen_ids = sort(collect(keys(ref[:gen]))) - load_ids = sort(collect(keys(ref[:load]))) - shunt_ids = sort(collect(keys(get(ref, :shunt, Dict())))) - - bus_to_idx = Dict(id => i for (i, id) in enumerate(bus_ids)) - branch_to_idx = Dict(id => i for (i, id) in enumerate(branch_ids)) - gen_to_idx = Dict(id => i for (i, id) in enumerate(gen_ids)) - load_to_idx = Dict(id => i for (i, id) in enumerate(load_ids)) - shunt_to_idx = Dict(id => i for (i, id) in enumerate(shunt_ids)) - - return IDMapping(bus_ids, branch_ids, gen_ids, load_ids, shunt_ids, - bus_to_idx, branch_to_idx, gen_to_idx, load_to_idx, shunt_to_idx) +function IDMapping(data::NamedTuple) + isempty(data.bus) && throw(ArgumentError("Network has no buses")) + bus_ids = sort([b.bus_i for b in data.bus]) + branch_ids = sort([br.index for br in data.branch]) + gen_ids = sort([g.index for g in data.gen]) + return IDMapping( + bus_ids, branch_ids, gen_ids, + Dict(id => i for (i, id) in enumerate(bus_ids)), + Dict(id => i for (i, id) in enumerate(branch_ids)), + Dict(id => i for (i, id) in enumerate(gen_ids)), + ) end """ - IDMapping(n::Int, m::Int, k::Int, n_load::Int; n_shunt::Int=0) + IDMapping(n::Int, m::Int, k::Int) -Create an identity mapping (1:n → 1:n, etc.) for direct/programmatic constructors. +Create identity mappings for direct programmatic constructors. """ -function IDMapping(n::Int, m::Int, k::Int, n_load::Int; n_shunt::Int=0) - bus_ids = collect(1:n) - branch_ids = collect(1:m) - gen_ids = collect(1:k) - load_ids = collect(1:n_load) - shunt_ids = collect(1:n_shunt) - - bus_to_idx = Dict(i => i for i in 1:n) - branch_to_idx = Dict(i => i for i in 1:m) - gen_to_idx = Dict(i => i for i in 1:k) - load_to_idx = Dict(i => i for i in 1:n_load) - shunt_to_idx = Dict(i => i for i in 1:n_shunt) - - return IDMapping(bus_ids, branch_ids, gen_ids, load_ids, shunt_ids, - bus_to_idx, branch_to_idx, gen_to_idx, load_to_idx, shunt_to_idx) -end - -function Base.show(io::IO, m::IDMapping) - print(io, "IDMapping($(length(m.bus_ids)) buses, $(length(m.branch_ids)) branches, ", - "$(length(m.gen_ids)) gens, $(length(m.load_ids)) loads, ", - "$(length(m.shunt_ids)) shunts)") +function IDMapping(n::Int, m::Int, k::Int) + return IDMapping( + collect(1:n), collect(1:m), collect(1:k), + Dict(i => i for i in 1:n), Dict(i => i for i in 1:m), Dict(i => i for i in 1:k), + ) end -# ============================================================================= -# Network Data Preparation Helper -# ============================================================================= - -""" - _prepare_network_data(net::Dict) → (pm_data, ref, id_map) - -Preprocess a PowerModels network dictionary exactly once. -Returns a deepcopy with standardized costs/thermal limits, the build_ref result, -and the IDMapping. -""" -function _prepare_network_data(net::Dict) - pm_data = deepcopy(net) - PM.standardize_cost_terms!(pm_data, order=2) - PM.calc_thermal_limits!(pm_data) - ref = PM.build_ref(pm_data)[:it][:pm][:nw][0] - id_map = IDMapping(ref) - return (pm_data, ref, id_map) +function Base.show(io::IO, mapping::IDMapping) + print(io, "IDMapping($(length(mapping.bus_ids)) buses, ", + "$(length(mapping.branch_ids)) branches, $(length(mapping.gen_ids)) gens)") end diff --git a/src/types/show.jl b/src/types/show.jl index ea2fad1..ad1bd71 100644 --- a/src/types/show.jl +++ b/src/types/show.jl @@ -253,6 +253,9 @@ function _problem_status_str(model::JuMP.Model) end end +_problem_status_str(model::ExaModels.AbstractExaModel) = "ExaModel" +_problem_status_str(model) = "unknown" + function _dc_cache_list(cache::DCSensitivityCache) [string(f) for f in _DC_CACHE_FIELDS if !isnothing(getfield(cache, f))] end diff --git a/test/common.jl b/test/common.jl index cb83cbf..a1b6613 100644 --- a/test/common.jl +++ b/test/common.jl @@ -20,9 +20,8 @@ # which defines its own load_test_case inline). # # Data loaders: -# load_test_case — parse MATPOWER case from PowerModels' bundled test data, -# returns a basic (sequentially-renumbered) network dict -# load_raw_case — same, but returns the raw (non-basic) network dict +# load_test_case — parse a MATPOWER case into a PowerIO.Network (via PowerDiff.parse_file) +# load_pm_case — parse a PowerModels dictionary for oracle comparisons only # # Programmatic networks: # create_2bus_network — minimal 2-bus DC network with known closed-form solution @@ -48,20 +47,31 @@ using ForwardDiff using Ipopt using JuMP: MOI -# PowerModels test data directory +# PowerModels test data directory and PowerDiff-owned PGLib artifact handle const PM_DATA_DIR = joinpath(dirname(pathof(PowerModels)), "..", "test", "data", "matpower") +const PD_PGLIB_DIR = PowerDiff.get_path(:pglib) + +# Build PowerDiff network tables (the NamedTuple that DCNetwork/ACNetwork consume, +# see PowerDiff._network_data) directly, for programmatic test networks. Values are +# taken as-is — already normalized — like the removed hand-built ParsedCase path. +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, vmax=1.1) = + (; 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) +pd_case(bus, gen, branch; name="case", baseMVA=100.0) = (; name, baseMVA, bus, gen, branch) """ load_test_case(case_name::String) -Load and prepare a PowerModels test case. -Returns a basic network dictionary or nothing if not found. +Load a PowerModels test fixture through `PowerDiff.parse_file`, returning a `PowerIO.Network`. """ function load_test_case(case_name::String) case_path = joinpath(PM_DATA_DIR, case_name) if isfile(case_path) - raw = PowerModels.parse_file(case_path) - return PowerModels.make_basic_network(raw) + return PowerDiff.parse_file(case_path) else @warn "Test case not found: $case_path" return nothing @@ -69,19 +79,40 @@ function load_test_case(case_name::String) end """ - load_raw_case(case_name::String) + load_pm_case(case_name::String; basic=false) -Load a raw (non-basic) PowerModels network. +Load a PowerModels dictionary for oracle comparisons. """ -function load_raw_case(case_name::String) +function load_pm_case(case_name::String; basic::Bool=false) case_path = joinpath(PM_DATA_DIR, case_name) if isfile(case_path) - return PowerModels.parse_file(case_path) + data = PowerModels.parse_file(case_path) + return basic ? PowerModels.make_basic_network(data) : data else return nothing end end +load_raw_case(case_name::String) = load_pm_case(case_name) + +""" + load_ac_pf_state(case_name::String) + +Solve an AC power flow with PowerModels, then wrap the resulting voltage vector +in PowerDiff's typed AC network representation. +""" +function load_ac_pf_state(case_name::String) + data = load_test_case(case_name) + pm_data = load_pm_case(case_name; basic=true) + if isnothing(data) || isnothing(pm_data) + return nothing + end + + PowerModels.compute_ac_pf!(pm_data) + v = PowerModels.calc_basic_bus_voltage(pm_data) + return ACPowerFlowState(ACNetwork(data), v) +end + """ create_2bus_network(; fmax=100.0, gmax=10.0, cl=10.0, cq=0.0, tau=0.0) diff --git a/test/mwe_unified.jl b/test/mwe_unified.jl index 30b3ad5..6181c22 100644 --- a/test/mwe_unified.jl +++ b/test/mwe_unified.jl @@ -22,10 +22,9 @@ using PowerDiff using PowerModels -# Load a test network +# Load a test network through PowerDiff's PowerIO parser (returns a PowerIO.Network) case_path = joinpath(dirname(pathof(PowerModels)), "..", "test", "data", "matpower", "case14.m") -data = PowerModels.parse_file(case_path) -net_data = PowerModels.make_basic_network(data) +net_data = PowerDiff.parse_file(case_path) # ============================================================================= # DC Power Flow Example (non-OPF) @@ -122,12 +121,15 @@ println("\n" * "=" ^ 60) println("=== AC Power Flow: Symbol-Based Sensitivities ===") println("=" ^ 60) -# Solve AC power flow -PowerModels.compute_ac_pf!(net_data) +# Solve AC power flow with PowerModels as an external oracle for the voltage vector, +# then wrap it in PowerDiff's typed AC state (the supported API since the dict path was removed) +pm_basic = PowerModels.make_basic_network(PowerModels.parse_file(case_path)) +PowerModels.compute_ac_pf!(pm_basic) +v = PowerModels.calc_basic_bus_voltage(pm_basic) # Create ACNetwork and ACPowerFlowState ac_net = ACNetwork(net_data) -state = ACPowerFlowState(net_data) +state = ACPowerFlowState(ac_net, v) println("AC Network: n=$(ac_net.n) buses, m=$(ac_net.m) branches") println("Voltage magnitudes: |v| = ", round.(abs.(state.v), digits=4)) diff --git a/test/runtests.jl b/test/runtests.jl index 7141e47..4040409 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -43,10 +43,11 @@ include("common.jl") @test_skip false else dc_net = DCNetwork(net) + nd = PowerDiff._network_data(net) - @test dc_net.n == length(net["bus"]) - @test dc_net.m == length(net["branch"]) - @test dc_net.k == length(net["gen"]) + @test dc_net.n == length(nd.bus) + @test dc_net.m == length(nd.branch) + @test dc_net.k == length(nd.gen) @test size(dc_net.A) == (dc_net.m, dc_net.n) @test size(dc_net.G_inc) == (dc_net.n, dc_net.k) @test length(dc_net.b) == dc_net.m @@ -56,6 +57,23 @@ include("common.jl") end end +# Regression: calc_demand_vector(::NamedTuple) must index by the sorted IDMapping, +# not by file order, so demand lands on the right bus when bus IDs are unsorted. +@testset "calc_demand_vector aligns with sorted IDMapping" begin + # Per-bus demand (loads already aggregated into bus pd); only bus_i and pd vary. + buses = [pd_bus(10, 1; pd=1.0), pd_bus(2, 1; pd=0.0), pd_bus(5, 1; pd=3.0)] + data = pd_case(buses, NamedTuple[], NamedTuple[]; name="unsorted") + + d = calc_demand_vector(data) + id_map = PowerDiff.IDMapping(data) + + # Sorted bus order is [2, 5, 10]; demand attaches per the IDMapping, not enumerate order. + @test d == [0.0, 3.0, 1.0] + @test d[id_map.bus_to_idx[10]] == 1.0 + @test d[id_map.bus_to_idx[5]] == 3.0 + @test d[id_map.bus_to_idx[2]] == 0.0 +end + # Verify DCOPFProblem builds, solves, and satisfies basic feasibility. # Bound tolerances are ~1e-6: Ipopt converges to ~1e-8 complementarity, # but bound projection adds O(1e-6) slack. @@ -312,8 +330,9 @@ end # Solve with our implementation # Use small τ for numerical stability in KKT system - dc_net = DCNetwork(net; tau=1e-3) - prob = DCOPFProblem(dc_net, calc_demand_vector(net)) + typed = load_test_case("case5.m") + dc_net = DCNetwork(typed; tau=1e-3) + prob = DCOPFProblem(dc_net, calc_demand_vector(typed)) sol = solve!(prob) # For LMPs, use the power balance duals directly (ν_bal) @@ -735,9 +754,7 @@ end end # AC power flow on case14 - pf_data = deepcopy(net) - PowerModels.compute_ac_pf!(pf_data) - state = ACPowerFlowState(pf_data) + state = load_ac_pf_state("case14.m") dvm_dp = calc_sensitivity(state, :vm, :p) dvm_dq = calc_sensitivity(state, :vm, :q) @@ -852,6 +869,9 @@ include("test_kkt_vjp_jvp.jl") include("test_acpf_jacobian.jl") include("test_acpf_va_flow.jl") include("test_parameter_transforms.jl") +include("test_parser_parity.jl") +include("test_non_matpower_parsers.jl") +include("test_ac_opf_exa_backend.jl") include("test_ac_opf_all_sens.jl") include("test_ac_topology_sens.jl") include("test_angle_diff_duals.jl") diff --git a/test/smoke_rts_gmlc.jl b/test/smoke_rts_gmlc.jl index 4a8bc82..2eb56e5 100644 --- a/test/smoke_rts_gmlc.jl +++ b/test/smoke_rts_gmlc.jl @@ -26,11 +26,13 @@ # - RTS_GMLC.m lives in PowerSystems/PowerSystemCaseBuilder, not PowerModels # - File path depends on installed package version # - AC OPF + ForwardDiff KKT Jacobian is slow (~30s+) -# - DC line workaround required (PowerModels limitation) +# - RTS_GMLC carries dcline records and piecewise linear generator costs, +# which are converted in the PowerDiff smoke path using Test using LinearAlgebra using PowerDiff +using PowerIO using PowerModels PowerModels.silence() @@ -38,22 +40,28 @@ PowerModels.silence() # ───────────────────────────────────────────────────────────────────────────── # Locate RTS_GMLC.m # ───────────────────────────────────────────────────────────────────────────── -const RTS_PATHS = [ - joinpath(dirname(dirname(pathof(PowerModels))), "test", "data", "matpower", "RTS_GMLC.m"), - # PowerSystems / PowerSystemCaseBuilder fallbacks - filter(isfile, [ - joinpath(d, "data", "matpower", "RTS_GMLC.m") - for d in readdir(joinpath(homedir(), ".julia", "packages", "PowerSystemCaseBuilder"); join=true) - if isdir(d) - ])..., - filter(isfile, [ - joinpath(d, "data", "matpower", "RTS_GMLC.m") - for d in readdir(joinpath(homedir(), ".julia", "packages", "PowerSystems"); join=true) - if isdir(d) - ])..., -] - -const RTS_PATH = let found = filter(isfile, RTS_PATHS) +function _candidate_rts_paths() + paths = String[ + joinpath(dirname(dirname(pathof(PowerModels))), "test", "data", "matpower", "RTS_GMLC.m"), + ] + for pkg in ("PowerSystemCaseBuilder", "PowerSystems") + root = joinpath(homedir(), ".julia", "packages", pkg) + isdir(root) || continue + for d in readdir(root; join=true) + isdir(d) && push!(paths, joinpath(d, "data", "matpower", "RTS_GMLC.m")) + end + end + for depot in DEPOT_PATH + artifacts = joinpath(depot, "artifacts") + isdir(artifacts) || continue + for (root, _, files) in walkdir(artifacts) + "RTS_GMLC.m" in files && push!(paths, joinpath(root, "RTS_GMLC.m")) + end + end + return paths +end + +const RTS_PATH = let found = filter(isfile, _candidate_rts_paths()) if isempty(found) error(""" RTS_GMLC.m not found. Install one of: @@ -64,18 +72,69 @@ const RTS_PATH = let found = filter(isfile, RTS_PATHS) first(found) end +function _linearize_piecewise_cost!(gen) + get(gen, "model", 2) == 1 || return false + + points = Float64.(get(gen, "cost", Float64[])) + ncost = Int(get(gen, "ncost", length(points) ÷ 2)) + ncost >= 1 && length(points) >= 2 * ncost || throw(ArgumentError( + "piecewise linear generator cost has inconsistent point data")) + + xs = points[1:2:(2 * ncost)] + ys = points[2:2:(2 * ncost)] + if ncost == 1 + slope = 0.0 + intercept = ys[1] + else + dx = xs[end] - xs[1] + !iszero(dx) || throw(ArgumentError( + "piecewise linear generator cost has duplicate endpoint output")) + slope = (ys[end] - ys[1]) / dx + intercept = ys[1] - slope * xs[1] + end + + gen["model"] = 2 + gen["ncost"] = 3 + gen["cost"] = [0.0, slope, intercept] + return true +end + +function _for_powerdiff_smoke(data) + clean = deepcopy(data) + if haskey(clean, "dcline") && !isempty(clean["dcline"]) + println(" Emptying $(length(clean["dcline"])) DC line(s) before PowerIO conversion") + empty!(clean["dcline"]) + end + if haskey(clean, "gen") + n_pwl = count(_linearize_piecewise_cost!, values(clean["gen"])) + n_pwl > 0 && println(" Linearizing $n_pwl piecewise linear generator cost curve(s) for PowerDiff smoke path") + end + return clean +end + +function _voltage_vector(pm_data, bus_ids) + bus_by_id = Dict(Int(bus["bus_i"]) => bus for bus in values(pm_data["bus"])) + return ComplexF64[ + bus_by_id[id]["vm"] * cis(deg2rad(bus_by_id[id]["va"])) + for id in bus_ids + ] +end + println("Using RTS_GMLC.m from: $RTS_PATH") # ───────────────────────────────────────────────────────────────────────────── # Load network # ───────────────────────────────────────────────────────────────────────────── raw = PowerModels.parse_file(RTS_PATH) +raw_for_powerdiff = _for_powerdiff_smoke(raw) +powerio_net = PowerIO.from_powermodels(raw_for_powerdiff) # Verify this is truly non-basic (bus IDs are not 1:n) bus_ids = sort([raw["bus"][k]["bus_i"] for k in keys(raw["bus"])]) @assert bus_ids != collect(1:length(bus_ids)) "RTS_GMLC should have non-sequential bus IDs" println("Bus IDs: $(bus_ids[1])..$(bus_ids[end]) ($(length(bus_ids)) buses)") println("Branches: $(length(raw["branch"])), Generators: $(length(raw["gen"])), Loads: $(length(raw["load"]))") +!isempty(get(raw, "dcline", Dict())) && println("DCLines excluded from PowerDiff smoke path: $(length(raw["dcline"]))") @testset "RTS_GMLC Smoke Tests" begin @@ -83,9 +142,9 @@ println("Branches: $(length(raw["branch"])), Generators: $(length(raw["gen"])), # DCNetwork # ================================================================= @testset "DCNetwork construction" begin - dc_net = DCNetwork(raw) + dc_net = DCNetwork(powerio_net) @test dc_net.n == length(bus_ids) - @test dc_net.m == length(raw["branch"]) + @test dc_net.m == length(raw_for_powerdiff["branch"]) # build_ref() filters inactive generators (158 total, 96 active) @test dc_net.k == length(dc_net.id_map.gen_ids) @test dc_net.k > 0 @@ -102,8 +161,8 @@ println("Branches: $(length(raw["branch"])), Generators: $(length(raw["gen"])), # DC Power Flow # ================================================================= @testset "DC power flow" begin - dc_net = DCNetwork(raw) - d = calc_demand_vector(raw) + dc_net = DCNetwork(powerio_net) + d = calc_demand_vector(powerio_net) @test length(d) == dc_net.n @test sum(d) > 0 # nonzero demand @@ -127,7 +186,7 @@ println("Branches: $(length(raw["branch"])), Generators: $(length(raw["gen"])), # ================================================================= local dc_prob # share across DC testsets @testset "DC OPF solve" begin - dc_prob = DCOPFProblem(raw) + dc_prob = DCOPFProblem(powerio_net) sol = solve!(dc_prob) @test sol.objective > 0 @test all(isfinite, sol.pg) @@ -165,18 +224,12 @@ println("Branches: $(length(raw["branch"])), Generators: $(length(raw["gen"])), end # ================================================================= - # ACNetwork (requires emptying dcline — PowerModels limitation) + # ACNetwork # ================================================================= - raw_ac = deepcopy(raw) - if !isempty(raw_ac["dcline"]) - println(" Emptying $(length(raw_ac["dcline"])) DC line(s) for AC compatibility") - empty!(raw_ac["dcline"]) - end - @testset "ACNetwork construction" begin - ac_net = ACNetwork(raw_ac) + ac_net = ACNetwork(powerio_net) @test ac_net.n == length(bus_ids) - @test ac_net.m == length(raw_ac["branch"]) + @test ac_net.m == length(raw_for_powerdiff["branch"]) @test ac_net.id_map.bus_ids == bus_ids end @@ -184,9 +237,10 @@ println("Branches: $(length(raw["branch"])), Generators: $(length(raw["gen"])), # AC Power Flow # ================================================================= @testset "AC power flow" begin - pf_data = deepcopy(raw_ac) + pf_data = deepcopy(raw_for_powerdiff) PowerModels.compute_ac_pf!(pf_data) - state = ACPowerFlowState(pf_data) + ac_net = ACNetwork(powerio_net) + state = ACPowerFlowState(ac_net, _voltage_vector(pf_data, ac_net.id_map.bus_ids)) @test all(isfinite, abs.(state.v)) @test all(v -> 0.8 < abs(v) < 1.2, state.v) # reasonable voltage magnitudes @@ -207,7 +261,7 @@ println("Branches: $(length(raw["branch"])), Generators: $(length(raw["gen"])), # AC OPF # ================================================================= @testset "AC OPF" begin - ac_prob = ACOPFProblem(raw_ac) + ac_prob = ACOPFProblem(powerio_net) sol = solve!(ac_prob) @test sol.objective > 0 @test all(isfinite, sol.vm) @@ -234,9 +288,9 @@ println("Branches: $(length(raw["branch"])), Generators: $(length(raw["gen"])), # Cross-validate: DC basic vs non-basic (SVD comparison) # ================================================================= @testset "DC basic vs non-basic SVD match" begin - basic = PowerModels.make_basic_network(deepcopy(raw)) - prob_nb = DCOPFProblem(raw) - prob_b = DCOPFProblem(basic) + basic = PowerModels.make_basic_network(deepcopy(raw_for_powerdiff)) + prob_nb = DCOPFProblem(powerio_net) + prob_b = DCOPFProblem(PowerIO.from_powermodels(basic)) sol_nb = solve!(prob_nb) sol_b = solve!(prob_b) diff --git a/test/test_ac_opf_all_sens.jl b/test/test_ac_opf_all_sens.jl index ff54aba..c434356 100644 --- a/test/test_ac_opf_all_sens.jl +++ b/test/test_ac_opf_all_sens.jl @@ -38,8 +38,7 @@ end @testset "AC OPF All Sensitivities" begin pm_path = joinpath(dirname(pathof(PowerModels)), "..", "test", "data", "matpower") file = joinpath(pm_path, "case5.m") - pm_data = PowerModels.parse_file(file) - pm_data = PowerModels.make_basic_network(pm_data) + pm_data = PowerDiff.parse_file(file) operands = [:va, :vm, :pg, :qg, :lmp, :qlmp] n, m, k = 5, 7, 5 # case5.m dimensions: 5 buses, 7 branches, 5 generators @@ -97,24 +96,15 @@ end ε = 1e-5 # Find buses with significant load (sorted for determinism) load_buses = Int[] - for (lid, load) in pm_data["load"] - bus = load["load_bus"] - if load["pd"] > 0.1 && !(bus in load_buses) - push!(load_buses, bus) - end + for bus in findall(>(0.1), prob.network.pd) + push!(load_buses, bus) end sort!(load_buses) for bus_idx in load_buses[1:min(2, length(load_buses))] - pm_pert = deepcopy(pm_data) - # Perturb all loads at this bus - for (lid, load) in pm_pert["load"] - if load["load_bus"] == bus_idx - load["pd"] += ε - end - end - - prob_pert = ACOPFProblem(pm_pert; silent=true) + net_pert = deepcopy(prob.network) + net_pert.pd[bus_idx] += ε + prob_pert = ACOPFProblem(net_pert; silent=true) sol_pert = solve!(prob_pert) for op in [:va, :vm, :pg, :qg] @@ -143,22 +133,14 @@ end ε = 1e-5 # Find buses with reactive load load_buses = Int[] - for (lid, load) in pm_data["load"] - bus = load["load_bus"] - if abs(load["qd"]) > 0.001 && !(bus in load_buses) - push!(load_buses, bus) - end + for bus in findall(>(0.001), abs.(prob.network.qd)) + push!(load_buses, bus) end for bus_idx in load_buses[1:min(2, length(load_buses))] - pm_pert = deepcopy(pm_data) - for (lid, load) in pm_pert["load"] - if load["load_bus"] == bus_idx - load["qd"] += ε - end - end - - prob_pert = ACOPFProblem(pm_pert; silent=true) + net_pert = deepcopy(prob.network) + net_pert.qd[bus_idx] += ε + prob_pert = ACOPFProblem(net_pert; silent=true) sol_pert = solve!(prob_pert) for op in [:va, :vm, :pg, :qg] @@ -186,10 +168,9 @@ end ε = 1e-5 for gen_idx in 1:min(2, k) - pm_pert = deepcopy(pm_data) - pm_pert["gen"][string(gen_idx)]["cost"][1] += ε - - prob_pert = ACOPFProblem(pm_pert; silent=true) + net_pert = deepcopy(prob.network) + net_pert.cq[gen_idx] += ε + prob_pert = ACOPFProblem(net_pert; silent=true) sol_pert = solve!(prob_pert) for op in [:va, :vm, :pg, :qg] @@ -217,10 +198,9 @@ end ε = 1e-5 for gen_idx in 1:min(2, k) - pm_pert = deepcopy(pm_data) - pm_pert["gen"][string(gen_idx)]["cost"][2] += ε - - prob_pert = ACOPFProblem(pm_pert; silent=true) + net_pert = deepcopy(prob.network) + net_pert.cl[gen_idx] += ε + prob_pert = ACOPFProblem(net_pert; silent=true) sol_pert = solve!(prob_pert) for op in [:va, :vm, :pg, :qg] @@ -248,10 +228,9 @@ end ε = 1e-5 for branch_idx in 1:min(2, m) - pm_pert = deepcopy(pm_data) - pm_pert["branch"][string(branch_idx)]["rate_a"] += ε - - prob_pert = ACOPFProblem(pm_pert; silent=true) + net_pert = deepcopy(prob.network) + net_pert.rate_a[branch_idx] += ε + prob_pert = ACOPFProblem(net_pert; silent=true) sol_pert = solve!(prob_pert) for op in [:va, :vm, :pg, :qg] @@ -300,7 +279,7 @@ end end # ========================================================================= - # Invalidation clears all cached data + # Invalidation clears derived cached data but preserves static KKT constants # ========================================================================= @testset "Invalidation clears cache" begin prob = ACOPFProblem(pm_data; silent=true) @@ -317,6 +296,6 @@ end @test isnothing(prob.cache.dz_dcq) @test isnothing(prob.cache.dz_dcl) @test isnothing(prob.cache.dz_dfmax) - @test isnothing(prob.cache.kkt_constants) + @test !isnothing(prob.cache.kkt_constants) end end diff --git a/test/test_ac_opf_exa_backend.jl b/test/test_ac_opf_exa_backend.jl new file mode 100644 index 0000000..6b7ad20 --- /dev/null +++ b/test/test_ac_opf_exa_backend.jl @@ -0,0 +1,113 @@ +using PowerDiff +using Ipopt +using JuMP: optimizer_with_attributes +using LinearAlgebra +using Test + +import PowerDiff: kkt, flatten_variables + +@testset "AC OPF Exa Backend" begin + pm_data = load_test_case("case5.m") + + if isnothing(pm_data) + @test_skip false + else + function _assert_small_kkt(prob, sol; tol=2e-3) + K = kkt(flatten_variables(sol, prob), prob) + @test norm(K, Inf) < tol + end + + function _assert_solution_parity(prob_exa, sol_exa, prob_jump, sol_jump; + primal_atol=1e-5, dual_atol=1e-3, price_atol=1e-5, obj_atol=1e-4, rtol=1e-5) + for field in (:va, :vm, :pg, :qg) + @test getfield(sol_exa, field) ≈ getfield(sol_jump, field) atol=primal_atol rtol=rtol + end + + for field in ( + :nu_p_bal, :nu_q_bal, :nu_ref_bus, + :nu_p_fr, :nu_p_to, :nu_q_fr, :nu_q_to, + :lam_angle_lb, :lam_angle_ub, + :mu_vm_lb, :mu_vm_ub, + :rho_pg_lb, :rho_pg_ub, :rho_qg_lb, :rho_qg_ub, + ) + @test getfield(sol_exa, field) ≈ getfield(sol_jump, field) atol=dual_atol rtol=rtol + end + + @test sol_exa.objective ≈ sol_jump.objective atol=obj_atol rtol=1e-8 + + for arc in keys(sol_jump.p) + @test sol_exa.p[arc] ≈ sol_jump.p[arc] atol=primal_atol rtol=rtol + @test sol_exa.q[arc] ≈ sol_jump.q[arc] atol=primal_atol rtol=rtol + end + + @test calc_lmp(sol_exa, prob_exa) ≈ calc_lmp(sol_jump, prob_jump) atol=price_atol rtol=rtol + @test calc_qlmp(sol_exa, prob_exa) ≈ calc_qlmp(sol_jump, prob_jump) atol=price_atol rtol=rtol + end + + @testset "Constructs and solves" begin + prob = ACOPFProblem(pm_data; backend=:exa, silent=true) + sol = solve!(prob) + + @test sol.objective > 0 + @test all(isfinite, sol.va) + @test all(isfinite, sol.vm) + @test all(isfinite, sol.pg) + @test all(isfinite, sol.qg) + _assert_small_kkt(prob, sol) + end + + @testset "Rejects custom optimizer" begin + optimizer = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0) + @test_throws ArgumentError ACOPFProblem(pm_data; backend=:exa, optimizer) + end + + @testset "Solve parity with JuMP backend" begin + prob_jump = ACOPFProblem(pm_data; backend=:jump, silent=true) + prob_exa = ACOPFProblem(pm_data; backend=:exa, silent=true) + + sol_jump = solve!(prob_jump) + sol_exa = solve!(prob_exa) + + _assert_solution_parity(prob_exa, sol_exa, prob_jump, sol_jump) + _assert_small_kkt(prob_exa, sol_exa) + end + + @testset "Switching update parity" begin + prob_jump = ACOPFProblem(pm_data; backend=:jump, silent=true) + prob_exa = ACOPFProblem(pm_data; backend=:exa, silent=true) + + sw = ones(prob_jump.network.m) + sw[1] = 0.85 + sw[min(2, end)] = 0.7 + + update_switching!(prob_jump, sw) + update_switching!(prob_exa, sw) + + sol_jump = solve!(prob_jump) + sol_exa = solve!(prob_exa) + + _assert_solution_parity(prob_exa, sol_exa, prob_jump, sol_jump; + primal_atol=2e-5, dual_atol=1e-3, obj_atol=1e-4, rtol=2e-5) + _assert_small_kkt(prob_exa, sol_exa; tol=3e-3) + end + + @testset "Sensitivity parity" begin + prob_jump = ACOPFProblem(pm_data; backend=:jump, silent=true) + prob_exa = ACOPFProblem(pm_data; backend=:exa, silent=true) + + for (op, param, atol, rtol) in [ + (:pg, :d, 1e-4, 1e-4), + (:va, :sw, 1e-4, 1e-4), + (:lmp, :d, 5e-4, 5e-4), + ] + S_jump = calc_sensitivity(prob_jump, op, param) + S_exa = calc_sensitivity(prob_exa, op, param) + + @test size(S_exa) == size(S_jump) + @test S_exa.row_to_id == S_jump.row_to_id + @test S_exa.col_to_id == S_jump.col_to_id + @test Matrix(S_exa) ≈ Matrix(S_jump) atol=atol rtol=rtol + end + end + end +end diff --git a/test/test_ac_opf_sens.jl b/test/test_ac_opf_sens.jl index 917ef7c..6e65a9c 100644 --- a/test/test_ac_opf_sens.jl +++ b/test/test_ac_opf_sens.jl @@ -28,8 +28,7 @@ using Test pm_path = joinpath(dirname(pathof(PowerModels)), "..", "test", "data", "matpower") file = joinpath(pm_path, "case5.m") - pm_data = PowerModels.parse_file(file) - pm_data = PowerModels.make_basic_network(pm_data) + pm_data = PowerDiff.parse_file(file) # Create and solve AC OPF @testset "ACOPFProblem construction and solving" begin @@ -48,8 +47,8 @@ using Test @test length(sol.qg) == 5 # Voltage magnitudes should be within limits - @test all(sol.vm .>= 0.9) - @test all(sol.vm .<= 1.1) + @test all(sol.vm .>= prob.network.vm_min .- 1e-6) + @test all(sol.vm .<= prob.network.vm_max .+ 1e-6) end @@ -285,9 +284,9 @@ using Test end @testset "Single ∂K/∂sw column includes angle-limit terms" begin - pm_tight = deepcopy(pm_data) - pm_tight["branch"]["1"]["angmax"] = 0.05 - prob = ACOPFProblem(pm_tight; silent=true) + net_tight = ACNetwork(pm_data) + net_tight.angmax[1] = 0.05 + prob = ACOPFProblem(net_tight; silent=true) sol = solve!(prob) @test abs(sol.lam_angle_lb[1]) + abs(sol.lam_angle_ub[1]) > 1e-3 diff --git a/test/test_ac_topology_sens.jl b/test/test_ac_topology_sens.jl index 73012c3..1e3c938 100644 --- a/test/test_ac_topology_sens.jl +++ b/test/test_ac_topology_sens.jl @@ -30,30 +30,8 @@ import PowerDiff: admittance_matrix # PQ Newton solver: pf_residual_pq / solve_pf_pq from common.jl @isdefined(pf_residual_pq) || include("common.jl") -function _branch_flows(v, Y, branch_data) - pf = zeros(length(branch_data)) - for (_, br) in branch_data - ℓ = br["index"] - f_bus = br["f_bus"] - t_bus = br["t_bus"] - Y_ft = Y[f_bus, t_bus] - I_ℓ = Y_ft * (v[f_bus] - v[t_bus]) - pf[ℓ] = real(v[f_bus] * conj(I_ℓ)) - end - return pf -end - -function _branch_currents_mag(v, Y, branch_data) - im_vec = zeros(length(branch_data)) - for (_, br) in branch_data - ℓ = br["index"] - f_bus = br["f_bus"] - t_bus = br["t_bus"] - Y_ft = Y[f_bus, t_bus] - im_vec[ℓ] = abs(Y_ft * (v[f_bus] - v[t_bus])) - end - return im_vec -end +_branch_flows(v, net) = real.(branch_power(net, v)) +_branch_currents_mag(v, net) = abs.(branch_current(net, v)) function _perturbed_voltage(state::ACPowerFlowState, param::Symbol, branch_idx::Int, epsilon::Float64, p_target, q_target) @@ -65,7 +43,7 @@ function _perturbed_voltage(state::ACPowerFlowState, param::Symbol, branch_idx:: end Y_pert = admittance_matrix(net_pert) v_pert = solve_pf_pq(Y_pert, state.v, p_target, q_target, state.idx_slack) - return v_pert, Y_pert + return v_pert, Y_pert, net_pert end @testset "AC PF Topology Sensitivities (:g, :b)" begin @@ -73,10 +51,7 @@ end @testset "Finite-difference verification — case5" begin file = joinpath(pm_path, "case5.m") - pm_data = PowerModels.parse_file(file) - net_data = PowerModels.make_basic_network(pm_data) - PowerModels.compute_ac_pf!(net_data) - state = ACPowerFlowState(net_data) + state = load_ac_pf_state("case5.m") n = state.n m = state.m @@ -110,8 +85,8 @@ end for e in test_branches @testset "∂/∂$(param)_$e" begin - v_p, Y_p = _perturbed_voltage(state, param, e, ε, p_base, q_base) - v_m, Y_m = _perturbed_voltage(state, param, e, -ε, p_base, q_base) + v_p, Y_p, net_p = _perturbed_voltage(state, param, e, ε, p_base, q_base) + v_m, Y_m, net_m = _perturbed_voltage(state, param, e, -ε, p_base, q_base) fd_vm = (abs.(v_p) - abs.(v_m)) / (2ε) if norm(fd_vm) > 1e-6 @@ -125,15 +100,15 @@ end @test rel_err_va < 1e-2 end - fd_f = (_branch_flows(v_p, Y_p, state.branch_data) - - _branch_flows(v_m, Y_m, state.branch_data)) / (2ε) + fd_f = (_branch_flows(v_p, net_p) - + _branch_flows(v_m, net_m)) / (2ε) if norm(fd_f) > 1e-6 rel_err_f = norm(Matrix(S_f)[:, e] - fd_f) / norm(fd_f) @test rel_err_f < 1e-2 end - fd_im = (_branch_currents_mag(v_p, Y_p, state.branch_data) - - _branch_currents_mag(v_m, Y_m, state.branch_data)) / (2ε) + fd_im = (_branch_currents_mag(v_p, net_p) - + _branch_currents_mag(v_m, net_m)) / (2ε) if norm(fd_im) > 1e-6 rel_err_im = norm(Matrix(S_im)[:, e] - fd_im) / norm(fd_im) @test rel_err_im < 1e-2 @@ -143,12 +118,47 @@ end end end + @testset "Transformer, phase shift, and parallel-line finite differences" begin + buses = [ + pd_bus(1, 3; vmax=1.1, vmin=0.9), + pd_bus(2, 1; vmax=1.1, vmin=0.9), + pd_bus(3, 1; vmax=1.1, vmin=0.9), + ] + gens = [ + pd_gen(1, 1; pg=0.5, qmax=1.0, qmin=-1.0, vg=1.0, pmax=2.0, pmin=0.0, cost=(1.0, 1.0, 0.0)), + ] + branches = [ + pd_branch(1, 1, 2; br_r=0.01, br_x=0.10, br_b=0.02, rate_a=2.0, rate_b=2.0, rate_c=2.0, tap=1.05, shift=0.12, angmin=-π / 3, angmax=π / 3), + pd_branch(2, 1, 2; br_r=0.02, br_x=0.20, br_b=0.01, rate_a=2.0, rate_b=2.0, rate_c=2.0, tap=1.00, shift=0.00, angmin=-π / 3, angmax=π / 3), + pd_branch(3, 2, 3; br_r=0.01, br_x=0.15, br_b=0.03, rate_a=2.0, rate_b=2.0, rate_c=2.0, tap=0.97, shift=-0.08, angmin=-π / 3, angmax=π / 3), + ] + net = ACNetwork(pd_case(buses, gens, branches; name="topology_fd")) + state = ACPowerFlowState(net, [1.01 + 0.02im, 0.98 - 0.04im, 1.02 + 0.01im]) + non_slack = [i for i in 1:state.n if i != state.idx_slack] + injections = state.v .* conj.(state.Y * state.v) + p_target = real.(injections)[non_slack] + q_target = imag.(injections)[non_slack] + ε = 1e-6 + + for param in (:g, :b) + S_vm = calc_sensitivity(state, :vm, param) + S_va = calc_sensitivity(state, :va, param) + S_f = calc_sensitivity(state, :f, param) + S_im = calc_sensitivity(state, :im, param) + for e in 1:net.m + v_p, _, net_p = _perturbed_voltage(state, param, e, ε, p_target, q_target) + v_m, _, net_m = _perturbed_voltage(state, param, e, -ε, p_target, q_target) + @test Matrix(S_vm)[:, e] ≈ (abs.(v_p) - abs.(v_m)) / (2ε) atol=1e-6 rtol=1e-4 + @test Matrix(S_va)[:, e] ≈ (angle.(v_p) - angle.(v_m)) / (2ε) atol=1e-6 rtol=1e-4 + @test Matrix(S_f)[:, e] ≈ (_branch_flows(v_p, net_p) - _branch_flows(v_m, net_m)) / (2ε) atol=1e-6 rtol=1e-4 + @test Matrix(S_im)[:, e] ≈ (_branch_currents_mag(v_p, net_p) - _branch_currents_mag(v_m, net_m)) / (2ε) atol=1e-6 rtol=1e-4 + end + end + end + @testset "Smoke tests — all 10 combinations" begin file = joinpath(pm_path, "case5.m") - pm_data = PowerModels.parse_file(file) - net_data = PowerModels.make_basic_network(pm_data) - PowerModels.compute_ac_pf!(net_data) - state = ACPowerFlowState(net_data) + state = load_ac_pf_state("case5.m") n = state.n m = state.m @@ -175,10 +185,7 @@ end @testset "Error: state without ACNetwork" begin file = joinpath(pm_path, "case5.m") - pm_data = PowerModels.parse_file(file) - net_data = PowerModels.make_basic_network(pm_data) - PowerModels.compute_ac_pf!(net_data) - full_state = ACPowerFlowState(net_data) + full_state = load_ac_pf_state("case5.m") raw_state = ACPowerFlowState(full_state.v, full_state.Y; idx_slack=full_state.idx_slack, branch_data=full_state.branch_data) @@ -189,9 +196,7 @@ end @testset "Non-basic network — col_to_id maps to branch IDs" begin file = joinpath(pm_path, "case5.m") - pm_data = PowerModels.parse_file(file) - PowerModels.compute_ac_pf!(pm_data) - state = ACPowerFlowState(pm_data) + state = load_ac_pf_state("case5.m") S = calc_sensitivity(state, :vm, :g) @test S.formulation == :acpf @@ -202,10 +207,7 @@ end @testset "Sensitivity metadata" begin file = joinpath(pm_path, "case14.m") - pm_data = PowerModels.parse_file(file) - net_data = PowerModels.make_basic_network(pm_data) - PowerModels.compute_ac_pf!(net_data) - state = ACPowerFlowState(net_data) + state = load_ac_pf_state("case14.m") for param in (:g, :b) S = calc_sensitivity(state, :vm, param) diff --git a/test/test_apf_integration.jl b/test/test_apf_integration.jl index efb4e9b..d00c75d 100644 --- a/test/test_apf_integration.jl +++ b/test/test_apf_integration.jl @@ -281,7 +281,7 @@ end @test_skip false else # case5.m has non-sequential bus IDs when not made basic - raw = PowerModels.parse_file(case_path) + raw = PowerDiff.parse_file(case_path) dc_net = DCNetwork(raw) # non-basic network apf_net = to_apf_network(dc_net) diff --git a/test/test_dc_opf_verification.jl b/test/test_dc_opf_verification.jl index 6297e89..9e2730c 100644 --- a/test/test_dc_opf_verification.jl +++ b/test/test_dc_opf_verification.jl @@ -234,9 +234,8 @@ end sol_base.pg[i] < dc_net.gmax[i] - 0.01, 1:dc_net.k) @test !isnothing(gen_idx) - net_pert = load_test_case("case5.m") - net_pert["gen"][string(gen_idx)]["cost"][2] += delta - dc_net_pert = DCNetwork(net_pert) + dc_net_pert = deepcopy(dc_net) + dc_net_pert.cl[gen_idx] += delta prob_pert = DCOPFProblem(dc_net_pert, d) sol_pert = solve!(prob_pert) lmp_pert = calc_lmp(sol_pert, dc_net_pert) @@ -260,9 +259,8 @@ end sol_base.pg[i] < dc_net.gmax[i] - 0.01, 1:dc_net.k) @test !isnothing(gen_idx) - net_pert = load_test_case("case5.m") - net_pert["gen"][string(gen_idx)]["cost"][1] += delta - dc_net_pert = DCNetwork(net_pert) + dc_net_pert = deepcopy(dc_net) + dc_net_pert.cq[gen_idx] += delta prob_pert = DCOPFProblem(dc_net_pert, d) sol_pert = solve!(prob_pert) lmp_pert = calc_lmp(sol_pert, dc_net_pert) diff --git a/test/test_jvp_vjp.jl b/test/test_jvp_vjp.jl index 5147587..8b001cc 100644 --- a/test/test_jvp_vjp.jl +++ b/test/test_jvp_vjp.jl @@ -22,8 +22,8 @@ # round-trip dict_to_vec/vec_to_dict, and error handling for invalid IDs. @testset "JVP / VJP" begin - raw = PowerModels.parse_file(joinpath(PM_DATA_DIR, "case5.m")) - basic = PowerModels.make_basic_network(deepcopy(raw)) + raw = PowerDiff._network_data(PowerDiff.parse_file(joinpath(PM_DATA_DIR, "case5.m"))) + basic = _make_basic_case(raw) # ================================================================= # Basic network JVP @@ -177,9 +177,7 @@ # Complex sensitivity (AC PF :v operand) # ================================================================= @testset "Complex sensitivity JVP/VJP" begin - pf_data = deepcopy(basic) - PowerModels.compute_ac_pf!(pf_data) - state = ACPowerFlowState(pf_data) + state = load_ac_pf_state("case5.m") S = calc_sensitivity(state, :v, :p) diff --git a/test/test_kkt_vjp_jvp.jl b/test/test_kkt_vjp_jvp.jl index 6c05600..5cb82cf 100644 --- a/test/test_kkt_vjp_jvp.jl +++ b/test/test_kkt_vjp_jvp.jl @@ -20,9 +20,8 @@ # match the materialized Sensitivity matrix path. @testset "KKT VJP/JVP" begin - basic = PowerModels.make_basic_network( - PowerModels.parse_file(joinpath(PM_DATA_DIR, "case5.m"))) - raw = PowerModels.parse_file(joinpath(PM_DATA_DIR, "case5.m")) + raw = PowerDiff._network_data(PowerDiff.parse_file(joinpath(PM_DATA_DIR, "case5.m"))) + basic = _make_basic_case(raw) # ================================================================= # DC OPF @@ -174,7 +173,7 @@ # AC OPF # ================================================================= @testset "AC OPF" begin - ac_data = PowerModels.parse_file(joinpath(PM_DATA_DIR, "case5.m")) + ac_data = PowerDiff.parse_file(joinpath(PM_DATA_DIR, "case5.m")) ac_prob = ACOPFProblem(ac_data) solve!(ac_prob) @@ -216,11 +215,7 @@ @testset "DC PF" begin net = DCNetwork(basic) n = net.n - d = zeros(n) - for (_, load) in basic["load"] - bus = load["load_bus"] - d[bus] += load["pd"] - end + d = calc_demand_vector(basic) state = DCPowerFlowState(net, d) @testset "VJP matches S' * w" begin @@ -282,14 +277,14 @@ # Slow path (no cache) — AC OPF # ================================================================= @testset "AC OPF slow path (no cache)" begin - ac_data2 = PowerModels.parse_file(joinpath(PM_DATA_DIR, "case5.m")) + ac_data2 = PowerDiff.parse_file(joinpath(PM_DATA_DIR, "case5.m")) ac_ref = ACOPFProblem(ac_data2); solve!(ac_ref) S = calc_sensitivity(ac_ref, :lmp, :d) w = randn(size(S, 1)) v = randn(size(S, 2)) # Fresh problem - ac_data3 = PowerModels.parse_file(joinpath(PM_DATA_DIR, "case5.m")) + ac_data3 = PowerDiff.parse_file(joinpath(PM_DATA_DIR, "case5.m")) ac_slow = ACOPFProblem(ac_data3); solve!(ac_slow) @test isnothing(ac_slow.cache.dz_dd) @@ -325,9 +320,7 @@ @test_throws ArgumentError jvp(prob, :lmp, :d, Dict(9999 => 1.0)) # ACPowerFlowState — not supported, gives helpful error - pf_data = deepcopy(basic) - PowerModels.compute_ac_pf!(pf_data) - ac_state = ACPowerFlowState(pf_data) + ac_state = load_ac_pf_state("case5.m") @test_throws ArgumentError vjp(ac_state, :vm, :p, randn(5)) @test_throws ArgumentError jvp(ac_state, :vm, :p, randn(5)) end diff --git a/test/test_non_matpower_parsers.jl b/test/test_non_matpower_parsers.jl new file mode 100644 index 0000000..417070e --- /dev/null +++ b/test/test_non_matpower_parsers.jl @@ -0,0 +1,121 @@ +using PowerIO + +function _assert_network_tables_compatible(actual, baseline; label, demand_rtol=1e-8) + @testset "$label tables" begin + @test length(actual.bus) == length(baseline.bus) + @test length(actual.gen) == length(baseline.gen) + @test length(actual.branch) == length(baseline.branch) + @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]) + end +end + +function _assert_constructs_and_solves(data; label) + @testset "$label constructors" begin + dc = DCNetwork(data) + @test dc.n == length(data.bus) + @test dc.k == length(data.gen) + @test dc.m == length(data.branch) + @test all(isfinite, dc.b) + + dc_prob = DCOPFProblem(dc) + dc_sol = solve!(dc_prob) + @test all(isfinite, dc_sol.va) + @test all(isfinite, dc_sol.pg) + @test all(isfinite, dc_sol.f) + + ac = ACNetwork(data) + @test size(admittance_matrix(ac)) == (ac.n, ac.n) + ac_prob = ACOPFProblem(ac; silent=true) + @test ac_prob.n_gen == length(data.gen) + end +end + +@testset "Non MATPOWER Parser Support" begin + matpower_net = PowerDiff.parse_file("pglib_opf_case14_ieee.m"; library=:pglib) + baseline = PowerDiff._network_data(matpower_net) + + @testset "PowerModels JSON" begin + text, warnings = PowerIO.to_format(matpower_net, "powermodels-json") + @test isempty(warnings) + + parsed = PowerDiff.parse_file(IOBuffer(text); from=:powermodels) + @test PowerIO.source_format(parsed) == "PowerModelsJson" + 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) + parsed_path = PowerDiff.parse_file(path; from="powermodels-json") + @test PowerIO.source_format(parsed_path) == "PowerModelsJson" + _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) + + parsed = PowerDiff.parse_file(IOBuffer(text); from=:egret) + @test PowerIO.source_format(parsed) == "EgretJson" + data = PowerDiff._network_data(parsed) + _assert_network_tables_compatible(data, baseline; label="Egret JSON") + _assert_constructs_and_solves(data; label="Egret JSON") + + @test_throws ArgumentError PowerDiff.parse_file(IOBuffer(text); filetype="json") + end + + @testset "PSS/E RAW" begin + text, warnings = PowerIO.to_format(matpower_net, "psse") + @test warnings isa AbstractVector + + parsed = PowerDiff.parse_file(IOBuffer(text); from=:psse) + @test PowerIO.source_format(parsed) == "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) + parsed_path = PowerDiff.parse_file(path) + @test PowerIO.source_format(parsed_path) == "Psse" + _assert_network_tables_compatible( + PowerDiff._network_data(parsed_path), baseline; label="PSS/E RAW path", demand_rtol=1e-5) + end + end +end + +@testset "Programmatic AC OPF constructor validation" begin + buses = [pd_bus(10, 1), pd_bus(20, 1)] + gens = [pd_gen(1, 10; pg=0.0, qmax=1.0, qmin=-1.0, pmax=2.0, pmin=0.0, cost=(1.0, 1.0, 0.0))] + branches = [pd_branch(1, 10, 20; br_r=0.01, br_x=0.1, rate_a=2.0)] + data = pd_case(buses, gens, branches; name="no_ref_bus") + + net = ACNetwork(data; idx_slack=2) + @test net.idx_slack == 2 + @test net.ref_bus_keys == [2] + prob = ACOPFProblem(net; silent=true) + @test prob.data.ref_bus_keys == [2] + + typed_ref_data = pd_case( + [pd_bus(10, 3), pd_bus(20, 1)], + gens, + branches; + name="typed_ref_bus", + ) + override_net = ACNetwork(typed_ref_data; idx_slack=2) + @test override_net.idx_slack == 2 + @test override_net.ref_bus_keys == [2] + @test ACOPFProblem(override_net; silent=true).data.ref_bus_keys == [2] + @test_throws ArgumentError ACNetwork(typed_ref_data; idx_slack=0) + @test_throws ArgumentError ACNetwork(typed_ref_data; idx_slack=3) + + y = ComplexF64[1 -1; -1 1] + raw_net = ACNetwork(y) + @test_throws ArgumentError ACOPFProblem(raw_net; silent=true) +end diff --git a/test/test_nonbasic.jl b/test/test_nonbasic.jl index 8e96400..ed2c0f5 100644 --- a/test/test_nonbasic.jl +++ b/test/test_nonbasic.jl @@ -21,9 +21,18 @@ # arbitrary element IDs. Uses case5.m with bus IDs [1,2,3,4,10] — bus 10 # maps to sequential index 5 via IDMapping. +# Renumber bus ids to a dense 1..n space, operating on PowerDiff network tables. +function _make_basic_case(data) + bus_map = Dict(id => i for (i, id) in enumerate(sort([b.bus_i for b in data.bus]))) + buses = [(; b..., bus_i=bus_map[b.bus_i]) for b in data.bus] + gens = [(; g..., gen_bus=bus_map[g.gen_bus]) for g in data.gen] + branches = [(; br..., f_bus=bus_map[br.f_bus], t_bus=bus_map[br.t_bus]) for br in data.branch] + return (; data.name, data.baseMVA, bus=buses, gen=gens, branch=branches) +end + @testset "Non-Basic Network Support" begin - raw = PowerModels.parse_file(joinpath(PM_DATA_DIR, "case5.m")) - basic = PowerModels.make_basic_network(deepcopy(raw)) + raw = PowerDiff._network_data(PowerDiff.parse_file(joinpath(PM_DATA_DIR, "case5.m"))) + basic = _make_basic_case(raw) # ================================================================= # IDMapping @@ -150,7 +159,7 @@ # Generator sensitivity: row_to_id should be gen IDs dpg_dd_nb = calc_sensitivity(prob_nb, :pg, :d) - @test dpg_dd_nb.row_to_id == sort(collect(parse(Int, k) for k in keys(raw["gen"]))) + @test dpg_dd_nb.row_to_id == sort([gen.index for gen in raw.gen]) # Branch sensitivity: row_to_id should be branch IDs df_dsw_nb = calc_sensitivity(prob_nb, :f, :sw) @@ -203,10 +212,7 @@ # AC Power Flow # ================================================================= @testset "AC power flow non-basic" begin - # Solve AC PF using PowerModels on non-basic network - pf_data_nb = deepcopy(raw) - PowerModels.compute_ac_pf!(pf_data_nb) - state_nb = ACPowerFlowState(pf_data_nb) + state_nb = load_ac_pf_state("case5.m") # Voltage sensitivity should be finite dvm_dp = calc_sensitivity(state_nb, :vm, :p) @@ -244,7 +250,7 @@ dpg_dsw = calc_sensitivity(prob_nb, :pg, :sw) @test all(isfinite, Matrix(dpg_dsw)) - @test dpg_dsw.row_to_id == sort(collect(parse(Int, k) for k in keys(raw["gen"]))) + @test dpg_dsw.row_to_id == sort([gen.index for gen in raw.gen]) end # ================================================================= @@ -351,37 +357,24 @@ end # ================================================================= - # Stored ref on network types + # Typed cached network data # ================================================================= - @testset "Stored ref field" begin + @testset "Typed cached network data" begin dc_nb = DCNetwork(raw) - @test !isnothing(dc_nb.ref) - @test haskey(dc_nb.ref, :bus) + @test length(dc_nb.demand) == dc_nb.n + @test length(dc_nb.pg_init) == dc_nb.n ac_nb = ACNetwork(raw) - @test !isnothing(ac_nb.ref) - @test haskey(ac_nb.ref, :bus) + @test length(ac_nb.pd) == ac_nb.n + @test length(ac_nb.gen_bus) == length(ac_nb.pmin) - # Programmatic constructor has ref=nothing + # Programmatic constructor carries typed zero vectors n, m, k = 2, 1, 1 A = sparse([1.0 -1.0]) G_inc = sparse(reshape([1.0, 0.0], 2, 1)) b = [-10.0] dc_prog = DCNetwork(n, m, k, A, G_inc, b) - @test isnothing(dc_prog.ref) - end - - # ================================================================= - # IDMapping shunt support - # ================================================================= - @testset "IDMapping shunt fields" begin - dc_nb = DCNetwork(raw) - id_map = dc_nb.id_map - - @test isa(id_map.shunt_ids, Vector{Int}) - @test isa(id_map.shunt_to_idx, Dict{Int,Int}) - # shunt_ids should be sorted - @test issorted(id_map.shunt_ids) + @test dc_prog.demand == zeros(n) end # ================================================================= @@ -445,103 +438,34 @@ @testset "IDMapping constructor validation" begin # Unsorted bus_ids should throw @test_throws ArgumentError IDMapping( - [3, 1, 2], [1, 2], [1], [1], Int[], - Dict(3=>1, 1=>2, 2=>3), Dict(1=>1, 2=>2), Dict(1=>1), Dict(1=>1), Dict{Int,Int}()) + [3, 1, 2], [1, 2], [1], + Dict(3=>1, 1=>2, 2=>3), Dict(1=>1, 2=>2), Dict(1=>1)) # Unsorted branch_ids should throw @test_throws ArgumentError IDMapping( - [1, 2, 3], [2, 1], [1], [1], Int[], - Dict(1=>1, 2=>2, 3=>3), Dict(2=>1, 1=>2), Dict(1=>1), Dict(1=>1), Dict{Int,Int}()) + [1, 2, 3], [2, 1], [1], + Dict(1=>1, 2=>2, 3=>3), Dict(2=>1, 1=>2), Dict(1=>1)) # Unsorted gen_ids should throw @test_throws ArgumentError IDMapping( - [1, 2], [1], [3, 1], [1], Int[], - Dict(1=>1, 2=>2), Dict(1=>1), Dict(3=>1, 1=>2), Dict(1=>1), Dict{Int,Int}()) - - # Unsorted load_ids should throw - @test_throws ArgumentError IDMapping( - [1, 2], [1], [1], [5, 2], Int[], - Dict(1=>1, 2=>2), Dict(1=>1), Dict(1=>1), Dict(5=>1, 2=>2), Dict{Int,Int}()) - - # Unsorted shunt_ids should throw - @test_throws ArgumentError IDMapping( - [1, 2], [1], [1], [1], [3, 1], - Dict(1=>1, 2=>2), Dict(1=>1), Dict(1=>1), Dict(1=>1), Dict(3=>1, 1=>2)) + [1, 2], [1], [3, 1], + Dict(1=>1, 2=>2), Dict(1=>1), Dict(3=>1, 1=>2)) # Length mismatch: bus_ids vs bus_to_idx @test_throws ArgumentError IDMapping( - [1, 2, 3], [1], [1], [1], Int[], - Dict(1=>1, 2=>2), Dict(1=>1), Dict(1=>1), Dict(1=>1), Dict{Int,Int}()) + [1, 2, 3], [1], [1], + Dict(1=>1, 2=>2), Dict(1=>1), Dict(1=>1)) # Length mismatch: branch_ids vs branch_to_idx @test_throws ArgumentError IDMapping( - [1, 2], [1, 2], [1], [1], Int[], - Dict(1=>1, 2=>2), Dict(1=>1), Dict(1=>1), Dict(1=>1), Dict{Int,Int}()) + [1, 2], [1, 2], [1], + Dict(1=>1, 2=>2), Dict(1=>1), Dict(1=>1)) # Valid construction should work id_map = IDMapping( - [1, 5, 10], [1, 2], [1], [1], Int[], - Dict(1=>1, 5=>2, 10=>3), Dict(1=>1, 2=>2), Dict(1=>1), Dict(1=>1), Dict{Int,Int}()) + [1, 5, 10], [1, 2], [1], + Dict(1=>1, 5=>2, 10=>3), Dict(1=>1, 2=>2), Dict(1=>1)) @test id_map.bus_ids == [1, 5, 10] @test id_map.bus_to_idx[10] == 3 end - - # ================================================================= - # _remap_ref_to_sequential - # ================================================================= - @testset "_remap_ref_to_sequential" begin - # Build ref from non-basic network - pm_data = deepcopy(raw) - PowerModels.standardize_cost_terms!(pm_data, order=2) - PowerModels.calc_thermal_limits!(pm_data) - ref = PowerModels.build_ref(pm_data)[:it][:pm][:nw][0] - id_map = IDMapping(ref) - - seq_ref = PowerDiff._remap_ref_to_sequential(ref, id_map) - - # Sequential ref should have keys 1:n for buses - n = length(id_map.bus_ids) - @test sort(collect(keys(seq_ref[:bus]))) == collect(1:n) - - # Sequential ref should have keys 1:m for branches - m = length(id_map.branch_ids) - @test sort(collect(keys(seq_ref[:branch]))) == collect(1:m) - - # Sequential ref should have keys 1:k for generators - k = length(id_map.gen_ids) - @test sort(collect(keys(seq_ref[:gen]))) == collect(1:k) - - # Bus-level fields should be remapped - for (seq_idx, bus_data) in seq_ref[:bus] - @test bus_data["bus_i"] == seq_idx - end - - # Branch f_bus / t_bus should be sequential - for (seq_idx, br_data) in seq_ref[:branch] - @test 1 <= br_data["f_bus"] <= n - @test 1 <= br_data["t_bus"] <= n - end - - # Generator gen_bus should be sequential - for (seq_idx, gen_data) in seq_ref[:gen] - @test 1 <= gen_data["gen_bus"] <= n - end - - # Arcs should all use sequential IDs - for (l, i, j) in seq_ref[:arcs] - @test 1 <= l <= m - @test 1 <= i <= n - @test 1 <= j <= n - end - - # bus_gens keys should be sequential bus IDs - for bus_idx in keys(seq_ref[:bus_gens]) - @test 1 <= bus_idx <= n - end - - # ref_buses should be remapped - for ref_bus in keys(seq_ref[:ref_buses]) - @test 1 <= ref_bus <= n - end - end end diff --git a/test/test_parameter_transforms.jl b/test/test_parameter_transforms.jl index 9e7b3c6..428fec4 100644 --- a/test/test_parameter_transforms.jl +++ b/test/test_parameter_transforms.jl @@ -26,13 +26,10 @@ using Test @testset "Parameter Transforms" begin pm_path = joinpath(dirname(pathof(PowerModels)), "..", "test", "data", "matpower") file = joinpath(pm_path, "case5.m") - pm_data = PowerModels.parse_file(file) - net_data = PowerModels.make_basic_network(pm_data) + net_data = PowerDiff.parse_file(file) @testset "AC PF: ∂/∂d = -∂/∂p (algebraic)" begin - pf_data = deepcopy(net_data) - PowerModels.compute_ac_pf!(pf_data) - state = ACPowerFlowState(pf_data) + state = load_ac_pf_state("case5.m") # Test all operands that work with :p for op in [:vm, :v, :im, :va, :f] @@ -45,9 +42,7 @@ using Test end @testset "AC PF: ∂/∂qd = -∂/∂q (algebraic)" begin - pf_data = deepcopy(net_data) - PowerModels.compute_ac_pf!(pf_data) - state = ACPowerFlowState(pf_data) + state = load_ac_pf_state("case5.m") for op in [:vm, :v, :im, :va, :f] S_q = calc_sensitivity(state, op, :q) @@ -58,9 +53,7 @@ using Test end @testset "AC PF: Jacobian blocks + demand transform is invalid" begin - pf_data = deepcopy(net_data) - PowerModels.compute_ac_pf!(pf_data) - state = ACPowerFlowState(pf_data) + state = load_ac_pf_state("case5.m") # :p and :q as operands don't have :p/:q as native params, # so the :d/:qd transforms cannot apply @@ -95,9 +88,7 @@ using Test end @testset "Introspection includes transform-derived symbols" begin - pf_data = deepcopy(net_data) - PowerModels.compute_ac_pf!(pf_data) - state = ACPowerFlowState(pf_data) + state = load_ac_pf_state("case5.m") params = parameter_symbols(state) @test :d in params diff --git a/test/test_parser_parity.jl b/test/test_parser_parity.jl new file mode 100644 index 0000000..180c7bb --- /dev/null +++ b/test/test_parser_parity.jl @@ -0,0 +1,174 @@ +import PowerIO + +const _INLINE_CASE = """ +function mpc = case_inline +mpc.version = '2'; +mpc.baseMVA = 100; +mpc.bus = [1 2 50 10 1 -2 1 1.0 0 230 1 1.1 0.9; 2 1 0 0 0 0 1 1.0 0 230 1 1.1 0.9]; +mpc.gen = [1 80 0 100 -100 1 100 1 150 0; 2 20 0 50 -50 1 100 0 50 0]; +mpc.branch = [1 2 0.01 0.1 0.02 0 0 0 0 0 0 -360 360; 1 2 0.02 0.2 0.01 100 100 100 1 0 1 -60 60]; +mpc.gencost = [2 0 0 3 0.01 2 3; 2 0 0 3 0.02 3 4]; +mpc.areas = [1 1]; +mpc.bus_name = ['one'; 'two']; +""" + +# `parse_file`/`parse_matpower` return a PowerIO.Network; `_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. +_inline_data() = PowerDiff._network_data(PowerDiff.parse_matpower(IOBuffer(_INLINE_CASE))) + +@testset "MATPOWER Parser Semantics" begin + @testset "Inline arrays and normalization" begin + data = _inline_data() + + @test data isa NamedTuple + @test data.name == "case_inline" + @test data.baseMVA == 100.0 + @test length(data.bus) == 2 + @test length(data.gen) == 1 + @test length(data.branch) == 1 + @test data.gen[1].index == 1 + @test data.branch[1].index == 2 + @test data.bus[1].bus_type == 3 + # 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.branch[1].tap == 1.0 + @test data.branch[1].rate_a > 0 + @test data.branch[1].angmin ≈ -π / 3 + @test data.branch[1].angmax ≈ π / 3 + @test data.gen[1].cost == (100.0, 200.0, 3.0) + end + + @testset "Multiline arrays and artifact path" begin + parsed = PowerDiff.parse_file("pglib_opf_case14_ieee.m"; library=:pglib) + @test parsed isa PowerIO.Network + nd = PowerDiff._network_data(parsed) + @test length(nd.bus) == 14 + @test length(nd.branch) == 20 + @test PowerDiff.get_path(:pglib) == PD_PGLIB_DIR + end + + @testset "Rejected inputs" begin + @test_throws ArgumentError PowerDiff.parse_file("case.json") + @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];") + @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") + @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") + @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.Network + @test_throws ArgumentError PowerDiff.parse_file(IOBuffer(_INLINE_CASE); backend=:native) + end +end + +# Field-for-field equality of two PowerDiff network tables; floats with ≈, ints with ==. +function _assert_netdata_equal(a, b, label) + @testset "$label" begin + @test a.baseMVA ≈ b.baseMVA + @test length(a.bus) == length(b.bus) + @test length(a.gen) == length(b.gen) + @test length(a.branch) == length(b.branch) + for (x, y) in zip(a.bus, b.bus) + @test x.bus_i == y.bus_i + @test x.bus_type == y.bus_type + @test x.pd ≈ y.pd && x.qd ≈ y.qd && x.gs ≈ y.gs && x.bs ≈ y.bs + @test x.vm ≈ y.vm && x.va ≈ y.va && x.vmin ≈ y.vmin && x.vmax ≈ y.vmax + end + for (x, y) in zip(a.gen, b.gen) + @test x.gen_bus == y.gen_bus + @test x.pg ≈ y.pg && x.qg ≈ y.qg + @test x.pmin ≈ y.pmin && x.pmax ≈ y.pmax && x.qmin ≈ y.qmin && x.qmax ≈ y.qmax + @test all(x.cost .≈ y.cost) + 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.rate_a ≈ y.rate_a + @test x.tap ≈ y.tap && x.shift ≈ y.shift && x.angmin ≈ y.angmin && x.angmax ≈ y.angmax + end + end +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 + end +end + +@testset "Typed AC Pi Model" begin + buses = [ + pd_bus(1, 3; vmax=1.1, vmin=0.9), + pd_bus(2, 1; vmax=1.1, vmin=0.9), + pd_bus(3, 1; vmax=1.1, vmin=0.9), + ] + gens = [ + pd_gen(1, 1; pg=0.5, qmax=1.0, qmin=-1.0, vg=1.0, pmax=2.0, pmin=0.0, cost=(1.0, 1.0, 0.0)), + ] + branches = [ + pd_branch(1, 1, 2; br_r=0.01, br_x=0.10, br_b=0.02, rate_a=2.0, rate_b=2.0, rate_c=2.0, tap=1.05, shift=0.12, angmin=-π / 3, angmax=π / 3), + pd_branch(2, 1, 2; br_r=0.02, br_x=0.20, br_b=0.01, rate_a=2.0, rate_b=2.0, rate_c=2.0, tap=1.00, shift=0.00, angmin=-π / 3, angmax=π / 3), + pd_branch(3, 2, 3; br_r=0.01, br_x=0.15, br_b=0.03, rate_a=2.0, rate_b=2.0, rate_c=2.0, tap=0.97, shift=-0.08, angmin=-π / 3, angmax=π / 3), + ] + data = pd_case(buses, gens, branches; name="pi_model") + net = ACNetwork(data) + v = [1.01 + 0.02im, 0.98 - 0.04im, 1.02 + 0.01im] + + rows = Int[] + cols = Int[] + vals = ComplexF64[] + expected_current = ComplexF64[] + for l in 1:net.m + y = net.g[l] + im * net.b[l] + tap = net.tap[l] * cis(net.shift[l]) + yff = (y + net.g_fr[l] + im * net.b_fr[l]) / abs2(tap) + yft = -y / conj(tap) + ytf = -y / tap + ytt = y + net.g_to[l] + im * net.b_to[l] + append!(rows, (net.f_bus[l], net.f_bus[l], net.t_bus[l], net.t_bus[l])) + append!(cols, (net.f_bus[l], net.t_bus[l], net.f_bus[l], net.t_bus[l])) + append!(vals, (yff, yft, ytf, ytt)) + push!(expected_current, yff * v[net.f_bus[l]] + yft * v[net.t_bus[l]]) + end + expected_y = sparse(rows, cols, vals, net.n, net.n) + + @test Matrix(admittance_matrix(net)) ≈ Matrix(expected_y) + @test branch_current(net, v) ≈ expected_current + @test branch_power(net, v) ≈ v[net.f_bus] .* conj.(expected_current) +end + +@testset "Status Filtering" begin + parsed = _inline_data() + @test length(parsed.gen) == 1 + @test length(parsed.branch) == 1 +end diff --git a/test/test_sensitivity_column.jl b/test/test_sensitivity_column.jl index 1517df5..b8dbc40 100644 --- a/test/test_sensitivity_column.jl +++ b/test/test_sensitivity_column.jl @@ -126,8 +126,7 @@ @testset "AC PF" begin pm_data = load_test_case("case5.m") @test !isnothing(pm_data) - ac_net = ACNetwork(pm_data) - ac_state = ACPowerFlowState(pm_data) + ac_state = load_ac_pf_state("case5.m") for (op, param) in [(:vm, :p), (:va, :q), (:im, :p), (:f, :q)] S = calc_sensitivity(ac_state, op, param) @@ -194,7 +193,7 @@ # Non-basic network (arbitrary element IDs) # ========================================================================= @testset "Non-basic network" begin - pm_data = load_raw_case("case5.m") + pm_data = load_test_case("case5.m") @test !isnothing(pm_data) # case5.m has bus IDs [1, 2, 3, 4, 10] diff --git a/test/test_sensitivity_coverage.jl b/test/test_sensitivity_coverage.jl index c2fd701..a994ca7 100644 --- a/test/test_sensitivity_coverage.jl +++ b/test/test_sensitivity_coverage.jl @@ -27,8 +27,7 @@ using Test # Load test case pm_path = joinpath(dirname(pathof(PowerModels)), "..", "test", "data", "matpower") file = joinpath(pm_path, "case5.m") - pm_data = PowerModels.parse_file(file) - net_data = PowerModels.make_basic_network(pm_data) + net_data = PowerDiff.parse_file(file) @testset "DC Power Flow — all 6 combinations" begin net = DCNetwork(net_data) @@ -98,9 +97,7 @@ using Test @testset "AC Power Flow — 24 native + 10 transform = 34 combinations" begin # Solve AC power flow first - pf_data = deepcopy(net_data) - PowerModels.compute_ac_pf!(pf_data) - state = ACPowerFlowState(pf_data) + state = load_ac_pf_state("case5.m") # 24 native combinations native_combos = [ diff --git a/test/test_update_switching.jl b/test/test_update_switching.jl index 625179c..f259dea 100644 --- a/test/test_update_switching.jl +++ b/test/test_update_switching.jl @@ -96,20 +96,18 @@ # ========================================================================= @testset "DC OPF: open lines do not enforce angle bounds" begin - dc_data = deepcopy(net_data) + dc_net = DCNetwork(net_data) # case5 has parallel lines 5 and 6 between the same buses - dc_data["branch"]["5"]["angmin"] = 1.0 - dc_data["branch"]["5"]["angmax"] = 2.0 - dc_data["branch"]["6"]["angmin"] = -2.0 - dc_data["branch"]["6"]["angmax"] = -1.0 - - dc_net = DCNetwork(dc_data) + dc_net.angmin[5] = 1.0 + dc_net.angmax[5] = 2.0 + dc_net.angmin[6] = -2.0 + dc_net.angmax[6] = -1.0 sw_new = copy(dc_net.sw) sw_new[5] = 0.0 sw_new[6] = 0.0 dc_net.sw .= sw_new - prob = DCOPFProblem(dc_net, calc_demand_vector(dc_data)) + prob = DCOPFProblem(dc_net, calc_demand_vector(net_data)) sol = solve!(prob) @test all(isfinite, sol.va) @@ -170,14 +168,14 @@ # ========================================================================= @testset "AC OPF: open lines do not enforce angle bounds" begin - ac_data = deepcopy(net_data) + ac_net = ACNetwork(net_data) # case5 has parallel lines 5 and 6 between the same buses - ac_data["branch"]["5"]["angmin"] = 1.0 - ac_data["branch"]["5"]["angmax"] = 2.0 - ac_data["branch"]["6"]["angmin"] = -2.0 - ac_data["branch"]["6"]["angmax"] = -1.0 + ac_net.angmin[5] = 1.0 + ac_net.angmax[5] = 2.0 + ac_net.angmin[6] = -2.0 + ac_net.angmax[6] = -1.0 - prob = ACOPFProblem(ac_data) + prob = ACOPFProblem(ac_net) sw_new = copy(prob.network.sw) sw_new[5] = 0.0 sw_new[6] = 0.0 diff --git a/test/unified/test_interface.jl b/test/unified/test_interface.jl index ee7d7d1..b06d49e 100644 --- a/test/unified/test_interface.jl +++ b/test/unified/test_interface.jl @@ -19,8 +19,8 @@ using Test @testset "Unified Architecture" begin # Load a test network case_path = joinpath(dirname(pathof(PowerModels)), "..", "test", "data", "matpower", "case5.m") - data = PowerModels.parse_file(case_path) - net_data = PowerModels.make_basic_network(data) + net_data = PowerDiff.parse_file(case_path) + pm_data = PowerModels.make_basic_network(PowerModels.parse_file(case_path)) @testset "Abstract Type Hierarchy" begin # Test that types inherit correctly @@ -140,11 +140,11 @@ using Test end @testset "ACNetwork" begin - # Construct from PowerModels network ac_net = ACNetwork(net_data) + nd = PowerDiff._network_data(net_data) @test ac_net isa AbstractPowerNetwork - @test ac_net.n == length(net_data["bus"]) - @test ac_net.m == length(net_data["branch"]) + @test ac_net.n == length(nd.bus) + @test ac_net.m == length(nd.branch) # Admittance matrix reconstruction Y = admittance_matrix(ac_net) @@ -159,10 +159,10 @@ using Test @testset "ACPowerFlowState" begin # Solve AC power flow - PowerModels.compute_ac_pf!(net_data) + PowerModels.compute_ac_pf!(pm_data) # Construct state from network - state = ACPowerFlowState(net_data) + state = ACPowerFlowState(ACNetwork(net_data), PowerModels.calc_basic_bus_voltage(pm_data)) @test state isa AbstractPowerFlowState @test length(state.v) == state.n @test size(state.Y) == (state.n, state.n) @@ -176,8 +176,8 @@ using Test @testset "AC Power Flow Sensitivity" begin # Solve AC power flow - PowerModels.compute_ac_pf!(net_data) - state = ACPowerFlowState(net_data) + PowerModels.compute_ac_pf!(pm_data) + state = ACPowerFlowState(ACNetwork(net_data), PowerModels.calc_basic_bus_voltage(pm_data)) dvm_dp = calc_sensitivity(state, :vm, :p) @test dvm_dp isa Sensitivity @@ -210,8 +210,8 @@ using Test @test operand_symbols(prob) == [:va, :pg, :f, :psh, :lmp] @test parameter_symbols(prob) == [:d, :sw, :cq, :cl, :fmax, :b] - PowerModels.compute_ac_pf!(net_data) - ac_state = ACPowerFlowState(net_data) + PowerModels.compute_ac_pf!(pm_data) + ac_state = ACPowerFlowState(ACNetwork(net_data), PowerModels.calc_basic_bus_voltage(pm_data)) @test Set(operand_symbols(ac_state)) == Set([:vm, :v, :im, :va, :f, :p, :q]) @test Set(parameter_symbols(ac_state)) == Set([:p, :q, :va, :vm, :d, :qd, :g, :b]) diff --git a/test/unified/test_sensitivity_verification.jl b/test/unified/test_sensitivity_verification.jl index 881fd83..d6d3905 100644 --- a/test/unified/test_sensitivity_verification.jl +++ b/test/unified/test_sensitivity_verification.jl @@ -27,8 +27,7 @@ using Test # Load a test network case_path = joinpath(dirname(pathof(PowerModels)), "..", "test", "data", "matpower", "case5.m") - data = PowerModels.parse_file(case_path) - net_data = PowerModels.make_basic_network(data) + net_data = PowerDiff.parse_file(case_path) @testset "DC Power Flow Switching Sensitivity" begin net = DCNetwork(net_data) @@ -254,8 +253,7 @@ using Test @testset "AC Voltage-Power Sensitivity" begin # Solve AC power flow - PowerModels.compute_ac_pf!(net_data) - state = ACPowerFlowState(net_data) + state = load_ac_pf_state("case5.m") # Get analytical sensitivity via new API dvm_dp = calc_sensitivity(state, :vm, :p)