Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 10 additions & 8 deletions .claude/CLAUDE.md
Original file line number Diff line number Diff line change
Expand Up @@ -65,9 +65,9 @@ Exported solver-model types and functions (see `src/PowerFlows.jl`):

**MCPB (mixed current-power-balance).** PQ buses use divided current balance (imag-first row order); PV buses use real-power balance + `|V|²` constraint with **only 2 vars/bus** (no Q state — the key difference from rectangular's 3). REF = (P_gen, Q_gen). System size is exactly 2n. Status: opt-in, NOT default; do not deprecate rectangular. Validated to polar parity. Performance: for NR/TR ≈ rectangular (no net win); for **LM, Mixed decisively beats Rectangular** (rectangular-LM fails to converge at ~10k buses; Mixed-LM converges like Polar-LM with the smallest 2n state). Jacobian kernels are called as concrete top-level functions, never stored as abstract `::Function` fields (that forces dynamic dispatch on the hot path).

**`validate_voltage_magnitudes`** exists for polar, rectangular, and mixed. For rect/mixed it checks squared bounds (`e²+f² ∈ [min²,max²]`) for PQ and PV (PV `(e,f)` are real state vars and `|V|²` can drift mid-iteration before the constraint row pins it); REF skipped. Toggle via `solver_settings[:validate_voltage_magnitudes]`.
**`validate_voltage_magnitudes`** exists for polar, rectangular, and mixed. For rect/mixed it checks squared bounds (`e²+f² ∈ [min²,max²]`) for PQ and PV (PV `(e,f)` are real state vars and `|V|²` can drift mid-iteration before the constraint row pins it); REF skipped. Toggle via the `validate_voltage_magnitudes` keyword of the formulation constructor (a `SolutionParameters` field; the old `solver_settings` Dict is gone).

**Linear-solver backends (PNM-owned).** PowerFlows no longer hand-rolls a KLU cache — KLU is **not** a direct dependency. Backends come from PNM (`PNM.KLULinSolveCache{Tv,Ti}`, `PNM.AAFactorCache` for AppleAccelerate) plus a PowerFlows-defined `PardisoLinSolveCache` in `ext/PowerFlowsPardisoExt.jl`. Selection = PNM preference default + per-solve kwarg (AC: `solver_settings[:linear_solver]`; DC: kwarg).
**Linear-solver backends (PNM-owned).** PowerFlows no longer hand-rolls a KLU cache — KLU is **not** a direct dependency. Backends come from PNM (`PNM.KLULinSolveCache{Tv,Ti}`, `PNM.AAFactorCache` for AppleAccelerate) plus a PowerFlows-defined `PardisoLinSolveCache` in `ext/PowerFlowsPardisoExt.jl`. Selection = PNM preference default + the `linear_solver` keyword (AC: a `SolutionParameters` field, resolved to a `String` at construction; DC: per-solve kwarg). `"Dense"` is rejected at resolution time.
- Index width gates the AppleAccelerate path: `INDEX_TYPE = @static Sys.isapple() ? Int64 : Int32` (drives `J_INDEX_TYPE`/`REC_INDEX_TYPE`). AA's Apple `libSparse` ABI needs Int64 `columnStarts`; KLU uses Int32 elsewhere. Any cache-type `Union` must list BOTH `KLULinSolveCache{Float64,Int32}` and `{…,Int64}` (PNM's DC ABA factorization is always Int64 — omitting it MethodErrors on Linux).
- The KLU and AA backends do NOT share generic functions: `PNM.solve!/full_factor!/symbolic_factor!/numeric_refactor!/tsolve!` are KLU-only (`.KLUWrapper`); AA's live in `PNM.AccelerateWrapper.*`. PowerFlows defines local dispatch over both cache types. `AAFactorCache` is Int-only with NO transpose solve, so voltage-stability factors (need Aᵀ\b) stay KLU-only.
- MKLPardiso is x86_64-only; `resolve_linear_solver_backend` rejects it when `Sys.ARCH !== :x86_64`; Pardiso tests gate on `Pardiso.mkl_is_available()`.
Expand All @@ -76,7 +76,7 @@ Exported solver-model types and functions (see `src/PowerFlows.jl`):

**LCC / network reduction.** PNM's zero-impedance reduction merges LCC-terminal buses; `data.lcc.arcs` stores REDUCED tuples. Post-processing (`get_lcc_names`, `arc_to_lcc`) must key with reduced tuples via `get_arc_tuple(PSY.get_arc(lcc), nrd)` where `nrd = PNM.get_network_reduction_data(data.power_network_matrix)` — keying with raw tuples is a KeyError. Degree-2 parity tests must build systems with `reduce_reactive_power_injectors=false` (the default drops susceptive-FA shunts).

