diff --git a/docs/Project.toml b/docs/Project.toml index 1b5ffbae1..faf126f3d 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -9,6 +9,9 @@ PowerSystemCaseBuilder = "f00506e0-b84f-492a-93c2-c0a9afc4364e" PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd" PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" +[sources] +PowerNetworkMatrices = {path = ".."} + [compat] Documenter = "1" julia = "^1.10" diff --git a/docs/make.jl b/docs/make.jl index e36e48de9..b398f49db 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -4,40 +4,42 @@ using Literate using DocumenterInterLinks links = InterLinks( + "Julia" => "https://docs.julialang.org/en/v1/objects.inv", "PowerSystems" => "https://sienna-platform.github.io/PowerSystems.jl/stable/", "PowerSystemCaseBuilder" => "https://sienna-platform.github.io/PowerSystemCaseBuilder.jl/stable/", ) include(joinpath(@__DIR__, "make_tutorials.jl")) -make_tutorials() +make_literate_folder("tutorials") +make_literate_folder("how_to_guides") pages = OrderedDict( "Welcome Page" => "index.md", "Tutorials" => Any[ - "Getting Started" => "tutorials/getting_started.md", - "Incidence, BA and ABA matrices" => "tutorials/tutorial_Incidence_BA_ABA_matrices.md", - "PTDF matrix" => "tutorials/tutorial_PTDF_matrix.md", - "VirtualPTDF matrix" => "tutorials/tutorial_VirtualPTDF_matrix.md", - "LODF matrix" => "tutorials/tutorial_LODF_matrix.md", - "VirtualLODF matrix" => "tutorials/tutorial_VirtualLODF_matrix.md", - "Industry DFAX values" => "tutorials/tutorial_DFAX.md", - "Radial Reduction" => "tutorials/tutorial_RadialReduction.md", - "Degree Two Reduction" => "tutorials/tutorial_DegreeTwoReduction.md", + "Introduction" => "tutorials/generated_introduction.md", + "Analysis at Scale" => "tutorials/generated_analysis_at_scale.md", ], "How-To Guides" => Any[ - "Compute Network Matrices" => "how_to_guides/compute_network_matrices.md", - "Choose a Linear Solver" => "how_to_guides/choose_linear_solver.md", + "Build Multiple Matrices" => "how_to_guides/generated_build_multiple_matrices.md", + "Choose a Linear Solver" => "how_to_guides/generated_choose_linear_solver.md", + "Reproduce Industry DFAX Values" => "how_to_guides/generated_reproduce_dfax_values.md", + "Define and Apply Contingencies" => "how_to_guides/generated_contingencies.md", + "Diagnose Network Connectivity" => "how_to_guides/generated_diagnose_connectivity.md", + ], + "Reference" => Any[ + "Matrix Overview and Indexing" => "reference/network_matrices_overview.md", + "Matrix Types" => "reference/matrix_types.md", + "Public API" => "reference/public.md", + "Internals" => "reference/internals.md", ], "Explanation" => Any[ - "Computational Considertaions" => "explanation/computational_considerations.md", "DC Power Flow Approximation" => "explanation/dc_power_flow_approximation.md", + "Computational Considerations" => "explanation/computational_considerations.md", "Network Reduction Theory" => "explanation/network_reduction_theory.md", "Flowgate Methodology" => "explanation/flowgate_methodology.md", - ], - "Reference" => Any[ - "Matrix Overview" => "reference/network_matrices_overview.md", - "Public API" => "reference/public.md", - "Internals" => "reference/internals.md", + "Concurrency and the KLU Lock" => "explanation/concurrency.md", + "Equivalent Representation of Reduced Branches" => "explanation/equivalent_branches.md", + "Slack Distribution and Reference-Bus Conventions" => "explanation/slack_conventions.md", ], ) diff --git a/docs/make_tutorials.jl b/docs/make_tutorials.jl index a00daf8a1..cadcba0cb 100644 --- a/docs/make_tutorials.jl +++ b/docs/make_tutorials.jl @@ -216,9 +216,9 @@ end # Download links: # - **Deployed / CI**: absolute URLs under `_DOCS_BASE_URL` when `_downloads_use_absolute_urls()` is true. # - **Local**: bare filenames (siblings of `generated_*.md` in `docs/src/tutorials/`). -function add_download_links(content, jl_file, ipynb_file) +function add_download_links(content, jl_file, ipynb_file, subdir = "tutorials") script_link, notebook_link = if _downloads_use_absolute_urls() - ("$_DOCS_BASE_URL/tutorials/$(jl_file)", "$_DOCS_BASE_URL/tutorials/$(ipynb_file)") + ("$_DOCS_BASE_URL/$subdir/$(jl_file)", "$_DOCS_BASE_URL/$subdir/$(ipynb_file)") else (jl_file, ipynb_file) end @@ -345,8 +345,8 @@ end # - If a markdown cell contains one or more image fragments, append exactly one # "view online" fallback note at the end of that cell. # - If the note already exists in the cell, no change is applied. -function add_image_links(nb::Dict, outputfile_base::AbstractString) - tutorial_url = "$_DOCS_BASE_URL/tutorials/$(outputfile_base)/" +function add_image_links(nb::Dict, outputfile_base::AbstractString, subdir = "tutorials") + tutorial_url = "$_DOCS_BASE_URL/$subdir/$(outputfile_base)/" msg = "_If image is not available when viewing in a Jupyter notebook, view the tutorial online [here]($tutorial_url)._" cells = get(nb, "cells", []) for (idx, cell) in enumerate(cells) @@ -409,29 +409,30 @@ end # Process tutorials with Literate ######################################################### -# Generate tutorial markdown + notebook artifacts from literate .jl sources. +# Generate markdown + notebook artifacts from the literate .jl sources in +# docs/src/ (e.g. "tutorials", "how_to_guides"). # # Pipeline: -# 1) discover tutorial .jl files (excluding helper files starting with "_") +# 1) discover .jl files in the folder (excluding helper files starting with "_") # 2) generate Documenter-flavored markdown with injected download links # 3) generate notebook with admonition conversion, setup preface, and image note -function make_tutorials() - tutorials_dir = abspath(joinpath(@__DIR__, "src", "tutorials")) +function make_literate_folder(subdir::String = "tutorials") + src_dir = abspath(joinpath(@__DIR__, "src", subdir)) # Exclude helper scripts that start with "_" - if isdir(tutorials_dir) - tutorial_files = + if isdir(src_dir) + source_files = filter( x -> endswith(x, ".jl") && !startswith(x, "_"), - readdir(tutorials_dir), + readdir(src_dir), ) - if !isempty(tutorial_files) - # Clean up old generated tutorial files - tutorial_outputdir = tutorials_dir - clean_old_generated_files(tutorial_outputdir) + if !isempty(source_files) + # Clean up old generated files + outputdir = src_dir + clean_old_generated_files(outputdir) - for file in tutorial_files + for file in source_files @show file - infile_path = joinpath(tutorials_dir, file) + infile_path = joinpath(src_dir, file) execute = if occursin("EXECUTE = TRUE", uppercase(readline(infile_path))) true @@ -443,7 +444,7 @@ function make_tutorials() # Generate markdown Literate.markdown(infile_path, - tutorial_outputdir; + outputdir; name = outputfile, credit = false, flavor = Literate.DocumenterFlavor(), @@ -453,6 +454,7 @@ function make_tutorials() insert_md(content), file, string(outputfile, ".ipynb"), + subdir, ) ), execute = execute) @@ -461,13 +463,17 @@ function make_tutorials() # preprocess_admonitions_for_notebook converts Documenter admonitions to blockquotes # so they render in Jupyter; markdown output keeps !!! style for Documenter. Literate.notebook(infile_path, - tutorial_outputdir; + outputdir; name = outputfile, credit = false, execute = false, preprocess = preprocess_admonitions_for_notebook, postprocess = nb -> - add_image_links(add_pkg_status_to_notebook(nb), outputfile)) + add_image_links( + add_pkg_status_to_notebook(nb), + outputfile, + subdir, + )) end end end diff --git a/docs/src/explanation/computational_considerations.md b/docs/src/explanation/computational_considerations.md index 347a8ee2c..8d223beaf 100644 --- a/docs/src/explanation/computational_considerations.md +++ b/docs/src/explanation/computational_considerations.md @@ -1,70 +1,60 @@ -## Computational Considerations - -### Matrix Construction - -All matrices in `PowerNetworkMatrices.jl` are derived from the `Ybus` matrix (i.e. building any matrix starts with building the `Ybus`). Additionally, all network reductions are applied to the `Ybus` matrix prior to computing the downstream matrices. This design choice is key for enabling high performance and code maintainability: Looping through the system objects is required only when building the `Ybus` (slow) and subsequent operations are built on fast matrix operations on (often sparse) matrices. In addition, network reductions are only defined for the `Ybus` but can be applied uniformly across all matrices. - -### Sparsity - -Power networks are sparse; most buses connect to only a few others. This sparsity is exploited for computational efficiency via sparse linear solvers: - - - Incidence and admittance matrices are very sparse. - - Common sensitivity matrices (e.g. PTDF and LODF) are dense. - -### Automatic Sparsification Tolerance - -The PTDF and LODF are obtained by solving against the reduced susceptance matrix `ABA = Aᵀ B A`, which is sparse. In theory, however, **the inverse of a sparse matrix is dense**: `ABA` is a grounded graph Laplacian, and `ABA⁻¹` has essentially no zeros even though `ABA` does. So the sensitivity matrices come out dense, and for a large network a single PTDF column has one entry per bus — tens of thousands of numbers, almost all of them negligible (a branch is effectively insensitive to an injection electrically far away). - -To recover sparsity we apply a **tolerance** and drop entries below it (`droptol!`). The size of that cutoff is the question this package answers automatically. Every PTDF/LODF constructor takes a `tol` keyword: - -```julia -tol::Union{Float64, AutoTolerance} = AutoTolerance() -``` - - - A **`Float64`** is an explicit, *absolute* cutoff: any entry with `|x| ≤ tol` is dropped. Use it to pin an exact result (`tol = eps()`) or to sparsify by a fixed number. This is the backward-compatible path. - - - An **[`AutoTolerance`](@ref)** (the default) chooses the cutoff from the data instead of a hand-tuned guess. It drops an entry of a row only when it is below the precision the input data can justify, *relative to that row's own peak*: - - ``` - drop entry j of row i when |row_i[j]| < α · max|row_i|, - α = clamp(safety · δ, 1e-6, 1e-2) - ``` - - where `δ` is the relative precision of the branch reactances (auto-discovered from their significant figures, or set explicitly via `data_precision`). Because the cutoff is *relative to each row's peak*, the achieved column sparsity is independent of the matrix scale and of how ill-conditioned `ABA` is. The 1-norm condition number of `ABA` is still estimated and logged as a diagnostic, but it never enters the cutoff. - -Sparsification only matters at scale, so `AutoTolerance` acts only where it pays off: - - - **On-demand (virtual) matrices at or above `AUTO_TOLERANCE_BUS_LIMIT` buses** are sparsified per requested row/column — this is what lets a column of a large system come back sparse instead of dense. - - **Small systems and the dense `PTDF`/`LODF` constructors are returned exactly** (`AutoTolerance` is a no-op there), preserving their dense type and numerical values. Pass a `Float64` `tol` to sparsify those explicitly. - -For very large studies, prefer the [`VirtualPTDF`](@ref)/[`VirtualLODF`](@ref) variants: they compute rows on demand and, with the default `AutoTolerance`, store each one sparsely. - -#### Accuracy and limitations - -Sparsification trades exactness for memory, and the relative per-row rule has consequences worth understanding before you rely on a sparsified matrix for a sensitive calculation. Each dropped entry is below `α · max|row|` (`α ≤ 1e-2`, and typically `α ≈ δ ≈ 5e-4`), so the dominant sensitivities are never touched — but the following hold: - - - **The error is one-signed, not zero-mean.** A dropped entry is set to exactly zero, never rounded, so each row's total mass strictly decreases. When you sum many small entries of a row (for example, aggregating a flow contribution across many buses), the truncation errors accumulate in the same direction instead of cancelling. The bias is bounded by `(number of dropped entries) · α · max|row|`. - - - **The cutoff is per row, so global invariants are not preserved.** Each row is sparsified against *its own* peak, with no coupling between rows or columns, so quantities depending on cross-row or cross-column structure — Kirchhoff's current law, a column sum, or the *difference* of two entries — are not conserved. The sharp case is two buses `j, k` both far from branch `i`: `PTDF[i,j]` and `PTDF[i,k]` may be similar in magnitude yet fall on opposite sides of the cutoff. Each *absolute* error stays below the threshold, but the *relative* error on the (tiny) difference `PTDF[i,j] − PTDF[i,k]` can approach 100%. - - **Auto-discovered precision assumes a power-of-10 base.** With `data_precision = :auto`, `δ` counts the significant figures of the branch reactances. Decimal significant-figure counts are invariant under multiplication by a power of ten (the conventional 100 MVA base) but **not** under an arbitrary impedance base: data converted by a non-power-of-10 base (e.g. `Z_base = kV²/MVA = 190.44 Ω`) reads more figures than it carries, so `:auto` *over-estimates* precision. The direction is safe — a smaller `α`, hence *less* aggressive dropping — but on such data prefer an explicit `data_precision`. - - **Contingency (Woodbury) corrections do not amplify the error.** In [`VirtualMODF`](@ref) the cutoff is applied to the *final* post-contingency row, after the exact Woodbury solve; the correction is computed from exact factorization solves, never sparsified rows. The error stays bounded by the cutoff however near-critical (near-islanding) the contingency is, even when the Woodbury update is severely ill-conditioned. Verified directly in the test suite. - -When you need an exact result — to preserve KCL, to difference two small sensitivities, or to validate against a reference — pass a `Float64` `tol`: use `tol = eps()` for an unsparsified matrix, or a deliberate fixed cutoff for a reproducible absolute tolerance. `AutoTolerance` is the memory-versus-accuracy lever; the explicit `Float64` paths remain available for when exactness matters more than size. - -### Matrix Sizes - -A system with $N_b$ buses and $N_a$ arcs has matrix dimensions: - - - Incidence: $N_a × N_b$ (sparse) - - Admittance: $N_b × N_b$ (sparse) - - PTDF: $N_a × N_b$ (dense) - - LODF: $N_a × N_a$ (dense) - -### Computational Complexity - -| Operation | Complexity | Notes | -|:----------------- |:-------------------- |:------------------------------ | -| Incidence Matrix | O($N_a$) | Simple topology scan | -| Admittance Matrix | O($N_a$) | Includes electrical parameters | -| PTDF | O($N_b^3$) | Requires matrix inversion | -| LODF | O($N_a \cdot N_b^2$) | Derived from PTDF | +# Computational Considerations + +## Sparsity + +Power networks are sparse — most buses connect to only a few others — and this is +exploited via sparse linear solvers, so the connectivity matrices ([`IncidenceMatrix`](@ref) and [`Ybus`](@ref)) +are very sparse. The sensitivity matrices ([`PTDF`](@ref), [`LODF`](@ref)), however, are obviously dense: e.g., one entry per bus in every +[`PTDF`](@ref) column, most of them negligible because a branch is nearly insensitive +to an injection electrically far away. + +## Sparsification and tolerance + +Those negligible entries can be dropped to recover sparsity. The **tolerance** is the +cutoff below which an entry is set to zero. The default [`AutoTolerance`](@ref) picks +it from the data as a *relative per-row* drop — an entry is dropped when +`|x| < α · max|row|` — which keeps large matrices sparse while leaving small systems +and the dense constructors exact. The exact rule (including the bus-count gate that +makes it a no-op on small systems and the `Float64` `tol` alternative) is in the +[`AutoTolerance`](@ref) docstring. For very large studies prefer the +[`VirtualPTDF`](@ref)/[`VirtualLODF`](@ref) variants, which compute rows on demand +and store each one sparsely. + +### Accuracy and limitations + +Sparsification trades exactness for memory. Because each dropped entry is below +`α · max|row|` (`α ≤ 1e-2`, typically `α ≈ 5e-4`), the dominant sensitivities are +never touched — but two properties are worth understanding before relying on a +sparsified matrix for a sensitive calculation: + + - **The error is one-signed, not zero-mean.** A dropped entry becomes exactly + zero, so a row's total mass strictly decreases. Summing many small entries of a + row — e.g. aggregating a transfer's flow contribution across a subsystem's buses — + accumulates that truncation in one direction instead of cancelling. The bias is + bounded by `(number of dropped entries) · α · max|row|`. + - **The cutoff is per row, so cross-row/column invariants are not preserved.** Each + row is sparsified against its own peak, so quantities that couple different rows or + columns — Kirchhoff's current law, a column sum, the *difference* of two entries — + are not conserved. Two buses `j, k` both far from branch `i` can have similar + `PTDF[i,j]`, `PTDF[i,k]` land on opposite sides of the cutoff: each absolute error + stays under threshold, yet the relative error on their tiny difference can approach + 100%. + +Contingency corrections do not compound this: in [`VirtualMODF`](@ref) the cutoff is +applied to the *final* post-contingency row, after the exact Woodbury solve, so the +error stays bounded by the cutoff however near-critical the contingency. + +When you need an exact result — to preserve KCL, to difference two small +sensitivities, or to validate against a reference — pass a `tol::Float64` +(`tol = eps()` for an unsparsified matrix, or a deliberate fixed cutoff). + +## Matrix sizes and complexity + +A system with $N_b$ buses and $N_a$ arcs: + +| Operation | Dimensions | Complexity | Notes | +|:----------------- |:-------------------- |:-------------------- |:------------------------------ | +| Incidence Matrix | $N_a × N_b$ (sparse) | $O(N_a)$ | Simple topology scan | +| Admittance Matrix | $N_b × N_b$ (sparse) | $O(N_a)$ | Includes electrical parameters | +| PTDF | $N_a × N_b$ (dense) | $O(N_b^3)$ | Requires matrix inversion | +| LODF | $N_a × N_a$ (dense) | $O(N_a \cdot N_b^2)$ | Derived from PTDF | diff --git a/docs/src/explanation/concurrency.md b/docs/src/explanation/concurrency.md new file mode 100644 index 000000000..8c86dc44f --- /dev/null +++ b/docs/src/explanation/concurrency.md @@ -0,0 +1,85 @@ +# Concurrency and the KLU lock + +The virtual matrices are designed to be **thread-safe to read**, but reading them +concurrently does **not** make them faster. This page explains why the solver +work serializes, what that means for multithreaded code, and why the per-arc +solve loop cannot be sped up by threading. + +## Concurrent `getindex` is safe but serialized + +Indexing a [`VirtualPTDF`](@ref), [`VirtualLODF`](@ref), or [`VirtualMODF`](@ref) +from several threads produces correct results with no data races. Two locking +layers guarantee it: + + - A **per-cache `ReentrantLock`** on each virtual matrix guards its row cache and + its solver scratch buffers. Cache lookups and inserts, and the factorization + solve itself, happen under this lock. + - A **process-wide `_LIBKLU_LOCK`** (`src/KLUWrapper/KLUWrapper.jl`) wraps *every* + libklu call in the entire process, across all matrices and all threads. + +The consequence is that the expensive part of a row computation — the KLU solve — +runs one at a time. Concurrent readers get correct answers, but they queue: the +throughput of `N` threads all missing the cache is essentially the throughput of +one. **Do not expect a parallel speedup from threading factor solves.** The value +of the locking is safety and correctness under concurrency, not scaling. + +The cache-miss path is written to hold locks as briefly as correctness allows. +The shared `cached_row_lookup` pattern (`src/row_cache.jl`) takes the cache lock +to test for a hit, runs the row computation, then takes the lock again to insert, +double-checking for a row a concurrent producer may have inserted in the +meantime. The compute itself is still serialized — through the per-cache solver +lock and `_LIBKLU_LOCK` — because that is where the libklu work lives. + +## Why libklu must be serialized + +`_LIBKLU_LOCK` is not conservative caution; it reflects a measured property of +the library. A design with a **pool of independent KLU caches** — distinct +`Numeric`, `Symbolic`, and `Common` objects per thread, so that each thread could +in principle solve without touching another's state — was implemented and then +removed. Empirically, per-thread objects did **not** prevent libklu state +corruption: distinct handles still interfered, producing intermittent wrong +results. The global lock is what makes concurrent use correct, so it stays, and +the per-thread pool was dropped because it added complexity with no throughput +benefit once the global lock was required anyway. What remains is one factor and +one cache per virtual matrix. + +The Apple Accelerate backend (`src/AccelerateWrapper/`) has no documented +cross-handle corruption issue analogous to libklu's, so its solves are guarded by +the per-cache `solver_lock` alone, without a process-wide lock. It is still +serialized per matrix, for the same buffer-safety reason. + +## Why the per-arc solve loop is inherently serial + +Building a PTDF or answering a set of contingency queries means one linear solve +per arc (or per contingency) against the factorized `ABA`. It is tempting to +parallelize that loop. It does not work, for two independent reasons: + + 1. **KLU cannot do concurrent solves**, even given per-thread workspaces. The + library serializes internally (this is the same property that forced + `_LIBKLU_LOCK`), so handing each thread its own scratch buffers does not let + the solves overlap — they still queue in the library. + 2. **The query pattern is incremental, not batched.** Sienna consumers ask for + rows one arc or one contingency at a time, as a study progresses. There is no + point at which a large batch of right-hand sides is available to solve + together, so a multi-RHS reformulation — the usual way to extract parallelism + from a factorized system — does not fit the access pattern. + +Because both the library and the workload resist it, threading the solve loop is +not a lever this package offers. The realized performance work went into making +each individual build cheaper (for example, the Ybus adjacency assembly), not +into running solves in parallel. + +## Practical guidance + + - **Reading virtual matrices from multiple threads is safe.** Correctness is + guaranteed; use it when a thread happens to need a row. + - **Do not thread hoping for solver speedup.** The solves serialize on + `_LIBKLU_LOCK` (KLU) or the per-cache lock (Accelerate); more threads means more + queueing, not more throughput. + - **Independent, non-solver work can still overlap.** The locks cover libklu + calls and each cache's buffers, not your surrounding logic. Parallelism is + worth pursuing above the solve, in how you organize a study, rather than inside + it. + +For the solver backends themselves, see the +[choose-a-linear-solver how-to](../how_to_guides/generated_choose_linear_solver.md). diff --git a/docs/src/explanation/dc_power_flow_approximation.md b/docs/src/explanation/dc_power_flow_approximation.md index 95097f330..d691489d4 100644 --- a/docs/src/explanation/dc_power_flow_approximation.md +++ b/docs/src/explanation/dc_power_flow_approximation.md @@ -1,24 +1,43 @@ -## The DC Power Flow Approximation +# The DC Power Flow Approximation -Many network matrices in PowerNetworkMatrices.jl rely on the DC power flow approximation. +Many network matrices in PowerNetworkMatrices.jl rely on the DC power flow +approximation. It linearizes the AC power flow equations into a purely +algebraic relationship between active-power injections and voltage angles, +which is what makes the sensitivity matrices ([`PTDF`](@ref), [`LODF`](@ref)) fast to build and +cheap to reason about. -### Assumptions: +### Assumptions - 1. **Voltage Magnitude**: All bus voltages are approximately 1.0 per unit - 2. **Small Angles**: Voltage angle differences are small (< 15°) - 3. **Resistance**: Line resistance is negligible compared to reactance - 4. **Active Power**: Only active power flows are considered + 1. **Voltage magnitude**: all bus voltages are approximately 1.0 per unit + 2. **Small angles**: voltage angle differences are small (< 15°), so + ``\sin(\theta_i - \theta_j) \approx \theta_i - \theta_j`` + 3. **Resistance**: line resistance is negligible compared to reactance + 4. **Active power**: only active power flows are considered -### When DC Approximation Works Well: +Under these assumptions the branch flow becomes a linear function of the angle +difference, ``P_{ij} \approx (\theta_i - \theta_j)/X_{ij}``, and the whole +network collapses to the susceptance-weighted graph Laplacian that the +[`BA_Matrix`](@ref) and [`ABA_Matrix`](@ref) encode. + +### When the DC approximation works well - Transmission systems (high voltage) - Normal operating conditions - Security and market analysis - Planning studies -### When to Be Cautious: +### When to be cautious - Distribution systems (high R/X ratios) - Large angle differences - Voltage-constrained systems - Detailed reactive power analysis + +## References + +For the DC power-flow approximation and its assumptions see B. Stott, +J. Jardim, and O. Alsaç, "DC Power Flow Revisited," *IEEE Transactions on +Power Systems*, vol. 24, no. 3, pp. 1290–1300, 2009; and A. J. Wood, +B. F. Wollenberg, and G. B. Sheblé, *Power Generation, Operation, and +Control*, 3rd ed., Wiley, 2013. See also +[https://en.wikipedia.org/wiki/Power-flow_study](https://en.wikipedia.org/wiki/Power-flow_study). diff --git a/docs/src/explanation/equivalent_branches.md b/docs/src/explanation/equivalent_branches.md new file mode 100644 index 000000000..5676de5c4 --- /dev/null +++ b/docs/src/explanation/equivalent_branches.md @@ -0,0 +1,117 @@ +# Equivalent representation of reduced branches + +When a network reduction collapses a group of branches into one, the reduced +network must carry a single **equivalent branch** in place of the group. This +page explains how the equivalent's electrical parameters and ratings are formed, +and — the subtle part — **why the rating policies differ** while the impedance +aggregation does not. + +The aggregated-branch types themselves (`BranchesParallel`, +`MixedBranchesParallel`, `BranchesSeries`, the internal +`ThreeWindingTransformerWinding`, and the resulting `EquivalentBranch`) are +documented by their docstrings in the +[full public API](../reference/public.md) under the internal (non-exported) symbols. + +## Two kinds of group + +Reductions produce two topological groupings: + + - **Parallel groups** (`AbstractBranchesParallel`): several branches sharing the + same pair of buses. Physically they are alternative paths between the same two + nodes. + - **Series chains** (`BranchesSeries`): a run of branches through intermediate + degree-two buses that carry no injection, so the chain behaves as a single + branch between its endpoints. + +A three-winding transformer contributes its own wye-to-star structure, handled +by the internal `ThreeWindingTransformerWinding`. + +## Impedance aggregation is physics — one answer + +The equivalent electrical parameters are not chosen; they are *derived* from the +requirement that the reduced branch present the same terminal behavior as the +group it replaces. PNM does this by building the group's equivalent admittance +(`populate_equivalent_ybus!`) and reading the physical parameters back off it +(`get_equivalent_physical_branch_parameters`, `src/common.jl`), yielding an +`EquivalentBranch` with series `r`/`x`, shunt `g`/`b` at each end, tap, +and phase shift. + +The combining rules follow directly: + + - **Parallel:** admittances add. The group's series susceptance is the sum of the + members', `b = Σ bᵢ` — more parallel paths means a stronger (lower-impedance) + connection. + - **Series:** impedances add, so susceptances combine reciprocally, + `b = 1 / Σ (1/bᵢ)` — a chain is weaker than its strongest link. + +There is a single correct answer here because the electrical behavior is fixed by +Kirchhoff's laws.[^circuits] The equivalent is exact for the linear (DC) model. + +## Rating aggregation is policy — several answers + +A **rating** is not an electrical quantity the way impedance is. It is a limit +imposed by the study, and "the limit of the group" is a genuinely ambiguous +question whose answer depends on what the study is protecting against. That is +why PNM exposes several rating strategies for a parallel group rather than one +(`src/BranchesParallel.jl`): + + - **[`get_sum_of_max_rating`](@ref) — nominal capacity.** `Σ Sᵢ`, treating every + circuit as independently loadable to its own thermal limit. This is the least + conservative aggregate; it assumes flow can be steered freely across the group + so that each circuit reaches its limit at once. It answers "how much could this + corridor carry in the best case?" + - **[`get_single_element_contingency_rating`](@ref) — N-1 security.** `Σ Sᵢ − maxᵢ Sᵢ`, the capacity that survives when the largest circuit in the group + trips. It answers "what can I still rely on after losing one element?" For a + group of one it is zero, correctly, because there is nothing left after the sole + circuit trips. + - **[`get_impedance_averaged_rating`](@ref) — realistic DC loading.** The + susceptance-weighted average `Σ fᵢ Sᵢ` with `fᵢ = bᵢ / Σ bₖ`. This reflects how + DC flow *actually* divides across a parallel group: current follows the path of + least impedance, so the low-impedance (high-susceptance) circuit carries the + larger share and reaches its limit first. Because flow cannot in fact be steered + arbitrarily, the sum-of-max is optimistic and this weighted figure is closer to + the binding constraint. It requires a finite, non-zero total susceptance and + throws an `ArgumentError` otherwise. + +The three are ordered from least to most physically constrained: +`sum_of_max ≥ impedance_averaged`, and the N-1 figure answers a different +(security) question entirely. Which one is "right" is a modeling decision, not a +computation the package can make for you — hence three named policies instead of +a silent default. + +### Series chains and emergency ratings + +Ratings propagate through the two group types differently, again for physical +reasons: + + - **Series chain rating** is the **weakest link**: `min` over the chain's + members, because a series path can carry no more than its most limited segment + (`get_equivalent_rating` on a `BranchesSeries`). When a member of the chain + is itself a parallel group, it contributes its *single-element-contingency* + rating — the conservative N-1 figure — rather than its optimistic sum. + - **Emergency ratings** follow the same shapes: a parallel group sums the + members' emergency ratings, a series chain takes the minimum (weakest link). + Where a branch has no distinct emergency (`rating_b`) value, its normal rating + is used as the post-contingency limit. + +### Availability + +A group is available only if all of its members are: parallel and series groups +both require every branch present and in service. Losing any one member makes the +equivalent unavailable, which keeps the reduced model consistent with outages of +the underlying branches. + +## Why this separation matters + +The clean split — **impedance is derived, rating is chosen** — is the key idea. +Impedance aggregation has a unique physical answer and PNM computes it once. +Rating aggregation encodes an operator's risk posture, so PNM refuses to pick for +you and instead names the policies (`sum_of_max`, `single_element_contingency`, +`impedance_averaged`) so a study selects the one matching its purpose: raw +capacity, N-1 security, or realistic DC loading. + +## References + +[^circuits]: Series and parallel combination of admittances/impedances is + elementary circuit theory; see + [Series and parallel circuits](https://en.wikipedia.org/wiki/Series_and_parallel_circuits). diff --git a/docs/src/explanation/flowgate_methodology.md b/docs/src/explanation/flowgate_methodology.md index 6b7a344ad..38d8a110c 100644 --- a/docs/src/explanation/flowgate_methodology.md +++ b/docs/src/explanation/flowgate_methodology.md @@ -1,14 +1,16 @@ # Flowgate Methodology -This page explains how the `VirtualMODF` matrix in PowerNetworkMatrices can be -used to evaluate *flowgates* — post-contingency distribution factors of +This page explains how the [`VirtualMODF`](@ref) matrix in PowerNetworkMatrices can +be used to evaluate *flowgates* — post-contingency distribution factors of monitored transmission elements. It describes the mathematics behind the Woodbury-based computation and shows how to query distribution factors using the current API. For a hands-on walkthrough that maps every industry DFAX flavor (GSF, LSF, LODF, OTDF, transfer DFAX, flowgate DFAX, and N-k DFAX) onto the matrices -this page describes, see the [Industry DFAX values](@ref) tutorial. +this page describes, see the +[Reproduce industry DFAX values](../how_to_guides/generated_reproduce_dfax_values.md) +how-to guide. ## Background @@ -18,10 +20,58 @@ source-to-sink transfer that appears as flow on the monitored element after the contingency occurs. In the DC power-flow model, this quantity can be expressed in closed form in -terms of the base-case PTDF and the LODF. `VirtualMODF` generalizes that +terms of the base-case [`PTDF`](@ref) and the [`LODF`](@ref). [`VirtualMODF`](@ref) generalizes that relationship by computing the full post-contingency PTDF row directly, which extends naturally to multi-element contingencies. +For the constructor signatures and the contingency/modification types named +below, see the [Public API Reference](../reference/public.md); for a task-oriented +walkthrough of building and querying a `VirtualMODF` — including the modification +type model — see the +[contingencies how-to](../how_to_guides/generated_contingencies.md). + +## Why a low-rank (Woodbury) update + +A contingency changes the network by removing (or scaling) a handful of +branches. In DC terms, that perturbs the reduced susceptance matrix ``ABA`` in a +way that is **low rank**: outaging ``k`` branches is a rank-``k`` modification, +because each branch contributes a single susceptance term ``b_c\,a_c a_c^\top`` +to ``ABA`` (with ``a_c`` the branch's incidence column). The post-contingency +susceptance matrix is therefore + +```math +\widetilde{ABA} \; = \; ABA \; - \; U\,\Sigma\,U^\top, +``` + +where ``U`` collects the incidence columns of the outaged branches and +``\Sigma`` their susceptance changes — a matrix of rank ``k \ll N``. + +The naïve way to get post-contingency sensitivities would be to rebuild and +re-factorize ``\widetilde{ABA}`` for *every* contingency, an ``O(N^3)`` +factorization each time. The Woodbury (Sherman–Morrison–Woodbury) matrix +identity[^woodbury] avoids this entirely. It expresses the inverse of a low-rank update in +terms of the *already-computed* factorization of the base ``ABA`` plus the +solution of a small ``k \times k`` system: + +```math +\widetilde{ABA}^{-1} \; = \; ABA^{-1} + \; + \; ABA^{-1} U \bigl(\Sigma^{-1} - U^\top ABA^{-1} U\bigr)^{-1} U^\top ABA^{-1}. +``` + +The base ``ABA`` is factorized once at construction. Each contingency then costs +only a few solves against that stored factorization to form the Woodbury factors, +and the expensive ``O(N^3)`` work is never repeated. This is the core reason +post-contingency analysis of many contingencies is tractable at scale: the +dominant cost is paid once and reused, and the per-contingency cost scales with +the (small) number of outaged elements, not the network size. + +The same structure is why there is **no dense `MODF` type** in the package. A +materialized post-contingency factor matrix would be one dense ``N_a \times N_b`` +matrix *per contingency* — the product of two already-large dimensions with a +third — which is prohibitive for any realistic contingency list. `VirtualMODF` +instead keeps only the base factorization and the small Woodbury factors, and +materializes individual post-contingency rows on demand. + ## How `VirtualMODF` computes post-contingency PTDF rows Given a base-case PTDF and a contingency described by a `NetworkModification`, @@ -47,7 +97,7 @@ of two entries of the post-contingency row: ``` For a single-element (N-1) contingency this is equivalent to the explicit -LODF expansion +LODF expansion[^lodf] ```math \mathrm{DF} \; = \; \mathrm{PTDF}[m, s] - \mathrm{PTDF}[m, k] @@ -59,9 +109,9 @@ without additional derivation. ## Describing a contingency -`VirtualMODF` queries are keyed by a `NetworkModification` (or by a -`ContingencySpec` or a `PSY.Outage` that resolves to one). A -`NetworkModification` can be built in several ways: +[`VirtualMODF`](@ref) queries are keyed by a [`NetworkModification`](@ref) (or by a +[`ContingencySpec`](@ref) or a [`Outage`](@extref PowerSystems.Outage) that resolves +to one). A [`NetworkModification`](@ref) can be built in several ways: ```julia using PowerSystems @@ -78,10 +128,12 @@ mod_branch = NetworkModification(vmodf, branch) mod_outage = NetworkModification(vmodf, sys, outage) ``` -When a `VirtualMODF` is constructed from a `PSY.System`, all `PSY.Outage` -supplemental attributes in the system are automatically resolved and -registered. Registered contingencies can be inspected with -`get_registered_contingencies(vmodf)` and queried directly by `PSY.Outage`. +When a [`VirtualMODF`](@ref) is constructed from a +[`System`](@extref PowerSystems.System), all +[`Outage`](@extref PowerSystems.Outage) supplemental attributes in the system are +automatically resolved and registered. Registered contingencies can be inspected +with [`get_registered_contingencies`](@ref) and queried directly by +[`Outage`](@extref PowerSystems.Outage). ## Querying post-contingency rows @@ -117,7 +169,8 @@ df = row[bus_lookup[source_bus]] - row[bus_lookup[sink_bus]] ## Caching and sparsification -`VirtualMODF` maintains two caches, both keyed by `NetworkModification`: +[`VirtualMODF`](@ref) maintains two caches, both keyed by +[`NetworkModification`](@ref): - A Woodbury-factor cache (one entry per contingency, populated on first query for that contingency). @@ -125,16 +178,23 @@ df = row[bus_lookup[source_bus]] - row[bus_lookup[sink_bus]] produced for each monitored arc. The maximum cache size is controlled by the `max_cache_size` keyword (MiB per contingency). -The `tol` keyword of the `VirtualMODF` constructor enables row-level -sparsification: entries whose magnitude is below `tol` are dropped from the -cached row. This reduces memory use and downstream arithmetic cost when many -rows are retained, at the expense of discarding small distribution-factor -contributions. The default `tol = eps()` keeps all entries. - -`clear_caches!(vmodf)` drops the Woodbury and row caches but retains the +The `tol` keyword of the `VirtualMODF` constructor +(`tol::Union{Float64, AutoTolerance}`, default `DEFAULT_AUTO_TOLERANCE`) controls +row-level sparsification: entries whose magnitude falls below the resolved cutoff +are dropped from the cached row. This reduces memory use and downstream +arithmetic cost when many rows are retained, at the expense of discarding small +distribution-factor contributions. Crucially, the cutoff is applied to the +*final* post-contingency row — after the exact Woodbury solve — so sparsification +never enters the correction itself, and the bound holds even when the Woodbury +update is severely ill-conditioned (a near-islanding contingency). Pass an +explicit `Float64` `tol` (e.g. `eps()`) when you need every entry retained; see +[Computational considerations](computational_considerations.md) for the +per-row [`AutoTolerance`](@ref) rule and its accuracy trade-offs. + +[`clear_caches!`](@ref) drops the Woodbury and row caches but retains the contingency registrations, so subsequent queries will simply recompute. Use -`clear_all_caches!(vmodf)` to also drop the registrations (after which the -`VirtualMODF` can no longer be queried). +[`clear_all_caches!`](@ref) to also drop the registrations (after which the +[`VirtualMODF`](@ref) can no longer be queried). ## Relationship to other matrix types @@ -151,8 +211,22 @@ contingency registrations, so subsequent queries will simply recompute. Use - The implementation assumes DC power flow (lossless, linearized). Voltage and stability limits that define some flowgate transfer capabilities must be handled externally. - - MOD-030 flowgate screening (OTDF thresholding, AFC and ATC arithmetic, - interconnection-wide congestion management procedures) is not provided by - this package. `VirtualMODF` computes the distribution factors that such a - layer would consume; the MOD-030 policy vocabulary is not part of the - current API. + - MOD-030 flowgate screening[^mod030] (OTDF thresholding, AFC and ATC + arithmetic, interconnection-wide congestion management procedures) is not + provided by this package. [`VirtualMODF`](@ref) computes the distribution + factors that such a layer would consume; the MOD-030 policy vocabulary is not + part of the current API. + +## References + +[^woodbury]: The Sherman–Morrison–Woodbury identity; see G. H. Golub and + C. F. Van Loan, *Matrix Computations*, 4th ed., Johns Hopkins, 2013, §2.1.4, + or [Woodbury matrix identity](https://en.wikipedia.org/wiki/Woodbury_matrix_identity). +[^mod030]: NERC Reliability Standard MOD-030, *Flowgate Methodology* (Available + Flowgate Capability). [https://www.nerc.com](https://www.nerc.com) +[^lodf]: The line outage distribution factor and its expression in terms of the + base-case PTDF are standard results; see A. J. Wood, B. F. Wollenberg, and + G. B. Sheblé, *Power Generation, Operation, and Control*, 3rd ed., Wiley, 2013, + and J. Guo, Y. Fu, Z. Li, and M. Shahidehpour, "Direct Calculation of Line + Outage Distribution Factors," *IEEE Transactions on Power Systems*, vol. 24, + no. 3, pp. 1633–1634, 2009. diff --git a/docs/src/explanation/network_reduction_theory.md b/docs/src/explanation/network_reduction_theory.md index 379e9e09d..75b9aedd9 100644 --- a/docs/src/explanation/network_reduction_theory.md +++ b/docs/src/explanation/network_reduction_theory.md @@ -1,34 +1,77 @@ # Network Reduction Theory -This document explains the theory and mathematics behind network reduction techniques in power systems. - -## Why Network Reduction? - -Power system networks can be very large, with thousands of buses and branches. Network reduction techniques simplify these networks while preserving essential characteristics for analysis. - -### Benefits of Reduction: - - 1. **Computational Efficiency**: Smaller matrices are faster to compute and invert - 2. **Focus**: Reduces complexity to focus on areas of interest - 3. **Scalability**: Makes large-scale studies tractable - 4. **Clarity**: Simplifies visualization and understanding - -### Preservation Goals: - -Network reduction aims to preserve: - - - Power flow relationships at retained buses - - Impedance relationships between key points - - Stability characteristics (when applicable) - - Essential topology features - -The type of reduction (e.g. `RadialReduction`, `DegreeTwoReduction`, `WardReduction`) determines the extent to which these characteristics are retained - -## Radial Branch Reduction - -### What is a Radial Branch? - -A radial branch is one that connects to a bus that has only one connection to the rest of the network. Think of it as a "dead-end" in the network. +This page explains the theory and mathematics behind network reduction. Power +systems can have thousands of buses and branches; reduction shrinks the network — +making large studies tractable and matrices faster to build and invert — while +preserving the power-flow relationships at the retained buses. Which +characteristics survive depends on the strategy ([`RadialReduction`](@ref), +[`DegreeTwoReduction`](@ref), [`WardReduction`](@ref)). + +## The graph and susceptance structure reduction operates on + +Network reduction is fundamentally a *graph* operation, and to understand why it +is well-defined it helps to see the three matrices that encode the network's +topology and electrical strength. These are the same building blocks the DC +sensitivity matrices are assembled from, and what every reduction manipulates +under the hood. (For constructors and accessors see the +[matrix type reference](../reference/matrix_types.md); the discussion here is about +*what they mean*.) + +### The incidence matrix: pure topology + +The [`IncidenceMatrix`](@ref) ``A`` is the oriented node–arc incidence matrix: one +row per arc, one column per bus, with a ``+1`` at the arc's *from* bus, a ``-1`` at +its *to* bus, and zeros elsewhere. It carries **topology only** — which bus +connects to which, and with what orientation — and nothing electrical. The +reference-bus column is dropped so the downstream susceptance matrix is +non-singular. Reading ``A`` column by column recovers each bus's degree, which is +exactly the quantity radial (degree 1) and degree-two reductions key off of. + +### The BA matrix: topology weighted by electrical strength + +The [`BA_Matrix`](@ref) is the product ``B A``, where ``B`` is the diagonal matrix +of branch susceptances (``b = 1/x`` under the DC approximation). Where ``A`` says +*which* buses a branch connects, ``BA`` scales each connection by *how electrically +strong* it is. Mapped onto bus angles it returns branch flows — the linear operator +behind ``P_{ij} = (\theta_i - \theta_j)/x_{ij}``. + +### The ABA matrix: the grounded graph Laplacian + +The [`ABA_Matrix`](@ref) is ``A^\top B A``, the reduced nodal susceptance matrix — +a **weighted graph Laplacian**[^laplacian] with the reference bus grounded out. Solving +``ABA\,\theta = P`` *is* the DC power flow, and inverting it produces the dense +[`PTDF`](@ref)/[`LODF`](@ref) sensitivities. Because it is a Laplacian, eliminating +a bus is a Kron elimination on this matrix — which is precisely why degree-two +reduction (below) is exact rather than approximate. + +### Why this makes reduction well-posed + +Every reduction is defined on this susceptance graph and then propagates uniformly +to the downstream matrices (see also +[Computational considerations](computational_considerations.md), on why reductions +are applied to the [`Ybus`](@ref) first). A radial bus is a degree-1 node; a +degree-two bus a degree-2 node; Ward reduction[^ward] Kron-eliminates the external +subgraph. Because all are operations on the incidence/susceptance structure, the +same reduction map applies to [`PTDF`](@ref), [`LODF`](@ref), and their virtual +variants without re-deriving anything per matrix. + +### A subtlety: the susceptance graph is not the admittance graph + +Connectivity and reduction do not always see the same network. +[`find_subnetworks`](@ref) walks the **admittance** graph of the [`Ybus`](@ref), +whereas `ABA` is built from the **susceptance** graph. These differ for branches +with ``r > 0`` and ``x = 0``: such a branch has finite admittance but zero +susceptance (``b = 1/x`` absent), so it appears in the [`Ybus`](@ref) graph but +*not* in [`BA_Matrix`](@ref). A network that is one connected island electrically +can therefore fragment into several components in the susceptance graph, leaving +blocks with no reference bus and a **singular `ABA`**. This is why zero-impedance +handling must resolve both endpoints of such a branch to a common node before +building the susceptance matrix. + +## Radial branch reduction + +A radial branch connects to a bus with only one connection to the rest of the +network — a "dead-end." ``` Main Network --- Bus A --- Bus B (radial) @@ -36,197 +79,103 @@ Main Network --- Bus A --- Bus B (radial) Bus C (radial) ``` -In this example, buses B and C are radial - they each connect to only one other bus. - -### Mathematical Basis - -Radial buses can be reduced because power flow to/from them is uniquely determined by the connection bus. The reduction involves: - - 1. **Identifying radial buses**: Buses with degree = 1 - 2. **Transferring loads**: Moving load/generation to the connection point - 3. **Removing the branch**: Eliminating the radial connection - -### Why Radial Reduction Works - -For a radial bus $k$ connected to bus $j$: - -$$P_k = \frac{V_j V_k}{X_{jk}} \sin(\theta_j - \theta_k)$$ - -Under DC approximation with $V_j \approx V_k \approx 1$: - -$$P_k \approx \frac{\theta_j - \theta_k}{X_{jk}}$$ - -Since $P_k$ is determined by the load at bus $k$, the angle $\theta_k$ is uniquely determined: - -$$\theta_k = \theta_j - P_k X_{jk}$$ +For a radial bus ``k`` connected to bus ``j``, the DC relation +``P_k \approx (\theta_j - \theta_k)/X_{jk}`` gives -The radial bus angle is completely determined by its parent bus, so it can be eliminated by transferring the load. - -### Implications of Radial Reduction - -**Preserved:** - - - Power flow at non-radial buses - - Voltage angles at non-radial buses - - System topology at the core network - -**Changed:** - - - Number of buses and branches (reduced) - - Detailed behavior at radial locations - - Local voltage profiles - -**When to Use:** - - - Radial connections are not of interest for analysis - - Focusing on transmission backbone - - Computational efficiency is important - -**When to Avoid:** - - - Studying distribution feeders (often radial) - - Detailed voltage analysis needed - - Radial buses have critical measurements +```math +\theta_k = \theta_j - P_k X_{jk}. +``` -## Degree Two Reduction (Kron Reduction) +The radial-bus angle is completely determined by its parent, so the bus can be +eliminated by transferring its load to the parent. Power flows and angles at the +core (non-radial) network are preserved exactly; only the eliminated locations lose +their explicit representation. Radial reduction is therefore for studies that +target the transmission backbone rather than distribution feeders. -### What is a Degree Two Bus? +## Degree-two (Kron) reduction -A degree two bus connects exactly two other buses, acting as a "pass-through" point. +A degree-two bus connects exactly two others, acting as a pass-through: ``` Bus A --- Bus B (degree 2) --- Bus C ``` -Bus B has degree 2 - it connects A and C but has no other connections. - -### Mathematical Basis: Kron Reduction - -Kron reduction is a systematic method to eliminate buses from the admittance matrix. Degree two reduction is a simple case of Kron reduction for eliminating degree two nodes. - -Given admittance matrix: - -$$\begin{bmatrix} I_r \\ I_e \end{bmatrix} = \begin{bmatrix} Y_{rr} & Y_{re} \\ Y_{er} & Y_{ee} \end{bmatrix} \begin{bmatrix} V_r \\ V_e \end{bmatrix}$$ - -where subscript $r$ denotes retained buses and $e$ denotes eliminated buses. - -If $I_e = 0$ (no injection at eliminated buses): - -$$Y_{ee} V_e = -Y_{er} V_r$$ -$$V_e = -Y_{ee}^{-1} Y_{er} V_r$$ - -Substituting back: - -$$I_r = (Y_{rr} - Y_{re} Y_{ee}^{-1} Y_{er}) V_r = Y_{reduced} V_r$$ - -The reduced admittance matrix is: - -$$Y_{reduced} = Y_{rr} - Y_{re} Y_{ee}^{-1} Y_{er}$$ - -### For a Single Degree Two Bus - -Consider bus $k$ connecting buses $i$ and $j$: +Kron reduction[^kron] eliminates buses from the admittance matrix. Partitioning into +retained (``r``) and eliminated (``e``) buses with no injection at the eliminated +buses (``I_e = 0``): +```math +Y_{reduced} = Y_{rr} - Y_{re} Y_{ee}^{-1} Y_{er}. ``` -Bus i ----(y_ik)---- Bus k ----(y_kj)---- Bus j -``` - -The equivalent admittance directly between $i$ and $j$ after eliminating $k$ is: - -$$y_{ij}^{new} = y_{ij}^{old} + \frac{y_{ik} \cdot y_{kj}}{y_{ik} + y_{kj} + y_{kk}}$$ - -If there's no shunt at bus $k$ ($y_{kk} = 0$): - -$$y_{ij}^{new} = y_{ij}^{old} + \frac{y_{ik} \cdot y_{kj}}{y_{ik} + y_{kj}}$$ - -This is analogous to combining series impedances in circuit theory. - -### Physical Interpretation - -Eliminating a degree two bus: - - - Combines the series impedances of the two connecting branches - - Creates an equivalent direct connection - - Preserves the overall impedance between endpoints - -### Why Degree Two Reduction Works - -The key insight is that a degree two bus with no injection ($P_k = 0$) serves only to pass power through. Its voltage angle is determined by: - -$$P_k = \frac{\theta_i - \theta_k}{X_{ik}} + \frac{\theta_j - \theta_k}{X_{kj}} = 0$$ - -Solving for $\theta_k$: - -$$\theta_k = \frac{X_{kj} \theta_i + X_{ik} \theta_j}{X_{ik} + X_{kj}}$$ - -The angle is a weighted average of its neighbors, so eliminating it and creating a direct equivalent connection preserves the flow relationship. - -### Implications of Degree Two Reduction - -**Preserved:** - - - Power flows at retained buses - - Voltage angles at retained buses - - Overall impedance relationships - - Equivalence for power flow studies -**Changed:** +For a single degree-two bus ``k`` between ``i`` and ``j`` (no shunt at ``k``) this +collapses to the series-impedance combination - - Number of buses (reduced) - - Detailed behavior at eliminated bus - - Branch topology (new equivalent branches created) - -**When to Use:** - - - Pass-through buses are not of analytical interest - - Focusing on injection/load buses - - Simplifying large transmission systems - - Reducing computational burden - -**When to Avoid:** - - - The degree two bus has measurements or controls - - Detailed branch flows needed at that location - - Bus has significant shunt elements - - Studying protection or relay settings - -## Combining Reductions - -Multiple reductions can be applied sequentially. In these cases, the order of the reduction matters as after each reduction, new candidate buses may appear. In general, it is recommended that radial reduction is applied before degree two reduction as applying radial reduction first may expose new degree two buses. - -## Practical Considerations - -### Load and Generation at Reduced Buses - -Eliminating buses without any connected injection components can be done cleanly. For eliminated buses that have load or generation attached, those devices must be mapped to retained buses. For buses with shunt admittances, this mapping process can affect the equivalent admittance matrix. `RadialReduction` and `DegreeTwoReduction` have options for specifying buses that should not be eliminated (e.g. due to the presence of connected injectors). - -### Validation - -Always verify reductions: - - - Compare power flows before and after - - Check that key quantities are preserved - - Validate on a small test system first - -## Limitations - -### What Reduction Cannot Do: - - - Preserve detailed voltage profiles everywhere - - Capture local dynamics at eliminated buses - - Represent phenomena requiring those buses (e.g., local stability) - - Maintain exact AC power flow at eliminated locations - -### When Full Network is Required: - - - Detailed state estimation - - Protection coordination - - Local voltage studies - - Distributed generation integration at eliminated buses - -## Further Reading +```math +y_{ij}^{new} = y_{ij}^{old} + \frac{y_{ik}\,y_{kj}}{y_{ik} + y_{kj}}. +``` -For practical application, see: +Equivalently, a pass-through bus with ``P_k = 0`` has an angle that is the +reactance-weighted average of its neighbors, +``\theta_k = (X_{kj}\theta_i + X_{ik}\theta_j)/(X_{ik} + X_{kj})``, so replacing it +with a direct equivalent branch preserves the flow relationship. Retained-bus flows +and overall impedances are preserved; the eliminated bus and its explicit branches +are replaced by one equivalent branch. Avoid it where the degree-two bus carries +measurements, controls, or a significant shunt. + +## Ward reduction + +Radial and degree-two reduction eliminate *individual* buses by their local degree. +[`WardReduction`](@ref) instead removes a whole **external subsystem** at once, +keeping only a chosen set of **study buses** and preserving how the external network +responds as seen from them. + +Partition the buses into three groups: the **study** buses to retain, the +**external** buses to eliminate, and the **boundary** buses — study buses that a +retained branch connects directly to the external area. Ward equivalencing is a Kron +(Schur-complement) elimination of the external block of the [`Ybus`](@ref): with the +admittance matrix partitioned into study/boundary (``s``) and external (``e``) parts, + +```math +Y_{eq} = Y_{ss} - Y_{se}\,Y_{ee}^{-1}\,Y_{es}. +``` - - [RadialReduction](@ref) - - [DegreeTwoReduction](@ref) +The correction term ``-Y_{se}Y_{ee}^{-1}Y_{es}`` is fully dense over the boundary +buses; it is realized as a set of **equivalent branches between boundary buses** plus +**equivalent shunt admittances** at them, so the reduced network reproduces the +driving-point and transfer behavior the study area would see from the original +external network. Unlike degree-two reduction — which is exact for the DC flows it +preserves — a Ward equivalent is exact only for the operating state its external +injections encode; it is a boundary-matched approximation for other states, which is +why the study area is chosen to contain everything of interest. + +Ward reduction needs a non-empty boundary: with no branch crossing between study and +external areas there is nothing to match against, and the external buses cannot be +folded onto the study set by impedance criteria (a degenerate case the implementation +flags rather than guessing an equivalent). + +## Combining reductions + +Reductions can be applied in sequence, and order matters — each pass can expose new +candidates for the next. Apply [`RadialReduction`](@ref) before +[`DegreeTwoReduction`](@ref): removing dead-ends often exposes new degree-two buses. + +Eliminating a bus with no injection is clean; a bus carrying load or generation has +those devices mapped to a retained bus, and shunt admittances affect the equivalent +admittance matrix. [`RadialReduction`](@ref) and [`DegreeTwoReduction`](@ref) accept +buses to protect from elimination (e.g. injector hosts). Reduction cannot preserve +detailed voltage profiles, local dynamics, or exact AC behavior at eliminated +locations — keep the full network for state estimation, protection coordination, or +local voltage studies. + +## References + +[^laplacian]: The matrix ``A^\top B A`` is a weighted graph Laplacian; see + [https://en.wikipedia.org/wiki/Laplacian_matrix](https://en.wikipedia.org/wiki/Laplacian_matrix). +[^kron]: For Kron reduction and its graph-theoretic interpretation see + F. Dörfler and F. Bullo, "Kron Reduction of Graphs With Applications to + Electrical Networks," *IEEE Transactions on Circuits and Systems I*, vol. 60, + no. 1, pp. 150–163, 2013. See also + [https://en.wikipedia.org/wiki/Kron_reduction](https://en.wikipedia.org/wiki/Kron_reduction). +[^ward]: J. B. Ward, "Equivalent Circuits for Power-Flow Studies," + *Transactions of the AIEE*, vol. 68, no. 1, pp. 373–382, 1949. diff --git a/docs/src/explanation/slack_conventions.md b/docs/src/explanation/slack_conventions.md new file mode 100644 index 000000000..34bf3f434 --- /dev/null +++ b/docs/src/explanation/slack_conventions.md @@ -0,0 +1,112 @@ +# Slack distribution & reference-bus conventions + +A PTDF answers "if I inject one unit of power at bus *j*, how much flows on branch +*i*?" But power is conserved: an injection somewhere must be balanced by a +withdrawal somewhere else. **Where that balancing withdrawal goes is the slack +convention**, and it changes the PTDF.[^ptdf] This page explains single versus +distributed slack, the role of the reference bus, and how the choice shows up in +the sensitivities. + +For the mechanics of configuring a distributed slack — the `dist_slack` keyword and +its per-matrix input type (`Dict{Int, Float64}` for [`PTDF`](@ref)/[`VirtualPTDF`](@ref) +versus `Vector{Float64}` for [`VirtualLODF`](@ref)/[`VirtualMODF`](@ref)) — see the +[matrix type reference](../reference/matrix_types.md). + +## Why a reference is needed at all + +The DC power flow solves `θ = ABA⁻¹ · p` for bus angles from injections. `ABA` is +a grounded graph Laplacian: without grounding it is singular (angles are only +defined up to a constant, and injections must sum to zero). Grounding it means +designating a **reference bus** whose angle is fixed and whose row/column is +removed from the system. That reference bus is also, by construction, the bus that +absorbs the network's power imbalance — the **single slack**. + +So every PTDF is implicitly *relative to a slack*. `PTDF[i, j]` is the flow on +branch *i* caused by injecting at bus *j* **and withdrawing at the slack**. There +is no such thing as a slack-free PTDF. + +## Single slack (the default) + +By default `dist_slack` is empty, and the matrix uses a **single reference bus** +as the slack. This is the standard convention and is what you want when the model +has one designated slack generator absorbing imbalance. + +The reference bus is excluded from the solve, so its own column in the PTDF is +zero — injecting at the slack and withdrawing at the slack moves nothing. + +**Effect of the reference-bus choice.** Moving the slack to a different bus shifts +each PTDF row by a constant: every sensitivity is measured against a different +balancing point. The choice is therefore not arbitrary for interpreting a single +`PTDF[i, j]` value. It does *not*, however, change physically meaningful flows: if +the injections you apply already sum to zero (a real transfer, generation minus +load), the resulting branch flows are independent of which bus was chosen as +slack. LODF values are likewise invariant to the slack choice. + +## Distributed slack + +A single slack is an idealization — real imbalance is picked up by many +generators according to participation factors, not one bus. A **distributed +slack** spreads the balancing withdrawal across several buses by weight. Instead +of the whole compensating withdrawal landing on one reference bus, each +participating bus *k* absorbs a share `wₖ / Σ w`. + +Concretely, the distributed-slack PTDF is the single-slack PTDF with each row's +weighted average subtracted: + +``` +PTDF_dist[i, j] = PTDF[i, j] − Σₖ (wₖ / Σ w) · PTDF[i, k] +``` + +so the sensitivities are now measured relative to the *weighted set* of slack +buses rather than one. The weights are normalized internally, so only their +ratios matter. + +Two structural requirements come with distributed slack (enforced by the +constructors): + + - There must be exactly **one reference bus** in the system when a non-empty + `dist_slack` is supplied — distributed slack redistributes the imbalance, but + the grounding of `ABA` is still a single reference. + - The weight vector must cover **every bus** (length equal to the bus count); + buses that do not participate simply carry weight zero. + +## The type differs by matrix + +The way you pass the weights depends on the matrix, and this catches people out: + +| Matrix | `dist_slack` type | +|:-------------------------------------------- |:--------------------------------------------------------- | +| [`PTDF`](@ref), [`VirtualPTDF`](@ref) | `Dict{Int, Float64}` — bus number → weight | +| [`VirtualLODF`](@ref), [`VirtualMODF`](@ref) | `Vector{Float64}` — one weight per bus, in bus-axis order | + +For [`PTDF`](@ref)/[`VirtualPTDF`](@ref) you supply a **dictionary keyed by bus +number**, which is convenient because you name only the participating buses; +internally it is normalized and expanded to a per-bus vector. For +[`VirtualLODF`](@ref)/[`VirtualMODF`](@ref) you supply the **positional vector** +directly. In both cases the default is empty, +meaning single-reference-bus slack. + +Why the LODF/MODF forms take a raw vector rather than a dict is an interface +detail, not a semantic difference — the underlying meaning (weighted balancing +withdrawal) is identical. Just match the type to the matrix you are building. + +## Choosing a convention + + - Use the **default single slack** for most analyses, and for validating against + references that assume one slack bus. + - Use a **distributed slack** when you want sensitivities consistent with how + generation actually rebalances — for economic dispatch, participation-factor + studies, or comparison with an AC solution whose losses and response are spread + across machines. + +The distributed slack does not make the DC model more or less exact; it changes +the *question* the PTDF answers, from "balanced at one bus" to "balanced across a +weighted fleet." + +## References + +[^ptdf]: Power transfer distribution factors, the reference/slack convention, and + distributed slack are standard sensitivity-analysis results; see A. J. Wood, + B. F. Wollenberg, and G. B. Sheblé, *Power Generation, Operation, and Control*, + 3rd ed., Wiley, 2013 (sensitivity factors), together with the PTDF/LODF + definitions in the Siemens PSS®E and PowerWorld Simulator documentation. diff --git a/docs/src/how_to_guides/build_multiple_matrices.jl b/docs/src/how_to_guides/build_multiple_matrices.jl new file mode 100644 index 000000000..9c0dc41dc --- /dev/null +++ b/docs/src/how_to_guides/build_multiple_matrices.jl @@ -0,0 +1,68 @@ +# # How to Build Multiple Matrices Without Repeating Work + +# Every matrix constructor that takes a [`System`](@extref PowerSystems.System) +# rebuilds the same intermediates from scratch — the [`Ybus`](@ref), the incidence +# matrix `A` ([`IncidenceMatrix`](@ref)), and the susceptance-weighted `BA` +# ([`BA_Matrix`](@ref)). This guide shows how to compute the shared pieces once and +# feed them to the constructors that accept pre-built matrices. + +using PowerNetworkMatrices +import PowerSystemCaseBuilder as PSB + +sys = PSB.build_system(PSB.PSITestSystems, "c_sys5"); + +# ## Build the shared intermediates once + +# The construction dependency chain is +# +# > `Ybus` → `IncidenceMatrix`, `BA_Matrix` → +# > `ABA_Matrix` / `PTDF`, and `PTDF` → `LODF`. +# +# The [`Ybus`](@ref) is the expensive shared root. Build it — and the incidence and +# BA matrices derived from it — exactly once: + +ybus = Ybus(sys) +A = IncidenceMatrix(ybus) +BA = BA_Matrix(ybus) + +# ## Reuse them across constructors + +# [`PTDF`](@ref) accepts the incidence and BA matrices directly, skipping its own +# [`Ybus`](@ref) build: + +ptdf = PTDF(A, BA) + +# [`LODF`](@ref) can be built straight from a [`PTDF`](@ref) you already have, +# reusing that work too — no second factorization of the network: + +lodf = LODF(A, ptdf) + +# Alternatively, the factorized [`ABA_Matrix`](@ref) route builds [`LODF`](@ref) +# from the same `A` and `BA`. All three inputs must share the same network +# reduction — which they do here, because they all descend from one `ybus`: + +aba = ABA_Matrix(ybus; factorize = true) +lodf_via_aba = LODF(A, aba, BA) + +# Virtual matrices likewise accept a pre-built [`Ybus`](@ref), so the lazy forms +# reuse the same root: + +vptdf = VirtualPTDF(ybus) + +# !!! note "Keep reductions consistent" +# +# Constructors that combine pre-built matrices (e.g. `LODF(A, ABA, BA)`) require +# every input to have been built with the **same** `network_reductions`. Because +# they all derive from a single [`Ybus`](@ref) here, they are automatically +# consistent. Pass `network_reductions` once, to the [`Ybus`](@ref) call, and +# everything downstream inherits it. See the [`NetworkReduction`](@ref) docstring +# for the keyword and its rules. + +# ## See also +# +# - [Matrix overview & indexing](@ref) — every matrix type, its axes, and how the +# shared intermediates fit together. +# - [How to Choose a Linear Solver](@ref) — the factorization cost that reuse +# avoids repeating. +# - [`NetworkReduction`](@ref) — supplying reductions via `network_reductions` +# to the shared [`Ybus`](@ref). diff --git a/docs/src/how_to_guides/choose_linear_solver.jl b/docs/src/how_to_guides/choose_linear_solver.jl new file mode 100644 index 000000000..c135c9de9 --- /dev/null +++ b/docs/src/how_to_guides/choose_linear_solver.jl @@ -0,0 +1,149 @@ +# # How to Choose a Linear Solver + +# This guide helps you select the appropriate linear solver for your network matrix computations. + +# ## Available Solvers + +# Pass the solver name as the `linear_solver` keyword to any matrix constructor +# ([`PTDF`](@ref), [`LODF`](@ref), [`ABA_Matrix`](@ref), …). `PowerNetworkMatrices.jl` supports four: +# +# 1. **`"KLU"`** - sparse +# [KLU](https://github.com/DrTimothyAldenDavis/SuiteSparse) factorization +# Always available (built-in `KLUWrapper` submodule); the default off Apple hardware. +# 2. **`"AppleAccelerateLU"`** - sparse LU via Apple's +# [Accelerate sparse solvers](https://developer.apple.com/documentation/accelerate/sparse_solvers). +# Always compiled in (built-in `AccelerateWrapper` submodule) but +# runtime-gated to macOS 15.5+ on Apple hardware, where it is the default. +# 3. **`"MKLPardiso"`** - Intel's +# [oneMKL PARDISO](https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/), +# wrapped by [`Pardiso.jl`](https://github.com/JuliaSparse/Pardiso.jl). A +# weak-dependency package extension: only loaded once you also add and import +# `Pardiso.jl`. +# 4. **`"Dense"`** - dense +# [LU](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.lu) +# from Julia's LinearAlgebra stdlib, for small or debugging cases. +# +# The default is platform-dependent: `AppleAccelerateLU` on macOS 15.5+ (Apple +# hardware), `KLU` everywhere else. KLU and Apple Accelerate are always-present +# submodules — only MKL Pardiso is an optional extension. + +# The examples below use a small test system loaded with `PowerSystemCaseBuilder`: + +using PowerNetworkMatrices +import PowerNetworkMatrices as PNM +import PowerSystemCaseBuilder as PSB + +sys = PSB.build_system(PSB.PSITestSystems, "c_sys5"); + +# ## Choosing the Right Solver + +# ### Use KLU When: +# +# - Working with typical power systems (most cases) +# - System size is medium to large (> 100 buses) +# - You want good performance without special dependencies +# - Running on any platform (Linux, macOS, Windows) +# +# Use [`PTDF`](@ref) with the KLU solver (the default off Apple hardware): + +ptdf_matrix = PTDF(sys) # platform default +# or explicitly: +ptdf_matrix = PTDF(sys; linear_solver = "KLU"); + +# ### Use Apple Accelerate When: +# +# - Running on Apple-silicon macOS 15.5 or newer +# - You want the platform-tuned sparse LU (it is the default there) +# +# Select it explicitly with: + +# ```julia +# ptdf_matrix = PTDF(sys; linear_solver = "AppleAccelerateLU"); +# ``` + +# ### Use Dense When: +# +# - System is very small (< 30 buses) +# - You're debugging or validating results +# - Matrix operations are simple and small-scale +# +# Specify the Dense solver explicitly: + +ptdf_matrix = PTDF(sys; linear_solver = "Dense"); + +# ### Use MKLPardiso When: +# +# - You have Intel processors +# - Running on Linux or Windows (not available on Apple silicon) +# - Maximum performance is critical +# - Working with very large systems (> 1000 buses) +# +# MKLPardiso lives in a weak-dependency package extension, so first add and +# import `Pardiso.jl` to load it, then request the solver: + +# ```julia +# using Pardiso # loads the MKLPardisoExt extension +# ptdf_matrix = PTDF(sys; linear_solver = "MKLPardiso") +# ``` + +# ## Performance Considerations + +# ### System Size + +# | Buses | Recommended Solver | +# |:------- |:------------------ | +# | < 30 | Dense or KLU | +# | 30-1000 | KLU | +# | > 1000 | KLU or MKLPardiso | + +# ### Platform Availability + +# | Solver | Linux | Windows | macOS | +# |:----------------- |:----- |:------- |:----------------- | +# | KLU | ✓ | ✓ | ✓ | +# | AppleAccelerateLU | ✗ | ✗ | ✓ (Apple, 15.5+) | +# | Dense | ✓ | ✓ | ✓ | +# | MKLPardiso | ✓ | ✓ | ✗ | +# +# `AppleAccelerateLU` needs no extra package — it is built in. Only +# `MKLPardiso` requires installing and importing `Pardiso.jl`. + +# ### Singular `ABA`: prefer KLU + +# The two default solvers differ in how they treat a singular `ABA` matrix. +# `AppleAccelerateLU` silently factorizes a singular matrix and returns garbage, +# whereas KLU raises. Prefer `"KLU"` whenever the `ABA` matrix may be singular — for +# example a full outage that isolates a bus, such as a 3-winding transformer's +# zero-injection star bus: + +# ```julia +# ptdf = PTDF(sys; linear_solver = "KLU") # safe when singularity is possible +# ``` + +# ### Persisting a preferred backend + +# The active sparse backend can be pinned across Julia sessions with `Preferences.jl` +# via the (non-exported) helpers in `src/linalg_settings.jl` — +# `PNM.set_linalg_backend_preference` / `PNM.get_linalg_backend_preference`, +# `PNM.set_linalg_backend_check` / `PNM.get_linalg_backend_check`, and +# `PNM.check_linalg_backend` (reports the active BLAS/LAPACK library and whether the +# requested backend is loaded). + +# ## Troubleshooting + +# ### MKLPardiso Not Available + +# If you get an error when using MKLPardiso: +# +# 1. Confirm you have added and imported `Pardiso.jl` (the extension only +# loads once `Pardiso` is available) +# 2. Verify you're on Linux or Windows (not macOS) +# 3. Check that you have Intel processors +# +# Fall back to KLU if MKLPardiso is unavailable — it is always present. + +# ## Related Topics +# +# - [How to Build Multiple Matrices Without Repeating Work](@ref) - use these solvers to build matrices +# - [Computational Considerations](@ref) - the reasoning behind the platform defaults +# - The [`AutoTolerance`](@ref) docstring - the orthogonal `tol` sparsification setting diff --git a/docs/src/how_to_guides/choose_linear_solver.md b/docs/src/how_to_guides/choose_linear_solver.md deleted file mode 100644 index d8ab04139..000000000 --- a/docs/src/how_to_guides/choose_linear_solver.md +++ /dev/null @@ -1,105 +0,0 @@ -# How to Choose a Linear Solver - -This guide helps you select the appropriate linear solver for your network matrix computations. - -## Available Solvers - -`PowerNetworkMatrices.jl` supports three linear solver methods: - - 1. **KLU** (default) - Sparse solver using KLU factorization - 2. **Dense** - Dense matrix operations - 3. **MKLPardiso** - Intel's MKL Pardiso solver (Intel systems only) - -## Choosing the Right Solver - -### Use KLU When: - - - Working with typical power systems (most cases) - - System size is medium to large (> 100 buses) - - You want good performance without special dependencies - - Running on any platform (Linux, macOS, Windows) - -Use [`PTDF`](@ref) with the default KLU solver: - -```julia -ptdf_matrix = PTDF(sys) # KLU is the default -# or explicitly: -ptdf_matrix = PTDF(sys; linear_solver = "KLU") -``` - -### Use Dense When: - - - System is very small (< 30 buses) - - You're debugging or validating results - - Matrix operations are simple and small-scale - -Specify the Dense solver explicitly: - -```julia -ptdf_matrix = PTDF(sys; linear_solver = "Dense") -``` - -### Use MKLPardiso When: - - - You have Intel processors - - Running on Linux or Windows (not available on macOS) - - Maximum performance is critical - - Working with very large systems (> 1000 buses) - -Specify the MKLPardiso solver: - -```julia -ptdf_matrix = PTDF(sys; linear_solver = "MKLPardiso") -``` - -## Performance Considerations - -### System Size - -| Buses | Recommended Solver | -|:------- |:------------------ | -| < 30 | Dense or KLU | -| 30-1000 | KLU | -| > 1000 | KLU or MKLPardiso | - -### Platform Availability - -| Solver | Linux | Windows | macOS | -|:---------- |:----- |:------- |:----- | -| KLU | ✓ | ✓ | ✓ | -| Dense | ✓ | ✓ | ✓ | -| MKLPardiso | ✓ | ✓ | ✗ | - -## Switching Solvers - -You can easily switch between solvers to compare performance: - -```julia -using BenchmarkTools - -# Benchmark KLU -@btime ptdf_klu = PTDF($sys; linear_solver = "KLU") - -# Benchmark Dense -@btime ptdf_dense = PTDF($sys; linear_solver = "Dense") - -# Benchmark MKLPardiso (if available) -@btime ptdf_mkl = PTDF($sys; linear_solver = "MKLPardiso") -``` - -## Troubleshooting - -### MKLPardiso Not Available - -If you get an error when using MKLPardiso: - - 1. Verify you're on Linux or Windows (not macOS) - 2. Check that you have Intel processors - 3. Ensure MKL dependencies are installed - -Fall back to KLU if MKLPardiso is unavailable. - -## Related Topics - - - [How to Compute Network Matrices](@ref) - Learn how to use these solvers - - [PTDF matrix](@ref) - Detailed walkthrough with solver comparisons diff --git a/docs/src/how_to_guides/compute_network_matrices.md b/docs/src/how_to_guides/compute_network_matrices.md deleted file mode 100644 index b103fcd00..000000000 --- a/docs/src/how_to_guides/compute_network_matrices.md +++ /dev/null @@ -1,161 +0,0 @@ -# How to Compute Network Matrices - -This guide shows you how to compute various network matrices for your power system. - -## Prerequisites - - - `PowerNetworkMatrices.jl` installed - - A power system model loaded (see [Quick Start Guide](@ref)) - -## Computing PTDF Matrix - -To compute the Power Transfer Distribution Factor matrix using [`PTDF`](@ref): - -```julia -using PowerNetworkMatrices -const PNM = PowerNetworkMatrices - -# Assuming you have a system loaded -ptdf_matrix = PNM.PTDF(sys) - -# Access the matrix data (in standard arcs × buses orientation) -matrix_data = PNM.get_ptdf_data(ptdf_matrix) -``` - -### Indexing PTDF Elements - -The PTDF matrix is indexed by **arc tuples** `(from_bus, to_bus)` and **bus numbers**: - -```julia -# Access PTDF element for arc (1, 2) and bus 3 -ptdf_matrix[(1, 2), 3] - -# Inspect available axes and lookup dictionaries -PNM.get_axes(ptdf_matrix) -PNM.get_lookup(ptdf_matrix) -``` - -## Computing LODF Matrix - -To compute the Line Outage Distribution Factor matrix using [`LODF`](@ref): - -```julia -lodf_matrix = PNM.LODF(sys) -matrix_data = PNM.get_lodf_data(lodf_matrix) -``` - -### Indexing LODF Elements - -The LODF matrix is indexed by **arc tuples** for both dimensions: - -```julia -# Access LODF element: flow change on arc (1, 4) due to outage of arc (2, 3) -lodf_matrix[(1, 4), (2, 3)] - -# Inspect available axes and lookup dictionaries -PNM.get_axes(lodf_matrix) -PNM.get_lookup(lodf_matrix) -``` - -## Computing Virtual PTDF Matrix - -For systems where you need virtual representation using [`VirtualPTDF`](@ref): - -```julia -vptdf_matrix = PNM.VirtualPTDF(sys) - -# Access element by arc tuple and bus number -vptdf_matrix[(1, 2), 3] -``` - -## Computing Virtual LODF Matrix - -Similarly for virtual LODF using [`VirtualLODF`](@ref): - -```julia -vlodf_matrix = PNM.VirtualLODF(sys) - -# Access element by arc tuples -vlodf_matrix[(1, 2), (3, 4)] -``` - -## Computing Virtual MODF Matrix - -For post-contingency / post-modification PTDF rows using [`VirtualMODF`](@ref): - -```julia -vmodf_matrix = PNM.VirtualMODF(sys) - -# Inspect contingencies auto-registered from PSY.Outage supplemental attributes -PNM.get_registered_contingencies(vmodf_matrix) - -# Access the post-modification PTDF row for monitored arc (1, 2) under a contingency -contingency = first(values(PNM.get_registered_contingencies(vmodf_matrix))) -vmodf_matrix[(1, 2), contingency] -``` - -## Computing Incidence and BA Matrices - -For the fundamental network topology matrices: - -Compute the incidence matrix using [`IncidenceMatrix`](@ref): - -```julia -incidence_matrix = PNM.IncidenceMatrix(sys) - -# Axes are (arc_tuples, bus_numbers) -PNM.get_axes(incidence_matrix) -``` - -Compute the BA matrix (Bus-Admittance) using [`BA_Matrix`](@ref): - -```julia -ba_matrix = PNM.BA_Matrix(sys) -``` - -Compute the ABA matrix using [`ABA_Matrix`](@ref): - -```julia -aba_matrix = PNM.ABA_Matrix(sys) -``` - -## Working with Pre-computed Matrices - -If you have already computed the incidence and BA matrices, you can use them to compute [`PTDF`](@ref): - -```julia -# Compute base matrices first -ba_matrix = PNM.BA_Matrix(sys) -a_matrix = PNM.IncidenceMatrix(sys) - -# Use them to compute PTDF -ptdf_matrix = PNM.PTDF(a_matrix, ba_matrix) -``` - -## Understanding Axes and Lookup Dictionaries - -All network matrices store `axes` and `lookup` fields that describe how rows and columns map to physical network elements: - - - **`axes`**: A tuple of vectors containing the identifiers for each dimension - - **`lookup`**: A tuple of dictionaries mapping identifiers to integer indices - -For matrices involving branches (IncidenceMatrix, PTDF, LODF), branches are represented as **arc tuples** `(from_bus_number, to_bus_number)` rather than branch name strings. This provides a compact, unambiguous identifier for each directed branch in the network. - -| Matrix | Dimension 1 (rows) | Dimension 2 (columns) | -|:----------------- |:------------------ |:--------------------- | -| `IncidenceMatrix` | Arc tuples | Bus numbers | -| `PTDF` | Arc tuples | Bus numbers | -| `LODF` | Arc tuples | Arc tuples | -| `Ybus` | Bus numbers | Bus numbers | -| `VirtualPTDF` | Arc tuples | Bus numbers | -| `VirtualLODF` | Arc tuples | Arc tuples | -| `VirtualMODF` | Arc tuples | Bus numbers | - -!!! note - - For backward compatibility, branch name strings can also be used to index PTDF and LODF matrices. This uses the `get_branch_multiplier` function internally to map names to arc tuples. Using arc tuples directly is recommended. - -## Next Steps - - - Learn about choosing linear solvers for optimal performance - - Understand the theory behind network matrices in the Explanation section diff --git a/docs/src/how_to_guides/contingencies.jl b/docs/src/how_to_guides/contingencies.jl new file mode 100644 index 000000000..bb9c6973d --- /dev/null +++ b/docs/src/how_to_guides/contingencies.jl @@ -0,0 +1,176 @@ +# # How to Define and Apply Contingencies + +# This guide shows how to compute post-contingency PTDF rows with +# [`VirtualMODF`](@ref) — the lazy Multiple Outage Distribution Factor matrix. +# You will attach outages to a system, let them auto-register, query monitored +# arcs under a contingency, and build manual modifications when you need full +# control. + +# !!! note +# There is **no dense `MODF` type**. Post-contingency factors are only +# available through [`VirtualMODF`](@ref), which computes rows on demand via the +# Woodbury identity. See [Flowgate Methodology](@ref) for the theory. + +# ## Prerequisites +# +# - `PowerNetworkMatrices.jl` and `PowerSystems.jl` installed +# - A power system model + +using PowerNetworkMatrices +import PowerNetworkMatrices as PNM +import PowerSystems as PSY +import PowerSystemCaseBuilder as PSB + +sys = PSB.build_system(PSB.PSITestSystems, "c_sys5"); + +# ## Attach Outages to the System + +# Contingencies are defined as [`PSY.Outage`](@extref PowerSystems.Outage) supplemental attributes on the +# components they trip. When a contingency only needs to *exist*, use a +# [`PSY.FixedForcedOutage`](@extref PowerSystems.FixedForcedOutage) with `outage_status = 1.0` (outaged) and attach it to +# each branch: + +for branch in PSY.get_components(PSY.ACTransmission, sys) + outage = PSY.FixedForcedOutage(; outage_status = 1.0) + PSY.add_supplemental_attribute!(sys, branch, outage) +end + +# ## Build the VirtualMODF + +# Registration is **automatic**: the constructor scans the system for [`PSY.Outage`](@extref PowerSystems.Outage) +# attributes and resolves each into a [`ContingencySpec`](@ref). There is no +# public `register_contingency` — construct the matrix from a system that already +# carries its outages. + +vmodf = VirtualMODF(sys) + +# Inspect what was registered with [`get_registered_contingencies`](@ref). It +# returns a `Dict{UUID, ContingencySpec}` keyed by the source outage's UUID: + +contingencies = get_registered_contingencies(vmodf) + +# ## Query a Monitored Arc Under a Contingency + +# Index the matrix as `vmodf[monitored_arc, spec]`. The monitored arc is an arc +# tuple `(from, to)` (or its integer index); the returned value is the full +# post-contingency PTDF row for that arc — one sensitivity per bus. + +# Pick a monitored arc and outage a *different* arc — monitoring an element that +# the contingency itself outages is undefined and raises. Here the spec is built +# straight from an arc with the convenience [`NetworkModification`](@ref) +# constructor: + +arcs = PNM.get_arc_axis(vmodf); +monitored_arc = arcs[1]; +ctg = NetworkModification(vmodf, arcs[2]); + +# The returned row carries one post-contingency sensitivity per bus: + +row = vmodf[monitored_arc, ctg] + +# The second index accepts three equivalent forms — the +# [`NetworkModification`](@ref) used above, a [`ContingencySpec`](@ref) from the +# registered set, or the original [`PSY.Outage`](@extref PowerSystems.Outage) (by +# its registered UUID). All resolve to the same [`NetworkModification`](@ref) and +# share the cached Woodbury factors, so repeated queries for one contingency across +# different monitored arcs reuse work: + +# ```julia +# spec = first(values(contingencies)) # a registered ContingencySpec +# vmodf[monitored_arc, spec] +# vmodf[monitored_arc, spec.modification] # its NetworkModification +# +# branch = first(PSY.get_components(PSY.ACTransmission, sys)) +# outage = first(PSY.get_supplemental_attributes(branch)) +# vmodf[monitored_arc, outage] # the PSY.Outage, by UUID +# ``` + +# ## The modification type model + +# Under the convenience constructor sit a few value types, layered from the smallest +# unit up to the solver-ready form. It helps to know them before dropping to the +# manual path: +# +# | Type | Represents | Scope | +# |:----------------------------- |:---------------------------------------------------------------------- |:-------------------------- | +# | [`ArcModification`](@ref) | A susceptance change on one aggregated arc, plus optional Ybus Pi-model deltas | One arc | +# | [`ShuntModification`](@ref) | A diagonal admittance change on one bus | One bus | +# | [`NetworkModification`](@ref) | A canonical, [`System`](@extref PowerSystems.System)-independent bundle of arc and shunt changes plus islanding status | Whole modification | +# | [`ContingencySpec`](@ref) | A [`NetworkModification`](@ref) tagged with the source [`PSY.Outage`](@extref PowerSystems.Outage) UUID | One registered contingency | +# +# [`NetworkModification`](@ref) is the canonical representation: once built it holds no +# reference to the source [`System`](@extref PowerSystems.System) and serves as the +# cache key inside [`VirtualMODF`](@ref) (its `label` is excluded from equality, so two +# physically identical modifications compare equal regardless of name). +# +# !!! note +# Partial (non-full-outage) susceptance changes are supported only on **direct and +# parallel** arcs. Series-reduced arcs and 3-winding transformer windings accept +# only a full outage of the equivalent; anything else raises an error. + +# ## Build a Modification Manually + +# The convenience constructor used above (`NetworkModification(matrix, arc)`, or +# `NetworkModification(matrix, branch)`) is the simplest path — it looks up the +# arc's susceptance and populates the deltas for you. When you want full control, +# assemble the low-level building blocks instead. An [`ArcModification`](@ref) is a +# susceptance change on one arc (`delta_b` negative for an outage); a +# [`ShuntModification`](@ref) is an admittance change on one bus. Both are indexed +# by their **integer** position in the matrix: + +# ```julia +# arc_index = PNM.get_arc_lookup(vmodf)[(1, 4)] +# arc_mod = ArcModification(arc_index, -5.0) # Δb removes the arc's susceptance +# +# bus_index = PNM.get_bus_lookup(vmodf)[3] +# shunt_mod = ShuntModification(bus_index, ComplexF32(-0.1im)) +# +# # Combine arc and shunt changes into one modification (label, arcs, shunts, islanding) +# custom = NetworkModification("arc_and_shunt", [arc_mod], [shunt_mod], false) +# vmodf[monitored_arc, custom] +# ``` + +# Prefer the convenience constructors over hand-built [`ArcModification`](@ref) +# values: they compute physically consistent `delta_b` and Pi-model deltas from the +# network data, which is otherwise your responsibility to get right. + +# ## One-Shot Post-Modification Rows from a VirtualPTDF + +# If you already hold a [`VirtualPTDF`](@ref) and want a single post-modification +# row without registering contingencies, use +# [`get_post_modification_ptdf_row`](@ref). It applies a [`NetworkModification`](@ref) +# through the same Woodbury correction: + +vptdf = VirtualPTDF(sys) +varcs = PNM.get_arc_axis(vptdf); +mod = NetworkModification(vptdf, varcs[2]); +row_oneshot = get_post_modification_ptdf_row(vptdf, varcs[1], mod) + +# Indexing is the equivalent form — it returns the same row: + +isapprox(vptdf[varcs[1], mod], row_oneshot) + +# This function does **no caching** — each call recomputes. When querying many +# monitored arcs for the *same* modification, precompute once with +# [`compute_woodbury_factors`](@ref) and reuse via [`apply_woodbury_correction`](@ref). + +# ## Contingencies and Network Reduction + +# If you build the [`VirtualMODF`](@ref) with `network_reductions`, any branch that a +# contingency outages or monitors must survive every reduction step. Outage and +# monitored-component buses are auto-protected from reduction. Declare monitored +# branches on the outage so their buses are kept: + +# ```julia +# monitored_line = PSY.get_component(PSY.ACTransmission, sys, "2") +# PSY.set_monitored_components!(outage, [monitored_line]) +# ``` + +# Querying a monitored arc that was reduced away raises a clear error rather than +# silently returning the base row. + +# ## See Also +# +# - [Public API Reference](@ref) — full docstrings for [`NetworkModification`](@ref), +# [`ContingencySpec`](@ref), [`compute_woodbury_factors`](@ref), and the rest +# - [Flowgate Methodology](@ref) — the Woodbury post-contingency theory diff --git a/docs/src/how_to_guides/diagnose_connectivity.jl b/docs/src/how_to_guides/diagnose_connectivity.jl new file mode 100644 index 000000000..31da2bddb --- /dev/null +++ b/docs/src/how_to_guides/diagnose_connectivity.jl @@ -0,0 +1,103 @@ +# # How to Diagnose a Disconnected Network + +# A singular `ABA` matrix or a failed DC power flow is frequently just a +# disconnected network: an island with no reference bus leaves `ABA` singular. +# Checking connectivity first localizes the problem before you dig into the numerics. +# This guide walks through the check, then deliberately breaks a network so you can see +# exactly what a fragmented result looks like — and how to get back. + +using PowerNetworkMatrices +import PowerSystems as PSY +import PowerSystemCaseBuilder as PSB + +sys = PSB.build_system(PSB.PSITestSystems, "c_sys5"); + +# ## Step 1 — Confirm a healthy network is connected + +# [`validate_connectivity`](@ref) returns `true` when the system forms a single +# connected component: + +validate_connectivity(sys) + +# [`find_subnetworks`](@ref) shows the decomposition behind that answer: a `Dict` +# mapping each island's reference bus to the set of bus numbers in it. A connected +# system yields a **single** entry: + +find_subnetworks(sys) + +# Both functions also accept an already-built [`AdjacencyMatrix`](@ref) or +# [`Ybus`](@ref), so a matrix you already have on hand is reused instead of rebuilt: + +adj = AdjacencyMatrix(sys) +validate_connectivity(adj) + +# ## Step 2 — Disconnect a bus and watch it split + +# To see a fragmented result on a real system, let's isolate one bus. A bus goes silent +# when every branch touching it is out of service, so we find bus `5`'s incident +# branches and mark them unavailable — `Ybus` (and therefore the connectivity check) +# only includes available branches: + +isolated_bus = 5 +incident = [ + br for br in PSY.get_components(PSY.ACBranch, sys) if + PSY.get_number(PSY.get_from(PSY.get_arc(br))) == isolated_bus || + PSY.get_number(PSY.get_to(PSY.get_arc(br))) == isolated_bus +] + +for br in incident + PSY.set_available!(br, false) +end + +# The network is now split. `validate_connectivity` reports it: + +validate_connectivity(sys) + +# ...and `find_subnetworks` returns **two** entries — the main island, and bus `5` +# stranded on its own: + +find_subnetworks(sys) + +# There is the diagnosis. That second island — the isolated `{5}` — has no reference +# bus of its own, which is exactly the block that would have made `ABA` singular. The +# bus set tells you precisely which buses to reconnect (or which island to study in +# isolation). Here it points straight back at the bus we broke. + +# ## Step 3 — Reconnect and recover + +# Restoring the branches we took out puts the network back together — bus `5` rejoins +# the main island and connectivity is whole again: + +for br in incident + PSY.set_available!(br, true) +end + +validate_connectivity(sys) + +# ...and the decomposition is back to a single island, identical to where we started: + +find_subnetworks(sys) + +# ## Choosing a traversal algorithm + +# The lower-level `find_subnetworks(M, bus_numbers; subnetwork_algorithm)` — which runs +# over a raw sparse connectivity matrix — lets you pick how the graph is walked: +# +# - [`iterative_union_find`](@ref) (the **default**) — an iterative union-find +# disjoint-set, safe on networks of any size. +# - [`depth_first_search`](@ref) — a recursive traversal. +# +# Both return the **same** island decomposition, so the choice is about performance, +# not correctness. Prefer the default union-find; it avoids the deep recursion that +# `depth_first_search` can hit on very large networks. The `subnetwork_algorithm` +# keyword also threads through the matrix constructors, so islands are detected the same +# way at build time as by an explicit [`find_subnetworks`](@ref) call: + +ABA_Matrix(sys; subnetwork_algorithm = depth_first_search); + +# ## See also +# +# - [Matrix overview & indexing](@ref) — the [`AdjacencyMatrix`](@ref) and +# [`Ybus`](@ref) graphs these checks traverse, and per-island axes. +# - [Network Reduction Theory](@ref) — how the susceptance graph can fragment into +# more islands than the admittance graph, and why that matters for `ABA`. diff --git a/docs/src/how_to_guides/reproduce_dfax_values.jl b/docs/src/how_to_guides/reproduce_dfax_values.jl new file mode 100644 index 000000000..3201cfcdf --- /dev/null +++ b/docs/src/how_to_guides/reproduce_dfax_values.jl @@ -0,0 +1,350 @@ +# # Reproduce industry DFAX values + +# To validate `PowerNetworkMatrices` against the industry's "DFAX" vocabulary +# (Distribution Factors, as used in PSS/E's DFAX activity and downstream by +# the NERC IDC for TLR and Congestion Management Procedures), +# follow the recipe below. Every flavor of DFAX — GSF, LSF, LODF, OTDF, transfer +# DFAX, flowgate DFAX, and multi-element (N-k) DFAX — is a special case of +# one unified formula. This guide reproduces each case with executable examples on +# the RTS-GMLC system and checks the results against their closed-form references. + +# For the Woodbury-identity derivation behind post-contingency PTDF rows, see +# the [Flowgate Methodology](@ref) explanation. + +# ## What is DFAX? + +# "DFAX" is shorthand for *distribution factor*. The term originated with +# PSS/E's DFAX activity, which writes a `.dfx` file consumed by the NERC +# Interchange Distribution Calculator (IDC). The IDC uses these distribution +# factors during Transmission Loading Relief (TLR) procedures and Congestion +# Management Procedures (CMP) — for example, applying a 5% threshold to +# decide whether a given source-to-sink transfer is a "significant" contributor to +# a monitored flowgate. + +# In practice DFAX is an umbrella term that covers several specific +# quantities. The table below maps each industry term to the +# `PowerNetworkMatrices` primitive that computes it: + +# | Industry term | What it answers | PNM primitive | +# |:------------------------- |:------------------------------------------------ |:---------------------------------------------------------------------- | +# | GSF / ISF | Δflow on `m` per 1 MW injection at bus `b` | `PTDF[m, b]` | +# | LSF | Same as GSF for loads (opposite sign) | `-PTDF[m, b]` | +# | LODF | Flow redistribution after one branch trips | `LODF[m, c]` | +# | OTDF | GSF with a contingency already in place | `PTDF[m,b] + LODF[m,c]·PTDF[c,b]`, or one entry of a `VirtualMODF` row | +# | Transfer DFAX (pre-cont.) | Fraction of a source→sink transfer reaching `m` | `PTDF[m,:]·(s_v − k_v)` | +# | Flowgate DFAX | Transfer DFAX on a (monitored, contingency) pair | `VirtualMODF[m, ctg]·(s_v − k_v)` | +# | Multi-element (N-k) DFAX | Same with multiple simultaneous outages | `VirtualMODF` (multi-arc `NetworkModification`) | + +# ### The Phase Shifter Factor (PSF) +# +# The **Phase Shifter Factor** belongs to the same vocabulary but is not a +# first-class primitive, because a phase shift is not a topology change — it enters +# the DC model as a pair of nodal injections. A shift ``\alpha`` on branch ``c`` +# from bus ``f`` to bus ``t``, with series susceptance ``b_c``, is equivalent to +# injecting ``+b_c\,\alpha`` at ``f`` and ``-b_c\,\alpha`` at ``t``. Since the +# `PTDF` already maps injections to flows, the sensitivity of a monitored arc ``m`` +# to that shift is read straight from two `PTDF` columns: +# +# ```math +# \mathrm{PSF}[m, c] \;=\; b_c \,\bigl(\mathrm{PTDF}[m, f] - \mathrm{PTDF}[m, t]\bigr), +# ``` +# +# where ``b_c`` is `PSY.get_series_susceptance` of the phase-shifting transformer. +# The rest of this guide covers the flow-based distribution factors. + +# ## The unified DFAX formula + +# In the DC power-flow model every flavor of DFAX is a special case of the +# same quantity. For a monitored arc ``m``, a source participation vector +# ``s_v``, a sink participation vector ``k_v``, and a (possibly empty) set of +# network modifications ``C``, + +# ```math +# \mathrm{DFAX}(m,\ s \to k,\ C) \;=\; \mathrm{PTDF}_C[m,\,:] \cdot (s_v - k_v), +# ``` + +# where ``\mathrm{PTDF}_C`` is the post-modification PTDF (equal to the base +# ``\mathrm{PTDF}`` when ``C = \emptyset``). The formula has two degrees of +# freedom — *who is shifting* (the source/sink vectors) and *what state the +# network is in* (the contingency ``C``). Each section below fixes one or +# both. + +# The reference (slack) bus is implicit in `PTDF`: a row of `PTDF` already +# encodes "inject at bus ``b``, absorb at the slack". So setting ``k_v = 0`` +# in the formula means "let the slack absorb the sink", and the GSF section +# below reduces to a single `PTDF` entry. + +# ## Setup + +# All subsequent sections build on this setup block. It loads the RTS-GMLC +# system and constructs the three matrices the rest of the tutorial uses. + +import PowerSystems as PSY +using PowerNetworkMatrices +import PowerNetworkMatrices as PNM +import PowerSystemCaseBuilder as PSB +using DataFrames + +sys = PSB.build_system(PSB.PSISystems, "RTS_GMLC_DA_sys"); + +ptdf = PTDF(sys); +lodf = LODF(sys); +vmodf = VirtualMODF(sys); + +# `VirtualMODF` is the most general object — it can compute post-modification +# PTDF rows under any contingency. We also build `PTDF` and `LODF` up front +# because the pre-contingency and single-element-outage sections below use +# them directly (faster than going through Woodbury when those special cases +# apply). + +# ## GSF and LSF (no contingency, point source) + +# The simplest special case sets ``k_v = 0`` (the slack absorbs the sink) and +# ``s_v = e_b`` (a unit vector at one bus). The unified formula collapses to +# a single `PTDF` entry — this is the **Generation Shift Factor**: + +m = (107, 203); # monitored arc AB1 (Area 1 → Area 2) +b = 101; # injection bus in Area 1 +gsf = ptdf[m, b] + +# The **Load Shift Factor** is the same quantity with the opposite sign +# (loads withdraw power instead of inject): + +lsf = -gsf + +# ### Subsystem-aggregated GSF + +# In practice analysts care about a *subsystem* of generators (for example, +# all generators in an area) rather than a single bus. Build a participation +# vector by weighting each generator's bus by its `Pmax` share within the +# subsystem, then dot the vector with the `PTDF` row: + +area1_gens = filter( + g -> PSY.get_name(PSY.get_area(PSY.get_bus(g))) == "1", + collect(PSY.get_available_components(PSY.Generator, sys)), +); + +total_pmax = sum(PSY.get_max_active_power, area1_gens); + +src_weights = Dict{Int, Float64}(); +for g in area1_gens + bn = PSY.get_number(PSY.get_bus(g)) + src_weights[bn] = get(src_weights, bn, 0.0) + + PSY.get_max_active_power(g) / total_pmax +end + +gsf_area1 = sum(w * ptdf[m, bn] for (bn, w) in src_weights) + +# `gsf_area1` is the fraction of an aggregate 1 MW dispatch increase across +# all Area 1 generators (split by `Pmax`) that lands on AB1. The slack still +# absorbs the corresponding withdrawal — this is a *one-sided* shift. + +# If the swing should be distributed across many buses instead of falling on +# the single reference bus, pass a `dist_slack` dictionary to the `PTDF` +# constructor (see the [`PTDF`](@ref) constructor). That is a different +# concept from subsystem aggregation: `dist_slack` redefines the reference, +# whereas the participation vector above defines the *source* of the +# transfer. + +# ## Transfer DFAX (no contingency, multi-bus source and sink) + +# When both source and sink are subsystems, the unified formula is the +# difference of two weighted `PTDF` row dot-products: + +# ```math +# \mathrm{TDF}(m,\ s \to k) \;=\; \mathrm{PTDF}[m,\,:] \cdot s_v +# \;-\; \mathrm{PTDF}[m,\,:] \cdot k_v. +# ``` + +# Build the sink vector from Area 2 loads, max-active-power weighted. We +# filter to `PSY.PowerLoad` because the abstract `ElectricLoad` type also +# covers shunt admittance components (`FixedAdmittance`) that don't carry a +# real-power weight: + +area2_loads = filter( + l -> PSY.get_name(PSY.get_area(PSY.get_bus(l))) == "2", + collect(PSY.get_available_components(PSY.PowerLoad, sys)), +); + +total_load = sum(PSY.get_max_active_power, area2_loads); + +snk_weights = Dict{Int, Float64}(); +for l in area2_loads + bn = PSY.get_number(PSY.get_bus(l)) + snk_weights[bn] = get(snk_weights, bn, 0.0) + + PSY.get_max_active_power(l) / total_load +end + +# The pre-contingency transfer DFAX for Area 1 → Area 2 on AB1 is then: + +tdf_pre = + sum(w * ptdf[m, bn] for (bn, w) in src_weights) - + sum(w * ptdf[m, bn] for (bn, w) in snk_weights) + +# `tdf_pre` answers: *if Area 1 ramps up by 1 MW (split by generator `Pmax`) +# and Area 2's load grows by 1 MW (split by load size), what fraction of +# that transfer shows up on AB1?* In market and TLR settings, this is the +# pre-contingency component of the flowgate impact. + +# ## OTDF (single contingency, point source) + +# The **Outage Transfer Distribution Factor** is the GSF you would observe +# if a specific outage were already in effect. For a single-element +# contingency on arc ``c``, OTDF has a closed-form expression in terms of +# `PTDF` and `LODF`: + +# ```math +# \mathrm{OTDF}(m, b, c) \;=\; \mathrm{PTDF}[m, b] + \mathrm{LODF}[m, c] \cdot \mathrm{PTDF}[c, b]. +# ``` + +# This is the unified formula with ``C = \{c\}`` and the slack absorbing the +# sink. `VirtualMODF` computes the same quantity through the Woodbury +# identity, which generalizes naturally to multi-element contingencies (see +# the N-k section below). For a single outage the two routes agree: + +c = (113, 215); # contingency: AB2 outage +otdf_closed = ptdf[m, b] + lodf[m, c] * ptdf[c, b] + +# + +ctg = NetworkModification(vmodf, c); +row_c = vmodf[m, ctg]; +bus_lookup = PNM.get_bus_lookup(vmodf); +otdf_vmodf = row_c[bus_lookup[b]] + +# + +isapprox(otdf_closed, otdf_vmodf; rtol = 1e-10) + +# The `isapprox` check is this guide's internal validation: `VirtualMODF` +# and the closed-form LODF expansion are the same calculation expressed two +# different ways. Whenever both apply, they agree to floating-point +# tolerance. + +# ## Flowgate DFAX (single contingency, source–sink transfer) + +# A *flowgate* in NERC parlance is the pair `(monitored facility, contingency)`. The flowgate DFAX is the unified formula with both +# nontrivial source/sink vectors and a nonempty ``C``: the source–sink +# subtraction from the transfer-DFAX section applied to the post-contingency +# row from the OTDF section. We reuse `row_c`, `src_weights`, and +# `snk_weights` already in scope: + +flowgate_dfax = + sum(w * row_c[bus_lookup[bn]] for (bn, w) in src_weights) - + sum(w * row_c[bus_lookup[bn]] for (bn, w) in snk_weights) + +# The NERC 5% rule treats a transfer as a "significant" contributor to a +# flowgate when the absolute DFAX exceeds 0.05. The check is one line: + +significant = abs(flowgate_dfax) >= 0.05 + +# When `significant == true`, the transfer is subject to curtailment or +# mitigation under the relevant TLR procedure. + +# ## N-k DFAX (multi-element contingency) + +# When the contingency `C` contains more than one element, the closed-form +# LODF expansion of the OTDF section no longer applies — there is no scalar +# `LODF[m, c]` when `c` is itself a set. The unified formula still applies, +# and `VirtualMODF` is built to handle it directly. Build the multi-element +# modification by merging the `arc_modifications` of each single-arc +# `NetworkModification` into one combined object: + +mod_ab2 = NetworkModification(vmodf, (113, 215)); # AB2 outage +mod_ab3 = NetworkModification(vmodf, (123, 217)); # AB3 outage + +ctg_n2 = NetworkModification( + "AB2_and_AB3_out", + vcat(collect(mod_ab2.arc_modifications), + collect(mod_ab3.arc_modifications)), +); + +row_n2 = vmodf[m, ctg_n2]; + +flowgate_dfax_n2 = + sum(w * row_n2[bus_lookup[bn]] for (bn, w) in src_weights) - + sum(w * row_n2[bus_lookup[bn]] for (bn, w) in snk_weights) + +# Removing two of the three parallel Area 1 → Area 2 paths forces a much +# larger fraction of any inter-area transfer onto AB1, so `flowgate_dfax_n2` +# is substantially larger than the N-1 value computed in the previous +# section. The same indexing call (`vmodf[m, ctg]`) handles N-1, N-2, and +# higher orders — that is the operational advantage of going through +# `VirtualMODF`. + +# ## Capstone: assembling a DFAX report + +# In production use, an analyst typically wants a *table* of distribution +# factors covering several transfers, several monitored facilities, and +# several contingencies — the kind of report that a PSS/E `.dfx` file plus +# IDC post-processing produces. Building that report is one nested loop +# around the unified formula. + +# For this example we use one transfer (Area 1 → Area 2 from the transfer +# DFAX section), three monitored arcs (the three parallel Area 1 → Area 2 +# paths), and three contingencies (each of the other two paths individually, +# plus the N-2 double-outage from the previous section). Each contingency +# carries the set of arc tuples it outages so that we can skip the +# ill-defined case of monitoring an outaged element: + +monitored = [(107, 203), (113, 215), (123, 217)]; # AB1, AB2, AB3 + +ctg_ab2 = NetworkModification(vmodf, (113, 215)); +ctg_ab3 = NetworkModification(vmodf, (123, 217)); +ctg_ab2_ab3 = NetworkModification( + "AB2_and_AB3_out", + vcat(collect(ctg_ab2.arc_modifications), + collect(ctg_ab3.arc_modifications)), +); + +contingencies = [ + ("AB2 out", Set([(113, 215)]), ctg_ab2), + ("AB3 out", Set([(123, 217)]), ctg_ab3), + ("AB2 & AB3 out", Set([(113, 215), (123, 217)]), ctg_ab2_ab3), +]; + +rows = NamedTuple[] +for mon in monitored + for (label, outaged, ctg_k) in contingencies + mon in outaged && continue + row = vmodf[mon, ctg_k] + df = + sum(w * row[bus_lookup[bn]] for (bn, w) in src_weights) - + sum(w * row[bus_lookup[bn]] for (bn, w) in snk_weights) + push!( + rows, + ( + monitored = mon, + contingency = label, + dfax = df, + significant = abs(df) >= 0.05, + ), + ) + end +end + +report = sort(DataFrame(rows), :dfax; by = abs, rev = true) + +# The sort by `abs(dfax)` puts the largest flowgate impacts at the top — the +# ones a TLR coordinator would investigate first. Filtering to +# `report[report.significant, :]` would keep only NERC-significant rows. + +# This table is the kind of output that drives downstream congestion and +# seams-coordination workflows; assembling it requires nothing beyond the +# matrices in this guide. + +# ## When to use which primitive + +# | Need | Reach for | +# |:--------------------------------------------------- |:--------------------------------------------- | +# | Many transfers, no contingencies | `PTDF` | +# | One contingency, all monitored branches | `LODF` + `PTDF` | +# | Specific flowgates `(monitored, contingency)` pairs | `VirtualMODF` | +# | Multi-element / N-k contingencies | `VirtualMODF` | +# | Memory-constrained or sparse usage | `VirtualPTDF` / `VirtualLODF` / `VirtualMODF` | + +# `VirtualMODF` is strictly more general than `LODF`-based OTDF arithmetic, +# but the closed-form route in the OTDF section is faster when you only have +# a single outage and many bus injections to evaluate. The decision is about +# which *direction* of the matrix you traverse most often, not about which +# one is "correct". diff --git a/docs/src/index.md b/docs/src/index.md index 873831ecd..ffde204a6 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -7,62 +7,44 @@ CurrentModule = PowerNetworkMatrices ## Overview `PowerNetworkMatrices.jl` is a [`Julia`](http://www.julialang.org) package for -the evaluation of network matrices given the system's data. The package allows to compute -the matrices according to different methods, providing a flexible and powerful tool. +building the network matrices used in DC/AC power flow, sensitivity, and +contingency analysis. Given a `PowerSystems.jl` [`System`](@extref PowerSystems.System), it produces the +linear-algebra layer of the [Sienna](https://www.nlr.gov/analysis/sienna.html) +power-systems platform: it reads the network data and returns matrix objects, it +does not own the data model. -`PowerNetworkMatrices.jl` is an active project under development, and we welcome your feedback, -suggestions, and bug reports. +`PowerNetworkMatrices.jl` is an active project under development, and we welcome +your feedback, suggestions, and bug reports. -## Documentation Structure +## Installation -PowerNetworkMatrices.jl strives to follow the [Diataxis documentation framework](https://diataxis.fr/), -which organizes documentation according to the different needs of users. The documentation is -structured into four main sections: - -### Tutorials - -**Learning-oriented guides to help you get started** - -Tutorials are hands-on lessons that take you through practical examples step-by-step. -They are designed to help you learn by doing, building understanding through practical experience. - - - Start here if you're new to PowerNetworkMatrices.jl - - Follow along with executable examples - - Build foundational knowledge - -### How-To Guides - -**Task-oriented guides for accomplishing specific goals** - -How-to guides provide direct instructions for solving specific problems or completing -particular tasks. They assume you have basic knowledge and want to accomplish something specific. - - - Use when you know what you want to do - - Get straight to the solution - - Focus on practical application - -### Explanation - -**Understanding-oriented discussion of key topics** - -Explanations provide background, context, and deeper understanding of concepts, design -decisions, and the theory behind the implementation. - - - Understand the "why" behind the features - - Learn about the mathematical foundations - - Explore conceptual relationships - -### Reference - -**Information-oriented technical descriptions** - -Reference documentation provides detailed, technical information about the API, functions, -and data structures. It's organized for quick lookup of specific details. - - - Look up function signatures and parameters - - Find available methods and options - - Access complete API documentation - -* * * +```text +] add PowerNetworkMatrices +``` -PowerNetworkMatrices has been developed as part of the Scalable Integrated Infrastructure Planning (SIIP) initiative at the U.S. Department of Energy's National Laboratory of the Rockies (formerly known as NREL) ([NLR](https://www.nrel.gov/)). +## Where to start + +The documentation follows the [Diátaxis](https://diataxis.fr/) framework, split +into four kinds of material: + + - **[Tutorials](tutorials/generated_introduction.md)** — learning-oriented + journeys, each answering one question end to end. Start with the + [Introduction](tutorials/generated_introduction.md), which screens a line + outage on a small network with a PTDF, an LODF, and a reduction; then + [Analysis at Scale](tutorials/generated_analysis_at_scale.md) does the same at + scale with virtual matrices and cache control. + - **[How-To Guides](how_to_guides/generated_build_multiple_matrices.md)** — + task recipes for a specific goal (build matrices efficiently, choose a solver, + reproduce industry DFAX values, define contingencies, diagnose connectivity). + - **[Reference](reference/network_matrices_overview.md)** — exhaustive + descriptions of every matrix type, accessor, and setting. Begin at the + [Matrix Overview and Indexing](reference/network_matrices_overview.md) hub. + - **[Explanation](explanation/dc_power_flow_approximation.md)** — the concepts + and trade-offs: the DC approximation, network-reduction theory, computational + considerations, concurrency, and slack conventions. + +## About + +`PowerNetworkMatrices.jl` has been developed as part of the Scalable Integrated +Infrastructure Planning (SIIP) initiative at the U.S. Department of Energy's +National Laboratory of the Rockies (formerly known as NREL) ([NLR](https://www.nlr.gov/)). diff --git a/docs/src/network_matrices.jl b/docs/src/network_matrices.jl deleted file mode 100644 index 2c6b330ee..000000000 --- a/docs/src/network_matrices.jl +++ /dev/null @@ -1,77 +0,0 @@ -# # Network Matrices - -# `PowerNetworkMatrices.jl` is able to build classic power systems modeling network matrices such as -# [Ybus](https://en.wikipedia.org/wiki/Nodal_admittance_matrix), [PTDF](https://www.powerworld.com/WebHelp/Content/MainDocumentation_HTML/Power_Transfer_Distribution_Factors.htm) and LODF. - -# Check section [Network Matrices](@ref net_mat) for more details. - -# ## Overview -# -# Network matrices are implemented in `PowerNetworkMatrices.jl` as arrays that support -# indexing by arc tuples `(from_bus_number, to_bus_number)` and bus numbers. -# The Ybus is stored as a SparseMatrix and the PTDF and LODF are stored as dense matrices. -# **Note**: Ybus is converted to a dense matrix for printing in the REPL. -# The network matrices code implements the dfs algorithm to find islands. - -using PowerSystems -import PowerSystems as PSY -DATA_DIR = "../../../data" #hide -system_data = System(joinpath(DATA_DIR, "matpower/case14.m")) - -# ## Ybus -# -# The Ybus can be calculated as follows: -ybus = Ybus(system_data) - -# The matrix can be indexed using directly the bus numbers. In this example buses are numbered -# 1-14. However, in large systems buses don't usually follow sequential numbering. You can access -# the entries of the Ybus with direct indexing or using the buses. - -ybus_entry = ybus[3, 3] -# -bus3 = get_component(Bus, system_data, "Bus 3 HV") -ybus_entry = ybus[bus3, bus3] - -# We recognize that many models require matrix operations. For those cases, you can access the -# sparse data as follows: - -sparse_array = get_data(ybus) - -# ## PTDF -# -# The PTDF matrix can be calculated as follows: -ptdf = PTDF(system_data) - -# The PTDF matrix is indexed by **arc tuples** `(from_bus, to_bus)` for branches and -# **bus numbers** for buses. An arc tuple represents a directed connection between two buses. -# For example, an arc `(1, 2)` represents a branch from bus 1 to bus 2. -# -# Elements can be accessed using arc tuples and bus numbers. The `axes` field lists the -# identifiers for each dimension and the `lookup` dictionaries map those identifiers to -# integer matrix indices: - -get_axes(ptdf) - -get_lookup(ptdf) - -# Access a PTDF entry using an arc tuple and bus number taken from the axes: - -first_arc = get_axes(ptdf)[2][1] -first_bus = get_axes(ptdf)[1][1] - -ptdf_entry = ptdf[first_arc, first_bus] - -# PTDF also takes a vector of distributed slacks, for now this feature requires passing a -# vector of weights with the same number of elements as buses in the system. For more details -# check the API entry for [`PTDF`](@ref). - -# ## LODF -# -# The LODF matrix is indexed by arc tuples for both dimensions. Each element `lodf[arc_i, arc_j]` -# represents the change in flow on `arc_i` when `arc_j` is taken out of service. - -lodf = LODF(system_data) - -# Both lookup dictionaries map arc tuples to matrix indices: - -get_lookup(lodf) diff --git a/docs/src/reference/matrix_types.md b/docs/src/reference/matrix_types.md new file mode 100644 index 000000000..6276a4de2 --- /dev/null +++ b/docs/src/reference/matrix_types.md @@ -0,0 +1,138 @@ +# Matrix type reference + +Unless noted otherwise, every type below is a subtype of the common supertype +`PowerNetworkMatrix` and shares the same indexing and accessor interface. + +## The `PowerNetworkMatrix` supertype + +```julia +abstract type PowerNetworkMatrix{T} <: AbstractArray{T, 2} end +``` + +Because it is an `AbstractArray{T,2}`, every concrete matrix supports `size`, +`axes`, and `getindex`. Indexing is overloaded so rows and columns are addressed +by *power-system identifiers* (bus numbers, arc tuples, branch names, +[`PSY`](@extref PowerSystems.System) components) rather than integer positions — +see the [overview hub](network_matrices_overview.md) for the accepted key types. +Concrete types carry a `data` field (dense `Matrix` or +[`SparseMatrixCSC`](@extref Julia SparseArrays.SparseMatrixCSC)), an `axes` tuple +of identifier vectors, a `lookup` tuple of `Dict`s, and a +[`NetworkReductionData`](@ref) describing any reduction applied at construction. + +Several types store their `data` **transposed** for efficiency +(`stores_transpose` is `true` for [`PTDF`](@ref), [`LODF`](@ref), and +[`BA_Matrix`](@ref)). Indexing hides this; use the type-specific data accessors +([`get_ptdf_data`](@ref), [`get_lodf_data`](@ref)) to obtain the standard +(non-transposed) orientation. + +## Distribution-factor matrices + +[`PTDF`](@ref) and [`LODF`](@ref) are dense distribution-factor matrices sharing +construction options (`linear_solver`, `tol`) and arc-tuple indexing. + + - **[`PTDF`](@ref)** — the Power Transfer Distribution Factor matrix. + `PTDF[arc, bus]` is the sensitivity of the flow on `arc` to a unit injection + at `bus`, under the DC approximation. + - **[`LODF`](@ref)** — the Line Outage Distribution Factor matrix. + `LODF[monitored, outaged]` is the fraction of `outaged`'s pre-outage flow that + redistributes onto `monitored`. Both dimensions are arcs; diagonal entries + are structurally `-1.0`, and are preserved when sparsifying. + +Network reductions reach both through `Ybus(...; network_reductions = [...])`. +`LODF(A, PTDF)` warns and densifies if the supplied [`PTDF`](@ref) was itself sparsified, since +that degrades LODF accuracy. Only [`PTDF`](@ref) supports HDF5 serialization (see +the [`to_hdf5`](@ref) / [`from_hdf5`](@ref) docstrings). + +## Virtual (on-demand) matrices + +Virtual matrices trade compute for memory: instead of materializing a dense +matrix they store the factorized system data and compute any single row on +demand, caching each row in an LRU `RowCache`. The cache is bounded by +`max_cache_size` (default 100 MiB) as both a byte budget and a maximum row count; +once either fills, the least-recently-used row is cleared. They expose the same +identifier-based indexing as their dense counterparts, are best for large systems +where only a subset of rows is needed, and are **not** serializable. + + - **[`VirtualPTDF`](@ref)** — lazy [`PTDF`](@ref); entries mean the same thing. + - **[`VirtualLODF`](@ref)** — lazy [`LODF`](@ref); both dimensions are arcs. + - **[`VirtualMODF`](@ref)** — the on-demand Modification (post-contingency) + Distribution Factor matrix: the [`PTDF`](@ref) row of a monitored arc *after* a + modification/contingency, via Woodbury updates. **There is no dense `MODF` + type.** Query it with the [`ContingencySpec`](@ref) / [`NetworkModification`](@ref) + types (see the [contingencies how-to](../how_to_guides/generated_contingencies.md)). + +!!! note "Distributed slack: `Dict` vs `Vector`" + + [`PTDF`](@ref)/[`VirtualPTDF`](@ref) take `dist_slack` as a + `Dict{Int, Float64}` (bus → weight); [`VirtualLODF`](@ref) and + [`VirtualMODF`](@ref) take it as a `Vector{Float64}` (one weight per bus, + ordered like the bus axis). Weights need not sum to one — they are normalized + internally. The empty default uses a single reference bus. For why distributing + the slack changes the factors and how to choose weights, see + [Slack distribution and reference-bus conventions](../explanation/slack_conventions.md). + +## Admittance and network-structure matrices + + - **[`Ybus`](@ref)** — the complex nodal admittance matrix (`YBUS_ELTYPE` is + `ComplexF64`). `Ybus[i, j]` is the mutual admittance between buses `i` and `j` + (off-diagonal) or the self-admittance of `i` (diagonal). It is the foundation + matrix — every DC and virtual matrix is built from it — and stays complex, so + it is factorized with KLU rather than the real-only backends. Asymmetry is + legitimate for phase-shifting transformers and must not be "corrected". Pass + `make_arc_admittance_matrices = true` to also build the two + [`ArcAdmittanceMatrix`](@ref) objects. + - **[`ArcAdmittanceMatrix`](@ref)** — per-arc admittance in one direction + (`:FromTo` or `:ToFrom`), for power-flow use. Not built standalone by typical + users; it is produced inside [`Ybus`](@ref). Query the direction with + `get_direction`. + - **[`IncidenceMatrix`](@ref)** — the bus-branch incidence matrix `A`. + `A[arc, bus]` is `+1` at the from-bus, `-1` at the to-bus, `0` otherwise + (exactly two nonzeros per arc row). Structural building block for the DC + matrices. + - **[`BA_Matrix`](@ref)** — the susceptance-weighted incidence matrix `B · A`, + where `B` is the diagonal matrix of branch series susceptances (`b = 1/x` under + the DC approximation) and `A` the [`IncidenceMatrix`](@ref). Axes match + [`IncidenceMatrix`](@ref); stored transposed. The reference-bus column is dropped + (one fewer column than the bus count). + - **[`ABA_Matrix`](@ref)** — the reduced bus-susceptance matrix `Aᵀ · B · A` + (`A` the [`IncidenceMatrix`](@ref), `B` the branch susceptance — see + [`BA_Matrix`](@ref)) with reference buses removed for invertibility — the + DC-power-flow system matrix. + Its `K` field optionally holds a KLU factorization: build it factorized + (`factorize = true`), or call [`factorize`](@ref) afterward (it returns a fresh + factorized copy; [`is_factorized`](@ref) checks). The stored factorization is + consumed by the low-level `LODF(A::IncidenceMatrix, ABA::ABA_Matrix, BA::BA_Matrix)` + constructor — it reads `ABA.K` directly, so passing a pre-factorized ABA there + avoids re-factorizing. The `PTDF` constructors do not take an `ABA_Matrix` and + factorize internally, so there is no factorization to hand them. + - **[`AdjacencyMatrix`](@ref)** — a symmetric bus-by-bus connectivity matrix + (`Int8`): nonzero where two buses share a branch, zero on the diagonal. Used by + [`validate_connectivity`](@ref) and [`find_subnetworks`](@ref). + +## Concrete type aliases + +`PowerflowMatrixTypes.jl` defines concrete aliases for the fully-parameterized +matrix types that appear in hot paths and downstream dispatch. Use them in method +signatures when you need to fix a concrete storage layout. All are exported, as is +`YBUS_ELTYPE`. + +| Alias | Underlying type | Meaning | +|:---------------------------- |:------------------------------------------------- |:--------------------------------------------------- | +| `DC_PTDF_Matrix` | `PTDF{…, Matrix{Float64}}` | Dense PTDF with bus/arc axes. | +| `DC_vPTDF_Matrix` | `VirtualPTDF{…, K} where {K}` | Virtual PTDF (factorization type `K` left free). | +| `DC_BA_Matrix` | `BA_Matrix{…}` | BA matrix with bus/arc axes. | +| `DC_ABA_Matrix_Factorized` | `ABA_Matrix{…, KLULinSolveCache{Float64, Int64}}` | ABA matrix carrying a KLU factorization. | +| `DC_ABA_Matrix_Unfactorized` | `ABA_Matrix{…, Nothing}` | ABA matrix with no factorization (`K === nothing`). | +| `AC_Ybus_Matrix` | `Ybus{…}` | Ybus with integer bus axes. | + +## See also + + - [Matrix overview and indexing hub](network_matrices_overview.md) — how + `A[row, col]` resolves and the per-type axis summary. + - [Full public API](public.md) — authoritative docstrings and signatures. + - [`NetworkReduction`](@ref) — the reduction spec types passed via + `network_reductions`, with the ordering rules. + - [How to Define and Apply Contingencies](@ref) — types used to query + [`VirtualMODF`](@ref). + - [`AutoTolerance`](@ref) and [How to Choose a Linear Solver](@ref) — the + sparsification `tol` and `linear_solver` backends. diff --git a/docs/src/reference/network_matrices_overview.md b/docs/src/reference/network_matrices_overview.md index b04d560a7..e0d5f0a39 100644 --- a/docs/src/reference/network_matrices_overview.md +++ b/docs/src/reference/network_matrices_overview.md @@ -1,143 +1,303 @@ -## Core Matrix Types - -### Incidence Matrix ([`IncidenceMatrix`](@ref)) - -The incidence matrix $A$ represents the bus-arc connectivity of the network, showing which arcs connect to which buses. It's a fundamental building block for other matrices. - -**Properties:** - - - Rectangular matrix (arcs × buses) - - Values indicate connection direction (+1 for "from" bus, -1 for "to" bus, or 0) - - Rows are indexed by arc tuples `(from_bus_number, to_bus_number)` - - Columns are indexed by bus numbers - -**Purpose:** -The incidence matrix captures pure network topology without electrical parameters, making it useful for graph-theoretic analysis. - -### Adjacency Matrix ([`AdjacencyMatrix`](@ref)) - -The adjacency matrix represents the directed connectivity between buses, showing which buses are connected to each other. It's a fundamental building block for other matrices. - -**Properties:** - - - Square matrix (buses × buses) - - Values indicate connection direction (+1 for "from" bus, -1 for "to" bus, or 0) - - Sparsity pattern matches the bus admittance matrix. - -**Purpose:** -The adjacency matrix captures pure network topology without electrical parameters, making it useful for graph-theoretic analysis. - -### Bus Admittance Matrix ([`Ybus`](@ref)) - -The bus admittance matrix combines network topology with electrical parameters (impedance/admittance). - -**Properties:** - - - Square matrix (buses × buses) - - Diagonal elements are self-admittances - - Off-diagonal elements are mutual admittances - - Symmetric for passive networks - - Both dimensions are indexed by bus numbers - -**Purpose:** -Forms the basis for power flow calculations and relates bus voltages to injected currents. - -**Note:** -The [`ArcAdmittanceMatrix`](@ref) is a rectangular matrix (arcs x buses) consisting of the off diagonal entries of the [`Ybus`](@ref). These matrices are useful when computing line flows from bus voltages and can be optionally created when building the [`Ybus`](@ref). - -### BA Matrix ([`BA_Matrix`](@ref)) - -The BA matrix represents the arc-bus incidence matrix weighted by arc susceptances. - -**Purpose:** -Serves as an intermediate calculation in deriving ABA, PTDF and LODF matrices. - -### ABA Matrix ([`ABA_Matrix`](@ref)) - -The ABA matrix is derived from $A \cdot B \cdot A^T$ where $A$ is the incidence matrix and $B$ contains branch admittances. - -**Purpose:** -Serves as an intermediate calculation in deriving PTDF and LODF matrices. - -## Sensitivity Matrices - -### Power Transfer Distribution Factors ([`PTDF`](@ref)) - -PTDF matrices answer the question: "If I inject 1 MW at bus $i$ and withdraw 1 MW at bus $j$, how much does the flow on branch $k$ change?" - -**Key Characteristics:** - - - Linearized approximation of power flow - - Valid for small perturbations around operating point - - Fast to compute and evaluate - - Widely used in market operations and security analysis - - Rows are indexed by arc tuples `(from_bus, to_bus)`, columns by bus numbers - -**Note:** -See [`VirtualPTDF`](@ref) for cases where it is not possible to compute or store the full PTDF matrix. - -### Line Outage Distribution Factors ([`LODF`](@ref)) - -LODF matrices answer: "If branch $m$ fails, how much does the flow redistribute to branch $k$?" - -**Key Characteristics:** - - - Predicts post-contingency flows - - Essential for N-1 security analysis - - Computed from PTDF matrix - - Helps identify critical lines - - Both dimensions are indexed by arc tuples `(from_bus, to_bus)` - -**Mathematical Relationship:** -LODF is derived from PTDF using matrix operations that model the effect of removing a branch from the network. - -See [`VirtualLODF`](@ref) for cases where it is not possible to compute or store the full LODF matrix. - -### Modification Distribution Factors ([`VirtualMODF`](@ref)) - -The MODF answers: "Under a registered contingency or network modification (one or more branch outages, partial-susceptance changes, or shunt changes), what is the resulting PTDF row for a given monitored arc?" +# Matrix overview & indexing + +This page summarizes every matrix type, shows the shared +construction and indexing pattern, documents how element indexing (`A[row, column]`) +resolves, and documents the accessor functions used to read data, axes, lookups, +reference buses, reduction data, and provenance + +All matrix types are concrete subtypes of the abstract supertype +`PowerNetworkMatrix{T} <: AbstractArray{T, 2}` (`src/PowerNetworkMatrix.jl`). Because +they are `AbstractArray{T,2}` subtypes, they support the standard array interface +(`size`, `axes`, `getindex`, iteration), but indexing is overloaded so that rows and +columns are addressed by domain identifiers (bus numbers and arc tuples) rather than +by raw integer positions. `LinearIndices` and `CartesianIndices` are intentionally +disabled. + +## Matrix taxonomy + +Every matrix stores two structural fields: + + - **`axes`**: a 2-tuple of vectors listing the identifiers (arc tuples and/or + bus numbers) for each dimension. + - **`lookup`**: a 2-tuple of dictionaries mapping those identifiers to integer + positions into the stored `data`. + +The storage form is one of three kinds: + + - **Dense** — the full matrix is materialized as a `Matrix{Float64}` (or, for + `Ybus`, a `SparseMatrixCSC`). + - **Sparse** — stored as a `SparseMatrixCSC`. + - **Virtual / lazy** — no full matrix is stored; rows are computed on demand + and cached in an LRU row cache. Use these past large-system limits instead of + the dense forms. Virtual matrices are not serializable. + +| Matrix | Rows | Columns | Storage | Represents | +|:----------------------------- |:----------- |:----------- |:---------------- |:----------------------------------------------------------- | +| [`IncidenceMatrix`](@ref) | arc tuples | bus numbers | sparse | signed bus–arc topology (`+1` from-bus, `-1` to-bus) | +| [`AdjacencyMatrix`](@ref) | bus numbers | bus numbers | sparse | signed bus–bus connectivity (Ybus sparsity pattern) | +| [`Ybus`](@ref) | bus numbers | bus numbers | sparse (complex) | nodal admittance (topology + electrical parameters) | +| [`ArcAdmittanceMatrix`](@ref) | arc tuples | bus numbers | sparse (complex) | off-diagonal Ybus entries; built as part of `Ybus` | +| [`BA_Matrix`](@ref) | bus numbers | arc tuples | sparse | ``B A``: incidence weighted by branch susceptance | +| [`ABA_Matrix`](@ref) | bus numbers | bus numbers | sparse | ``A^\top B A`` DC susceptance matrix; optionally factorized | +| [`PTDF`](@ref) | arc tuples | bus numbers | dense | power transfer distribution factors | +| [`LODF`](@ref) | arc tuples | arc tuples | dense | line outage distribution factors | +| [`VirtualPTDF`](@ref) | arc tuples | bus numbers | virtual | lazy per-row PTDF | +| [`VirtualLODF`](@ref) | arc tuples | arc tuples | virtual | lazy per-row LODF | +| [`VirtualMODF`](@ref) | arc tuples | bus numbers | virtual | post-modification / post-contingency PTDF rows | + +Notes on the taxonomy: + + - `Ybus` is complex-valued (`ComplexF64`); all other numeric matrices are + real (`Float64`). `IncidenceMatrix` / `AdjacencyMatrix` store signed `Int8` + topology. + - `PTDF` and `LODF` store their data **transposed** internally; `getindex` + and [`get_ptdf_data`](@ref) / [`get_lodf_data`](@ref) hide + this so callers always see the standard `(row, column)` orientation. + - There is **no dense `MODF` type** — only [`VirtualMODF`](@ref). Contrast + with PTDF/LODF, which have both dense and virtual forms. + - `ArcAdmittanceMatrix` is produced as a byproduct of building `Ybus` (via a + construction keyword) rather than being independently constructed by typical + users. + +Full constructor signatures, keyword arguments, and concrete type aliases +(`DC_PTDF_Matrix`, `DC_ABA_Matrix_Factorized`, `AC_Ybus_Matrix`, …) are on the +[Matrix type reference](matrix_types.md). + +## Constructing matrices + +Every matrix type is a constructor that takes the +[`System`](@extref PowerSystems.System) and returns the matrix object. The call is +identical across types — only the name changes: + +```julia +using PowerNetworkMatrices +import PowerSystemCaseBuilder as PSB + +sys = PSB.build_system(PSB.PSITestSystems, "c_sys5") + +ptdf = PTDF(sys) +lodf = LODF(sys) +ybus = Ybus(sys) +aba = ABA_Matrix(sys) +``` + +The shared build-time keywords — `network_reductions`, `tol`, `linear_solver`, +`dist_slack` — work on every constructor that accepts them; each has its own how-to. +The lazy [`VirtualPTDF`](@ref) / [`VirtualLODF`](@ref) / [`VirtualMODF`](@ref) forms +build and index exactly like their materialized counterparts — swap the type name; +they compute rows on demand and cache them instead of storing the whole matrix. + +Some constructors also accept **already-built matrices** instead of a +[`System`](@extref PowerSystems.System), so shared intermediates (`Ybus`, incidence, +BA) are computed once and reused. See +[How to Build Multiple Matrices Without Repeating Work](@ref). + +## Arc-tuple indexing + +Matrices that involve branches identify each branch by an **arc tuple** — a +`Tuple{Int, Int}` of the form `(from_bus_number, to_bus_number)` giving the +directed connection between two buses. Arc tuples, rather than branch-name +strings, are the canonical branch identifier because they: + + - identify a network element compactly and unambiguously; + - survive network reductions, where named branches may be merged or eliminated + but the surviving equivalent arc keeps a well-defined endpoint pair; + - match the mathematical formulation, in which a branch is defined by its two + endpoint buses. + +## How `A[row, column]` resolves + +Indexing is fully overloaded on `PowerNetworkMatrix` (`src/PowerNetworkMatrix.jl`). +`A[row, column]` calls `to_index(A, row, column)`, which maps each supplied +identifier to an integer position through the per-dimension `lookup` dictionary +(via the internal `lookup_index` helper), then reads the underlying `data`. + +The accepted element types for `row` and `column`, and how each resolves: + +| Index value | Resolves via | Supported | +|:-------------------------------------- |:------------------------------------------------------------------------------------------------------------------------------------- |:------------------------------ | +| `Int` (bus number) | direct `lookup[i]` | ✅ | +| arc tuple `(from, to)::Tuple{Int,Int}` | direct `lookup[i]` | ✅ | +| `PSY.ACBus` | `lookup_index` specialization → `Base.to_index(bus) = get_number(bus)` | ✅ | +| `PSY.Arc` | `lookup_index` specialization → `Base.to_index(arc) = get_arc_tuple(arc)` | ✅ | +| branch-name `String` | dedicated `getindex` on `PTDF` / `LODF` / `VirtualPTDF` (maps name → arc via reduction data, applies parallel/aggregation multiplier) | ✅ (PTDF/LODF/VirtualPTDF only) | +| `Colon` (`:`) | returns the whole row/column | ✅ | +| `PowerNetworkMatrixKey` | `A[k]` splats `k.I` back into `A[k.I...]` | ✅ | +| raw `Int` position pair | dense positional fast path (`A.data[…]`) | ✅ | +| `PSY.ACBranch` | — | ❌ raises `KeyError` | + +!!! warning "Branch objects are not directly indexable" + + A [`ACBranch`](@extref PowerSystems.ACBranch) component **cannot** be passed as + an index — doing so raises a `KeyError`. Although `Base.to_index(::PSY.ACBranch)` + is defined (returning the branch's arc tuple), the matrix `getindex` path routes + only [`ACBus`](@extref PowerSystems.ACBus) and [`Arc`](@extref PowerSystems.Arc) + through `Base.to_index`; branch components are not converted. Index a branch by + its **arc tuple** (the branch's [`Arc`](@extref PowerSystems.Arc)), or, for + [`PTDF`](@ref)/[`LODF`](@ref)/[`VirtualPTDF`](@ref), by its **name string**. -**Key Characteristics:** +Only `PTDF`, `LODF`, and `VirtualPTDF` accept branch-name `String` indices; +`LODF` requires a `String` for both dimensions when using names. Name indexing +maps the name to an arc tuple through the network reduction data and multiplies +by the appropriate parallel/aggregation factor, so it is retained for backward +compatibility but is slower and less direct than arc-tuple indexing. - - On-demand computation of post-modification PTDF rows via the Woodbury matrix identity over a base PTDF. - - Caches Woodbury factors per modification and post-modification PTDF rows per `(monitored_arc, modification)` pair. - - Auto-registers `PSY.Outage` supplemental attributes attached to the source system. - - Rows are indexed by monitored arc tuples or arc indices; the second index is a [`ContingencySpec`](@ref), a [`NetworkModification`](@ref), or a `PSY.Outage`. +### Examples -**Use it when:** you need post-contingency or post-modification flow sensitivities for one or a few contingencies on a large system, and constructing the full LODF or recomputing PTDF per contingency is too expensive. +```julia +using PowerNetworkMatrices +import PowerSystems as PSY +import PowerSystemCaseBuilder as PSB -## Arc-Based Indexing +sys = PSB.build_system(PSB.PSITestSystems, "c_sys5") +ptdf = PTDF(sys) -All matrices that involve branches use **arc tuples** as identifiers instead of branch name strings. An arc tuple is a `Tuple{Int, Int}` of the form `(from_bus_number, to_bus_number)`, representing the directed connection between two buses. +# By bus number and arc tuple (canonical): +ptdf[(2, 3), 1] -This indexing approach: +# By PSY component objects: +bus1 = first(b for b in PSY.get_components(PSY.ACBus, sys) if PSY.get_number(b) == 1) +branch = first(PSY.get_components(PSY.ACBranch, sys)) +ptdf[PSY.get_arc(branch), bus1] # PSY.Arc row, PSY.ACBus column - - Provides compact, unambiguous identification of network elements - - Naturally handles network reductions where branches may be merged or eliminated - - Aligns with the mathematical formulation where branches are defined by their endpoint buses +# By branch name (PTDF/LODF/VirtualPTDF only): +ptdf[PSY.get_name(branch), 1] -Each matrix stores `axes` and `lookup` fields: +# NOT allowed — raises KeyError: +# ptdf[branch, 1] # a PSY.ACBranch object - - **`axes`**: A tuple of vectors listing the arc tuples and/or bus numbers for each dimension - - **`lookup`**: A tuple of dictionaries that map arc tuples or bus numbers to integer matrix indices +# Whole row / column with a Colon: +ptdf[:, 1] # column for bus 1 +``` -**Indexing summary by matrix type:** +The `Ybus` accepts bus numbers or `PSY.ACBus` objects on both dimensions: -| Matrix | Dimension 1 (rows) | Dimension 2 (columns) | -|:----------------- |:------------------ |:--------------------- | -| `IncidenceMatrix` | Arc tuples | Bus numbers | -| `BA_Matrix` | Arc tuples | Bus numbers | -| `PTDF` | Arc tuples | Bus numbers | -| `LODF` | Arc tuples | Arc tuples | -| `Ybus` | Bus numbers | Bus numbers | -| `VirtualPTDF` | Arc tuples | Bus numbers | -| `VirtualLODF` | Arc tuples | Arc tuples | -| `VirtualMODF` | Arc tuples | Bus numbers | +```julia +ybus = Ybus(sys) +ybus[3, 3] +ybus[PSY.get_number(bus1), PSY.get_number(bus1)] +``` -!!! note +!!! note "Reduced arcs are not indexable" - For backward compatibility, branch name strings can also be used to index PTDF and LODF matrices. This internally maps the branch name to its corresponding arc tuple via the network reduction data. Using arc tuples directly is recommended for new code. - -!!! note + When network reductions (e.g. `RadialReduction`, `DegreeTwoReduction`) are + applied, eliminated branches are absent from the matrix. Indexing with an arc + tuple that was reduced away raises an error. Inspect the surviving + identifiers with `PNM.get_axes(A)` (see + [Accessors: axes, lookups, and data](@ref)). + +## Accessors: axes, lookups, and data + +The functions below read structural and numeric data from any matrix — the backing +array, axes, lookup dictionaries, reference buses, reduction data, and system +provenance. Exported accessors are documented in full on the +[Full public API](public.md); the internal helpers are reached through the module +prefix (commonly aliased `PNM.`) and are documented here. + + - **Exported:** [`get_ptdf_data`](@ref), [`get_lodf_data`](@ref), + [`get_partial_lodf_row`](@ref), [`get_network_reduction_data`](@ref), + [`get_system_uuid`](@ref). + - **Internal (call via `PNM.`):** `get_data`, `get_axes`, `get_lookup`, + `get_bus_axis`, `get_arc_axis`, `get_bus_lookup`, `get_arc_lookup`, + `get_ref_bus`, `get_ref_bus_position`. + +### Data extraction + + - **`PNM.get_data(mat)`** — the raw backing array (`mat.data`) exactly as stored: + a complex [`SparseMatrixCSC`](@extref Julia SparseArrays.SparseMatrixCSC) for + [`Ybus`](@ref), the internally **transposed** dense matrix for + [`PTDF`](@ref)/[`LODF`](@ref). + - **[`get_ptdf_data`](@ref)** / **[`get_lodf_data`](@ref)** — the matrix in + standard (non-transposed) orientation, via a lazy `transpose` (not a copy). For + a [`VirtualLODF`](@ref), [`get_lodf_data`](@ref) instead returns the LRU cache + contents as a `Dict{Int, Vector{Float64}}` of already-computed rows. + - **[`get_partial_lodf_row`](@ref)** — the LODF row for a **partial** susceptance + change `delta_b` on one arc (full outage: `delta_b = -arc_susceptance`). The + entry point for partial outages that a plain [`VirtualLODF`](@ref) row (which + assumes a full outage) does not cover. + +```julia +import PowerNetworkMatrices as PNM + +PNM.get_axes(ptdf) # (bus-number vector, arc-tuple vector) +PNM.get_lookup(ptdf) # (bus lookup Dict, arc lookup Dict) +PNM.get_data(ybus) # raw SparseMatrixCSC +``` + +### Axes and lookups + +`PNM.get_axes(mat)` returns `mat.axes` and `PNM.get_lookup(mat)` returns +`mat.lookup`, each a 2-tuple ordered `(dimension 1, dimension 2)`. The axis vector +lists identifiers (bus numbers as `Int`, arcs as `Tuple{Int,Int}`) in position +order; the matching lookup maps each identifier back to its integer position in +`data`. These are the authoritative way to enumerate valid indices — especially +after a reduction, where some arcs/buses are no longer present. Defined for every +matrix type. + +The dimension-specific accessors `PNM.get_bus_axis` / `get_arc_axis` / +`get_bus_lookup` / `get_arc_lookup` select the correct dimension without the caller +knowing which index (1 or 2) is the bus or arc dimension for a given matrix type. +They are defined only for the dimensions a matrix actually has: + +| Matrix | `get_bus_axis` | `get_arc_axis` | +|:--------------------- |:-----------------:|:-----------------:| +| `IncidenceMatrix` | `axes[2]` | `axes[1]` | +| `AdjacencyMatrix` | `axes[1]` | — (both dims bus) | +| `Ybus` | `axes[1]` | — (both dims bus) | +| `ArcAdmittanceMatrix` | `axes[2]` | `axes[1]` | +| `BA_Matrix` | `axes[1]` | `axes[2]` | +| `ABA_Matrix` | `axes[1]` | — (both dims bus) | +| `PTDF` | `axes[1]` | `axes[2]` | +| `LODF` | — (both dims arc) | `axes[1]` | +| `VirtualPTDF` | `axes[1]` | `axes[2]` | +| `VirtualLODF` | — (both dims arc) | `axes[1]` | +| `VirtualMODF` | `axes[2]` | `axes[1]` | + +### Reference buses + +`PNM.get_ref_bus(mat)` returns the sorted reference (slack) bus numbers — one per +electrical island — and `PNM.get_ref_bus_position(mat)` their integer positions in +the bus dimension. Together they identify the slack bus(es) held fixed when the +matrix was built, which matters for interpreting [`PTDF`](@ref) columns and for +reduction/contingency math. Defined for the distribution-factor, incidence, +adjacency, BA/ABA, and arc-admittance matrices. + +### Reduction data + +**[`get_network_reduction_data`](@ref)** returns the [`NetworkReductionData`](@ref) +for the matrix — which buses/arcs were merged or eliminated and how they map back +(empty when no reduction was applied). This is the object queried by the reduction +accessors (`get_bus_reduction_map`, `get_removed_buses`, `get_reductions`, …); its +fields and accessors are documented on the [`NetworkReductionData`](@ref) docstring. + +!!! note "Serialization drops reduction data" - When network reductions are applied (e.g. `RadialReduction`, `DegreeTwoReduction`), some branches are eliminated from the network. Attempting to index a matrix with an arc tuple that was reduced will result in an error. Use `get_axes` to inspect the available arc tuples after reduction. + A [`PTDF`](@ref) loaded via [`from_hdf5`](@ref) carries an + **empty** [`NetworkReductionData`](@ref). See the [`to_hdf5`](@ref) docstring. + +### System provenance + +**[`get_system_uuid`](@ref)** returns the UUID of the +[`System`](@extref PowerSystems.System) the matrix was built from, or `nothing` for +types that do not track origin. [`VirtualPTDF`](@ref) and [`VirtualMODF`](@ref) +store it; it backs the consistency check that a matrix and a system passed together +share a source. + +## Reference map + +This overview is the entry point. Detailed reference lives on the sibling pages: + + - [Matrix type reference](matrix_types.md) — constructor signatures, keyword + arguments, and concrete type aliases for every matrix type. + - [How to Diagnose a Disconnected Network](@ref) — testing whether the network is + connected and enumerating electrical islands. + - [How to Define and Apply Contingencies](@ref) — `ArcModification`, + `ShuntModification`, `NetworkModification`, `ContingencySpec`, and the Woodbury + tooling. + - [`AutoTolerance`](@ref) and [How to Choose a Linear Solver](@ref) — the + sparsification `tol` and the linear-solver backends. + - [`to_hdf5`](@ref) / [`from_hdf5`](@ref) — HDF5 persistence (PTDF only). + - [Full public API](public.md) — the curated autodocs for every exported + symbol, including the reduction specs, `NetworkReductionData`, and the + aggregated-branch rating functions. + - [Internals](internals.md) — the KLU and Accelerate solver submodules. diff --git a/docs/src/reference/public.md b/docs/src/reference/public.md index cb5cbc1f2..f80c8b804 100644 --- a/docs/src/reference/public.md +++ b/docs/src/reference/public.md @@ -1,6 +1,170 @@ # Public API Reference +```@meta +CurrentModule = PowerNetworkMatrices +``` + +## Matrix types + +Dense, virtual, and structural network matrices, plus the concrete type aliases. +See [Matrix Overview](network_matrices_overview.md) for the taxonomy and indexing +rules. + +```@autodocs +Modules = [PowerNetworkMatrices] +Public = true +Private = false +Order = [:type, :constant, :function, :macro] +Filter = t -> t in ( + PTDF, LODF, VirtualPTDF, VirtualLODF, VirtualMODF, + Ybus, ArcAdmittanceMatrix, BA_Matrix, ABA_Matrix, + IncidenceMatrix, AdjacencyMatrix, + DC_PTDF_Matrix, DC_vPTDF_Matrix, DC_BA_Matrix, + DC_ABA_Matrix_Factorized, DC_ABA_Matrix_Unfactorized, + AC_Ybus_Matrix, YBUS_ELTYPE, +) +``` + +## Network reductions + +Reduction specifications, the reduction-data container, and its accessors. See the +[`NetworkReduction`](@ref) docstring for the `network_reductions` keyword and rules. + +```@autodocs +Modules = [PowerNetworkMatrices] +Public = true +Private = false +Order = [:type, :constant, :function, :macro] +Filter = t -> t in ( + NetworkReduction, RadialReduction, DegreeTwoReduction, WardReduction, + NetworkReductionData, + get_bus_reduction_map, get_network_reduction_data, + get_reductions, get_ward_reduction, +) +``` + +## Aggregated-branch ratings + +Rating-aggregation strategies for equivalent branches produced by network +reduction. + +```@autodocs +Modules = [PowerNetworkMatrices] +Public = true +Private = false +Order = [:type, :constant, :function, :macro] +Filter = t -> t in ( + get_sum_of_max_rating, + get_single_element_contingency_rating, + get_impedance_averaged_rating, +) +``` + +## Contingencies & modifications + +Modification and contingency specification types, Ybus-delta application, and the +Woodbury-based post-contingency PTDF update. See the +[contingencies how-to](../how_to_guides/generated_contingencies.md) for the type model +and worked examples. + +```@autodocs +Modules = [PowerNetworkMatrices] +Public = true +Private = false +Order = [:type, :constant, :function, :macro] +Filter = t -> t in ( + ArcModification, ShuntModification, NetworkModification, ContingencySpec, + apply_ybus_modification, compute_ybus_delta, + apply_woodbury_correction, compute_woodbury_factors, + get_post_modification_ptdf_row, get_registered_contingencies, +) +``` + +## Solvers & tolerance + +Sparsification tolerance and factorization controls. See +[How to Choose a Linear Solver](../how_to_guides/generated_choose_linear_solver.md) for +the solver backends. + +```@autodocs +Modules = [PowerNetworkMatrices] +Public = true +Private = false +Order = [:type, :constant, :function, :macro] +Filter = t -> t in ( + AutoTolerance, discover_data_precision, + factorize, is_factorized, +) +``` + +## Serialization + +HDF5 persistence for `PTDF` (only). + ```@autodocs Modules = [PowerNetworkMatrices] Public = true +Private = false +Order = [:type, :constant, :function, :macro] +Filter = t -> t in ( + to_hdf5, from_hdf5, get_system_uuid, +) +``` + +## Data accessors + +Extract the underlying numeric data from computed matrices. + +```@autodocs +Modules = [PowerNetworkMatrices] +Public = true +Private = false +Order = [:type, :constant, :function, :macro] +Filter = t -> t in ( + get_ptdf_data, get_lodf_data, get_partial_lodf_row, +) +``` + +## Cache management + +Control and reset the virtual-matrix row caches. + +```@autodocs +Modules = [PowerNetworkMatrices] +Public = true +Private = false +Order = [:type, :constant, :function, :macro] +Filter = t -> t in ( + clear_caches!, clear_all_caches!, +) +``` + +## Connectivity + +Island detection and connectivity validation. + +```@autodocs +Modules = [PowerNetworkMatrices] +Public = true +Private = false +Order = [:type, :constant, :function, :macro] +Filter = t -> t in ( + find_subnetworks, validate_connectivity, + depth_first_search, iterative_union_find, +) +``` + +## Internal (non-exported) symbols + +The symbols below are **not** exported and are **not** part of the supported API. +They are documented only so the manual covers every docstring shipped in the main +module; they may change without notice. Internal solver submodules +(`KLUWrapper`, `AccelerateWrapper`) are documented on the +[Internals](internals.md) page. + +```@autodocs +Modules = [PowerNetworkMatrices] +Public = false +Private = true +Order = [:type, :constant, :function, :macro] ``` diff --git a/docs/src/tutorials/analysis_at_scale.jl b/docs/src/tutorials/analysis_at_scale.jl new file mode 100644 index 000000000..29cec5571 --- /dev/null +++ b/docs/src/tutorials/analysis_at_scale.jl @@ -0,0 +1,137 @@ +# # Analysis at Scale + +# > Across every single-line outage, which contingencies drive a surviving line closest to — or past — its limit? +# +# A dense [`PTDF`](@ref) or [`LODF`](@ref) is an `O(N²)` array of `Float64`, which on a real interconnection is tens of gigabytes, and a screen reads each row exactly once. This tutorial makes use of **virtual** matrices that compute rows on demand. +# +# !!! note +# We use a small network (73 buses) for demonstration purposes. + +using PowerNetworkMatrices +import PowerNetworkMatrices as PNM +import PowerSystems as PSY +import PowerSystemCaseBuilder as PSB +using LinearAlgebra: dot +using DataFrames + +sys = PSB.build_system(PSB.PSISystems, "RTS_GMLC_DA_sys"); + +# ## Step 1 — Set up the study: base-case flows and limits + +# A post-contingency flow requires the flow each line carries at present and the limit it must stay under. + +# **Base-case flows.** Under the DC approximation a line's flow is its [`PTDF`](@ref) row dotted with the vector of net bus injections, that is, generation minus load. We build the injection vector from the system, accumulating per bus because a bus can host several generators and loads, ordered to match the matrix's bus axis. + +vptdf = VirtualPTDF(sys) +bus_lookup = PNM.get_bus_lookup(vptdf) + +injection = zeros(Float64, length(bus_lookup)) +for gen in PSY.get_components( + d -> !isa(d, Union{PSY.ElectricLoad, PSY.SynchronousCondenser}), + PSY.StaticInjection, sys) + PSY.get_available(gen) || continue + injection[bus_lookup[PSY.get_number(PSY.get_bus(gen))]] += PSY.get_active_power(gen) +end +for load in PSY.get_components(d -> !isa(d, PSY.FixedAdmittance), PSY.ElectricLoad, sys) + PSY.get_available(load) || continue + injection[bus_lookup[PSY.get_number(PSY.get_bus(load))]] -= PSY.get_active_power(load) +end + +# The system defaults to its per-unit *system base*, so these injections and the ratings read below are already on the same `100`-MVA base and are directly comparable. The injections do not need to sum to zero because the reference bus balances the difference. + +# The base flow on every line follows. Each `vptdf[arc, :]` computes that line's [`PTDF`](@ref) row on first access and caches it, so one pass touches each row once. + +arcs = vptdf.axes[1] +base_flow = Dict(arc => dot(vptdf[arc, :], injection) for arc in arcs); + +# !!! note +# Building *every* base flow this way touches the whole [`PTDF`](@ref), one row at a +# time. When only the base flows are needed, the sparse `ABA`/`BA` DC solve +# ([`ABA_Matrix`](@ref), [`BA_Matrix`](@ref)) obtains them in a single factorization; +# the row-at-a-time route here is what the N-1 screen below requires anyway. + +# **Line limits.** Each branch's rating is read and keyed by arc. Parallel branches share an arc, so their ratings are summed into the combined corridor limit. + +line_rating = Dict{Tuple{Int, Int}, Float64}() +for branch in PSY.get_components(PSY.ACTransmission, sys) + arc = PSY.get_arc(branch) + key = (PSY.get_number(PSY.get_from(arc)), PSY.get_number(PSY.get_to(arc))) + line_rating[key] = get(line_rating, key, 0.0) + PSY.get_rating(branch) +end + +# ## Step 2 — Screen every contingency on a VirtualLODF + +# The [`LODF`](@ref) gives the redistribution: when `outaged` trips, line `monitored` picks up `LODF[monitored, outaged]` of the outaged line's pre-trip flow. The post-contingency flow on `monitored` is therefore `base_flow[monitored] + LODF[monitored, outaged] · base_flow[outaged]`, and its **loading** is that quantity over its rating. + +# [`VirtualLODF`](@ref) is constructed and indexed exactly like the dense [`LODF`](@ref), with the same constructor and the same `[monitored, outaged]` indexing, but never forms the whole matrix. `max_cache_size` caps the row cache in MiB. + +vlodf = VirtualLODF(sys; max_cache_size = 100) + +# A **row**, meaning one monitored line's factors against every outage, is the unit the cache stores, so the sweep proceeds row by row: compute each monitored line's row once, then score it against every outage. For each outage we keep the single worst-loaded survivor. + +lines = vlodf.axes[1] +outage_col = vlodf.lookup[2] + +worst = Dict{Tuple{Int, Int}, @NamedTuple{line::Tuple{Int, Int}, loading::Float64}}() +for monitored in lines + factors = vlodf[monitored, :] + rating = line_rating[monitored] + f_m = base_flow[monitored] + for outaged in lines + outaged == monitored && continue + post = f_m + factors[outage_col[outaged]] * base_flow[outaged] + loading = abs(post) / rating + if !haskey(worst, outaged) || loading > worst[outaged].loading + worst[outaged] = (line = monitored, loading = loading) + end + end +end + +# Ranking the outages by the worst loading they cause gives the screen's result. + +screen = sort( + DataFrame(; + outaged_line = [o for o in keys(worst)], + worst_monitored = [w.line for w in values(worst)], + loading = [round(w.loading; digits = 3) for w in values(worst)], + ), + :loading; rev = true, +) + +# The day-ahead schedule is **not** N-1 secure: most single-line outages here push some surviving line over its rating. + +count(>(1.0), screen.loading) + +# The worst contingency drives a line to roughly `1.3×` its limit. The top pair is reciprocal — tripping `(107, 108)` overloads `(107, 203)` and vice versa — because they form a tightly coupled corridor in which each inherits essentially the entire flow of the other. The numbers are identical to those of a dense [`LODF`](@ref); the difference is that no dense matrix was ever stored. + +# The sweep visited every row, so all of them are now cached. + +length(vlodf.cache) + +# On RTS-GMLC that is approximately 108 short rows, a negligible memory footprint, which is why nothing was evicted. At realistic scale the situation is the opposite: the rows do not all fit, and `max_cache_size` is a hard ceiling. Once the cache is full the **least-recently-used** row is dropped. A full screen still completes, trading a bounded memory footprint for the recomputation of an evicted row should the screen return to it. That trade is what allows an N-1 screen to run on a grid whose dense [`LODF`](@ref) would not fit. + +# ## Step 3 — Reuse rows across operating points + +# This screen is not run once. It reruns at every operating point, but the [`LODF`](@ref) is a property of the **topology** rather than the dispatch: the factors do not change from hour to hour, only the base flows they multiply. A study that reruns the screen therefore wants its rows to remain resident rather than be recomputed on each pass. + +# In practice a defined set of facilities is monitored every cycle, here the inter-area tie corridors. Declaring them up front as `persistent_arcs` holds those rows in the cache and makes them **exempt from eviction**, so no amount of churn from the rest of a screen can force them to be re-solved. + +tie_lines = [(107, 203), (113, 215), (123, 217)] +vlodf_watch = VirtualLODF(sys; persistent_arcs = tie_lines, max_cache_size = 100) + +# Across repeated operating points only the base flows are recomputed, which is cheap, and the tie-line factors are read straight from the pinned rows. [`VirtualPTDF`](@ref) takes the same keyword. Pinned rows still count against the budget, so the constructor errors if the pinned set alone would exceed `max_cache_size`. + +# ## Step 4 — Reclaim the memory + +# For [`VirtualPTDF`](@ref) and [`VirtualLODF`](@ref), the row cache is emptied in place to free it. Pinned rows are removed too, with a warning. + +empty!(vlodf.cache) + +# ## Where to go next +# +# - [Reproduce industry DFAX values](@ref) — richer distribution-factor reports +# (transfer, flowgate, N-k) on this same system. +# - [How to Define and Apply Contingencies](@ref) — [`VirtualMODF`](@ref) for +# multi-element post-contingency factors, the next step past single-line screening. +# - [Computational Considerations](@ref) — the sparsity and complexity behind why the +# dense matrices are the ones to avoid at scale. diff --git a/docs/src/tutorials/getting_started.md b/docs/src/tutorials/getting_started.md deleted file mode 100644 index 78ff3342c..000000000 --- a/docs/src/tutorials/getting_started.md +++ /dev/null @@ -1,61 +0,0 @@ -# Quick Start Guide - -!!! note - - `PowerSystemCaseBuilder.jl` is a helper library that makes it easier to reproduce examples in the documentation and tutorials. Normally you would pass your local files to create the system data instead of calling the function `build_system`. - For more details visit [PowerSystemCaseBuilder Documentation](https://sienna-platform.github.io/PowerSystemCaseBuilder.jl/stable) - -For more details about loading data and adding more dynamic components check the -[Creating a System with Dynamic devices](https://sienna-platform.github.io/PowerSystems.jl/stable/tutorials/add_dynamic_data/) -section of the documentation in `PowerSystems.jl`. - -## Loading data - -Data can be loaded from a pss/e raw file and a pss/e dyr file. - -```@repl quick_start_guide -using PowerNetworkMatrices -using PowerSystemCaseBuilder - -import PowerNetworkMatrices as PNM -import PowerSystemCaseBuilder as PSB - -sys = PSB.build_system(PSB.PSITestSystems, "c_sys5") -``` - -## Computation of the PTDF matrix - -Once system data is loaded, network matrices can be evaluated. The following -example shows how the PTDF matrix is computed. - -The function `PTDF` is called for the evaluation of the matrix and other data. These -are stored in a structure of type `PTDF`. - -```@repl quick_start_guide -# evaluate the PTDF structure containing the matrix and other data. -ptdf_matrix = PNM.PTDF(sys); - -# show the PTDF matrix. -PNM.get_data(ptdf_matrix) -``` - -As it can be seen, the PTDF matrix is stored internally in transposed form for computational efficiency. -The function `get_ptdf_data` returns the data in the standard orientation (arcs × buses). - -The matrix axes are indexed by arc tuples `(from_bus_number, to_bus_number)` and bus numbers. -You can inspect the axes and lookup dictionaries as follows: - -```@repl quick_start_guide -# axes and lookup dictionaries describe the arc tuples and bus numbers for each dimension -PNM.get_axes(ptdf_matrix) -PNM.get_lookup(ptdf_matrix) -``` - -Elements can be accessed using arc tuples and bus numbers directly. The example below picks -the first arc and first bus from the matrix axes so it works for any system: - -```@repl quick_start_guide -some_arc = PNM.get_axes(ptdf_matrix)[2][1] -some_bus = PNM.get_axes(ptdf_matrix)[1][1] -ptdf_matrix[some_arc, some_bus] -``` diff --git a/docs/src/tutorials/introduction.jl b/docs/src/tutorials/introduction.jl new file mode 100644 index 000000000..536c39e72 --- /dev/null +++ b/docs/src/tutorials/introduction.jl @@ -0,0 +1,103 @@ +# # Introduction + +# This tutorial uses a [`PTDF`](@ref) to examine current power flow, a [`LODF`](@ref) to examine how a tripped line's power redirects, and network reductions to make the analysis faster to repeat. Together this will answer: +# +# > If a key transmission line trips, which other lines are most at risk of overloading? + +using PowerNetworkMatrices +import PowerNetworkMatrices as PNM +import PowerSystemCaseBuilder as PSB +using Logging + +# !!! note +# +# `PowerSystemCaseBuilder.jl` only supplies the ready-made example systems used +# throughout this documentation. To build a +# [`System`](@extref PowerSystems.System) from your own data, see the +# [PowerSystems.jl documentation](https://sienna-platform.github.io/PowerSystems.jl/stable). + +# ## Step 1 — Load the network + +# Network matrices are built from a [`PowerSystems.System`](@extref PowerSystems.System). + +sys = with_logger(NullLogger()) do + PSB.build_system(PSB.PSSEParsingTestSystems, "psse_14_network_reduction_test_system"); +end + +# We will focus on the line between buses `103` and `104`. + +# ## Step 2 — Determine important lines with the PTDF + +# A [`PTDF`](@ref) (Power Transfer Distribution Factor) matrix is indexed by an **arc tuple** `(from_bus, to_bus)` and a **bus number**, answering what proportion of power injected at a given bus flows through a given branch (when withdrawn at the reference bus). One column of the [`PTDF`](@ref) corresponds to one bus's influence on every branch. + +ptdf = PTDF(sys) +ptdf[(103, 104), 103] + +# Compare a bus that barely touches this line: + +ptdf[(103, 104), 102] + +# Since power is withdrawn at reference buses, power injected at a reference bus does not flow. + +ref_buses = PNM.get_ref_bus(ptdf) +@assert ptdf[(103, 104), only(ref_buses)] == 0.0 +@assert ptdf[(102, 103), only(ref_buses)] == 0.0 + +# ## Step 3 — Determine post-contingency power flow with the LODF + +# The [`LODF`](@ref) (Line Outage Distribution Factor) is indexed by two arc tuples, the **monitored** arc and the **outaged** arc, answering how flow is redistributed in an outage. + +lodf = LODF(sys) + +# We can filter for monitored arcs that absorb flow past a given threshold. + +outaged = (103, 104) +monitored = [(arc, lodf[arc, outaged]) for arc in lodf.axes[1] if arc != outaged] +filter!(pair -> abs(pair[2]) > 0.05, monitored) +sort!(monitored; by = pair -> -abs(pair[2])) +monitored + +# Branch `(102, 103)` inherits the entire flow of the outaged line, a factor of `-1.0`, because it is the series partner on the far side of bus `103`. The two parallel paths `101–115–102` and `101–117–118–104` each absorb roughly `65%`, and branch `(102, 104)` takes the remaining `35%`. Every other branch is unaffected. + +# A negative factor means the redistributed flow runs *against* the monitored branch's `(from, to)` orientation. Individual factors can also be read directly. + +lodf[(102, 103), (103, 104)] + +# ## Step 4 — Reduce the network to repeat the analysis cheaply + +# Reliability studies are rerun across every credible outage and every operating point, so the matrices should be as small as possible. Networks typically contain buses that do not affect a study of this kind: dead-end (radial) buses and pass-through (degree-two) buses. PNM can reduce them away, shrinking the matrices while leaving the surviving sensitivities unchanged. + +# Reductions are supplied to any constructor through the `network_reductions` keyword. Here we combine [`RadialReduction`](@ref), which drops dangling buses, with [`DegreeTwoReduction`](@ref), which fuses pass-through chains. + +reductions = NetworkReduction[RadialReduction(), DegreeTwoReduction()] +ptdf_reduced = PTDF(sys; network_reductions = reductions) +size(ptdf_reduced), size(ptdf) + +# The reduced matrix has fewer branch rows and fewer bus columns, but the sensitivity from Step 2 is unchanged. + +ptdf_reduced[(103, 104), 103], ptdf[(103, 104), 103] + +# The same holds for the [`LODF`](@ref) screen. + +lodf_reduced = LODF(sys; network_reductions = reductions) +lodf_reduced[(102, 103), (103, 104)], lodf[(102, 103), (103, 104)] + +# A reduction therefore gives a smaller and faster pair of matrices that answer identically for every element that survives it. + +# !!! warning +# When screening a *reduced* network against ratings, a degree-two merge fuses several +# branches into one equivalent arc whose limit is an aggregate. Use +# [`get_single_element_contingency_rating`](@ref) and the related aggregated-rating +# accessors rather than a single branch's `PSY.get_rating`. + +# Note that every lookup above used an **arc tuple** rather than a branch name. Arc tuples are an unambiguous identifier and they survive reductions, whereas named branches may be merged away, so they are the recommended identifier to use in code that runs both before and after a reduction. + +# ## Where to go next +# +# - [Analysis at Scale](@ref) — the second tutorial: screening *many* contingencies +# on a large network with the memory-light virtual matrices and cache control. +# - [Matrix overview & indexing](@ref) — the reference for every matrix type, its +# axes, indexing, and accessors. +# - [The DC Power Flow Approximation](@ref) — the theory behind these sensitivities. +# - [How to Define and Apply Contingencies](@ref) — post-contingency factors beyond +# the single-outage [`LODF`](@ref). diff --git a/docs/src/tutorials/tutorial_DFAX.md b/docs/src/tutorials/tutorial_DFAX.md deleted file mode 100644 index 2f37f64ae..000000000 --- a/docs/src/tutorials/tutorial_DFAX.md +++ /dev/null @@ -1,360 +0,0 @@ -# Industry DFAX values - -In this tutorial we map the industry's "DFAX" vocabulary (Distribution -Factors, as used in PSS/E's DFAX activity and downstream by the NERC IDC for -TLR and Congestion Management Procedures) onto the matrices provided by -`PowerNetworkMatrices`. Every flavor of DFAX — GSF, LSF, LODF, OTDF, transfer -DFAX, flowgate DFAX, and multi-element (N-k) DFAX — is a special case of one -unified formula. The tutorial walks through each case with executable -examples on the RTS-GMLC system. - -We assume you have already worked through the [PTDF matrix](@ref) and -[LODF matrix](@ref) tutorials. For the Woodbury-identity derivation behind -post-contingency PTDF rows, see the [Flowgate Methodology](@ref) explanation. - -## What is DFAX? - -"DFAX" is shorthand for *distribution factor*. The term originated with -PSS/E's DFAX activity, which writes a `.dfx` file consumed by the NERC -Interchange Distribution Calculator (IDC). The IDC uses these distribution -factors during Transmission Loading Relief (TLR) procedures and Congestion -Management Procedures (CMP) — for example, applying a 5% threshold to decide -whether a given source-to-sink transfer is a "significant" contributor to a -monitored flowgate. - -In practice DFAX is an umbrella term that covers several specific -quantities. The table below maps each industry term to the -`PowerNetworkMatrices` primitive that computes it: - -| Industry term | What it answers | PNM primitive | -|:------------------------- |:------------------------------------------------ |:---------------------------------------------------------------------- | -| GSF / ISF | Δflow on `m` per 1 MW injection at bus `b` | `PTDF[m, b]` | -| LSF | Same as GSF for loads (opposite sign) | `-PTDF[m, b]` | -| LODF | Flow redistribution after one branch trips | `LODF[m, c]` | -| OTDF | GSF with a contingency already in place | `PTDF[m,b] + LODF[m,c]·PTDF[c,b]`, or one entry of a `VirtualMODF` row | -| Transfer DFAX (pre-cont.) | Fraction of a source→sink transfer reaching `m` | `PTDF[m,:]·(s_v − k_v)` | -| Flowgate DFAX | Transfer DFAX on a (monitored, contingency) pair | `VirtualMODF[m, ctg]·(s_v − k_v)` | -| Multi-element (N-k) DFAX | Same with multiple simultaneous outages | `VirtualMODF` (multi-arc `NetworkModification`) | - -The Phase Shifter Factor (PSF) is part of the broader DFAX vocabulary but is -not a first-class primitive in `PowerNetworkMatrices`. Users who need it can -build it manually through `NetworkModification` and `Ybus`; the present -tutorial covers only flow-based distribution factors. - -## The unified DFAX formula - -In the DC power-flow model every flavor of DFAX is a special case of the -same quantity. For a monitored arc ``m``, a source participation vector -``s_v``, a sink participation vector ``k_v``, and a (possibly empty) set of -network modifications ``C``, - -```math -\mathrm{DFAX}(m,\ s \to k,\ C) \;=\; \mathrm{PTDF}_C[m,\,:] \cdot (s_v - k_v), -``` - -where ``\mathrm{PTDF}_C`` is the post-modification PTDF (equal to the base -``\mathrm{PTDF}`` when ``C = \emptyset``). The formula has two degrees of -freedom — *who is shifting* (the source/sink vectors) and *what state the -network is in* (the contingency ``C``). Each section below fixes one or -both. - -The reference (slack) bus is implicit in `PTDF`: a row of `PTDF` already -encodes "inject at bus ``b``, absorb at the slack". So setting ``k_v = 0`` -in the formula means "let the slack absorb the sink", and the GSF section -below reduces to a single `PTDF` entry. - -## Setup - -All subsequent sections build on this setup block. It loads the RTS-GMLC -system and constructs the three matrices the rest of the tutorial uses. - -```@repl tutorial_DFAX -using PowerSystems -using PowerNetworkMatrices -using PowerSystemCaseBuilder -using DataFrames - -const PSY = PowerSystems -const PNM = PowerNetworkMatrices -const PSB = PowerSystemCaseBuilder - -sys = PSB.build_system(PSB.PSISystems, "RTS_GMLC_DA_sys"); - -ptdf = PTDF(sys); -lodf = LODF(sys); -vmodf = VirtualMODF(sys); -``` - -`VirtualMODF` is the most general object — it can compute post-modification -PTDF rows under any contingency. We also build `PTDF` and `LODF` up front -because the pre-contingency and single-element-outage sections below use -them directly (faster than going through Woodbury when those special cases -apply). - -## GSF and LSF (no contingency, point source) - -The simplest special case sets ``k_v = 0`` (the slack absorbs the sink) and -``s_v = e_b`` (a unit vector at one bus). The unified formula collapses to -a single `PTDF` entry — this is the **Generation Shift Factor**: - -```@repl tutorial_DFAX -m = (107, 203); # monitored arc AB1 (Area 1 → Area 2) -b = 101; # injection bus in Area 1 -gsf = ptdf[m, b] -``` - -The **Load Shift Factor** is the same quantity with the opposite sign -(loads withdraw power instead of inject): - -```@repl tutorial_DFAX -lsf = -gsf -``` - -### Subsystem-aggregated GSF - -In practice analysts care about a *subsystem* of generators (for example, -all generators in an area) rather than a single bus. Build a participation -vector by weighting each generator's bus by its `Pmax` share within the -subsystem, then dot the vector with the `PTDF` row: - -```@repl tutorial_DFAX -area1_gens = filter( - g -> PSY.get_name(PSY.get_area(PSY.get_bus(g))) == "1", - collect(PSY.get_available_components(PSY.Generator, sys)), -); - -total_pmax = sum(PSY.get_max_active_power, area1_gens); - -src_weights = Dict{Int, Float64}(); -for g in area1_gens - bn = PSY.get_number(PSY.get_bus(g)) - src_weights[bn] = get(src_weights, bn, 0.0) + - PSY.get_max_active_power(g) / total_pmax -end - -gsf_area1 = sum(w * ptdf[m, bn] for (bn, w) in src_weights) -``` - -`gsf_area1` is the fraction of an aggregate 1 MW dispatch increase across -all Area 1 generators (split by `Pmax`) that lands on AB1. The slack still -absorbs the corresponding withdrawal — this is a *one-sided* shift. - -If the swing should be distributed across many buses instead of falling on -the single reference bus, pass a `dist_slack` dictionary to the `PTDF` -constructor (see the [PTDF matrix](@ref) tutorial). That is a different -concept from subsystem aggregation: `dist_slack` redefines the reference, -whereas the participation vector above defines the *source* of the -transfer. - -## Transfer DFAX (no contingency, multi-bus source and sink) - -When both source and sink are subsystems, the unified formula is the -difference of two weighted `PTDF` row dot-products: - -```math -\mathrm{TDF}(m,\ s \to k) \;=\; \mathrm{PTDF}[m,\,:] \cdot s_v - \;-\; \mathrm{PTDF}[m,\,:] \cdot k_v. -``` - -Build the sink vector from Area 2 loads, max-active-power weighted. We -filter to `PSY.PowerLoad` because the abstract `ElectricLoad` type also -covers shunt admittance components (`FixedAdmittance`) that don't carry a -real-power weight: - -```@repl tutorial_DFAX -area2_loads = filter( - l -> PSY.get_name(PSY.get_area(PSY.get_bus(l))) == "2", - collect(PSY.get_available_components(PSY.PowerLoad, sys)), -); - -total_load = sum(PSY.get_max_active_power, area2_loads); - -snk_weights = Dict{Int, Float64}(); -for l in area2_loads - bn = PSY.get_number(PSY.get_bus(l)) - snk_weights[bn] = get(snk_weights, bn, 0.0) + - PSY.get_max_active_power(l) / total_load -end -``` - -The pre-contingency transfer DFAX for Area 1 → Area 2 on AB1 is then: - -```@repl tutorial_DFAX -tdf_pre = - sum(w * ptdf[m, bn] for (bn, w) in src_weights) - - sum(w * ptdf[m, bn] for (bn, w) in snk_weights) -``` - -`tdf_pre` answers: *if Area 1 ramps up by 1 MW (split by generator `Pmax`) -and Area 2's load grows by 1 MW (split by load size), what fraction of -that transfer shows up on AB1?* In market and TLR settings, this is the -pre-contingency component of the flowgate impact. - -## OTDF (single contingency, point source) - -The **Outage Transfer Distribution Factor** is the GSF you would observe -if a specific outage were already in effect. For a single-element -contingency on arc ``c``, OTDF has a closed-form expression in terms of -`PTDF` and `LODF`: - -```math -\mathrm{OTDF}(m, b, c) \;=\; \mathrm{PTDF}[m, b] + \mathrm{LODF}[m, c] \cdot \mathrm{PTDF}[c, b]. -``` - -This is the unified formula with ``C = \{c\}`` and the slack absorbing the -sink. `VirtualMODF` computes the same quantity through the Woodbury -identity, which generalizes naturally to multi-element contingencies (see -the N-k section below). For a single outage the two routes agree: - -```@repl tutorial_DFAX -c = (113, 215); # contingency: AB2 outage -otdf_closed = ptdf[m, b] + lodf[m, c] * ptdf[c, b] - -ctg = NetworkModification(vmodf, c); -row_c = vmodf[m, ctg]; -bus_lookup = PNM.get_bus_lookup(vmodf); -otdf_vmodf = row_c[bus_lookup[b]] - -isapprox(otdf_closed, otdf_vmodf; rtol = 1e-10) -``` - -The `isapprox` check is the tutorial's internal validation: `VirtualMODF` -and the closed-form LODF expansion are the same calculation expressed two -different ways. Whenever both apply, they agree to floating-point -tolerance. - -## Flowgate DFAX (single contingency, source–sink transfer) - -A *flowgate* in NERC parlance is the pair `(monitored facility, contingency)`. The flowgate DFAX is the unified formula with both -nontrivial source/sink vectors and a nonempty ``C``: the source–sink -subtraction from the transfer-DFAX section applied to the post-contingency -row from the OTDF section. We reuse `row_c`, `src_weights`, and -`snk_weights` already in scope: - -```@repl tutorial_DFAX -flowgate_dfax = - sum(w * row_c[bus_lookup[bn]] for (bn, w) in src_weights) - - sum(w * row_c[bus_lookup[bn]] for (bn, w) in snk_weights) -``` - -The NERC 5% rule treats a transfer as a "significant" contributor to a -flowgate when the absolute DFAX exceeds 0.05. The check is one line: - -```@repl tutorial_DFAX -significant = abs(flowgate_dfax) >= 0.05 -``` - -When `significant == true`, the transfer is subject to curtailment or -mitigation under the relevant TLR / CMP procedure. - -## N-k DFAX (multi-element contingency) - -When the contingency `C` contains more than one element, the closed-form -LODF expansion of the OTDF section no longer applies — there is no scalar -`LODF[m, c]` when `c` is itself a set. The unified formula still applies, -and `VirtualMODF` is built to handle it directly. Build the multi-element -modification by merging the `arc_modifications` of each single-arc -`NetworkModification` into one combined object: - -```@repl tutorial_DFAX -mod_ab2 = NetworkModification(vmodf, (113, 215)); # AB2 outage -mod_ab3 = NetworkModification(vmodf, (123, 217)); # AB3 outage - -ctg_n2 = NetworkModification( - "AB2_and_AB3_out", - vcat(collect(mod_ab2.arc_modifications), - collect(mod_ab3.arc_modifications)), -); - -row_n2 = vmodf[m, ctg_n2]; - -flowgate_dfax_n2 = - sum(w * row_n2[bus_lookup[bn]] for (bn, w) in src_weights) - - sum(w * row_n2[bus_lookup[bn]] for (bn, w) in snk_weights) -``` - -Removing two of the three parallel Area 1 → Area 2 paths forces a much -larger fraction of any inter-area transfer onto AB1, so `flowgate_dfax_n2` -is substantially larger than the N-1 value computed in the previous -section. The same indexing call (`vmodf[m, ctg]`) handles N-1, N-2, and -higher orders — that is the operational advantage of going through -`VirtualMODF`. - -## Capstone: assembling a DFAX report - -In production use, an analyst typically wants a *table* of distribution -factors covering several transfers, several monitored facilities, and -several contingencies — the kind of report that a PSS/E `.dfx` file plus -IDC post-processing produces. Building that report is one nested loop -around the unified formula. - -For this example we use one transfer (Area 1 → Area 2 from the transfer -DFAX section), three monitored arcs (the three parallel Area 1 → Area 2 -paths), and three contingencies (each of the other two paths individually, -plus the N-2 double-outage from the previous section). Each contingency -carries the set of arc tuples it outages so that we can skip the -ill-defined case of monitoring an outaged element: - -```@repl tutorial_DFAX -monitored = [(107, 203), (113, 215), (123, 217)]; # AB1, AB2, AB3 - -ctg_ab2 = NetworkModification(vmodf, (113, 215)); -ctg_ab3 = NetworkModification(vmodf, (123, 217)); -ctg_ab2_ab3 = NetworkModification( - "AB2_and_AB3_out", - vcat(collect(ctg_ab2.arc_modifications), - collect(ctg_ab3.arc_modifications)), -); - -contingencies = [ - ("AB2 out", Set([(113, 215)]), ctg_ab2), - ("AB3 out", Set([(123, 217)]), ctg_ab3), - ("AB2 & AB3 out", Set([(113, 215), (123, 217)]), ctg_ab2_ab3), -]; - -rows = NamedTuple[] -for mon in monitored - for (label, outaged, ctg_k) in contingencies - mon in outaged && continue - row = vmodf[mon, ctg_k] - df = - sum(w * row[bus_lookup[bn]] for (bn, w) in src_weights) - - sum(w * row[bus_lookup[bn]] for (bn, w) in snk_weights) - push!( - rows, - ( - monitored = mon, - contingency = label, - dfax = df, - significant = abs(df) >= 0.05, - ), - ) - end -end - -report = sort(DataFrame(rows), :dfax; by = abs, rev = true) -``` - -The sort by `abs(dfax)` puts the largest flowgate impacts at the top — the -ones a TLR coordinator would investigate first. Filtering to -`report[report.significant, :]` would keep only NERC-significant rows. - -This table is the kind of output that drives downstream congestion and -seams-coordination workflows; assembling it requires nothing beyond the -matrices in this tutorial. - -## When to use which primitive - -| Need | Reach for | -|:--------------------------------------------------- |:--------------------------------------------- | -| Many transfers, no contingencies | `PTDF` | -| One contingency, all monitored branches | `LODF` + `PTDF` | -| Specific flowgates `(monitored, contingency)` pairs | `VirtualMODF` | -| Multi-element / N-k contingencies | `VirtualMODF` | -| Memory-constrained or sparse usage | `VirtualPTDF` / `VirtualLODF` / `VirtualMODF` | - -`VirtualMODF` is strictly more general than `LODF`-based OTDF arithmetic, -but the closed-form route in the OTDF section is faster when you only have -a single outage and many bus injections to evaluate. The decision is about -which *direction* of the matrix you traverse most often, not about which -one is "correct". diff --git a/docs/src/tutorials/tutorial_DegreeTwoReduction.md b/docs/src/tutorials/tutorial_DegreeTwoReduction.md deleted file mode 100644 index 2829fe680..000000000 --- a/docs/src/tutorials/tutorial_DegreeTwoReduction.md +++ /dev/null @@ -1,178 +0,0 @@ -# DegreeTwoReduction - -In this tutorial the `DegreeTwoReduction` network reduction algorithm is presented. This reduction eliminates buses with exactly two connections by combining the incident branches into a single equivalent branch while preserving the electrical characteristics of the network. - -Before diving into this tutorial we encourage the user to load `PowerNetworkMatrices`, hit the `?` key in the REPL terminal and look for the documentation of `DegreeTwoReduction`. - -## Understanding Degree-Two Buses - -Degree-two buses are nodes in the network topology that have exactly two connections. These intermediate buses can be eliminated by replacing the two incident branches with a single equivalent branch, simplifying the network while maintaining its electrical behavior. The reduction is performed recursively, identifying and eliminating chains of degree-two nodes to maximize network simplification. - -## Basic Usage of DegreeTwoReduction - -The `DegreeTwoReduction` can be applied when constructing various network matrices. The most common use case is with the `Ybus` matrix: - -```@repl tutorial_DegreeTwoReduction -using PowerNetworkMatrices -using PowerSystemCaseBuilder - -import PowerNetworkMatrices as PNM -import PowerSystemCaseBuilder as PSB - -# Load a test system -sys = PSB.build_system(PSSEParsingTestSystems, "psse_14_network_reduction_test_system") - -# Create Ybus with degree-two reduction -ybus = Ybus(sys; network_reductions = NetworkReduction[DegreeTwoReduction()]); -``` - -## Accessing Reduction Information - -After applying the reduction, you can access information about which buses were eliminated and how branches were combined: - -```@repl tutorial_DegreeTwoReduction -# Get the network reduction data -reduction_data = get_network_reduction_data(ybus); - -# View the series branch mapping -# This shows how multiple branches were combined into composite branches -PNM.get_series_branch_map(reduction_data) - -# View the removed buses -PNM.get_removed_buses(reduction_data) - -# View the removed arcs (series branches that were combined) -PNM.get_removed_arcs(reduction_data) -``` - -## Configuration Options - -The `DegreeTwoReduction` provides several configuration options: - -### Protecting Specific Buses - -You can protect certain buses from reduction even if they have degree two: - -```@repl tutorial_DegreeTwoReduction -# Create degree-two reduction that protects specific buses -reduction = DegreeTwoReduction(; irreducible_buses = [115]); - -# Apply to system (if these buses exist in the system) -ybus_protected = Ybus(sys; network_reductions = NetworkReduction[reduction]); -reduction_data_protected = get_network_reduction_data(ybus_protected); -# Compare with unprotected case: -PNM.get_removed_buses(reduction_data) -PNM.get_removed_buses(reduction_data_protected) -``` - -### Handling Reactive Power Injectors - -By default, `DegreeTwoReduction` reduces buses with reactive power injections. You can change this behavior: - -```@repl tutorial_DegreeTwoReduction -# Create reduction that preserves buses with reactive power injections -reduction = DegreeTwoReduction(; reduce_reactive_power_injectors = false); - -# Apply to system -ybus_preserve_reactive = Ybus(sys; network_reductions = NetworkReduction[reduction]); -``` - -## Combining with Other Network Matrices - -The `DegreeTwoReduction` can be applied to other network matrix types as well: - -```@repl tutorial_DegreeTwoReduction -# Apply to PTDF matrix -ptdf = PTDF(sys; network_reductions = NetworkReduction[DegreeTwoReduction()]); - -# Apply to LODF matrix -lodf = LODF(sys; network_reductions = NetworkReduction[DegreeTwoReduction()]); - -# Apply to BA Matrix -ba_matrix = BA_Matrix(sys; network_reductions = NetworkReduction[DegreeTwoReduction()]); - -# Apply to ABA Matrix -aba_matrix = ABA_Matrix(sys; network_reductions = NetworkReduction[DegreeTwoReduction()]); -``` - -## Benefits of Degree-Two Reduction - -Using `DegreeTwoReduction` provides several advantages: - - 1. **Smaller Matrices**: Eliminates intermediate buses from network matrices - 2. **Faster Computations**: Reduced matrix dimensions lead to faster operations - 3. **Simplified Topology**: Creates a more direct representation of the network - 4. **Preserved Accuracy**: Maintains exact electrical equivalence for the reduced network - -## Example: Comparing Matrix Sizes - -```@repl tutorial_DegreeTwoReduction -# Create Ybus without reduction -ybus_full = Ybus(sys); - -# Create Ybus with degree-two reduction -ybus_reduced = Ybus(sys; network_reductions = NetworkReduction[DegreeTwoReduction()]); - -# Compare sizes -size(ybus_full) -size(ybus_reduced) -``` - -## Understanding Series Branch Chains - -When degree-two buses are eliminated, the reduction algorithm identifies chains of series-connected branches. For example: - -``` -Bus A --- Branch 1 --- Bus B --- Branch 2 --- Bus C -``` - -If Bus B has degree two, it can be eliminated, and Branches 1 and 2 are combined into a single equivalent branch: - -``` -Bus A --- Equivalent Branch --- Bus C -``` - -The equivalent branch's electrical parameters (impedance, admittance) are calculated to preserve the overall electrical behavior. - -## Combining Multiple Reductions - -`DegreeTwoReduction` can be combined with other network reduction algorithms like `RadialReduction`: - -```@repl tutorial_DegreeTwoReduction -# Apply both radial and degree-two reductions -reductions = [RadialReduction(), DegreeTwoReduction()]; -ybus_multi = Ybus(sys; network_reductions = reductions); - -# Get combined reduction data -multi_reduction_data = get_network_reduction_data(ybus_multi); -``` - -## Order of Reductions - -When combining multiple reductions, the order can affect the final result: - -```@repl tutorial_DegreeTwoReduction -# First apply radial, then degree-two -reductions1 = [RadialReduction(), DegreeTwoReduction()]; -ybus1 = Ybus(sys; network_reductions = reductions1) - -# First apply degree-two, then radial -reductions2 = [DegreeTwoReduction(), RadialReduction()]; -ybus2 = Ybus(sys; network_reductions = reductions2) - -# Compare results -size(ybus1) -size(ybus2) -``` - -In this case, the result is the same, however this is not guaranteed. -In general, applying `RadialReduction` first is recommended, as it can create new degree-two buses that can then be eliminated by `DegreeTwoReduction`. - -## Important Notes - - - **Topology Preservation**: The reduction maintains essential network connectivity - - **Reference Bus Protection**: Reference (slack) buses are automatically protected from elimination - - **Parallel Paths**: The algorithm handles parallel branches correctly - - **Three-Winding Transformers**: Special handling for three-winding transformer connections - - **Reversibility**: The reduction maintains detailed mapping information for result interpretation - - **Electrical Equivalence**: Equivalent branches are computed to maintain exact electrical behavior diff --git a/docs/src/tutorials/tutorial_Incidence_BA_ABA_matrices.md b/docs/src/tutorials/tutorial_Incidence_BA_ABA_matrices.md deleted file mode 100644 index cb721e003..000000000 --- a/docs/src/tutorials/tutorial_Incidence_BA_ABA_matrices.md +++ /dev/null @@ -1,139 +0,0 @@ -# Incidence, BA and ABA matrices - -In this tutorial the `IncidenceMatrix`, `BA_matrix` and `ABA_matrix` are presented. -The methods used for their evaluation, as well as how data is stored is shown in -the following subsections. - -The matrices here presented are the building blocks for the computation of the PTDF and LODF matrices. - -## IncidenceMatrix - -The `PowerNetworkMatrices` package defines the structure `IncidenceMatrix`, which -store the Incidence Matrix of the considered system as well as the most relevant network data. - -At first, the `System` data is loaded. - -```@repl tutorial_Incidence_BA_ABA_matrices -using PowerNetworkMatrices -using PowerSystemCaseBuilder - -import PowerNetworkMatrices as PNM -import PowerSystemCaseBuilder as PSB - -sys = PSB.build_system(PSB.PSITestSystems, "c_sys5"); -``` - -Then the Incidence Matrix is computed as follows: - -```@repl tutorial_Incidence_BA_ABA_matrices -incidence_matrix = PNM.IncidenceMatrix(sys) -``` - -The `incidence_matrix` variable is a structure of type `IncidenceMatrix`. Getter functions are available to access additional -information: - -```@repl tutorial_Incidence_BA_ABA_matrices -# axis names: row and column names. -# row names: tuples of the arcs (from bus number, to bus number) -# column names: names of the buses -PNM.get_axes(incidence_matrix) - -# data: the incidence matrix data -PNM.get_data(incidence_matrix) - -# lookup: dictionary linking the arc tuples and bus numbers with the row -# and column numbers, respectively. -PNM.get_lookup(incidence_matrix) - -# ref_bus_positions: set containing the positions of the reference buses. -# this represents the positions where to add the column of zeros. Please refer to the -# example in the BA matrix for more details. -PNM.get_ref_bus_position(incidence_matrix) -``` - -Note that the number of columns is lower than the actual number of system buses since -the column related to the reference bus is discarded. - -## BA_Matrix - -The `BA_Matrix` is a structure containing the matrix coming from the product of the -`IncidenceMatrix` and the diagonal matrix containing the impedance of the system's branches ("B" matrix). - -The `BA_Matrix` is computed as follows: - -```@repl tutorial_Incidence_BA_ABA_matrices -ba_matrix = PNM.BA_Matrix(sys) -``` - -Note that the axes order matches the `IncidenceMatrix` (arcs x buses), but for computational considerations the raw data is the transposed matrix. - -The matrix data can similarly be accessed with getter functions: - -```@repl tutorial_Incidence_BA_ABA_matrices -PNM.get_data(ba_matrix) -``` - -Note that the number of columns is lower than the actual number of system buses since -the column related to the reference bus is discarded. - -To add the column of zeros related to the reference bus, it is necessary to use the -information from `get_ref_bus_position`. - -```@repl tutorial_Incidence_BA_ABA_matrices -# assumes a single reference bus -ref_bus_position = first(PNM.get_ref_bus_position(ba_matrix)) -new_ba_matrix = hcat( - ba_matrix.data[:, 1:(ref_bus_position - 1)], - zeros(size(ba_matrix, 1), 1), - ba_matrix.data[:, ref_bus_position:end], -) -``` - -However, trying to change the data field with a matrix of different dimension -will result in an error. - -```@repl tutorial_Incidence_BA_ABA_matrices -ba_matrix.data = new_ba_matrix -``` - -## ABA_Matrix - -The `ABA_Matrix` is a structure containing the matrix coming from the product of the -`IncidenceMatrix` and the `BA_Matrix`. -It features the same fields as the `IncidenceMatrix` and the `BA_Matrix`, plus the `K` one. -The field `ABA_Matrix.K` stores the LU factorization matrices (using the -methods contained in the package `KLU`). - -To evaluate the `ABA_Matrix`, the following command is sufficient: - -```@repl tutorial_Incidence_BA_ABA_matrices -aba_matrix = ABA_Matrix(sys); -``` - -By default the LU factorization matrices are not computed, leaving the `K` field empty: - -```@repl tutorial_Incidence_BA_ABA_matrices -isnothing(aba_matrix.K) -``` - -In case these are wanted, the keyword `factorize` must be true. - -```@repl tutorial_Incidence_BA_ABA_matrices -aba_matrix_with_LU = ABA_Matrix(sys; factorize = true); - -aba_matrix_with_LU.K -``` - -If the `ABA_Matrix` is already computed but the LU factorization was not performed, this can be done by considering the following command: - -```@repl tutorial_Incidence_BA_ABA_matrices -aba_matrix.K -aba_matrix = factorize(aba_matrix); -aba_matrix.K -``` - -The following command can then be used to check if the `ABA_Matrix` contains the LU factorization matrices: - -```@repl tutorial_Incidence_BA_ABA_matrices -is_factorized(aba_matrix) -``` diff --git a/docs/src/tutorials/tutorial_LODF_matrix.md b/docs/src/tutorials/tutorial_LODF_matrix.md deleted file mode 100644 index 634a7b311..000000000 --- a/docs/src/tutorials/tutorial_LODF_matrix.md +++ /dev/null @@ -1,119 +0,0 @@ -# LODF matrix - -In this tutorial the methods for computing the Line Outage Distribution Factor (`LODF`) are presented. -Before diving into this tutorial we encourage the user to load `PowerNetworkMatrices`, hit the `?` key in the REPL terminal and look for the documentation of the different `LODF` methods available. - -## Evaluation of the `LODF` matrix - -As for the `PTDF` matrix, the `LODF` one can be evaluated according to three different approaches: - - - `Dense`: considers functions for dense matrix multiplication and inversion - - `KLU`: considers functions for sparse matrix multiplication and inversion (**default**) - - `MKLPardiso`: uses Intel's MKL Pardiso solver for sparse matrix operations (only available on Intel-based systems running Linux or Windows) - -The evaluation of the `LODF` matrix can be easily performed starting from importing the system's data and then by simply calling the `LODF` method. - -```@repl tutorial_PTDF_matrix -using PowerNetworkMatrices -using PowerSystemCaseBuilder - -import PowerNetworkMatrices as PNM -import PowerSystemCaseBuilder as PSB - -# get the System data -sys = PSB.build_system(PSB.PSITestSystems, "c_sys5"); - -# compute the LODF matrix -lodf_1 = LODF(sys); - -lodf_2 = LODF(sys; linear_solver = "Dense"); - -# show matrix -get_lodf_data(lodf_1) -``` - -## Indexing the `LODF` matrix - -The `LODF` matrix is indexed by **arc tuples** `(from_bus_number, to_bus_number)` for both dimensions. Both the row (selected line) and column (outage line) use arc tuples as identifiers. - -```@repl tutorial_PTDF_matrix -# inspect the axes and lookup dictionaries (both dimensions are arc tuples) -get_axes(lodf_1) -get_lookup(lodf_1) -``` - -Elements of the `LODF` matrix can be accessed by arc tuples: - -```@repl tutorial_PTDF_matrix -# access LODF element: flow change on arc (1, 4) due to outage of arc (2, 3) -lodf_1[(1, 4), (2, 3)] -``` - -!!! note - - For backward compatibility, branch name strings can also be used to index the `LODF` matrix (e.g. `lodf_1["branch_name_1", "branch_name_2"]`). This uses `get_branch_multiplier` internally to map branch names to their corresponding arc tuples. Using arc tuples directly is recommended. - -## Computing `LODF` from pre-computed matrices - -Advanced users might be interested in computing the `LODF` matrix starting from either the `IncidenceMatrix` and `PTDF` structures (`CASE 1`), or by the information related to `IncidenceMatrix`, `BA_Matrix` and `ABA_Matrix` (`CASE 2`). - -```@repl tutorial_PTDF_matrix -# CASE 1 - -# get the Incidence and PTDF matrix -a = IncidenceMatrix(sys); -ptdf = PTDF(sys); - -# compute LODF matrix with the two network matrices -lodf_3 = LODF(a, ptdf); - -# CASE 2 - -# get the BA and ABA matrices (ABA matrix must include LU factorization -# matrices) -ba = BA_Matrix(sys); -aba = ABA_Matrix(sys; factorize = true); - -# compute LODF matrix with the three network matrices -lodf_4 = LODF(a, aba, ba); -``` - -**NOTE:** whenever the method `LODF(sys::System)` is used, the methods previously defined for `CASE 1` are executed in sequence. Therefore the method `LODF(a::IncidenceMatrix, ptdf::PTDF)` is the default one when evaluating the `LODF` matrix from the `System` data directly. - -## Available methods for the computation of the `LODF` matrix - -For those methods that either require the evaluation of the `PTDF` matrix, or that execute this evaluation internally, three different approaches can be used. - -As for the `PTDF` matrix, here too the optional argument `linear_solver` can be specified with either `KLU` (for sparse matrix calculation), `Dense` (for dense matrix calculation), or `MKLPardiso` (for Intel MKL Pardiso sparse solver). - -```@repl tutorial_PTDF_matrix -lodf_dense = LODF(sys; linear_solver = "Dense"); - -lodf_klu = LODF(sys; linear_solver = "KLU"); -``` - -**NOTE (1):** by default the "KLU" method is selected, which appeared to require significant less time and memory with respect to "Dense". -Please note that regardless of which method (`KLU`, `Dense`, or `MKLPardiso`) is used, the resulting `LODF` matrix is stored as a dense one. - -**NOTE (2):** for the moment, the method `LODF(a::IncidenceMatrix, aba::ABA_Matrix, ba::BA_Matrix)` will take `KLU` as `linear_solver` option. - -**Note on MKLPardiso**: The `MKLPardiso` solver option is only available on Intel-based systems running Linux or Windows. On other platforms (e.g., ARM-based systems or macOS), use `KLU` or `Dense` instead. - -## "Sparse" `LODF` matrix - -The `LODF` matrix can be computed in a "sparse" fashion by defining the input argument `tol`. If this argument is defined, then elements of the `LODF` matrix whose absolute values are below the set threshold are dropped. In addition, the matrix will be stored as a sparse one of type `SparseArrays.SparseMatrixCSC{Float64, Int}` type instead of `Matrix{Float64}` one. - -By considering an "extreme" value of 0.4 as `tol`, the `LODF` matrix can be computed as follows: - -```@repl tutorial_PTDF_matrix -lodf_sparse = LODF(sys; tol = 0.4); -get_lodf_data(lodf_sparse) -``` - -Note that in practice much smaller values of `tol` are typically used (e.g., 1e-5). - -**NOTE (1):** elements whose absolute values exceed the `tol` argument are removed from the `LODF` matrix *after* this has been computed. - -**NOTE (2):** the `tol` argument does not refer to the "sparsification" tolerance of the `PTDF` matrix that is computed in the `LODF` method. - -**NOTE (3):** in case the method `LODF(a::IncidenceMatrix, ptdf::PTDF)` is considered, an error will be thrown whenever the `tol` argument in the `PTDF` structure used as input is different than `1e-15`. diff --git a/docs/src/tutorials/tutorial_PTDF_matrix.md b/docs/src/tutorials/tutorial_PTDF_matrix.md deleted file mode 100644 index 507bbba3e..000000000 --- a/docs/src/tutorials/tutorial_PTDF_matrix.md +++ /dev/null @@ -1,203 +0,0 @@ -# PTDF matrix - -In this tutorial the methods for computing the Power Transfer Distribution Factors (`PTDF`) are presented. -Before diving into this tutorial we encourage the user to load `PowerNetworkMatrices`, hit the `?` key in the REPL terminal and look for the documention of the different `PTDF` methods available. - -## Evaluation of the `PTDF` matrix - -The `PTDF` matrix can be evaluated according to three different approaches: - - - `Dense`: considers functions for dense matrix multiplication and inversion - - `KLU`: considers functions for sparse matrix multiplication and inversion (**default**) - - `MKLPardiso`: uses Intel's MKL Pardiso solver for sparse matrix operations (only available on Intel-based systems running Linux or Windows) - -The evaluation of the `PTDF` matrix can be easily performed starting from importing the system's data and then by simply calling the `PTDF` method. - -```@repl tutorial_PTDF_matrix -using PowerSystems -using PowerNetworkMatrices -using PowerSystemCaseBuilder - -const PSY = PowerSystems -const PNM = PowerNetworkMatrices -const PSB = PowerSystemCaseBuilder - -sys = PSB.build_system(PSB.PSITestSystems, "c_sys5"); - -ptdf_1 = PTDF(sys); - -get_ptdf_data(ptdf_1) -``` - -Note that while the `PTDF` stores the transpose of the matrix data, the function `get_ptdf_data` returns the data in the standard orientation (arcs × buses). - -## Indexing the `PTDF` matrix - -The `PTDF` matrix is indexed by **arc tuples** `(from_bus_number, to_bus_number)` for the row dimension and **bus numbers** for the column dimension. - -```@repl tutorial_PTDF_matrix -# inspect the axes and lookup dictionaries -get_axes(ptdf_1) -get_lookup(ptdf_1) -``` - -Elements of the `PTDF` matrix can be accessed by arc tuple and bus number: - -```@repl tutorial_PTDF_matrix -# access PTDF element for arc (1, 2) and bus 3 -ptdf_1[(1, 2), 3] -``` - -!!! note - - For backward compatibility, branch name strings can also be used to index the `PTDF` matrix (e.g. `ptdf_1["branch_name", 3]`). This uses `get_branch_multiplier` internally to map the branch name to its corresponding arc tuple. Using arc tuples directly is recommended. - -## Computing `PTDF` from pre-computed matrices - -Advanced users might be interested in computing the `PTDF` matrix starting from either the data contained in the `IncidenceMatrix` and `BA_matrix` structures. - -```@repl tutorial_PTDF_matrix -# evaluate the BA_matrix and Incidence_Matrix -ba_matrix = BA_Matrix(sys); -a_matrix = IncidenceMatrix(sys); - -# get the PTDF matrix starting from the values of the -# previously computed matrices -ptdf_2 = PTDF(a_matrix, ba_matrix); -get_ptdf_data(ptdf_2) -``` - -## Available methods for the computation of the `PTDF` matrix - -As previously mentioned, the `PTDF` matrix can be evaluated considering different approaches. The method can be selected by specifying the field `linear_solver` in the `PTDF` function. - -```@repl tutorial_PTDF_matrix -ptdf_dense = PTDF(sys; linear_solver = "Dense"); -get_ptdf_data(ptdf_dense) - -ptdf_klu = PTDF(sys; linear_solver = "KLU"); -get_ptdf_data(ptdf_klu) - -# For MKLPardiso (if available) -# ptdf_mkl = PTDF(sys; linear_solver = "MKLPardiso"); -# get_ptdf_data(ptdf_mkl) -``` - -By default the "KLU" method is selected, which appeared to require significant less time and memory with respect to "Dense". -Please note that regardless of which method (`KLU`, `Dense`, or `MKLPardiso`) is used, the resulting `PTDF` matrix is stored as a dense one. - -**Note on MKLPardiso**: The `MKLPardiso` solver option is only available on Intel-based systems running Linux or Windows. On other platforms (e.g., ARM-based systems or macOS), use `KLU` or `Dense` instead. - -## Evaluating the `PTDF` matrix considering distributed slack bus - -Whenever needed, the `PTDF` matrix can be computed with a distributed slack bus. To do so, a vector of type `Dict{Int64, Float64}` needs to be defined, specifying the weights for each bus of the system. These weights identify how the load on the slack bus is redistributed accross the system. - -```@repl tutorial_PTDF_matrix -# consider equal distribution accross each bus for this example -buscount = length(PSY.get_available_components(PSY.ACBus, sys)); -dist_slack = 1 / buscount * ones(buscount); -dist_slack_dict = Dict(i => dist_slack[i] / sum(dist_slack) for i in 1:buscount); -``` - -Once the dictionary of the weights is defined, the `PTDF` matrix can be computed by defining the input argument `dist_slack` (empty array `Dict{Int64, Float64}()` by default): - -```@repl tutorial_PTDF_matrix -ptdf_distr = PTDF(sys; dist_slack = dist_slack_dict); -``` - -The difference between a the matrix computed with and without the `dist_slack` field defined can be seen as follows: - -```@repl tutorial_PTDF_matrix -# with no distributed slack bus -get_ptdf_data(ptdf_klu) -# with distributed slack bus -get_ptdf_data(ptdf_distr) -``` - -## "Sparse" `PTDF` matrix - -The `PTDF` matrix can be computed in a "sparse" fashion by defining the input argument `tol`. If this argument is defined, then elements of the `PTDF` matrix whose absolute values are below the set threshold are dropped. In addition, the matrix will be stored as a sparse one of type `SparseArrays.SparseMatrixCSC{Float64, Int}` instead of `Matrix{Float64}`. - -By considering an "extreme" value of 0.2 as `tol`, the `PTDF` matrix can be computed as follows: - -```@repl tutorial_PTDF_matrix -ptdf_sparse = PTDF(sys; tol = 0.2); -get_ptdf_data(ptdf_sparse) -``` - -**NOTE:** In practice, much smaller values are typically used for `tol`(e.g., 1e-5). - -## Network Reductions - -The `PTDF` matrix can be computed with network reductions applied to simplify the system topology. Network reductions eliminate certain buses and branches while preserving the electrical characteristics of the network. This can significantly reduce computation time and memory usage for large systems. - -Two types of network reductions are supported: - - - `RadialReduction`: Eliminates radial (leaf) buses that have only one connection - - `DegreeTwoReduction`: Eliminates degree-two buses (buses with exactly two connections) by combining their incident branches - -For detailed information about these reductions, see the [RadialReduction](@ref) and [DegreeTwoReduction](@ref) tutorials. - -### Using Network Reductions with PTDF - -To apply network reductions, pass a vector of `NetworkReduction` objects to the `network_reductions` keyword argument: - -```@repl tutorial_PTDF_matrix -# Apply radial reduction -ptdf_radial = PTDF(sys; network_reductions = NetworkReduction[RadialReduction()]); - -# Apply degree-two reduction -ptdf_degree_two = PTDF(sys; network_reductions = NetworkReduction[DegreeTwoReduction()]); - -# Combine multiple reductions (order matters - RadialReduction first is recommended) -ptdf_combined = PTDF( - sys; - network_reductions = NetworkReduction[RadialReduction(), DegreeTwoReduction()], -); -``` - -### Protecting Specific Buses from Reduction - -Both reduction types allow you to specify buses that should not be eliminated using the `irreducible_buses` parameter: - -```@repl tutorial_PTDF_matrix -# Protect specific buses from radial reduction -reduction = RadialReduction(; irreducible_buses = [1, 2]) -ptdf_protected = PTDF(sys; network_reductions = NetworkReduction[reduction]); -``` - -### DegreeTwoReduction Options - -The `DegreeTwoReduction` has an additional option to control whether buses with reactive power injections are reduced: - -```@repl tutorial_PTDF_matrix -# Preserve buses with reactive power injections -reduction = DegreeTwoReduction(; reduce_reactive_power_injectors = false) -ptdf_preserve_reactive = PTDF(sys; network_reductions = NetworkReduction[reduction]); -``` - -### Accessing Reduction Information - -After computing the PTDF matrix with reductions, you can access information about what was reduced: - -```@repl tutorial_PTDF_matrix -ptdf_reduced = PTDF(sys; network_reductions = [RadialReduction(), DegreeTwoReduction()]); - -# Get the reduction data -reduction_data = get_network_reduction_data(ptdf_reduced) -``` - -### Combining Reductions with Other Options - -Network reductions can be combined with other PTDF options like distributed slack and sparsification: - -```@repl tutorial_PTDF_matrix -ptdf_full_options = PTDF(sys; - linear_solver = "KLU", - dist_slack = dist_slack_dict, - tol = 1e-5, - network_reductions = [RadialReduction(), DegreeTwoReduction()], -); -``` - -**NOTE**: The reference (slack) bus is automatically protected from elimination during reductions. diff --git a/docs/src/tutorials/tutorial_RadialReduction.md b/docs/src/tutorials/tutorial_RadialReduction.md deleted file mode 100644 index 0dc624b46..000000000 --- a/docs/src/tutorials/tutorial_RadialReduction.md +++ /dev/null @@ -1,118 +0,0 @@ -# RadialReduction - -In this tutorial the `RadialReduction` network reduction algorithm is presented. This reduction eliminates radial (dangling) buses and their associated branches from the power network while preserving the electrical behavior of the core network. - -Before diving into this tutorial we encourage the user to load `PowerNetworkMatrices`, hit the `?` key in the REPL terminal and look for the documentation of `RadialReduction`. - -## Understanding Radial Branches - -Radial buses are leaf nodes in the network topology with only one connection. These buses do not affect the electrical behavior of the rest of the network and for certain applications can be safely eliminated to simplify network matrices and improve computational efficiency. - -## Basic Usage of RadialReduction - -The `RadialReduction` can be applied when constructing various network matrices. The most common use case is with the `Ybus` matrix: - -```@repl tutorial_RadialReduction -using PowerNetworkMatrices -using PowerSystemCaseBuilder - -import PowerNetworkMatrices as PNM -import PowerSystemCaseBuilder as PSB - -# Load a test system -sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"); - -# Create Ybus with radial reduction -ybus = Ybus(sys; network_reductions = NetworkReduction[RadialReduction()]) -``` - -## Accessing Reduction Information - -After applying the reduction, you can access information about which buses and branches were eliminated: - -```@repl tutorial_RadialReduction -# Get the network reduction data -reduction_data = get_network_reduction_data(ybus); - -# View the bus reduction mapping -# This shows which buses were reduced to which parent buses -get_bus_reduction_map(reduction_data) - -# View the reverse bus search mapping -# This maps each reduced bus to its ultimate parent -PNM.get_reverse_bus_search_map(reduction_data) - -# View the removed arcs (branches) -PNM.get_removed_arcs(reduction_data) -``` - -## Protecting Specific Buses from Reduction - -In some cases, you may want to preserve certain buses even if they are radial. This can be done using the `irreducible_buses` parameter: - -```@repl tutorial_RadialReduction -# Create radial reduction that protects buses 8 and 14 -reduction = RadialReduction(; irreducible_buses = [8, 14]); - -# Apply to system (buses must exist in the system) -ybus_protected = Ybus(sys; network_reductions = NetworkReduction[reduction]); - -# Bus 8 was radial, but preserved from reduction -reduction_data = get_network_reduction_data(ybus_protected); -get_bus_reduction_map(reduction_data) -``` - -## Combining with Other Network Matrices - -The `RadialReduction` can be applied to other network matrix types as well: - -```@repl tutorial_RadialReduction -# Apply to PTDF matrix -ptdf = PTDF(sys; network_reductions = NetworkReduction[RadialReduction()]); - -# Apply to LODF matrix -lodf = LODF(sys; network_reductions = NetworkReduction[RadialReduction()]); - -# Apply to Incidence Matrix -incidence = IncidenceMatrix(sys; network_reductions = NetworkReduction[RadialReduction()]); -``` - -## Benefits of Radial Reduction - -Using `RadialReduction` provides several advantages: - - 1. **Smaller Matrices**: Eliminates unnecessary rows and columns from network matrices - 2. **Faster Computations**: Reduced matrix dimensions lead to faster linear algebra operations - 3. **Better Conditioning**: Removing radial elements can improve numerical properties - 4. **Memory Efficiency**: Reduces storage requirements for large network models - -## Example: Comparing Matrix Sizes - -```@repl tutorial_RadialReduction -# Create Ybus without reduction -ybus_full = Ybus(sys); - -# Create Ybus with radial reduction -ybus_reduced = Ybus(sys; network_reductions = NetworkReduction[RadialReduction()]); - -# Compare sizes -size(ybus_full) -size(ybus_reduced) -``` - -## Combining Multiple Reductions - -`RadialReduction` can be combined with other network reduction algorithms like `DegreeTwoReduction`: - -```@repl tutorial_RadialReduction -# Apply both radial and degree-two reductions -reductions = [RadialReduction(), DegreeTwoReduction()]; -ybus_multi = Ybus(sys; network_reductions = reductions); -``` - -## Important Notes - - - **Reference Bus Protection**: Reference (slack) buses are automatically protected from elimination, regardless of their connectivity - - **Order Matters**: When combining multiple reductions, they are applied in the order specified in the vector - - **Reversibility**: The reduction maintains mapping information (`bus_reduction_map` and `reverse_bus_search_map`) that can be used for result interpretation - - **Electrical Equivalence**: The reduced network maintains the same electrical behavior as the original network for all non-eliminated elements diff --git a/docs/src/tutorials/tutorial_VirtualLODF_matrix.md b/docs/src/tutorials/tutorial_VirtualLODF_matrix.md deleted file mode 100644 index d073f2493..000000000 --- a/docs/src/tutorials/tutorial_VirtualLODF_matrix.md +++ /dev/null @@ -1,98 +0,0 @@ -# VirtualLODF - -The `VirtualLODF` structure follows the same philosofy as the `VirtualPTDF`: it contains rows of the original `LODF` matrix, evaluated and cached on demand. - -Refer to the different arguments of the `VirtualLODF` methods by looking at the "Public API Reference" page. - -## How the `VirtualLODF` works - -The `VirtualLODF` structure retains many of the similarities of the `VirtualPTDF`. However, its computation is more complex and requires some additional data. - -Starting from the system data, the `IncidenceMatrix`, `BA_Matrix` and `ABA_Matrix` (with relative LU factorization matrices) are evaluated. The `ABA_Matrix` and `BA_Matrix` are used for the computation of the diagonal elements of the `PTDF` matrix, and this vector is stored in the `VirtualLODF` structure together with the other structures mentioned above. - -Once the `VirtualLODF` is initialized, each row of the matrix can be evaluated separately and on user request. The algorithmic procedure is the following: - - 1. Define the `VirtualPTDF` structure - 2. Call any element of the matrix to define and store the relative row as well as showing the selected element - -Regarding point 2, if the row has been stored previously then the desired element is just loaded from the cache and shown. - -The flowchart below shows how the `VirtualLODF` is structured and how it works. Examples will be presented in the following sections. - -```@raw html - -``` - -## Initialize `VirtualLODF` and compute/access row/element - -As for the `LODF` matrix, at first the `System` data must be loaded. The "RTS-GMLC" systems is considered as example: - -```@repl tutorial_VirtualPTDF_matrix -using PowerNetworkMatrices -using PowerSystemCaseBuilder - -import PowerNetworkMatrices as PNM -import PowerSystemCaseBuilder as PSB - -sys = PSB.build_system(PSB.PSISystems, "RTS_GMLC_DA_sys"); -``` - -At this point the `VirtualLODF` is initialized with the following simple command: - -```@repl tutorial_VirtualPTDF_matrix -v_lodf = VirtualLODF(sys) -``` - -Now, an element of the matrix can be computed by using the arc tuples as indices: - -```@repl tutorial_VirtualPTDF_matrix -el_v_lodf = v_lodf[(221, 222), (202, 206)] -``` - -This element represent the portion flowing on arc (202, 206) now diverted on arc (221, 222) as a consequence of its outage. - -Alternatively, the value can be indexed by row and column numbers directly. In this case the row and column numbers are mapped by the dictonaries contained in the `lookup` field. - -```@repl tutorial_VirtualPTDF_matrix -row_number = v_lodf.lookup[1][(221, 222)] -col_number = v_lodf.lookup[2][(202, 206)] -el_C31_2_105_bis = v_lodf[row_number, col_number] -``` - -**NOTE**: this example was made for the sake of completeness and considering the actual arc tuples is recommended. - -As previously mentioned, in order to evaluate a single element of the `VirtualLODF`, the entire row related to the selected branch must be considered. For this reason it is cached for later calls. -This is evident by looking at the following example: - -```@repl tutorial_VirtualPTDF_matrix -sys_2k = PSB.build_system(PSB.PSYTestSystems, "tamu_ACTIVSg2000_sys"); - -v_lodf_2k = VirtualLODF(sys_2k); - -# evaluate PTDF row related to arc (5270, 5474) -@time v_lodf_2k[(5270, 5474), (2118, 2113)] - -# call same element after the row has been stored -@time v_lodf_2k[(5270, 5474), (2118, 2113)] -``` - -## "Sparse" `VirtualPTDF` - -Sparsification of each row can be achieved in the same fashion as for the `LODF` matrix, by removing those elements whose absolute values is below a certain threshold. - -As for the example show for the `LODF` matrix, here to a very high values of 0.4 is considered for the `tol` field. Again, this value is considered just for the sake of this example. - -```@repl tutorial_VirtualPTDF_matrix -# smaller system for the next examples -sys_2 = PSB.build_system(PSB.PSITestSystems, "c_sys5"); - -v_lodf_dense = VirtualLODF(sys_2); -v_lodf_sparse = VirtualLODF(sys_2; tol = 0.4); -``` - -Let's now evaluate a row and compare the results: - -```@repl tutorial_VirtualPTDF_matrix -original_row = [v_lodf_dense[(1, 2), j] for j in v_lodf_dense.axes[2]] -sparse_row = [v_lodf_sparse[(1, 2), j] for j in v_lodf_sparse.axes[2]] -``` diff --git a/docs/src/tutorials/tutorial_VirtualPTDF_matrix.md b/docs/src/tutorials/tutorial_VirtualPTDF_matrix.md deleted file mode 100644 index 81e80678d..000000000 --- a/docs/src/tutorials/tutorial_VirtualPTDF_matrix.md +++ /dev/null @@ -1,207 +0,0 @@ -# VirtualPTDF - -Contrary to the traditional `PTDF` matrix, the `VirtualPTDF` is a structure containing rows of the original matrix, related to specific system arcs. -The different rows of the `PTDF` matrix are cached in the `VirtualPTDF` structure as they are evaluated. This allows to keep just the portion of the original matrix which is of interest to the user, avoiding the unnecessary computation of the whole matrix. - -Refer to the different arguments of the `VirtualPTDF` methods by looking at the "Public API Reference" page. - -## How the `VirtualPTDF` works - -The `VirtualPTDF` is a structure containing everything needed to compute any row of the PTDF matrix and store it. To do so, the `VirtualPTDF` must store the BA matrix (coming from the `BA_Matrix` struct) and the inverse of the ABA matrix (coming from `ABA_MAtrix` struct). In particular, `KLU` is used to get the LU factorization matrices of the ABA matrix and these ones are stored, avoid the inversion. - -Once the `VirtualPTDF` is initialized, each row of the PTDF matrix can be evaluated separately. The algorithmic procedure is the following: - - 1. Define the `VirtualPTDF` structure - 2. Call any element of the matrix to define and store the relative row as well as showing the selected element - -Regarding point 2, if the row has been stored previously then the desired element is just loaded from the cache and shown. - -The flowchart below shows how the `VirtualPTDF` is structured and how it works. Examples will be presented in the following sections. - -```@raw html - -``` - -## Initialize `VirtualPTDF` and compute/access row/element - -As for the `PTDF` matrix, at first the `System` data must be loaded. The "RTS-GMLC" systems is considered as example: - -```@repl tutorial_VirtualPTDF_matrix -using PowerNetworkMatrices -using PowerSystems -using PowerSystemCaseBuilder - -import PowerSystems as PSY -import PowerNetworkMatrices as PNM -import PowerSystemCaseBuilder as PSB - -sys = PSB.build_system(PSB.PSISystems, "RTS_GMLC_DA_sys"); -``` - -At this point the `VirtualPTDF` is initialized with the following simple command: - -```@repl tutorial_VirtualPTDF_matrix -v_ptdf = VirtualPTDF(sys); -``` - -Now, an element of the matrix can be computed by calling the arc tuple and bus number: - -```@repl tutorial_VirtualPTDF_matrix -el_C31_105 = v_ptdf[(318, 321), 105] -``` - -Alternatively, the value can be indexed by row and column numbers directly. In this case the row and column numbers are mapped by the dictonaries contained in the `lookup` field. - -```@repl tutorial_VirtualPTDF_matrix -row_number = v_ptdf.lookup[1][(318, 321)] -col_number = v_ptdf.lookup[2][105] -el_C31_105_bis = v_ptdf[row_number, col_number] -``` - -**NOTE**: this example was made for the sake of completeness and considering the actual arc tuple and bus number is recommended. - -As previously mentioned, in order to evaluate a single element of the `VirtualPTDF`, the entire row related to the selected arc must be considered. For this reason it is cached in the `VirtualPTDF` structure for later calls. -This is evident by looking at the following example: - -```@repl tutorial_VirtualPTDF_matrix -sys_2k = PSB.build_system(PSB.PSYTestSystems, "tamu_ACTIVSg2000_sys"); - -v_ptdf_2k = VirtualPTDF(sys_2k); - -# evaluate PTDF row related to arc (5270, 5474) -@time v_ptdf_2k[(5270, 5474), 8155] - -# call same element after the row has been stored -@time v_ptdf_2k[(5270, 5474), 8155] -``` - -## `VirtualPTDF` with distributed slack bus - -As for the `PTDF` matrix, here too each row can be evaluated considering distributed slack buses. -A vector of type `Vector{Float64}` is defined, specifying the weights for each bus of the system. - -```@repl tutorial_VirtualPTDF_matrix -# smaller system for the next examples -sys_2 = PSB.build_system(PSB.PSITestSystems, "c_sys5"); - -# consider equal distribution accross each bus for this example -buscount = length(PSY.get_available_components(PSY.ACBus, sys_2)); -dist_slack = 1 / buscount * ones(buscount); -dis_slack_dict = Dict(i => dist_slack[i] / sum(dist_slack) for i in 1:buscount) -``` - -Now initialize the `VirtualPTDF` by defining the `dist_slack` field with the vector of weights previously computed: - -```@repl tutorial_VirtualPTDF_matrix -v_ptdf_distr = VirtualPTDF(sys_2; dist_slack = dis_slack_dict); -v_ptdf_orig = VirtualPTDF(sys_2); -``` - -Now check the difference with the same row evaluated without considering distributed slack bus. - -```@repl tutorial_VirtualPTDF_matrix -# get the first arc tuple from the axes -first_arc = v_ptdf_distr.axes[1][1] -row_distr = [v_ptdf_distr[first_arc, j] for j in v_ptdf_distr.axes[2]] -row_original = [v_ptdf_orig[first_arc, j] for j in v_ptdf_orig.axes[2]] -``` - -## "Sparse" `VirtualPTDF` - -Sparsification of each row can be achieved in the same fashion as for the `PTDF` matrix, by removing those elements whose absolute values is below a certain threshold. - -As for the example show for the `PTDF` matrix, here to a very high values of 0.2 is considered for the `tol` field. Again, this value is considered just for the sake of this example. - -```@repl tutorial_VirtualPTDF_matrix -v_ptdf_dense = VirtualPTDF(sys_2); -v_ptdf_sparse = VirtualPTDF(sys_2; tol = 0.2); -``` - -Let's now evaluate the same row as before and compare the results: - -```@repl tutorial_VirtualPTDF_matrix -first_arc = v_ptdf_dense.axes[1][1] -original_row = [v_ptdf_dense[first_arc, j] for j in v_ptdf_dense.axes[2]] -sparse_row = [v_ptdf_sparse[first_arc, j] for j in v_ptdf_sparse.axes[2]] -``` - -## Network Reductions - -The `VirtualPTDF` can be computed with network reductions applied to simplify the system topology. Network reductions eliminate certain buses and branches while preserving the electrical characteristics of the network. This can significantly reduce computation time and memory usage for large systems. - -Two types of network reductions are supported: - - - `RadialReduction`: Eliminates radial (leaf) buses that have only one connection - - `DegreeTwoReduction`: Eliminates degree-two buses (buses with exactly two connections) by combining their incident branches - -For detailed information about these reductions, see the [RadialReduction](@ref) and [DegreeTwoReduction](@ref) tutorials. - -### Using Network Reductions with VirtualPTDF - -To apply network reductions, pass a vector of `NetworkReduction` objects to the `network_reductions` keyword argument: - -```@repl tutorial_VirtualPTDF_matrix -# Apply radial reduction -v_ptdf_radial = VirtualPTDF(sys_2; network_reductions = NetworkReduction[RadialReduction()]); - -# Apply degree-two reduction -v_ptdf_degree_two = - VirtualPTDF(sys_2; network_reductions = NetworkReduction[DegreeTwoReduction()]); - -# Combine multiple reductions (order matters - RadialReduction first is recommended) -v_ptdf_combined = - VirtualPTDF( - sys_2; - network_reductions = NetworkReduction[RadialReduction(), DegreeTwoReduction()], - ); -``` - -### Protecting Specific Buses from Reduction - -Both reduction types allow you to specify buses that should not be eliminated using the `irreducible_buses` parameter: - -```@repl tutorial_VirtualPTDF_matrix -# Protect specific buses from radial reduction -reduction = RadialReduction(; irreducible_buses = [1, 2]) -v_ptdf_protected = VirtualPTDF(sys_2; network_reductions = NetworkReduction[reduction]); -``` - -### DegreeTwoReduction Options - -The `DegreeTwoReduction` has an additional option to control whether buses with reactive power injections are reduced: - -```@repl tutorial_VirtualPTDF_matrix -# Preserve buses with reactive power injections -reduction = DegreeTwoReduction(; reduce_reactive_power_injectors = false) -v_ptdf_preserve_reactive = - VirtualPTDF(sys_2; network_reductions = NetworkReduction[reduction]); -``` - -### Accessing Reduction Information - -After initializing the VirtualPTDF with reductions, you can access information about what was reduced: - -```@repl tutorial_VirtualPTDF_matrix -v_ptdf_reduced = - VirtualPTDF( - sys_2; - network_reductions = NetworkReduction[RadialReduction(), DegreeTwoReduction()], - ); - -# Get the reduction data -reduction_data = get_network_reduction_data(v_ptdf_reduced) -``` - -### Combining Reductions with Other Options - -Network reductions can be combined with other VirtualPTDF options like distributed slack and sparsification: - -```@repl tutorial_VirtualPTDF_matrix -v_ptdf_full_options = VirtualPTDF(sys_2; - dist_slack = dis_slack_dict, - tol = 1e-5, - network_reductions = NetworkReduction[RadialReduction(), DegreeTwoReduction()], -); -``` - -**NOTE**: The `VirtualPTDF` only supports `KLU` and `AppleAccelerate` linear solvers when using network reductions. The reference (slack) bus is automatically protected from elimination during reductions. diff --git a/src/AdjacencyMatrix.jl b/src/AdjacencyMatrix.jl index 28efbd06a..dcefb8a5b 100644 --- a/src/AdjacencyMatrix.jl +++ b/src/AdjacencyMatrix.jl @@ -7,7 +7,7 @@ This matrix describes the directed connectivity between buses, where non-zero en electrical connections through transmission lines, transformers, or other network elements. The matrix is indexed using bus numbers, which do not need to be sequential. Each element -`A[i,j]` is non-zero if there is a direct electrical connection between bus `i` and bus `j`. +``A[i,j]`` is non-zero if there is a direct electrical connection between bus ``i`` and bus ``j``. Diagonal elements are typically zero since self-loops are not meaningful in power network topology. # Fields diff --git a/src/BA_ABA_matrices.jl b/src/BA_ABA_matrices.jl index 6036e7fef..181431c4a 100644 --- a/src/BA_ABA_matrices.jl +++ b/src/BA_ABA_matrices.jl @@ -2,7 +2,10 @@ Structure containing the BA matrix and related network topology data. The BA matrix represents the branch-bus incidence matrix weighted by branch susceptances, -computed as the product of the incidence matrix A and the susceptance matrix B. +computed as the product ``B A`` of the incidence matrix ``A`` (the [`IncidenceMatrix`](@ref)) +and the susceptance matrix ``B`` — the diagonal matrix of branch series susceptances +(``b = 1/x`` under the DC approximation, or ``b = 1/(a x)`` for a branch with off-nominal +tap ratio ``a``). # Fields - `data::SparseArrays.SparseMatrixCSC{Float64, Int}`: @@ -45,7 +48,7 @@ stores_transpose(::BA_Matrix) = true """ BA_Matrix(sys::PSY.System; network_reductions::Vector{NetworkReduction} = Vector{NetworkReduction}(), kwargs...) -Construct a BA_Matrix from a PowerSystems.System by first building the underlying Ybus matrix +Construct a `BA_Matrix` from a `PSY.System` by first building the underlying [`Ybus`](@ref) matrix and then computing the branch-bus incidence matrix weighted by branch susceptances. # Arguments @@ -58,14 +61,14 @@ and then computing the branch-bus incidence matrix weighted by branch susceptanc Whether to include constant impedance loads as shunt admittances in the network model - `subnetwork_algorithm=iterative_union_find`: Algorithm used for identifying electrical islands and connected components -- Additional keyword arguments are passed to the underlying `Ybus` constructor +- Additional keyword arguments are passed to the underlying [`Ybus`](@ref) constructor # Returns - `BA_Matrix`: The constructed BA matrix structure containing the transposed branch-bus incidence matrix weighted by susceptances, along with network topology information # Notes -- This constructor creates a `Ybus` matrix internally and then converts it to a `BA_Matrix` +- This constructor creates a [`Ybus`](@ref) matrix internally and then converts it to a [`BA_Matrix`](@ref) - Network reductions can significantly improve computational efficiency for large systems - The resulting matrix supports DC power flow calculations and sensitivity analysis """ @@ -85,7 +88,7 @@ end """ BA_Matrix(ybus::Ybus) -Construct a BA_Matrix from a Ybus matrix. +Construct a `BA_Matrix` from a [`Ybus`](@ref) matrix. # Arguments - `ybus::Ybus`: The Ybus matrix from which to construct the BA matrix @@ -175,9 +178,10 @@ end """ Structure containing the ABA matrix and related power system analysis data. -The ABA matrix represents the bus susceptance matrix computed as A^T * B * A, where A is the -incidence matrix and B is the branch susceptance matrix. This matrix is fundamental for DC -power flow analysis, sensitivity calculations, and linear power system studies. +The ABA matrix represents the bus susceptance matrix computed as ``A^\\top B A``, where ``A`` is the +incidence matrix (the [`IncidenceMatrix`](@ref)) and ``B`` the branch susceptance matrix (see +[`BA_Matrix`](@ref)). This matrix is fundamental for DC power flow analysis, sensitivity +calculations, and linear power system studies. # Fields - `data::SparseArrays.SparseMatrixCSC{Float64, Int}`: @@ -197,14 +201,14 @@ power flow analysis, sensitivity calculations, and linear power system studies. Container for network reduction information applied during matrix construction # Mathematical Properties -- **Matrix Form**: ABA = A^T * B * A (bus susceptance matrix) -- **Dimensions**: (n_buses - n_ref) × (n_buses - n_ref) +- **Matrix Form**: ``\\mathrm{ABA} = A^\\top B A`` (bus susceptance matrix) +- **Dimensions**: `(n_buses - n_ref) × (n_buses - n_ref)` - **Symmetry**: Positive definite symmetric matrix (for connected networks) - **Sparsity**: Inherits sparsity pattern from network topology # Notes - Reference buses are excluded from the matrix to ensure invertibility -- Factorization enables efficient solving of linear systems Ax = b +- Factorization enables efficient solving of linear systems ``\\mathrm{ABA}\\, \\theta = P`` - Used primarily for DC power flow analysis and power system sensitivity studies - Supports various network reduction techniques for computational efficiency """ @@ -233,8 +237,8 @@ get_bus_lookup(M::ABA_Matrix) = M.lookup[1] """ ABA_Matrix(sys::PSY.System; factorize::Bool = false, network_reductions::Vector{NetworkReduction} = NetworkReduction[], kwargs...) -Construct an ABA_Matrix from a PowerSystems.System by computing A^T * B * A where A is the -incidence matrix and B is the branch susceptance matrix. The resulting matrix is fundamental +Construct an `ABA_Matrix` from a `PSY.System` by computing ``A^\\top B A``, where ``A`` is the +incidence matrix and ``B`` the branch susceptance matrix. The resulting matrix is fundamental for DC power flow analysis and power system sensitivity studies. # Arguments @@ -249,7 +253,7 @@ for DC power flow analysis and power system sensitivity studies. Whether to include constant impedance loads as shunt admittances in the network model - `subnetwork_algorithm=iterative_union_find`: Algorithm used for identifying electrical islands and connected components -- Additional keyword arguments are passed to the underlying `Ybus` constructor +- Additional keyword arguments are passed to the underlying [`Ybus`](@ref) constructor # Returns - `ABA_Matrix`: The constructed ABA matrix structure containing: @@ -261,7 +265,7 @@ for DC power flow analysis and power system sensitivity studies. 1. **Ybus Construction**: Creates admittance matrix from system data 2. **Incidence Matrix**: Computes bus-branch incidence matrix A 3. **BA Matrix**: Forms branch susceptance weighted incidence matrix -4. **ABA Computation**: Calculates A^T * B * A (bus susceptance matrix) +4. **ABA Computation**: Calculates ``A^\\top B A`` (bus susceptance matrix) 5. **Reference Bus Removal**: Excludes reference buses for invertibility 6. **Optional Factorization**: Performs KLU decomposition if requested @@ -287,8 +291,8 @@ end """ ABA_Matrix(ybus::Ybus; factorize::Bool = false) -Construct an ABA_Matrix from a Ybus matrix by computing A^T * B * A where A is the -incidence matrix and B is the branch susceptance matrix. The resulting matrix is fundamental +Construct an `ABA_Matrix` from a [`Ybus`](@ref) matrix by computing ``A^\\top B A``, where ``A`` is +the incidence matrix and ``B`` the branch susceptance matrix. The resulting matrix is fundamental for DC power flow analysis and power system sensitivity studies. Network reductions can be passed via the computed Ybus matrix. @@ -308,7 +312,7 @@ via the computed Ybus matrix. # Mathematical Process 1. **Incidence Matrix**: Computes bus-branch incidence matrix A (from Ybus matrix) 2. **BA Matrix**: Forms branch susceptance weighted incidence matrix -3. **ABA Computation**: Calculates A^T * B * A (bus susceptance matrix) +3. **ABA Computation**: Calculates ``A^\\top B A`` (bus susceptance matrix) 4. **Reference Bus Removal**: Excludes reference buses for invertibility 5. **Optional Factorization**: Performs KLU decomposition if requested diff --git a/src/BranchesParallel.jl b/src/BranchesParallel.jl index 3e90a7838..e1494ed7a 100644 --- a/src/BranchesParallel.jl +++ b/src/BranchesParallel.jl @@ -1,8 +1,35 @@ +""" + AbstractBranchesParallel <: PSY.ACTransmission + +Internal supertype for a parallel group: two or more branches spanning the **same** +bus pair (a double / multi circuit), produced when a network reduction folds +parallel branches into one equivalent arc. Subtypes are [`BranchesParallel`](@ref) +(homogeneous member types) and [`MixedBranchesParallel`](@ref) (mixed member types), +both `PSY.ACTransmission` so a group can stand in for a real branch. Not exported. + +The equivalent series susceptance of a parallel group is the sum of its members' +series susceptances. Three rating aggregations are available, each answering a +different operational question: [`get_sum_of_max_rating`](@ref), +[`get_single_element_contingency_rating`](@ref), and +[`get_impedance_averaged_rating`](@ref). +""" abstract type AbstractBranchesParallel <: PSY.ACTransmission end -# `arc_key` is the group's canonical arc in original bus numbers (the seed member's -# orientation). It is remapped with `nr` on read, so orientation no longer depends on the -# order of `branches`. +""" + BranchesParallel{T<:PSY.ACTransmission} <: AbstractBranchesParallel + +Homogeneous parallel group: every member has the same concrete branch type `T` (the +inner constructor errors if `T` is not concrete — use [`MixedBranchesParallel`](@ref) +for mixed types). Not exported; produced by network reduction, not constructed by +users. `BranchesParallel(branches)` infers `arc_key` from the first member. + +# Fields +- `branches::Vector{T}`: the parallel member branches. +- `arc_key::Tuple{Int, Int}`: the group's canonical arc in original bus numbers (the + seed member's orientation), remapped through the [`NetworkReductionData`](@ref) on + read so orientation does not depend on the order of `branches`. +- `equivalent_ybus`: cached 2×2 equivalent admittance block; `nothing` until populated. +""" mutable struct BranchesParallel{T <: PSY.ACTransmission} <: AbstractBranchesParallel branches::Vector{T} arc_key::Tuple{Int, Int} @@ -27,6 +54,19 @@ function BranchesParallel(branches::Vector{T}) where {T <: PSY.ACTransmission} return BranchesParallel{T}(branches, get_arc_tuple(first(branches)), nothing) end +""" + MixedBranchesParallel <: AbstractBranchesParallel + +Heterogeneous parallel group: members of differing concrete branch types held under +the abstract element type `PSY.ACTransmission`. The counterpart to the homogeneous +[`BranchesParallel`](@ref). Not exported; produced by network reduction. + +# Fields +- `branches::Vector{PSY.ACTransmission}`: the parallel member branches. +- `arc_key::Tuple{Int, Int}`: canonical arc in original bus numbers (see + [`BranchesParallel`](@ref)). +- `equivalent_ybus`: cached 2×2 equivalent admittance block; `nothing` until populated. +""" mutable struct MixedBranchesParallel <: AbstractBranchesParallel branches::Vector{PSY.ACTransmission} arc_key::Tuple{Int, Int} diff --git a/src/BranchesSeries.jl b/src/BranchesSeries.jl index b12d55274..aa8733e24 100644 --- a/src/BranchesSeries.jl +++ b/src/BranchesSeries.jl @@ -1,3 +1,25 @@ +""" + BranchesSeries <: PSY.ACTransmission + +A chain of branches connected in series through eliminated degree-2 buses, as +produced by [`DegreeTwoReduction`](@ref). Members are bucketed by concrete type; a +member may itself be a parallel group, so a series chain can nest parallel blocks. +Subtypes `PSY.ACTransmission` so the chain can stand in for a real branch. Not +exported. `BranchesSeries()` builds an empty chain; `add_branch!(bs, branch, orientation)` +appends a segment with its `:FromTo` / `:ToFrom` orientation. + +The equivalent series susceptance is the reciprocal of the sum of member reciprocal +susceptances. The chain rating is set by its weakest link — the minimum member +rating, where a nested parallel member contributes its N-1 +[`get_single_element_contingency_rating`](@ref). + +# Fields +- `branches::Dict{DataType, Vector{<:PSY.ACTransmission}}`: members bucketed by concrete type. +- `needs_insertion_order::Bool`: `true` when the chain mixes types and needs `insertion_order`. +- `insertion_order::Vector{Tuple{DataType, Int}}`: physical ordering along the chain. +- `segment_orientations::Vector{Symbol}`: per-segment `:FromTo` / `:ToFrom` orientation. +- `equivalent_ybus`: cached 2×2 equivalent admittance block; `nothing` until populated. +""" mutable struct BranchesSeries <: PSY.ACTransmission branches::Dict{DataType, Vector{<:PSY.ACTransmission}} needs_insertion_order::Bool diff --git a/src/IncidenceMatrix.jl b/src/IncidenceMatrix.jl index e1a3e1da5..fa98389fe 100644 --- a/src/IncidenceMatrix.jl +++ b/src/IncidenceMatrix.jl @@ -4,15 +4,15 @@ Structure containing the network incidence matrix and related topology data. The incidence matrix A represents the bus-branch connectivity of the power network, where each row corresponds to a branch and each column corresponds to a bus. Elements are: - +1 for the "from" bus of a branch -- -1 for the "to" bus of a branch +- -1 for the "to" bus of a branch - 0 for buses not connected to the branch # Fields - `data::SparseArrays.SparseMatrixCSC{Int8, Int}`: - The incidence matrix data with dimensions (n_branches × n_buses). Values are {-1, 0, +1} + The incidence matrix data with dimensions `(n_branches × n_buses)`. Values are {-1, 0, +1} representing the directed connectivity between branches and buses - `axes::Ax`: - Tuple containing (arc_identifiers, bus_numbers) where arcs are branch endpoint pairs + Tuple containing `(arc_identifiers, bus_numbers)` where arcs are branch endpoint pairs and buses are the network bus numbers - `lookup::L <: NTuple{2, Dict}`: Tuple of dictionaries providing fast lookup from arc/bus identifiers to matrix indices @@ -22,14 +22,14 @@ each row corresponds to a branch and each column corresponds to a bus. Elements Container for network reduction information applied during matrix construction # Mathematical Properties -- **Matrix Dimensions**: (n_branches × n_buses) +- **Matrix Dimensions**: `(n_branches × n_buses)` - **Element Values**: {-1, 0, +1} representing directed branch-bus connectivity - **Row Sum**: Each row sums to zero (conservation at branch level) -- **Rank**: Rank is (n_buses - n_islands) for connected networks +- **Rank**: Rank is `(n_buses - n_islands)` for connected networks - **Sparsity**: Very sparse with exactly 2 non-zero elements per branch row # Applications -- **Power Flow**: Forms the basis for DC power flow equations: P = A^T * f +- **Power Flow**: Forms the basis for DC power flow equations: ``P = A^\\top f`` - **Sensitivity Analysis**: Used in PTDF and LODF calculations - **Network Analysis**: Identifies connected components and network structure - **Topology Processing**: Enables network reduction and equivalencing algorithms @@ -38,7 +38,7 @@ each row corresponds to a branch and each column corresponds to a bus. Elements - Each branch contributes exactly one row with two non-zero entries (+1, -1) - Reference buses are preserved in the matrix but identified separately - Supports various network reduction techniques for computational efficiency -- Essential building block for BA_Matrix and ABA_Matrix constructions +- Essential building block for [`BA_Matrix`](@ref) and [`ABA_Matrix`](@ref) constructions """ struct IncidenceMatrix{Ax <: NTuple{2, Vector}, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{Int8} @@ -103,7 +103,7 @@ and creating the bus-branch connectivity matrix fundamental to power system anal Whether to include constant impedance loads as shunt admittances in the network model - `subnetwork_algorithm=iterative_union_find`: Algorithm used for identifying electrical islands and connected components -- Additional keyword arguments are passed to the underlying `Ybus` constructor +- Additional keyword arguments are passed to the underlying [`Ybus`](@ref) constructor # Returns - `IncidenceMatrix`: The constructed incidence matrix structure containing: @@ -119,8 +119,8 @@ and creating the bus-branch connectivity matrix fundamental to power system anal 5. **Network Reductions**: Applies specified reduction algorithms if provided # Applications -- **Foundation Matrix**: Essential for constructing BA_Matrix and ABA_Matrix -- **DC Power Flow**: Enables linearized power flow analysis through P = A^T * f +- **Foundation Matrix**: Essential for constructing [`BA_Matrix`](@ref) and [`ABA_Matrix`](@ref) +- **DC Power Flow**: Enables linearized power flow analysis through ``P = A^\\top f`` - **Sensitivity Analysis**: Required for PTDF, LODF, and other sensitivity calculations - **Network Analysis**: Supports topology processing and network equivalencing @@ -168,8 +168,8 @@ structure already captured in the Ybus matrix. 5. **Metadata Transfer**: Preserves reference bus positions and network reduction information # Mathematical Properties -- **Matrix Form**: A[i,j] = +1 if branch i originates at bus j, -1 if it terminates at bus j, 0 otherwise -- **Dimensions**: (n_branches × n_buses) including all network branches and buses +- **Matrix Form**: ``A[i,j] = +1`` if branch ``i`` originates at bus ``j``, ``-1`` if it terminates at bus ``j``, ``0`` otherwise +- **Dimensions**: `(n_branches × n_buses)` including all network branches and buses - **Sparsity**: Exactly 2 non-zero entries per branch row (except for isolated buses) - **Consistency**: Maintains the same network topology and reduction state as the source Ybus @@ -177,7 +177,7 @@ structure already captured in the Ybus matrix. - This constructor is more efficient when a Ybus matrix is already available - Preserves all network reduction information from the source matrix - Isolated buses are handled explicitly to maintain network completeness -- Essential for creating downstream matrices like BA_Matrix and ABA_Matrix from existing Ybus +- Essential for creating downstream matrices like [`BA_Matrix`](@ref) and [`ABA_Matrix`](@ref) from an existing Ybus """ function IncidenceMatrix(ybus::Ybus) nr = ybus.network_reduction_data diff --git a/src/KLUWrapper/KLUWrapper.jl b/src/KLUWrapper/KLUWrapper.jl index 7d50b5914..8bbfa269c 100644 --- a/src/KLUWrapper/KLUWrapper.jl +++ b/src/KLUWrapper/KLUWrapper.jl @@ -12,7 +12,7 @@ designed for the access patterns of `PowerNetworkMatrices`: This module is intentionally lighter than `KLU.jl`: it owns no Julia-side copies of the matrix values, exposes the symbolic/numeric split directly, and -binds only the SuiteSparse_long (`klu_l_*`, `klu_zl_*`) entry points used by +binds only the `SuiteSparse_long` (`klu_l_*`, `klu_zl_*`) entry points used by the package. """ module KLUWrapper diff --git a/src/NetworkReduction.jl b/src/NetworkReduction.jl index f7b9dff7d..24459dd48 100644 --- a/src/NetworkReduction.jl +++ b/src/NetworkReduction.jl @@ -2,13 +2,49 @@ NetworkReduction Abstract base type for all network reduction algorithms used in power network analysis. -Network reductions are mathematical transformations that eliminate buses and branches +Network reductions are mathematical transformations that eliminate buses and branches while preserving the electrical behavior of the remaining network elements. Concrete implementations include: - [`RadialReduction`](@ref): Eliminates radial (dangling) buses and branches - [`DegreeTwoReduction`](@ref): Eliminates buses with exactly two connections - [`WardReduction`](@ref): Reduces external buses while preserving study bus behavior + +# Applying reductions + +Reductions are applied only through the `network_reductions` keyword, a +`Vector{NetworkReduction}` accepted by every matrix constructor that builds from a +`System` — [`PTDF`](@ref), [`Ybus`](@ref), [`BA_Matrix`](@ref), [`ABA_Matrix`](@ref), +[`LODF`](@ref), [`VirtualPTDF`](@ref), [`VirtualLODF`](@ref), and [`VirtualMODF`](@ref). +The specs are applied in vector order; the default (an empty vector) applies no +reduction. + +```julia +ptdf = PTDF(sys; network_reductions = NetworkReduction[RadialReduction(), DegreeTwoReduction()]) +``` + +Write the vector with the `NetworkReduction[...]` element-type prefix: a bare +`[RadialReduction()]` infers the narrower `Vector{RadialReduction}`, which the keyword +(typed `Vector{NetworkReduction}`) will not accept. The prefix is unnecessary only when +the vector already holds two or more different spec types. + +# Ordering and validation rules + +The applied specs are validated at construction (a violation throws or warns): + +- each reduction **type** may appear at most once; +- [`WardReduction`](@ref) must be **last** when present; +- `ZeroImpedanceBranchReduction` may not be listed — it is auto-applied during + [`Ybus`](@ref) construction; +- a [`DegreeTwoReduction`](@ref) placed before a [`RadialReduction`](@ref) warns, since + running radial first usually exposes more degree-two buses for the second pass. + +# Reading back what changed + +The applied reductions are recorded on the matrix; retrieve the record with +[`get_network_reduction_data`](@ref) and inspect it through the +[`NetworkReductionData`](@ref) accessors (e.g. `get_removed_buses`, `get_removed_arcs`, +and the `keys(get_bus_reduction_map(nrd))` bus-survival check). """ abstract type NetworkReduction end diff --git a/src/NetworkReductionData.jl b/src/NetworkReductionData.jl index f7c825636..e69685a9a 100644 --- a/src/NetworkReductionData.jl +++ b/src/NetworkReductionData.jl @@ -131,9 +131,15 @@ network reduction algorithms. - `all_branch_maps_by_type::BranchMapsByType`: Branch mappings organized by component type - `reductions::ReductionContainer`: Container tracking applied reduction algorithms - `name_to_arc_map::Dict{Type, DataStructures.SortedDict{String, Tuple{Tuple{Int, Int}, String}}}`: Lazily filled with the call to [`populate_branch_maps_by_type!`](@ref), maps string names to their corresponding arcs and the map where the arc can be found. -- `component_to_reduction_name_map::Dict{Type, Dict{String, String}}`: Lazily filled with the call to [`populate_branch_maps_by_type!`](@ref), maps component names to the names of the reduction entries used in name_to_arc_map. +- `component_to_reduction_name_map::Dict{Type, Dict{String, String}}`: Lazily filled with the call to [`populate_branch_maps_by_type!`](@ref), maps component names to the names of the reduction entries used in `name_to_arc_map`. - `filters_applied::Dict{Type, Function}`: Filters applied when populating branch maps by type - `direct_branch_name_map::Dict{String, Tuple{Int, Int}}`: Lazily filled, maps branch names to their corresponding arc tuples for direct branches + +Each field has a like-named `get_*` accessor (e.g. `get_irreducible_buses`, +`get_removed_buses`); [`get_bus_reduction_map`](@ref) and [`get_reductions`](@ref) are +exported, the rest are internal. To test whether a bus survived a reduction, use +`bus in keys(get_bus_reduction_map(nrd))` — radial and degree-two survivors appear as +keys, and Ward survivors are the `study_buses`. """ @kwdef mutable struct NetworkReductionData irreducible_buses::Set{Int} = Set{Int}() # Buses that are not reduced in the network reduction @@ -500,7 +506,7 @@ function Base.empty!(rb::NetworkReductionData) end """ - get_retained_branches_names(network_reduction_data::NetworkReductionData) + get_retained_branches_names(network_reduction_data::NetworkReductionData) Gets the branch names that are retained after network reduction. This method only returns the branch names from non-three winding transformer branches that have a one-to-one correspondence with @@ -521,7 +527,7 @@ function get_retained_branches_names(network_reduction_data::NetworkReductionDat end """ - get_ac_transmission_types(network_reduction_data::NetworkReductionData) + get_ac_transmission_types(network_reduction_data::NetworkReductionData) Gets the concrete types of all AC transmission branches included in an instance of NetworkReductionData diff --git a/src/Ybus.jl b/src/Ybus.jl index c8e6a0be1..b4dfe589b 100644 --- a/src/Ybus.jl +++ b/src/Ybus.jl @@ -8,7 +8,7 @@ electrical parameters needed for power flow calculations and network analysis. # Fields - `data::SparseArrays.SparseMatrixCSC{YBUS_ELTYPE, Int}`: Sparse Y-bus matrix with complex admittance values - `adjacency_data::SparseArrays.SparseMatrixCSC{Int8, Int}`: Network connectivity information -- `axes::Ax`: Tuple of bus axis vectors for indexing (bus_numbers, bus_numbers) +- `axes::Ax`: Tuple of bus axis vectors for indexing `(bus_numbers, bus_numbers)` - `lookup::L`: Tuple of lookup dictionaries mapping bus numbers to matrix indices - `subnetwork_axes::Dict{Int, Ax}`: Bus axes for each electrical island/subnetwork - `arc_subnetwork_axis::Dict{Int, Vector{Tuple{Int, Int}}}`: Arc axes for each subnetwork @@ -25,7 +25,7 @@ electrical parameters needed for power flow calculations and network analysis. # Usage The Y-bus is fundamental for: -- Power flow analysis: V = Y⁻¹I +- Power flow analysis: ``V = Y^{-1} I`` - Short circuit calculations - Network impedance analysis - Sensitivity analysis (PTDF/LODF) @@ -101,7 +101,7 @@ end Build a Y-bus matrix from the system and return its default network reduction data. This function constructs a Y-bus matrix with no network reductions applied and returns -the resulting `NetworkReductionData`, which contains the basic bus and branch mappings +the resulting [`NetworkReductionData`](@ref), which contains the basic bus and branch mappings for the system without any reduction algorithms. # Arguments @@ -300,8 +300,8 @@ connecting to a virtual star bus. Each available winding is mapped separately. - `br::PSY.ThreeWindingTransformer`: Three-winding transformer to add # Implementation Details -- Only adds arcs for available windings (checked via PSY.get_available_*) -- Maintains transformer3W_map and reverse_transformer3W_map +- Only adds arcs for available windings (checked via `PSY.get_available_*`) +- Maintains `transformer3W_map` and `reverse_transformer3W_map` - Each winding is numbered (1=primary, 2=secondary, 3=tertiary) """ function add_to_branch_maps!( @@ -1146,17 +1146,17 @@ end Generate unique arc axis from from-bus and to-bus index vectors. -Creates a vector of unique (from_bus, to_bus) tuples representing the arcs (branches) +Creates a vector of unique `(from_bus, to_bus)` tuples representing the arcs (branches) in the system. Used for constructing arc admittance matrices and organizing network topology data. # Arguments -- `fb::Vector{Int}`: Vector of from-bus indices into bus_axis -- `tb::Vector{Int}`: Vector of to-bus indices into bus_axis +- `fb::Vector{Int}`: Vector of from-bus indices into `bus_axis` +- `tb::Vector{Int}`: Vector of to-bus indices into `bus_axis` - `bus_axis::Vector{Int}`: Vector of bus numbers # Returns -- `Vector{Tuple{Int, Int}}`: Unique arcs as (from_bus_number, to_bus_number) tuples +- `Vector{Tuple{Int, Int}}`: Unique arcs as `(from_bus_number, to_bus_number)` tuples # Examples ```julia @@ -1190,7 +1190,7 @@ corresponding arc list for matrix indexing. - `ybus::Ybus`: Y-bus matrix containing subnetwork information # Returns -- `Dict{Int, Tuple{Vector{Int}, Vector{Tuple{Int, Int}}}}`: Dictionary mapping reference bus numbers to (bus_axis, arc_axis) tuples for each subnetwork +- `Dict{Int, Tuple{Vector{Int}, Vector{Tuple{Int, Int}}}}`: Dictionary mapping reference bus numbers to `(bus_axis, arc_axis)` tuples for each subnetwork # Implementation Details - Combines bus axes from `ybus.subnetwork_axes` with arc axes from `ybus.arc_subnetwork_axis` @@ -1226,7 +1226,7 @@ list for matrix indexing. - `ybus::Ybus`: Y-bus matrix containing subnetwork information # Returns -- `Dict{Int, Tuple{Vector{Tuple{Int, Int}}, Vector{Int}}}`: Dictionary mapping reference bus numbers to (arc_axis, bus_axis) tuples for each subnetwork +- `Dict{Int, Tuple{Vector{Tuple{Int, Int}}, Vector{Int}}}`: Dictionary mapping reference bus numbers to `(arc_axis, bus_axis)` tuples for each subnetwork # Implementation Details - Swaps order compared to `make_bus_arc_subnetwork_axes` (arc first, bus second) diff --git a/src/apply_zero_impedance_reduction.jl b/src/apply_zero_impedance_reduction.jl index 7780d51b4..884dedf21 100644 --- a/src/apply_zero_impedance_reduction.jl +++ b/src/apply_zero_impedance_reduction.jl @@ -59,7 +59,12 @@ function _is_zero_impedance_arc( # Transformer-bearing arcs are excluded from zero-impedance bus merging. _any_transformer(parallel_br) && return false return any( - _is_zero_impedance_branch(br, susceptance_threshold, min_x_eps, resistance_tolerance) + _is_zero_impedance_branch( + br, + susceptance_threshold, + min_x_eps, + resistance_tolerance, + ) for br in parallel_br ) end diff --git a/src/auto_tolerance.jl b/src/auto_tolerance.jl index ca3b29e12..37aeb04f0 100644 --- a/src/auto_tolerance.jl +++ b/src/auto_tolerance.jl @@ -3,25 +3,43 @@ Request automatic, condition-aware sparsification of a PTDF/LODF matrix. The matrix is sparsified with a *relative* per-row cutoff: an entry is dropped when - - |entry| < α · max|row|, α = clamp(safety · δ, $(string(1e-6)), $(string(1e-2))) - -where `δ` is the relative precision of the branch data. Because the cutoff is +```math +|\\mathrm{entry}| < \\alpha \\cdot \\max|\\mathrm{row}|, +\\qquad \\alpha = \\mathrm{clamp}(\\mathrm{safety} \\cdot \\delta, \\, 10^{-6}, \\, 10^{-2}) +``` +where ``\\delta`` is the relative precision of the branch data. Because the cutoff is relative to each row's own peak, columns of large, ill-conditioned systems stay sparse regardless of the conditioning of `ABA`; the 1-norm condition estimate of `ABA` is still computed and logged as a diagnostic, but never multiplies the cutoff. -- `data_precision`: relative precision `δ` of the branch parameters. `:auto` - (default) discovers it from the branch reactances (see +- `data_precision`: relative precision ``\\delta`` of the branch parameters. `:auto` + (default) discovers it from the branch susceptances (see [`discover_data_precision`](@ref)); a `Float64` sets it explicitly (e.g. `1e-3` for reactances good to 0.1%). - `safety`: aggressiveness multiplier on `δ`; `> 1` sparsifies more, `< 1` less. - `quantile`: only used when `data_precision = :auto`; which quantile of the per-branch significant-figure counts to adopt. -A plain `Float64` `tol` is still accepted by every constructor and applies a -fixed *absolute* cutoff (backward compatible / exact tolerance). +# Where the cutoff applies + +An `AutoTolerance` is a **no-op below `AUTO_TOLERANCE_BUS_LIMIT`** (2000 buses): small +systems and the test cases are returned exactly, and its relative drop is reserved for +the large virtual matrices ([`VirtualPTDF`](@ref) / [`VirtualLODF`](@ref) / +[`VirtualMODF`](@ref)). On the **dense** [`PTDF`](@ref) / [`LODF`](@ref) path it is also +a no-op, preserving the `Matrix{Float64}` element type. + +A plain `Float64` `tol` is accepted by every constructor and applies a fixed *absolute* +cutoff (``|\\mathrm{entry}| < \\mathrm{tol}``) at any system size, dense or virtual — the backward-compatible, +exact-tolerance path. + +# Examples + +```julia +PTDF(sys; tol = 1e-5) # fixed absolute cutoff, any size +PTDF(sys; tol = AutoTolerance(; safety = 5.0)) # sparsify large virtual matrices harder +PTDF(sys; tol = AutoTolerance(; data_precision = 1e-3)) # pin precision instead of discovering it +``` """ struct AutoTolerance{D <: Union{Float64, Symbol}} data_precision::D @@ -114,14 +132,14 @@ end """ discover_data_precision(susceptances; q = 0.5, maxsig = 10, rtol = 1e-9) -> Float64 -Estimate relative data precision from branch susceptances `b_k`. Recovers the -reactances `x_k = 1/b_k` (the susceptance hides the original precision; the +Estimate relative data precision from branch susceptances ``b_k``. Recovers the +reactances ``x_k = 1/b_k`` (the susceptance hides the original precision; the reciprocal does not), counts the significant figures of each, and returns -`0.5·10^(-(s-1))` at the `q`-quantile of those counts, clamped to `[eps, 1e-2]`. -`maxsig` is coupled to `rtol`: rounding to `s` figures carries a relative error of -`0.5·10^(1-s)`, so `rtol = 1e-9` first accepts `s = 10` and no real data resolves +``0.5 \\cdot 10^{-(s-1)}`` at the ``q``-quantile of those counts, clamped to `[eps, 1e-2]`. +`maxsig` is coupled to `rtol`: rounding to ``s`` figures carries a relative error of +``0.5 \\cdot 10^{1-s}``, so `rtol = 1e-9` first accepts ``s = 10`` and no real data resolves further. Full-precision data (e.g. computed equivalent branches) hits this `maxsig` -cap and yields the tightest precision `5e-10`. +cap and yields the tightest precision ``5 \\times 10^{-10}``. """ function discover_data_precision( susceptances::AbstractVector{Float64}; diff --git a/src/common.jl b/src/common.jl index 03ec26dc3..a2824aa1e 100644 --- a/src/common.jl +++ b/src/common.jl @@ -214,7 +214,8 @@ end """ Validates that the user bus input is consistent with the ybus axes and the prior reductions. -Is used to check `irreducible_buses` for `Radial` and `DegreeTwo` reductions and `study_buses` for `WardReduction`. +Is used to check `irreducible_buses` for [`RadialReduction`](@ref) and +[`DegreeTwoReduction`](@ref), and `study_buses` for [`WardReduction`](@ref). """ function validate_buses(A::PowerNetworkMatrix, buses::Set{Int}) reverse_bus_search_map = get_network_reduction_data(A).reverse_bus_search_map @@ -259,8 +260,8 @@ reference bus positions. BA matrix. NOTE: -- evaluates A with "calculate_A_matrix", or extract A.data (if A::IncidenceMatrix) -- evaluates BA with "calculate_BA_matrix", or extract BA.data (if A::BA_Matrix) +- evaluates A with `calculate_A_matrix`, or extracts `A.data` (if `A::IncidenceMatrix`) +- evaluates BA with `calculate_BA_matrix`, or extracts `BA.data` (if `BA::BA_Matrix`) """ function calculate_ABA_matrix( A::SparseArrays.SparseMatrixCSC{Int8, Int}, @@ -378,6 +379,16 @@ function populate_equivalent_ybus!( return end +""" + get_equivalent_physical_branch_parameters(segment, nr::NetworkReductionData) -> EquivalentBranch + +Physical branch parameters equivalent to a reduced group of branches. Lazily builds the +group's 2×2 equivalent Ybus block (via `populate_equivalent_ybus!`) if not already +cached, then decomposes that block into an [`EquivalentBranch`](@ref) — series `r`/`x`, +from/to shunt `g`/`b`, and the transformer `tap`/`shift` — so a series chain or set of +parallel branches can be represented as a single Pi-model branch. `segment` is a +[`BranchesSeries`](@ref) or [`AbstractBranchesParallel`](@ref) reduction group. +""" function get_equivalent_physical_branch_parameters( segment::Union{AbstractBranchesParallel, BranchesSeries}, nr::NetworkReductionData, @@ -488,7 +499,7 @@ end -> Tuple{Symbol, Union{Tuple{Int, Int}, Nothing}} Classify a branch component by looking up which reverse map it belongs to in the -`NetworkReductionData`. Returns `(tag, arc_tuple)` where `tag` is one of: +[`NetworkReductionData`](@ref). Returns `(tag, arc_tuple)` where `tag` is one of: - `:direct` -- branch is the sole branch on its arc - `:parallel` -- branch is one of several parallel branches on its arc - `:series` -- branch is part of a series chain on its arc @@ -587,15 +598,18 @@ end Compute the change in equivalent arc susceptance when multiple components are simultaneously tripped from a series chain. -For a series chain with segments of susceptance b₁, b₂, ..., bₙ, the equivalent -susceptance is: b_eq = 1 / (1/b₁ + 1/b₂ + ... + 1/bₙ). +For a series chain with segments of susceptance ``b_1, b_2, \\ldots, b_n``, the +equivalent susceptance is +```math +b_\\mathrm{eq} = \\frac{1}{1/b_1 + 1/b_2 + \\cdots + 1/b_n}. +``` Segments can be individual branches or `BranchesParallel` groups. When a tripped component is inside a parallel group, only that branch's susceptance is removed from the group — the rest of the parallel group remains in the series chain. -Returns Δb = b_new - b_old (always negative for outages). -If all segments are fully tripped, returns -b_eq (full arc outage). +Returns ``\\Delta b = b_\\mathrm{new} - b_\\mathrm{old}`` (always negative for outages). +If all segments are fully tripped, returns ``-b_\\mathrm{eq}`` (full arc outage). """ function _compute_series_outage_delta_b( series_chain::BranchesSeries, diff --git a/src/connectivity_checks.jl b/src/connectivity_checks.jl index 3ffe5e3da..40f5ba64f 100644 --- a/src/connectivity_checks.jl +++ b/src/connectivity_checks.jl @@ -245,7 +245,8 @@ a the ABA or Adjacency Matrix. - `bus_numbers::Vector{Int}`: vector containing the indices of the system's buses. - `subnetwork_algorithm::Function`: - algorithm for computing subnetworks. Valid options are iterative_union_find (default) and depth_first_search + algorithm for computing subnetworks. Valid options are [`iterative_union_find`](@ref) + (default) and [`depth_first_search`](@ref) """ function find_subnetworks( M::SparseArrays.SparseMatrixCSC, diff --git a/src/degree_two_reduction.jl b/src/degree_two_reduction.jl index 04662fb9e..3ff42ba53 100644 --- a/src/degree_two_reduction.jl +++ b/src/degree_two_reduction.jl @@ -252,7 +252,7 @@ end reduced_indices::Set{Int}, irreducible_indices::Set{Int}) -Recursively build a chain in one direction from current_node, avoiding prev_node. +Recursively build a chain in one direction from `current_node`, avoiding `prev_node`. """ function _get_partial_chain_recursive!( current_chain::Vector{Int}, diff --git a/src/linalg_settings.jl b/src/linalg_settings.jl index 4d2efb265..edd250fdc 100644 --- a/src/linalg_settings.jl +++ b/src/linalg_settings.jl @@ -91,7 +91,7 @@ set_linalg_backend_preference(linalglib::Symbol) = get_linalg_backend_preference() = Preferences.@load_preference("linalg_backend") -"Set a preference whether to run check_linalg_backend at the package loading time." +"Set a preference whether to run `check_linalg_backend` at the package loading time." set_linalg_backend_check(check::Bool) = Preferences.@set_preferences!("linalg_backend_check" => check) @@ -99,7 +99,7 @@ get_linalg_backend_check() = Preferences.@load_preference("linalg_backend_check" function check_lbt_library() lb_msg(lib) = """The $(lib) library is being used for Julia's BLAS and LAPACK routines, - including dense linear algebra operations such as `BLAS.gemm``.""" + including dense linear algebra operations such as `BLAS.gemm`.""" blas_config = lowercase(string(LinearAlgebra.BLAS.get_config())) if contains(blas_config, "mkl") diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl index c1f274dd6..e4739913e 100644 --- a/src/lodf_calculations.jl +++ b/src/lodf_calculations.jl @@ -2,9 +2,9 @@ Structure containing the Line Outage Distribution Factor (LODF) matrix and related power system data. The LODF matrix contains sensitivity coefficients that quantify how the outage of one transmission -line affects the power flows on all other lines in the system. Each element LODF[i,j] represents -the change in flow on line i when line j is taken out of service, normalized by the pre-outage -flow on line j. +line affects the power flows on all other lines in the system. Each element ``\\mathrm{LODF}[i,j]`` +represents the change in flow on line ``i`` when line ``j`` is taken out of service, normalized by +the pre-outage flow on line ``j``. # Fields - `data::M <: AbstractArray{Float64, 2}`: @@ -22,8 +22,8 @@ flow on line j. Container for network reduction information applied during matrix construction # Mathematical Properties -- **Matrix Form**: LODF[i,j] = ∂f_i/∂P_j where f_i is flow on line i, P_j is injection change due to line j outage -- **Dimensions**: (n_branches × n_branches) for all transmission lines in the system +- **Matrix Form**: ``\\mathrm{LODF}[i,j] = \\partial f_i / \\partial P_j`` where ``f_i`` is flow on line ``i`` and ``P_j`` is the injection change due to the outage of line ``j`` +- **Dimensions**: `(n_branches × n_branches)` for all transmission lines in the system - **Diagonal Elements**: Always -1 (100% flow reduction on the outaged line itself) - **Symmetry**: Generally non-symmetric matrix reflecting directional flow sensitivities - **Physical Meaning**: Values represent fraction of pre-outage flow that redistributes to other lines @@ -203,7 +203,8 @@ end Computes the LODF matrix using the internal Apple Accelerate backend (`AccelerateWrapper`). Available only on macOS. Shape mirrors `_calculate_LODF_matrix_KLU(a, ptdf)` exactly: factor the diagonal "demand" - matrix `diag(1 - PTDF·A)` and solve in place against `a · ptdf`. + matrix ``\\mathrm{diag}(1 - A \\, \\mathrm{PTDF})`` and solve in place against + ``A \\, \\mathrm{PTDF}``. # Arguments - `a::SparseArrays.SparseMatrixCSC{Int8, Int}`: Incidence Matrix @@ -221,8 +222,35 @@ end end end +# Numeric/default tol: original PTDF-based route, unchanged behavior. +function _lodf_from_system( + tol::Float64, + A::IncidenceMatrix, + BA::BA_Matrix, + Ymatrix::Ybus, + linear_solver::String, +) + # Keep the intermediate PTDF dense (tol = eps()); the from-PTDF LODF needs an + # unsparsified PTDF for accuracy, and only the LODF itself is sparsified. + ptdf = PTDF(A, BA; linear_solver = linear_solver, tol = eps()) + return LODF(A, ptdf; linear_solver = linear_solver, tol = tol) +end + +# AutoTolerance: build a factorized ABA so conditioning is available, then use +# the KLU-only ABA/BA constructor. +function _lodf_from_system( + spec::AutoTolerance, + A::IncidenceMatrix, + BA::BA_Matrix, + Ymatrix::Ybus, + ::String, +) + ABA = ABA_Matrix(Ymatrix; factorize = true) + return LODF(A, ABA, BA; tol = spec) +end + """ - LODF(sys::PSY.System; linear_solver::String = _default_linear_solver(), tol::Float64 = eps(), network_reductions::Vector{NetworkReduction} = NetworkReduction[], kwargs...) + LODF(sys::PSY.System; linear_solver::String = _default_linear_solver(), tol::Union{Float64, AutoTolerance} = DEFAULT_AUTO_TOLERANCE, network_reductions::Vector{NetworkReduction} = NetworkReduction[], kwargs...) Construct a Line Outage Distribution Factor (LODF) matrix from a PowerSystems.System by computing the sensitivity of line flows to single line outages. This is the primary constructor for LODF @@ -233,9 +261,11 @@ analysis starting from system data. # Keyword Arguments - `linear_solver::String = _default_linear_solver()`: - Linear solver algorithm for matrix computations. Options: "KLU", "Dense", "MKLPardiso" -- `tol::Float64 = eps()`: - Sparsification tolerance for dropping small matrix elements to reduce memory usage + Linear solver algorithm for matrix computations. Options: "KLU", "AppleAccelerateLU", "Dense", "MKLPardiso" +- `tol::Union{Float64, AutoTolerance} = DEFAULT_AUTO_TOLERANCE`: + Sparsification tolerance for dropping small matrix elements to reduce memory usage. + A `Float64` applies a fixed absolute cutoff at any size; an [`AutoTolerance`](@ref) + (the default) applies a relative per-row cutoff on large virtual matrices only. - `network_reductions::Vector{NetworkReduction} = NetworkReduction[]`: Vector of network reduction algorithms to apply before matrix construction - `include_constant_impedance_loads::Bool=true`: @@ -259,16 +289,18 @@ analysis starting from system data. 6. **Sparsification**: Applies tolerance threshold to reduce matrix density # Linear Solver Options -- **"KLU"**: Sparse LU factorization (default, recommended for most cases) +- **"KLU"**: Sparse LU factorization (default off Apple hardware, recommended for most cases) +- **"AppleAccelerateLU"**: Apple Accelerate sparse LU (default on macOS 15.5+ Apple hardware) - **"Dense"**: Dense matrix operations (faster for small systems) - **"MKLPardiso"**: Intel MKL Pardiso solver (requires MKL, best for very large systems) # Mathematical Foundation -The LODF matrix is computed using the relationship: -``` -LODF = (A * PTDF) / (1 - diag(A * PTDF)) +With ``H = A \\, \\mathrm{PTDF}``, the sensitivity of monitored line ``\\ell`` to the outage +of line ``e`` is +```math +\\mathrm{LODF}[\\ell, e] = \\frac{H[\\ell, e]}{1 - H[e, e]} ``` -where A is the incidence matrix and PTDF is the power transfer distribution factor matrix. +where ``A`` is the incidence matrix and ``\\mathrm{PTDF}`` is the power transfer distribution factor matrix. # Notes - Sparsification with `tol > eps()` can significantly reduce memory usage @@ -277,33 +309,6 @@ where A is the incidence matrix and PTDF is the power transfer distribution fact - Diagonal elements are always -1.0 representing complete flow loss on outaged lines - For very large systems, consider using "MKLPardiso" solver with appropriate chunk size """ -# Numeric/default tol: original PTDF-based route, unchanged behavior. -function _lodf_from_system( - tol::Float64, - A::IncidenceMatrix, - BA::BA_Matrix, - Ymatrix::Ybus, - linear_solver::String, -) - # Keep the intermediate PTDF dense (tol = eps()); the from-PTDF LODF needs an - # unsparsified PTDF for accuracy, and only the LODF itself is sparsified. - ptdf = PTDF(A, BA; tol = eps()) - return LODF(A, ptdf; linear_solver = linear_solver, tol = tol) -end - -# AutoTolerance: build a factorized ABA so conditioning is available, then use -# the KLU-only ABA/BA constructor. -function _lodf_from_system( - spec::AutoTolerance, - A::IncidenceMatrix, - BA::BA_Matrix, - Ymatrix::Ybus, - ::String, -) - ABA = ABA_Matrix(Ymatrix; factorize = true) - return LODF(A, ABA, BA; tol = spec) -end - function LODF( sys::PSY.System; linear_solver::String = _default_linear_solver(), @@ -320,7 +325,7 @@ function LODF( end """ - LODF(A::IncidenceMatrix, PTDFm::PTDF; linear_solver::String = _default_linear_solver(), tol::Float64 = eps()) + LODF(A::IncidenceMatrix, PTDFm::PTDF; linear_solver::String = _default_linear_solver(), tol::Union{Float64, AutoTolerance} = DEFAULT_AUTO_TOLERANCE) Construct a Line Outage Distribution Factor (LODF) matrix from existing incidence and PTDF matrices. This constructor is more efficient when the prerequisite matrices are already available. @@ -331,22 +336,23 @@ This constructor is more efficient when the prerequisite matrices are already av # Keyword Arguments - `linear_solver::String = _default_linear_solver()`: - Linear solver algorithm for matrix computations. Options: "KLU", "Dense", "MKLPardiso" -- `tol::Float64 = eps()`: + Linear solver algorithm for matrix computations. Options: "KLU", "AppleAccelerateLU", "Dense", "MKLPardiso" +- `tol::Union{Float64, AutoTolerance} = DEFAULT_AUTO_TOLERANCE`: Sparsification tolerance for the LODF matrix (not applied to input PTDF) # Returns - `LODF`: The constructed LODF matrix structure with line outage sensitivity coefficients # Mathematical Computation -The LODF matrix is computed using the formula: -``` -LODF = (A * PTDF) / (1 - diag(A * PTDF)) +With ``H = A \\, \\mathrm{PTDF}``, the sensitivity of monitored line ``\\ell`` to the outage +of line ``e`` is +```math +\\mathrm{LODF}[\\ell, e] = \\frac{H[\\ell, e]}{1 - H[e, e]} ``` where: -- A is the incidence matrix representing bus-branch connectivity -- PTDF contains power transfer distribution factors -- The denominator (1 - diagonal terms) accounts for the outaged line's own flow +- ``A`` is the incidence matrix (the [`IncidenceMatrix`](@ref)) representing bus-branch connectivity +- ``\\mathrm{PTDF}`` contains power transfer distribution factors +- The denominator ``1 - H[e,e]`` accounts for the outaged line's own flow # Important Notes - **PTDF Sparsification**: The input PTDF matrix should be non-sparsified (constructed with default tolerance) to avoid accuracy issues @@ -417,7 +423,7 @@ function LODF( end """ - LODF(A::IncidenceMatrix, ABA::ABA_Matrix, BA::BA_Matrix; linear_solver::String = "KLU", tol::Float64 = eps()) + LODF(A::IncidenceMatrix, ABA::ABA_Matrix, BA::BA_Matrix; linear_solver::String = "KLU", tol::Union{Float64, AutoTolerance} = DEFAULT_AUTO_TOLERANCE) Construct a Line Outage Distribution Factor (LODF) matrix from incidence, ABA, and BA matrices. This constructor provides direct control over the underlying matrix computations and is most @@ -425,29 +431,29 @@ efficient when the prerequisite matrices with factorization are already availabl # Arguments - `A::IncidenceMatrix`: The incidence matrix containing bus-branch connectivity information -- `ABA::ABA_Matrix`: The bus susceptance matrix (A^T * B * A), preferably with KLU factorization -- `BA::BA_Matrix`: The branch susceptance weighted incidence matrix (B * A) +- `ABA::ABA_Matrix`: The bus susceptance matrix ``A^\\top B A``, preferably with KLU factorization +- `BA::BA_Matrix`: The branch susceptance weighted incidence matrix ``B A`` # Keyword Arguments - `linear_solver::String = "KLU"`: This constructor is intentionally KLU-only because `ABA.K` is always a KLU factorization. The keyword is kept for API consistency; passing any other value will error. -- `tol::Float64 = eps()`: +- `tol::Union{Float64, AutoTolerance} = DEFAULT_AUTO_TOLERANCE`: Sparsification tolerance for dropping small matrix elements # Returns - `LODF`: The constructed LODF matrix structure with line outage sensitivity coefficients # Mathematical Computation -This method computes LODF using the factorized form: -``` -LODF = (A * ABA^(-1) * BA) / (1 - diag(A * ABA^(-1) * BA)) +This method computes LODF using the factorized form ``H = A\\, \\mathrm{ABA}^{-1} \\mathrm{BA}``, +```math +\\mathrm{LODF}[\\ell, e] = \\frac{H[\\ell, e]}{1 - H[e, e]} ``` where: -- A is the incidence matrix -- ABA^(-1) uses the factorized form from the ABA matrix (requires `ABA.K` to be factorized) -- BA is the susceptance-weighted incidence matrix +- ``A`` is the incidence matrix (the [`IncidenceMatrix`](@ref)) +- ``\\mathrm{ABA}^{-1}`` uses the factorized form from the [`ABA_Matrix`](@ref) (requires `ABA.K` to be factorized) +- ``\\mathrm{BA}`` is the susceptance-weighted incidence matrix (the [`BA_Matrix`](@ref)) # Requirements and Limitations - **Factorization Required**: The ABA matrix should be pre-factorized (contains KLU factorization) for efficiency diff --git a/src/modf_definitions.jl b/src/modf_definitions.jl index d28598891..e3d2cfb68 100644 --- a/src/modf_definitions.jl +++ b/src/modf_definitions.jl @@ -110,7 +110,7 @@ No dependency on `PSY.System` after construction. Modification vectors are converted to tuples at construction time to guarantee immutability. This is required because `NetworkModification` is used as a `Dict` -key (via custom `hash`/`==`) in VirtualMODF caches; mutable fields would +key (via custom `hash`/`==`) in [`VirtualMODF`](@ref) caches; mutable fields would silently corrupt lookups if modified after insertion. """ struct NetworkModification @@ -185,12 +185,14 @@ end WoodburyFactors Cached Woodbury intermediates shared across monitored arcs for one contingency. -Computed from van Dijk et al. Eq. 29: - B_m⁻¹ = B_r⁻¹ - B_r⁻¹ U (A⁻¹ + U⊤ B_r⁻¹ U)⁻¹ U⊤ B_r⁻¹ +Computed from van Dijk et al. Eq. 29, +```math +B_m^{-1} = B_r^{-1} - B_r^{-1} U (A^{-1} + U^\\top B_r^{-1} U)^{-1} U^\\top B_r^{-1}. +``` # Fields -- `Z::Matrix{Float64}`: B⁻¹U matrix (n_bus × M), one column per modified arc -- `W_inv::Matrix{Float64}`: Pre-inverted W = (A⁻¹ + U⊤B⁻¹U)⁻¹ (M × M). For M ≤ 2, computed analytically; for M > 2, computed via LU factorization. +- `Z::Matrix{Float64}`: ``B^{-1}U`` matrix (n_bus × M), one column per modified arc +- `W_inv::Matrix{Float64}`: Pre-inverted ``W = (A^{-1} + U^\\top B^{-1} U)^{-1}`` (M × M). For M ≤ 2, computed analytically; for M > 2, computed via LU factorization. - `arc_indices::Vector{Int}`: Arc indices of modified arcs - `delta_b::Vector{Float64}`: Susceptance changes per modified arc - `is_islanding::Bool`: Whether this contingency islands the network diff --git a/src/network_modification.jl b/src/network_modification.jl index 319b74427..de87fdcd4 100644 --- a/src/network_modification.jl +++ b/src/network_modification.jl @@ -105,7 +105,7 @@ end """ _compute_arc_ybus_delta(nr, arc_tuple, delta_b) -> NTuple{4, YBUS_ELTYPE} -Compute the Pi-model Ybus delta `(ΔY11, ΔY12, ΔY21, ΔY22)` for an arc modification by +Compute the Pi-model Ybus delta ``(\\Delta Y_{11}, \\Delta Y_{12}, \\Delta Y_{21}, \\Delta Y_{22})`` for an arc modification by dispatching to the per-map handler that owns `arc_tuple`. """ function _compute_arc_ybus_delta( @@ -140,8 +140,8 @@ end """ $(TYPEDSIGNATURES) -Construct a full arc outage `NetworkModification` by bus-pair tuple. -Looks up arc susceptance from the matrix and sets `Δb = -b_arc`. +Construct a full arc outage [`NetworkModification`](@ref) by bus-pair tuple. +Looks up arc susceptance from the matrix and sets ``\\Delta b = -b_\\mathrm{arc}``. """ function NetworkModification(mat::PowerNetworkMatrix, arc::Tuple{Int, Int}) arc_lookup = get_arc_lookup(mat) @@ -158,7 +158,7 @@ end """ $(TYPEDSIGNATURES) -Construct a `NetworkModification` from a branch component using network +Construct a [`NetworkModification`](@ref) from a branch component using network reduction reverse maps to classify the branch as direct, parallel, or series. """ function NetworkModification(mat::PowerNetworkMatrix, branch::PSY.ACTransmission) @@ -175,7 +175,7 @@ end """ $(TYPEDSIGNATURES) -Construct a `NetworkModification` from a `ThreeWindingTransformer` component. +Construct a [`NetworkModification`](@ref) from a `PSY.ThreeWindingTransformer` component. Automatically decomposes the transformer into its three winding arcs and classifies each one. For a partial outage (single winding trip), use a `ThreeWindingTransformerWinding` instead. @@ -242,7 +242,7 @@ end """ $(TYPEDSIGNATURES) -Construct a `NetworkModification` from a `PSY.Outage` supplemental attribute. +Construct a [`NetworkModification`](@ref) from a `PSY.Outage` supplemental attribute. Resolves the outage's associated `ACTransmission` components through the system, classifies each by the matrix's network reduction maps, and builds the modification. Handles multi-component outages with series-chain grouping. @@ -609,7 +609,7 @@ end """ compute_ybus_delta(ybus::Ybus, mod::NetworkModification) -> SparseMatrixCSC{YBUS_ELTYPE, Int} -Compute the sparse ΔYbus matrix from a canonical `NetworkModification`. +Compute the sparse ``\\Delta Y_\\mathrm{bus}`` matrix from a canonical `NetworkModification`. Combines arc modifications (branch outages producing Pi-model deltas) and shunt modifications (diagonal admittance changes) into a single sparse delta. @@ -657,8 +657,8 @@ end """ apply_ybus_modification(ybus::Ybus, mod::NetworkModification) -> SparseMatrixCSC -Apply a canonical NetworkModification to a Ybus, returning the modified sparse matrix. -Convenience wrapper around `compute_ybus_delta`. +Apply a canonical [`NetworkModification`](@ref) to a [`Ybus`](@ref), returning the modified +sparse matrix. Convenience wrapper around [`compute_ybus_delta`](@ref). """ function apply_ybus_modification( ybus::Ybus, diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl index ce9da613f..6128da57f 100644 --- a/src/ptdf_calculations.jl +++ b/src/ptdf_calculations.jl @@ -2,15 +2,16 @@ Structure containing the Power Transfer Distribution Factor (PTDF) matrix and related power system data. The PTDF matrix contains sensitivity coefficients that quantify how power injections at buses -affect the power flows on transmission lines. Each element PTDF[i,j] represents the incremental -change in flow on line i due to a unit power injection at bus j, under DC power flow assumptions. +affect the power flows on transmission lines. Each element ``\\mathrm{PTDF}[i,j]`` represents the +incremental change in flow on line ``i`` due to a unit power injection at bus ``j``, under DC power +flow assumptions. # Fields - `data::M <: AbstractArray{Float64, 2}`: The PTDF matrix data stored in transposed form for computational efficiency. Element (i,j) represents the sensitivity of line j flow to bus i injection - `axes::Ax`: - Tuple containing (bus_numbers, branch_identifiers) for matrix dimensions + Tuple containing `(bus_numbers, branch_identifiers)` for matrix dimensions - `lookup::L <: NTuple{2, Dict}`: Tuple of dictionaries providing fast lookup from bus/branch identifiers to matrix indices - `subnetwork_axes::Dict{Int, Ax}`: @@ -21,9 +22,9 @@ change in flow on line i due to a unit power injection at bus j, under DC power Container for network reduction information applied during matrix construction # Mathematical Properties -- **Matrix Form**: PTDF[i,j] = ∂f_i/∂P_j where f_i is flow on line i, P_j is injection at bus j -- **Dimensions**: (n_buses × n_arcs) for all buses and impedance arcs -- **Linear Superposition**: Total flow = Σ(PTDF[i,j] × P_j) for all injections P_j +- **Matrix Form**: ``\\mathrm{PTDF}[i,j] = \\partial f_i / \\partial P_j`` where ``f_i`` is flow on line ``i`` and ``P_j`` is injection at bus ``j`` +- **Dimensions**: `(n_buses × n_arcs)` for all buses and impedance arcs +- **Linear Superposition**: ``f_i = \\sum_j \\mathrm{PTDF}[i,j] \\, P_j`` over all injections ``P_j`` - **Physical Meaning**: Values represent the fraction of bus injection that flows through each line - **Reference Bus**: Rows corresponding to reference buses are typically zero @@ -285,7 +286,7 @@ end end """ - PTDF(sys::PSY.System; dist_slack::Dict{Int, Float64} = Dict{Int, Float64}(), linear_solver = _default_linear_solver(), tol::Float64 = eps(), network_reductions::Vector{NetworkReduction} = NetworkReduction[], kwargs...) + PTDF(sys::PSY.System; dist_slack::Dict{Int, Float64} = Dict{Int, Float64}(), linear_solver = _default_linear_solver(), tol::Union{Float64, AutoTolerance} = DEFAULT_AUTO_TOLERANCE, network_reductions::Vector{NetworkReduction} = NetworkReduction[], kwargs...) Construct a Power Transfer Distribution Factor (PTDF) matrix from a PowerSystems.System by computing the sensitivity of transmission line flows to bus power injections. This is the primary constructor @@ -299,9 +300,11 @@ for PTDF analysis starting from system data. Dictionary mapping bus numbers to distributed slack weights for realistic slack modeling. Empty dictionary uses single slack bus (default behavior) - `linear_solver::String = _default_linear_solver()`: - Linear solver algorithm for matrix computations. Options: "KLU", "Dense", "MKLPardiso" -- `tol::Float64 = eps()`: - Sparsification tolerance for dropping small matrix elements to reduce memory usage + Linear solver algorithm for matrix computations. Options: "KLU", "AppleAccelerateLU", "Dense", "MKLPardiso" +- `tol::Union{Float64, AutoTolerance} = DEFAULT_AUTO_TOLERANCE`: + Sparsification tolerance for dropping small matrix elements to reduce memory usage. + A `Float64` applies a fixed absolute cutoff at any size; an [`AutoTolerance`](@ref) + (the default) applies a relative per-row cutoff on large virtual matrices only. - `network_reductions::Vector{NetworkReduction} = NetworkReduction[]`: Vector of network reduction algorithms to apply before matrix construction - `include_constant_impedance_loads::Bool=true`: @@ -320,7 +323,7 @@ for PTDF analysis starting from system data. 1. **Ybus Construction**: Creates system admittance matrix with specified reductions 2. **Incidence Matrix**: Builds bus-branch connectivity matrix A 3. **BA Matrix**: Computes branch susceptance weighted incidence matrix -4. **PTDF Computation**: Calculates power transfer distribution factors using A^T × B^(-1) × A +4. **PTDF Computation**: Calculates power transfer distribution factors by solving ``(A^\\top B A)\\, X = A^\\top B`` 5. **Distributed Slack**: Applies distributed slack correction if specified 6. **Sparsification**: Removes small elements based on tolerance threshold @@ -331,16 +334,19 @@ for PTDF analysis starting from system data. - **Physical Meaning**: Distributed slack better represents generator response to load changes # Linear Solver Options -- **"KLU"**: Sparse LU factorization (default, recommended for most cases) +- **"KLU"**: Sparse LU factorization (default off Apple hardware, recommended for most cases) +- **"AppleAccelerateLU"**: Apple Accelerate sparse LU (default on macOS 15.5+ Apple hardware) - **"Dense"**: Dense matrix operations (faster for small systems, higher memory usage) - **"MKLPardiso"**: Intel MKL Pardiso solver (requires MKL library, best for very large systems) # Mathematical Foundation -The PTDF matrix is computed as: +The PTDF matrix is computed as +```math +\\mathrm{PTDF} = B A (A^\\top B A)^{-1} ``` -PTDF = (A^T × B × A)^(-1) × A^T × B -``` -where A is the incidence matrix and B is the susceptance matrix. +where ``A`` is the incidence matrix (the [`IncidenceMatrix`](@ref)) and ``B`` the branch +susceptance matrix (see [`BA_Matrix`](@ref)). The `data` field holds the transpose +``(A^\\top B A)^{-1} A^\\top B``; use [`get_ptdf_data`](@ref) for the orientation above. # Notes - Results are valid under DC power flow assumptions (linear approximation) @@ -363,7 +369,7 @@ function PTDF(sys::PSY.System; end """ - PTDF(ybus::Ybus; dist_slack::Dict{Int, Float64} = Dict{Int, Float64}(), linear_solver = _default_linear_solver(), tol::Float64 = eps(), network_reductions::Vector{NetworkReduction} = NetworkReduction[], kwargs...) + PTDF(ybus::Ybus; dist_slack::Dict{Int, Float64} = Dict{Int, Float64}(), linear_solver = _default_linear_solver(), tol::Union{Float64, AutoTolerance} = DEFAULT_AUTO_TOLERANCE, network_reductions::Vector{NetworkReduction} = NetworkReduction[], kwargs...) Construct a Power Transfer Distribution Factor (PTDF) matrix from existing Ybus matrix. This constructor is more efficient when the prerequisite matrices are already available and provides @@ -377,9 +383,11 @@ direct control over the underlying matrix computations. Dictionary mapping bus numbers to distributed slack weights for realistic slack modeling. Empty dictionary uses single slack bus (default behavior) - `linear_solver::String = _default_linear_solver()`: - Linear solver algorithm for matrix computations. Options: "KLU", "Dense", "MKLPardiso" -- `tol::Float64 = eps()`: - Sparsification tolerance for dropping small matrix elements to reduce memory usage + Linear solver algorithm for matrix computations. Options: "KLU", "AppleAccelerateLU", "Dense", "MKLPardiso" +- `tol::Union{Float64, AutoTolerance} = DEFAULT_AUTO_TOLERANCE`: + Sparsification tolerance for dropping small matrix elements to reduce memory usage. + A `Float64` applies a fixed absolute cutoff at any size; an [`AutoTolerance`](@ref) + (the default) applies a relative per-row cutoff on large virtual matrices only. # Returns - `PTDF`: The constructed PTDF matrix structure containing: @@ -390,7 +398,7 @@ direct control over the underlying matrix computations. # Construction Process 1. **Incidence Matrix**: Builds bus-branch connectivity matrix A (from Ybus matrix) 2. **BA Matrix**: Computes branch susceptance weighted incidence matrix -3. **PTDF Computation**: Calculates power transfer distribution factors using A^T × B^(-1) × A +3. **PTDF Computation**: Calculates power transfer distribution factors by solving ``(A^\\top B A)\\, X = A^\\top B`` 4. **Distributed Slack**: Applies distributed slack correction if specified 5. **Sparsification**: Removes small elements based on tolerance threshold @@ -401,16 +409,19 @@ direct control over the underlying matrix computations. - **Physical Meaning**: Distributed slack better represents generator response to load changes # Linear Solver Options -- **"KLU"**: Sparse LU factorization (default, recommended for most cases) +- **"KLU"**: Sparse LU factorization (default off Apple hardware, recommended for most cases) +- **"AppleAccelerateLU"**: Apple Accelerate sparse LU (default on macOS 15.5+ Apple hardware) - **"Dense"**: Dense matrix operations (faster for small systems, higher memory usage) - **"MKLPardiso"**: Intel MKL Pardiso solver (requires MKL library, best for very large systems) # Mathematical Foundation -The PTDF matrix is computed as: -``` -PTDF = (A^T × B × A)^(-1) × A^T × B +The PTDF matrix is computed as +```math +\\mathrm{PTDF} = B A (A^\\top B A)^{-1} ``` -where A is the incidence matrix and B is the susceptance matrix. +where ``A`` is the incidence matrix (the [`IncidenceMatrix`](@ref)) and ``B`` the branch +susceptance matrix (see [`BA_Matrix`](@ref)). The `data` field holds the transpose +``(A^\\top B A)^{-1} A^\\top B``; use [`get_ptdf_data`](@ref) for the orientation above. # Notes - Results are valid under DC power flow assumptions (linear approximation) @@ -436,7 +447,7 @@ function PTDF(ybus::Ybus; end """ - PTDF(A::IncidenceMatrix, BA::BA_Matrix; dist_slack::Dict{Int, Float64} = Dict{Int, Float64}(), linear_solver = _default_linear_solver(), tol::Float64 = eps()) + PTDF(A::IncidenceMatrix, BA::BA_Matrix; dist_slack::Dict{Int, Float64} = Dict{Int, Float64}(), linear_solver = _default_linear_solver(), tol::Union{Float64, AutoTolerance} = DEFAULT_AUTO_TOLERANCE) Construct a Power Transfer Distribution Factor (PTDF) matrix from existing incidence and BA matrices. This constructor is more efficient when the prerequisite matrices are already available and provides @@ -451,22 +462,26 @@ direct control over the underlying matrix computations. Dictionary mapping bus numbers to distributed slack participation factors. Empty dictionary uses single slack bus (reference bus from matrices) - `linear_solver::String = _default_linear_solver()`: - Linear solver algorithm for matrix computations. Options: "KLU", "Dense", "MKLPardiso" -- `tol::Float64 = eps()`: - Sparsification tolerance for dropping small matrix elements to reduce memory usage + Linear solver algorithm for matrix computations. Options: "KLU", "AppleAccelerateLU", "Dense", "MKLPardiso" +- `tol::Union{Float64, AutoTolerance} = DEFAULT_AUTO_TOLERANCE`: + Sparsification tolerance for dropping small matrix elements to reduce memory usage. + A `Float64` applies a fixed absolute cutoff at any size; an [`AutoTolerance`](@ref) + (the default) applies a relative per-row cutoff on large virtual matrices only. # Returns - `PTDF`: The constructed PTDF matrix structure with injection-to-flow sensitivity coefficients # Mathematical Computation -The PTDF matrix is computed using the relationship: -``` -PTDF = (A^T × B × A)^(-1) × A^T × B +The PTDF matrix is computed using the relationship +```math +\\mathrm{PTDF} = B A (A^\\top B A)^{-1} ``` where: -- A is the incidence matrix representing bus-branch connectivity -- B is the diagonal susceptance matrix (embedded in BA matrix) +- ``A`` is the incidence matrix representing bus-branch connectivity +- ``B`` is the diagonal susceptance matrix (embedded in BA matrix) - The computation involves solving the ABA linear system for efficiency +- The `data` field holds the transpose ``(A^\\top B A)^{-1} A^\\top B``; use + [`get_ptdf_data`](@ref) for the orientation above # Distributed Slack Handling - **Single Slack**: Uses reference bus identified from input matrices diff --git a/src/radial_reduction.jl b/src/radial_reduction.jl index c50fab131..4b710d0b7 100644 --- a/src/radial_reduction.jl +++ b/src/radial_reduction.jl @@ -56,7 +56,7 @@ with only one connection that do not affect the electrical behavior of the core - `A::SparseArrays.SparseMatrixCSC{Int8, Int}`: The incidence matrix data representing bus-branch connectivity structure - `arc_map::Dict{Tuple{Int, Int}, Int}`: - Dictionary mapping branch endpoint pairs (from_bus, to_bus) to matrix row indices + Dictionary mapping branch endpoint pairs `(from_bus, to_bus)` to matrix row indices - `bus_map::Dict{Int, Int}`: Dictionary mapping bus numbers to matrix column indices - `ref_bus_positions::Set{Int}`: diff --git a/src/row_cache.jl b/src/row_cache.jl index 54dd8e5dc..e75732c0d 100644 --- a/src/row_cache.jl +++ b/src/row_cache.jl @@ -137,8 +137,8 @@ end """ Deletes a row from the stored matrix in cache not belonging to the -persistent_cache_keys set. Uses LRU (Least Recently Used) eviction strategy -based on access_order tracking. +`persistent_cache_keys` set. Uses LRU (Least Recently Used) eviction strategy +based on `access_order` tracking. """ function purge_one!(cache::RowCache) # Use LRU eviction: find oldest non-persistent key @@ -200,7 +200,7 @@ end `cutoff` is the resolved `SparsificationCutoff` stored on the matrix: an `AbsoluteCutoff` drops below a fixed value, a `RelativeCutoff` drops below -`fraction · max|row|` so columns of large cases stay sparse. +``\\mathrm{fraction} \\cdot \\max|\\mathrm{row}|`` so columns of large cases stay sparse. """ function cached_row_lookup( compute_row, diff --git a/src/serialization.jl b/src/serialization.jl index 5871fa544..855aecf6e 100644 --- a/src/serialization.jl +++ b/src/serialization.jl @@ -1,6 +1,13 @@ """ Serialize the PTDF to an HDF5 file. +Only the dense [`PTDF`](@ref) type can be serialized; there is no HDF5 path for +[`LODF`](@ref), [`Ybus`](@ref), the DC susceptance matrices, or any virtual matrix. The file +stores the matrix data, `tol`, axes, lookups, and per-subnetwork axes, but **not** +the [`NetworkReductionData`](@ref): a PTDF built with `network_reductions` loses +that context on a round-trip and is rehydrated with an empty reduction. Keep the +construction code if you need the reduction metadata. + # Arguments - `ptdf::PTDF`: matrix - `filename::AbstractString`: File to create @@ -45,7 +52,12 @@ function to_hdf5( end """ -Deserialize a PTDF from an HDF5 file. +Deserialize a PTDF from an HDF5 file. The convenience constructor `PTDF(filename)` +calls this. + +The returned [`PTDF`](@ref) reproduces the data, axes, lookups, and `tol` of the +serialized matrix, but its [`NetworkReductionData`](@ref) is always empty — the +reduction context is not persisted (see [`to_hdf5`](@ref)). # Arguments - `::Type{PTDF}`: diff --git a/src/virtual_lodf_calculations.jl b/src/virtual_lodf_calculations.jl index 93357fc44..a93dde7d9 100644 --- a/src/virtual_lodf_calculations.jl +++ b/src/virtual_lodf_calculations.jl @@ -11,7 +11,7 @@ The VirtualLODF struct is indexed using branch names. # Thread-safety -Concurrent `getindex` (and `get_partial_lodf_row`) is safe but serialized: +Concurrent `getindex` (and [`get_partial_lodf_row`](@ref)) is safe but serialized: every libklu solve runs under `_LIBKLU_LOCK` (process-wide) and the per-cache `solver_lock`, and the row cache is guarded by `cache_lock`. Multi-threaded callers can issue requests concurrently; the libklu work runs one at a time. @@ -28,7 +28,7 @@ callers can issue requests concurrently; the libklu work runs one at a time. Vector contiaining the element-wise reciprocal of the diagonal elements coming from multuiplying the PTDF matrix with th Incidence matrix - `PTDF_A_diag::Vector{Float64}`: - Raw diagonal elements of the PTDF·A product (H[e,e] values), before + Raw diagonal elements of the ``\\mathrm{PTDF} \\, A`` product (``H[e,e]`` values), before tolerance clamping. Used for partial susceptance change computations. - `arc_susceptances::Vector{Float64}`: Effective susceptance for each arc, extracted from the BA matrix. @@ -119,9 +119,9 @@ end """ _get_PTDF_A_diag(K, BA, A, ref_bus_positions) -> Vector{Float64} -Compute `diag(PTDF · A)`. Each row of `A` has exactly two nonzeros (+1 at the -from-bus, -1 at the to-bus), so the per-arc dot product reduces to two indexed -reads into the solved PTDF row after a one-time transpose of `A`. +Compute ``\\mathrm{diag}(\\mathrm{PTDF} \\, A)``. Each row of ``A`` has exactly two nonzeros +(``+1`` at the from-bus, ``-1`` at the to-bus), so the per-arc dot product reduces to two +indexed reads into the solved PTDF row after a one-time transpose of ``A``. """ function _get_PTDF_A_diag( K, @@ -501,21 +501,23 @@ Compute the partial LODF column for a susceptance change `delta_b` on arc `arc_i Concurrent callers serialize on `vlodf.solver_lock` and `_LIBKLU_LOCK`. Uses the Sherman-Morrison (matrix inversion lemma) formula derived from DC power flow -sensitivity analysis. For a change Δb in the susceptance of arc e, the change in flow -on monitoring arc ℓ per unit pre-change flow on arc e is: - - partial_LODF[ℓ, e] = α · (b_ℓ / b_e) · H[ℓ,e] / (1 - α · H[e,e]) - +sensitivity analysis. For a change ``\\Delta b`` in the susceptance of arc ``e``, the change in flow +on monitoring arc ``\\ell`` per unit pre-change flow on arc ``e`` is +```math +\\mathrm{partialLODF}[\\ell, e] = \\alpha \\, \\frac{b_\\ell}{b_e} \\, \\frac{H[\\ell, e]}{1 - \\alpha \\, H[e, e]} +``` where: -- α = -Δb / b_e (positive for outage/decrease, negative for increase) -- H[ℓ, e] = (A · (ABA)⁻¹ · BA)[ℓ, e] = b_e · C[e, ℓ] (computed via KLU solve) -- b_ℓ = susceptance of monitoring arc ℓ -- H[e,e] = PTDF_A_diag[e] - -When `delta_b = -b_e` (full outage), α = 1 and this reduces to the standard LODF column: - LODF[ℓ, e] = b_ℓ · C[e, ℓ] / (1 - H[e,e]) +- ``\\alpha = -\\Delta b / b_e`` (positive for outage/decrease, negative for increase) +- ``H[\\ell, e] = (A \\, (\\mathrm{ABA})^{-1} \\, \\mathrm{BA})[\\ell, e] = b_e \\, C[e, \\ell]`` (computed via KLU solve) +- ``b_\\ell`` is the susceptance of monitoring arc ``\\ell`` +- ``H[e, e]`` is `PTDF_A_diag[e]` + +When `delta_b = -b_e` (full outage), ``\\alpha = 1`` and this reduces to the standard LODF column +```math +\\mathrm{LODF}[\\ell, e] = \\frac{b_\\ell \\, C[e, \\ell]}{1 - H[e, e]}. +``` When `delta_b = 0`, returns zeros (no change). -The self-element (ℓ = e) is overridden to -1.0 for full outage per standard LODF convention. +The self-element (``\\ell = e``) is overridden to -1.0 for full outage per standard LODF convention. """ function _getindex_partial( vlodf::VirtualLODF, diff --git a/src/virtual_modf_calculations.jl b/src/virtual_modf_calculations.jl index c4c93bfe2..3f9ef881c 100644 --- a/src/virtual_modf_calculations.jl +++ b/src/virtual_modf_calculations.jl @@ -30,7 +30,7 @@ cache and skips the recomputation. - `A::SparseArrays.SparseMatrixCSC{Int8, Int}`: Incidence matrix. - `PTDF_A_diag::Vector{Float64}`: - Diagonal of `PTDF·A` (H[e,e] values). Lazily populated on the first + Diagonal of ``\\mathrm{PTDF} \\, A`` (``H[e,e]`` values). Lazily populated on the first read of `vmodf.PTDF_A_diag`; empty until then. - `arc_susceptances::Vector{Float64}`: Effective susceptance for each arc. @@ -41,7 +41,7 @@ cache and skips the recomputation. - `dist_slack::Vector{Float64}`: Distributed slack bus weights. - `axes::Ax`: - Tuple of (arc_axis, bus_axis). + Tuple of `(arc_axis, bus_axis)`. - `lookup::L`: Tuple of lookup dictionaries for indexing. - `valid_ix::Vector{Int}`: @@ -276,7 +276,7 @@ auto-applied during `Ybus` construction. A `Float64` applies a fixed absolute cutoff; an [`AutoTolerance`](@ref) (the default) applies a relative per-row cutoff so requested columns stay sparse on large systems. -- `max_cache_size::Int`: Max cache size in MiB per contingency (default: MAX_CACHE_SIZE_MiB) +- `max_cache_size::Int`: Max cache size in MiB per contingency (default: `MAX_CACHE_SIZE_MiB`) - `network_reductions::Vector{NetworkReduction}`: Network reductions to apply - `automatically_register_outages::Bool`: Register all system Outage attributes (default: true) """ @@ -487,8 +487,12 @@ end Compute the post-modification PTDF row for a monitored arc under the given modification. Gets or computes Woodbury factors, then applies the Woodbury correction. -For N-1 contingencies, the result satisfies: - post_ptdf[mon, :] = pre_ptdf[mon, :] + LODF[mon, e] * pre_ptdf[e, :] +For N-1 contingencies, the result satisfies +```math +\\mathrm{post}[\\mathrm{mon}, :] = \\mathrm{pre}[\\mathrm{mon}, :] + \\mathrm{LODF}[\\mathrm{mon}, e] \\, \\mathrm{pre}[e, :] +``` +where ``\\mathrm{pre}`` and ``\\mathrm{post}`` are the base and post-contingency PTDF and +``e`` is the outaged arc. """ function _compute_modf_entry( vmodf::VirtualMODF, @@ -624,11 +628,11 @@ end clear_all_caches!(vmodf::VirtualMODF) Clear all caches including contingency registrations. After calling this function, -the `VirtualMODF` object is effectively empty and cannot be queried — it has -no registered contingencies. To restore functionality, a new `VirtualMODF` must +the [`VirtualMODF`](@ref) object is effectively empty and cannot be queried — it has +no registered contingencies. To restore functionality, a new [`VirtualMODF`](@ref) must be constructed from a `PSY.System`. -Use `clear_caches!` instead to preserve contingency registrations while +Use [`clear_caches!`](@ref) instead to preserve contingency registrations while freeing computation cache memory. """ function clear_all_caches!(vmodf::VirtualMODF) diff --git a/src/virtual_ptdf_calculations.jl b/src/virtual_ptdf_calculations.jl index 91e2f1f67..75379f533 100644 --- a/src/virtual_ptdf_calculations.jl +++ b/src/virtual_ptdf_calculations.jl @@ -130,10 +130,12 @@ struct with an empty cache. - `linear_solver::String = _default_linear_solver()`: Linear solver to use for factorization. Options: "KLU", "AppleAccelerateLU". Defaults to "AppleAccelerateLU" on macOS 15.5+ and "KLU" elsewhere. -- `tol::Float64 = eps()`: - Tolerance related to sparsification and values to drop. +- `tol::Union{Float64, AutoTolerance} = DEFAULT_AUTO_TOLERANCE`: + Tolerance related to sparsification and values to drop. A `Float64` applies a + fixed absolute cutoff; an [`AutoTolerance`](@ref) (the default) applies a + relative per-row cutoff so requested rows stay sparse on large systems. - `max_cache_size::Int`: - max cache size in MiB (initialized as MAX_CACHE_SIZE_MiB). + max cache size in MiB (initialized as `MAX_CACHE_SIZE_MiB`). - `persistent_arcs::Vector{Tuple{Int, Int}} = Vector{Tuple{Int, Int}}()`: arcs to be evaluated as soon as the VirtualPTDF is created (initialized as empty vector of tuples). - `network_reduction::NetworkReduction`: @@ -208,10 +210,12 @@ The return is a VirtualPTDF struct with an empty cache. - `linear_solver::String = _default_linear_solver()`: Linear solver to use for factorization. Options: "KLU", "AppleAccelerateLU". Defaults to "AppleAccelerateLU" on macOS 15.5+ and "KLU" elsewhere. -- `tol::Float64 = eps()`: - Tolerance related to sparsification and values to drop. +- `tol::Union{Float64, AutoTolerance} = DEFAULT_AUTO_TOLERANCE`: + Tolerance related to sparsification and values to drop. A `Float64` applies a + fixed absolute cutoff; an [`AutoTolerance`](@ref) (the default) applies a + relative per-row cutoff so requested rows stay sparse on large systems. - `max_cache_size::Int`: - max cache size in MiB (initialized as MAX_CACHE_SIZE_MiB). + max cache size in MiB (initialized as `MAX_CACHE_SIZE_MiB`). - `persistent_arcs::Vector{Tuple{Int, Int}} = Vector{Tuple{Int, Int}}()`: arcs to be evaluated as soon as the VirtualPTDF is created (initialized as empty vector of tuples). """ diff --git a/src/virtual_ptdf_modification.jl b/src/virtual_ptdf_modification.jl index 4009b87fc..5172cabd8 100644 --- a/src/virtual_ptdf_modification.jl +++ b/src/virtual_ptdf_modification.jl @@ -1,7 +1,5 @@ -""" -Public API for computing post-modification PTDF rows from a `VirtualPTDF` -and a `NetworkModification`, using the Woodbury matrix identity. -""" +# Public API for computing post-modification PTDF rows from a VirtualPTDF and a +# NetworkModification, using the Woodbury matrix identity. """ compute_woodbury_factors(vptdf, mod) -> WoodburyFactors @@ -61,10 +59,10 @@ end One-shot convenience function: compute the post-modification PTDF row for a monitored arc under a network modification. Internally calls -`compute_woodbury_factors` then `apply_woodbury_correction`. +[`compute_woodbury_factors`](@ref) then [`apply_woodbury_correction`](@ref). No caching — each call recomputes. Use the two-step API -(`compute_woodbury_factors` + `apply_woodbury_correction`) when querying +([`compute_woodbury_factors`](@ref) + [`apply_woodbury_correction`](@ref)) when querying multiple monitored arcs for the same modification. $(TYPEDSIGNATURES) diff --git a/src/woodbury_kernel.jl b/src/woodbury_kernel.jl index f80431691..c2ccfa958 100644 --- a/src/woodbury_kernel.jl +++ b/src/woodbury_kernel.jl @@ -1,17 +1,15 @@ -""" -Shared Woodbury matrix identity kernel for computing post-modification -network sensitivity factors. Used by both VirtualPTDF and VirtualMODF. - -Implements van Dijk et al. Eq. 29: - B_m⁻¹ = B_r⁻¹ - B_r⁻¹ U (A⁻¹ + U⊤ B_r⁻¹ U)⁻¹ U⊤ B_r⁻¹ -""" +# Shared Woodbury matrix identity kernel for computing post-modification network +# sensitivity factors, used by both VirtualPTDF and VirtualMODF. Implements van Dijk +# et al. Eq. 29: +# +# B_m^-1 = B_r^-1 - B_r^-1 U (A^-1 + U' B_r^-1 U)^-1 U' B_r^-1 """ _invert_woodbury_W(W_mat, ::Val{M}) -> (W_inv::Matrix{Float64}, is_islanding::Bool) -Invert the M×M Woodbury W matrix. Dispatches on `Val{M}` so the compiler -can specialize each case. Analytical formulas for M=1 and M=2 avoid LU -factorization overhead. Falls back to LU for M > 2. +Invert the ``M \\times M`` Woodbury ``W`` matrix. Dispatches on `Val{M}` so the compiler +can specialize each case. Analytical formulas for ``M = 1`` and ``M = 2`` avoid LU +factorization overhead. Falls back to LU for ``M > 2``. """ function _invert_woodbury_W( W_mat::Matrix{Float64},