**Perf NR-cache reuse (polar).** `PolarNRCache` (slot `data.polar_nr_cache::RefValue{Union{Nothing,AbstractNRCache}}`) reuses residual/Jacobian/symbolic factorization across Q-limit retries and time steps; LCC or a changed subnetwork/slack forces a rebuild. `data.solver_cache::RefValue{Union{Nothing,SolverCache}}` holds the analogous per-solve cache — `DCSolverCache` (DC/PTDF) or `FastDecoupledCache` (FDNR factor-once B′/B″); the getters dispatch on the concrete subtype, so a cross-use is a loud `MethodError`, not a silent mis-read (no sentinel tag).
**Perf NR-cache reuse (polar).** `PolarNRCache` (slot `data.polar_nr_cache::RefValue{Union{Nothing,AbstractNRCache}}`) reuses residual/Jacobian/symbolic factorization across Q-limit retries and time steps; LCC or a changed subnetwork/slack forces a rebuild. Residuals, Jacobians and `HomotopyHessian` never store `data` — it is passed explicitly (`R(data, x, t)`, `J(data, t)`, `ACPowerFlowJacobian(data, residual, t)`), because the cache hangs off `data` and a stored back-reference would form a cycle; `test_nr_cache_reuse.jl` guards this. `data.solver_cache::RefValue{Union{Nothing,SolverCache}}` holds the analogous per-solve cache — `DCSolverCache` (DC/PTDF) or `FastDecoupledCache` (FDNR factor-once B′/B″); the getters dispatch on the concrete subtype, so a cross-use is a loud `MethodError`, not a silent mis-read (no sentinel tag).

**Benchmark measurement trap.** Repeated `_ac_power_flow`/`solve_power_flow!` on the same `data` warm-starts to 0-iteration convergence (lazy early-return). Perturb injections per rep or you measure nothing. Use iteration count (not wall-clock) as the robust metric; the wall-clock timer is noisy. Background heavy compute (10k benchmark, full perf suite) — never block synchronously in a subagent.

Expand All @@ -88,24 +88,26 @@ Exported solver-model types and functions (see `src/PowerFlows.jl`):

## Commands (verified against this clone)

This package uses **ReTest** and a `test/Project.toml` env (deps incl. PowerSystemCaseBuilder, ReTest, Pardiso, Aqua). Read the `sienna-test-environment` skill for the shared rules; PowerFlows specifics:
This package uses **ParallelTestRunner** and a `test/Project.toml` env (deps incl. PowerSystemCaseBuilder, ParallelTestRunner, Pardiso, Aqua) — one worker process per `test_*.jl` file, sharing nothing but `test/includes.jl`'s preamble. Read the `sienna-test-environment` skill for the shared rules; PowerFlows specifics:

```sh
# Compile-check between edits (package env, fast):
julia --project -e 'using PowerFlows'

# One-time per clone: make --project=test resolve PowerFlows to the WORKING TREE
# (else it can resolve the registered copy in ~/.julia/packages and run stale source,
# and new test/test_*.jl files are invisible to the glob in test/PowerFlowsTests.jl):
# and new test/test_*.jl files are invisible to the glob in test/runtests.jl):
julia --project=test -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd()))'
# Verify (must print the working-tree path, not ~/.julia/packages/...):
julia --project=test -e 'import Pkg; println(Base.find_package("PowerFlows"))'

# Run full suite:
julia --project=test test/runtests.jl

# Run a filtered subset via ReTest:
julia --project=test -e 'using PowerFlows; include("test/PowerFlowsTests.jl"); using .PowerFlowsTests, ReTest; retest(PowerFlowsTests, r"<regex>")'
# Run a subset filtered by FILE name (startswith), cap parallelism, or list discoverable tests:
julia --project=test test/runtests.jl test_dc_power_flow
julia --project=test test/runtests.jl --jobs=4
julia --project=test test/runtests.jl --list

# Docs:
julia --project=docs docs/make.jl
Expand All @@ -114,7 +116,7 @@ julia --project=docs docs/make.jl
julia --project=scripts/formatter -e 'include("scripts/formatter/formatter_code.jl")'
```

ReTest runs the whole suite and reports failures at the end (does not abort on first failure). Note `runtests.jl` aborts the whole run at the first exception outside a `@test`; "suite green" means the run REACHED the final `Main.PowerFlowsTests | <N>` summary with no Error column. Under recent PSY/IS, `PSY.System("file.raw")` may not parse PSS/E raw — use the PowerSystemCaseBuilder `PowerFlowFileParser` path for raw inputs in tests.
Each test file runs as its own testset in its own worker process; the runner reports pass/fail per file and does not abort the whole run on one file's failure. Under recent PSY/IS, `PSY.System("file.raw")` may not parse PSS/E raw — use the PowerSystemCaseBuilder `PowerFlowFileParser` path for raw inputs in tests.

## Auto-generated files / do-not-edit

Expand Down
18 changes: 9 additions & 9 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -37,18 +37,18 @@ PowerFlowsPardisoExt = "Pardiso"
[sources]
# Branch pins until the OpenAPI packages are registered and the psy6-line branches merge.
InfrastructureSystems = {rev = "IS4", url = "https://github.com/Sienna-Platform/InfrastructureSystems.jl.git"}
PowerSystems = {url = "https://github.com/Sienna-Platform/PowerSystems.jl.git", rev = "psy6"}
PowerSystems = {url = "https://github.com/Sienna-Platform/PowerSystems.jl.git", rev = "mb/remote-control"}
# Pinned to the branch rather than the registry: registry PNM v0.24 requires PowerSystems
# 5.11, which this line (5.10.0) does not satisfy.
PowerNetworkMatrices = {url = "https://github.com/Sienna-Platform/PowerNetworkMatrices.jl.git", rev = "psy6"}
PowerNetworkMatrices = {url = "https://github.com/Sienna-Platform/PowerNetworkMatrices.jl.git", rev = "jd/psy6_correctness"}

InfrastructureCoreOpenAPIModels = {url = "https://github.com/Sienna-Platform/PowerOpenAPIModels.git", rev = "main", subdir = "InfrastructureCoreOpenAPIModels.jl"}
InfrastructureTimeSeriesOpenAPIModels = {url = "https://github.com/Sienna-Platform/PowerOpenAPIModels.git", rev = "main", subdir = "InfrastructureTimeSeriesOpenAPIModels.jl"}
PowerCoreOpenAPIModels = {url = "https://github.com/Sienna-Platform/PowerOpenAPIModels.git", rev = "main", subdir = "PowerCoreOpenAPIModels.jl"}
PowerDynamicsOpenAPIModels = {url = "https://github.com/Sienna-Platform/PowerOpenAPIModels.git", rev = "main", subdir = "PowerDynamicsOpenAPIModels.jl"}
PowerInvestmentsOpenAPIModels = {url = "https://github.com/Sienna-Platform/PowerOpenAPIModels.git", rev = "main", subdir = "PowerInvestmentsOpenAPIModels.jl"}
PowerOpenAPIModels = {url = "https://github.com/Sienna-Platform/PowerOpenAPIModels.git", rev = "main", subdir = "PowerOpenAPIModels.jl"}
PowerOperationsOpenAPIModels = {url = "https://github.com/Sienna-Platform/PowerOpenAPIModels.git", rev = "main", subdir = "PowerOperationsOpenAPIModels.jl"}
InfrastructureCoreOpenAPIModels = {url = "https://github.com/Sienna-Platform/PowerOpenAPIModels.git", rev = "mb/remote-control", subdir = "InfrastructureCoreOpenAPIModels.jl"}
InfrastructureTimeSeriesOpenAPIModels = {url = "https://github.com/Sienna-Platform/PowerOpenAPIModels.git", rev = "mb/remote-control", subdir = "InfrastructureTimeSeriesOpenAPIModels.jl"}
PowerCoreOpenAPIModels = {url = "https://github.com/Sienna-Platform/PowerOpenAPIModels.git", rev = "mb/remote-control", subdir = "PowerCoreOpenAPIModels.jl"}
PowerDynamicsOpenAPIModels = {url = "https://github.com/Sienna-Platform/PowerOpenAPIModels.git", rev = "mb/remote-control", subdir = "PowerDynamicsOpenAPIModels.jl"}
PowerInvestmentsOpenAPIModels = {url = "https://github.com/Sienna-Platform/PowerOpenAPIModels.git", rev = "mb/remote-control", subdir = "PowerInvestmentsOpenAPIModels.jl"}
PowerOpenAPIModels = {url = "https://github.com/Sienna-Platform/PowerOpenAPIModels.git", rev = "mb/remote-control", subdir = "PowerOpenAPIModels.jl"}
PowerOperationsOpenAPIModels = {url = "https://github.com/Sienna-Platform/PowerOpenAPIModels.git", rev = "mb/remote-control", subdir = "PowerOperationsOpenAPIModels.jl"}

[compat]
InfrastructureCoreOpenAPIModels = "0.1"
Expand Down
8 changes: 4 additions & 4 deletions docs/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -25,13 +25,13 @@ PowerFlows = {path = ".."}
InfrastructureSystems = {rev = "IS4", url = "https://github.com/Sienna-Platform/InfrastructureSystems.jl.git"}
# Pinned to the branch rather than the registry: registry PNM v0.24 requires PowerSystems
# 5.11, which this line (5.10.0) does not satisfy.
PowerNetworkMatrices = {url = "https://github.com/Sienna-Platform/PowerNetworkMatrices.jl.git", rev = "psy6"}
PowerSystemCaseBuilder = {url = "https://github.com/Sienna-Platform/PowerSystemCaseBuilder.jl.git", rev = "psy6"}
PowerSystems = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerSystems.jl.git"}
PowerNetworkMatrices = {url = "https://github.com/Sienna-Platform/PowerNetworkMatrices.jl.git", rev = "jd/psy6_correctness"}
PowerSystemCaseBuilder = {url = "https://github.com/Sienna-Platform/PowerSystemCaseBuilder.jl.git", rev = "mb/remote-control"}
PowerSystems = {url = "https://github.com/Sienna-Platform/PowerSystems.jl.git", rev = "mb/remote-control"}
# PSCB's own [sources] pin for this is ignored once PSCB is a dependency rather than the root
# project; the docs env needs its own pin so Pkg can resolve PSCB's unregistered parser dep.
PowerTableDataParser = {url = "https://github.com/NLR-Sienna/PowerTableDataParser.jl.git", rev = "psy6"}
PowerFlowFileParser = {url = "https://github.com/Sienna-Platform/PowerFlowFileParser.jl.git", rev = "psy6"}
PowerFlowFileParser = {url = "https://github.com/Sienna-Platform/PowerFlowFileParser.jl.git", rev = "mb/remote-control"}

[compat]
Documenter = "^1.0"
2 changes: 2 additions & 0 deletions docs/src/tutorials/discrete_control_14bus.jl
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,7 @@ facts = FACTSControlDevice(;
voltage_setpoint = 1.0,
max_shunt_current = 100.0,
shunt_control_type = FACTSShuntControlType.STATCOM,
input_basis = CU,
)
add_component!(sys, facts)

Expand Down Expand Up @@ -176,6 +177,7 @@ facts_tight = FACTSControlDevice(;
voltage_setpoint = 1.0,
max_shunt_current = 5.0,
shunt_control_type = FACTSShuntControlType.STATCOM,
input_basis = CU,
)
add_component!(sys3, facts_tight)

Expand Down
22 changes: 14 additions & 8 deletions scripts/benchmarks/discrete_control_scaling.jl
Original file line number Diff line number Diff line change
Expand Up @@ -49,25 +49,30 @@ function _add_feeder!(sys::PSY.System, ref::PSY.ACBus, k::Int)
PSY.add_component!(sys,
PSY.PowerLoad(; name = "load_$k", available = true, bus = b_load,
active_power = 0.5, reactive_power = 0.25, base_power = 100.0,
max_active_power = 100.0, max_reactive_power = 100.0))
max_active_power = 100.0, max_reactive_power = 100.0,
input_basis = PSY.CU))
PSY.add_component!(sys,
PSY.PowerLoad(; name = "shload_$k", available = true, bus = b_sh,
active_power = 0.05, reactive_power = 0.025, base_power = 100.0,
max_active_power = 100.0, max_reactive_power = 100.0))
max_active_power = 100.0, max_reactive_power = 100.0,
input_basis = PSY.CU))
PSY.add_component!(sys,
PSY.Line(; name = "line_$k", available = true, active_power_flow = 0.0,
reactive_power_flow = 0.0, arc = PSY.Arc(; from = ref, to = b_sh),
r = 1e-2, x = 1e-2, b = (from = 0.0, to = 0.0),
rating = 10.0, angle_limits = (min = -pi / 2, max = pi / 2)))
rating = 10.0, angle_limits = (min = -pi / 2, max = pi / 2),
input_basis = PSY.CU))
PSY.add_component!(sys,
PSY.TwoWindingTransformer(; name = "tap_$k",
circuit = PSY.TransformerCircuit(; available = true,
arc = PSY.Arc(; from = ref, to = b_load),
r = 0.01, x = 0.10, tap = 1.0, rating = 1.0, base_power = 100.0,
control_objective = PSY.TransformerControlObjective.VOLTAGE)))
control_objective = PSY.TransformerControlObjective.VOLTAGE,
input_basis = PSY.CU),
input_basis = PSY.CU))
PSY.add_component!(sys,
PSY.SwitchedAdmittance(; name = "shunt_$k", available = true, bus = b_sh,
Y = 0.0 + 0.0im, initial_status = [0], number_of_steps = [4],
number_engaged = [0], number_of_steps = [4],
Y_increase = [0.0 + 0.05im], admittance_limits = (min = 0.9, max = 1.1)))
return nothing
end
Expand All @@ -80,7 +85,8 @@ function build_controlled_system(K::Int)
PSY.add_component!(sys, ref)
PSY.add_component!(sys,
PSY.Source(; name = "source", available = true, bus = ref,
active_power = 0.0, reactive_power = 0.0, R_th = 0.0, X_th = 1e-5))
active_power = 0.0, reactive_power = 0.0, R_th = 0.0, X_th = 1e-5,
input_basis = PSY.CU))
for k in 1:K
_add_feeder!(sys, ref, k)
end
Expand All @@ -91,8 +97,8 @@ end
# 0-iteration warm-start early return).
function _perturb_loads!(sys::PSY.System, rng)
for ld in PSY.get_components(PSY.PowerLoad, sys)
base = PSY.get_active_power(ld)
PSY.set_active_power!(ld, base * (1.0 + 0.1 * (rand(rng) - 0.5)))
base = PSY.get_active_power(ld, PSY.SU)
PSY.set_active_power!(ld, base * (1.0 + 0.1 * (rand(rng) - 0.5)) * PSY.SU)
end
return nothing
end
Expand Down
2 changes: 1 addition & 1 deletion scripts/benchmarks/formulation_solver_comparison.jl
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@ function run_system(group, name, build_kwargs, extra_settings)
merge(Dict{Symbol, Any}(:validate_voltage_magnitudes => false),
extra_settings)
end
pf = F{S}(; correct_bustypes = true, solver_settings = settings)
pf = F{S}(; correct_bustypes = true, solution_parameters = SolutionParameters(; settings...))
bench(pf, sys, "$fname / $sname")
end
return
Expand Down
10 changes: 5 additions & 5 deletions scripts/benchmarks/method_comparison.jl
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ const PSY = PowerSystems
const PF = PowerFlows

# ─────────────────────────────────────────────────────────────────────────────
# Custom logger that intercepts the @info messages emitted by the solver
# Custom logger that intercepts the @debug convergence messages emitted by the solver
# and extracts iteration count and residual norms.
# ─────────────────────────────────────────────────────────────────────────────
mutable struct MetricsCapture
Expand Down Expand Up @@ -183,7 +183,7 @@ function run_trial(sys, solver_type, solver_settings, x_solved, n, bus_types, K;
# Build PF object
pf = ACPowerFlow{solver_type}(;
correct_bustypes = true,
solver_settings = solver_settings,
solution_parameters = PF.SolutionParameters(; solver_settings...),
)
data = quietly(() -> PF.PowerFlowData(pf, sys))

Expand All @@ -209,10 +209,10 @@ function run_trial(sys, solver_type, solver_settings, x_solved, n, bus_types, K;
# consistently across all methods (including homotopy).
residual_obj = PF.ACPowerFlowResidual(data, 1)
if haskey(kwargs, :x0)
residual_obj(kwargs[:x0], 1)
residual_obj(data, kwargs[:x0], 1)
else
x0_default = PF.calculate_x0(data, 1)
residual_obj(x0_default, 1)
residual_obj(data, x0_default, 1)
end
init_res_L2 = norm(residual_obj.Rv, 2)
init_res_Linf = norm(residual_obj.Rv, Inf)
Expand All @@ -237,7 +237,7 @@ function run_trial(sys, solver_type, solver_settings, x_solved, n, bus_types, K;

# Compute final power flow residual directly — don't rely on log capture,
# which may be missing (homotopy) or absent on non-convergence.
residual_obj(x_final, 1)
residual_obj(data, x_final, 1)
final_res_L2 = norm(residual_obj.Rv, 2)
final_res_Linf = norm(residual_obj.Rv, Inf)

Expand Down
Loading