diff --git a/Project.toml b/Project.toml index a7927567..07152b5e 100644 --- a/Project.toml +++ b/Project.toml @@ -23,11 +23,6 @@ Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" [extensions] PowerFlowsPardisoExt = "Pardiso" -[sources] -# Track PNM main until a release with the KLU/AppleAccelerate/Pardiso solver -# wrappers is registered. CI resolves PNM from here. -PowerNetworkMatrices = {url = "https://github.com/NREL-Sienna/PowerNetworkMatrices.jl", rev = "jd/handle_exceptions_modf"} - [compat] DataFrames = "1" Pardiso = "1" @@ -38,7 +33,7 @@ JSON3 = "1" LineSearches = "7.4.0" LinearAlgebra = "1" Logging = "1" -PowerNetworkMatrices = "^0.21" +PowerNetworkMatrices = "^0.22" PowerSystems = "^5.10" SparseArrays = "1" StaticArrays = "^1.9" diff --git a/docs/make.jl b/docs/make.jl index 34e42dcb..5ce81af2 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -20,6 +20,7 @@ pages = OrderedDict( "Mixed Current-Power Balance Formulation" => "explanation/mixed_cpb_formulation.md", "Levenberg-Marquardt vs Gauss-Seidel" => "explanation/lm_vs_gauss_seidel.md", "LCC Model Implementation" => "explanation/lcc_model.md", + "Discrete Control via λ-Continuation" => "explanation/discrete_control.md", ], "Reference" => Any[ "Public API Reference" => "reference/api/public.md", diff --git a/docs/src/explanation/discrete_control.md b/docs/src/explanation/discrete_control.md new file mode 100644 index 00000000..e722643d --- /dev/null +++ b/docs/src/explanation/discrete_control.md @@ -0,0 +1,233 @@ +# Discrete Control Devices via λ-Continuation + +PowerFlows.jl supports two families of voltage-controlling discrete devices: +**tap-changing transformers** and **switched (stepping) shunts**. This page +explains what those devices are, why an outer-loop continuation strategy is +used instead of embedding the control law in the Jacobian, how the sigmoid +control law and its homotopy ramp work, and what happens when the continuous +solution is snapped to discrete settings. + +## What and why + +Commercial power flow tools (PSS/E, PSLF) simulate discrete control devices as +part of steady-state solution: a tap-changing transformer adjusts its turns +ratio to regulate the voltage at a designated bus; a switched shunt switches +reactive admittance blocks in or out for the same purpose. Without modeling +these devices, the solved bus voltages can differ materially from what a +commercial tool reports, because the network seen by the solver is wrong. + +The natural mathematical formulation (Agarwal, Pandey, Jereminov & Pileggi, +arXiv 1811.02000, §IV.C; Pandey PhD thesis 2018, §5.2.8.1) embeds a +differentiable sigmoid approximation of the discrete control law directly in +the residual and Jacobian, so the Newton iteration solves for device parameters +simultaneously with bus voltages. PowerFlows.jl deliberately does **not** use +that implicit embedding for the initial implementation. Instead it wraps the +existing inner solvers in an **outer-loop λ-continuation**: + +- The inner solver ([`NewtonRaphsonACPowerFlow`](@ref), + [`TrustRegionACPowerFlow`](@ref)) is called without any modification. +- Between outer iterations only `data` is mutated: Y-bus `nzval` entries for + tap devices; the reactive constant-impedance withdrawal matrix for shunt + devices. +- The outer loop is **formulation-agnostic**: it works identically for + [`ACPolarPowerFlow`](@ref), [`ACRectangularPowerFlow`](@ref), and + [`ACMixedPowerFlow`](@ref) because it calls `_solve_with_q_limits!`, the + existing Q-limit loop, as a black box. + +The dispatch point is `_ac_power_flow` in `solve_ac_power_flow.jl`: when +`data.controlled_devices` is non-empty it delegates to +`_control_continuation!`; otherwise the existing path is taken unchanged +(regression invariant). + +A reserved interface method `stamp_control!` exists on every device type but +is never called by the outer loop; it is an error stub that marks where a +future per-formulation implicit embedding would dispatch. Similarly, +`ControlledPhaseShifter` and `ControlledFACTS` are typed seams in the +hierarchy whose interface methods all call an internal `_seam_err` — they +occupy the right dispatch slots but are not implemented. + +## Device abstraction + +Two abstract families sit under `AbstractControlledDevice`: + +- `AbstractBranchControl` — devices that mutate the branch's 2×2 Y-bus block. + The only implemented leaf is `ControlledTap` (voltage-controlling + `TapTransformer`). +- `AbstractShuntControl` — devices that mutate the bus reactive + constant-impedance withdrawal. The only implemented leaf is + `ControlledSwitchedShunt` (voltage-controlling `SwitchedAdmittance`). + +The runtime container is `ControlledDeviceSet`, which holds +`Vector{ControlledTap}` and `Vector{ControlledSwitchedShunt}` as concretely +typed fields. All outer-loop traversal iterates `for d in set.taps` and `for +d in set.shunts`, so dispatch is monomorphic and the hot-path kernels are +allocation-free. + +For `ControlledTap`, `apply_parameter!` rewrites four `nzval` entries of the +sparse Y-bus in place using cached linear offsets (`nz_offsets::NTuple{4,Int}`) +resolved once at device-set construction. The delta update +`nzval[k] += Y_new − Y_old` preserves any parallel-branch contributions already +in the shared slot; `d.current` (not the lossy `nzval`) is the authoritative +source for the old parameter. For `ControlledSwitchedShunt`, +`apply_parameter!` is a scalar write to +`data.bus_reactive_power_constant_impedance_withdrawals[bus_ix, ts]`. + +## The sigmoid control law + +The continuous target for each device is a sigmoid function of the controlled +bus voltage magnitude $|V|$: + +$$\sigma(\ell, h, S, x, x_\text{set}) = \frac{h - \ell}{1 + e^{S(x - x_\text{set})}} + \ell$$ + +For a **tap transformer** controlled on its primary side (equation 46 of the +thesis), the sigmoid maps a low voltage to a high tap ratio and a high voltage +to a low tap ratio: $\ell = p_\min$, $h = p_\max$. When the controlled bus is +on the secondary side (equation 47), the limits are swapped +($\ell = p_\max$, $h = p_\min$), making the sigmoid increasing rather than +decreasing. The result is clamped to $[p_\min, p_\max]$. + +For a **switched shunt** (equation 9), low voltage calls for more susceptance +and high voltage for less: $\ell = b_\min$, $h = b_\max$; the sigmoid is +decreasing in $|V|$. + +The steepness $S$ controls how closely the smooth sigmoid approximates the +step function. It starts at `INITIAL_CONTROL_STEEPNESS = 100` and is ramped +by `CONTROL_STEEPNESS_GROWTH = 2.0` after each settling phase, up to +`MAX_CONTROL_STEEPNESS = 5000` (values from paper equation 10 / thesis). The +ramp happens only after the devices have settled at the current $S$, so the +solver is never asked to handle a stiff sigmoid before the network state is +compatible with it. + +### Plant-sign correction + +The sigmoid formulas above assume a particular orientation between the +controlled device and the controlled bus. Real networks can wire either end of +a transformer to either bus, so the closed-loop gain may be positive or negative +with the nominal sigmoid orientation. `_plant_sign` measures the sensitivity +$\partial|V|/\partial p$ by a small perturbation at the current parameter +value. If the resulting plant gain is positive, it means the nominal sigmoid +orientation (decreasing for primary-controlled tap) would produce positive +feedback; the `flip` flag is set and the sigmoid limits are swapped +(`_control_target` passes `(hi, lo)` instead of `(lo, hi)`), restoring +negative feedback regardless of device wiring. + +## The continuation engine + +The outer loop `_control_continuation!` runs up to +`MAX_CONTROL_OUTER_ITERATIONS = 20` passes. Each pass: + +1. Computes the sigmoid target $p^*$ for every device given the current + controlled-bus voltage and steepness $S$. +2. Applies an **adaptive under-relaxation** step + $p \leftarrow p + \omega(p^* - p)$. +3. Applies each update via `_continuation_to!`, the incremental robust + applicator. + +**Under-relaxation.** The factor $\omega$ is chosen so the local iteration-map +slope magnitude is at most `CONTROL_CONTRACTION = 0.7`. The closed-loop slope +is bounded by $|h - \ell| \cdot S/4 \cdot |\partial|V|/\partial p|$ (the +maximum sigmoid derivative times the plant gain), giving + +$$\omega \leq \frac{1 + \theta}{1 + |g'|}$$ + +where $\theta = 0.7$ is the contraction target and $g'$ is the closed-loop +gain bound. This guarantees monotone convergence at every steepness without +re-measuring the plant gain each iteration. The cap `CONTROL_RELAXATION_MAX = +0.8` prevents $\omega$ from being set too high when $g'$ is near zero. + +**Incremental applicator.** `_continuation_to!` walks the parameter from its +current value to the relaxed target in sub-steps. The first sub-step is +`MIN_LAMBDA_STEP = 1e-3` of the total interval; successful NR calls allow the +step to grow (factor `CONTROL_STEP_GROWTH = 1.5`, capped at +`MAX_LAMBDA_STEP = 1.0`); a failed NR call halves the step. If the step falls +below `MIN_LAMBDA_STEP` the parameter stays at the last converged point and +the applicator returns the reached value. + +**Settling and steepness ramp.** A pass is considered settled when every device +moves by less than `CONTROL_PARAM_TOL = 1e-5`. Steepness $S$ is advanced only +after settling, not after every iteration. If `MAX_CONTROL_OUTER_ITERATIONS` +is reached before full-steepness convergence, an honest `@warn` is emitted +(regulation may be loose); the solver does not claim failure on account of +steepness alone. + +**Oscillation safeguard.** If the direction of a device's parameter update +reverses sign on more than `CONTROL_OSCILLATION_LIMIT = 3` consecutive outer +passes, the device is frozen at its current parameter and a `@warn` is emitted. +The frozen device returns `Inf` as its gap, so the outer loop does not count +it as settled while other devices continue to converge. This mirrors the +PSS/E behavior of locking oscillating adjustments after a fixed iteration +count. + +## Snap and feasibility restoration + +After the continuous outer loop converges (or reaches its iteration limit), +`snap_and_restore!` discretizes every device: + +- **Tap.** The continuous parameter is snapped to the nearest value in the + pre-computed `levels` vector (a uniform grid of `NTP` tap positions from + `p_min` to `p_max`). +- **Shunt.** The target susceptance is placed using a block-greedy algorithm: + blocks are processed largest-first; each block's step count is floored to + avoid overshooting; a single ±1 bounded refinement pass corrects + under-committed blocks. Both `block_order` and `block_n` are pre-allocated + fields, so `snap_to_discrete` makes no heap allocations. + +After snapping all devices, the inner solver is called on the discretized +network. If it converges, the procedure is complete. If not, each device is +individually restored toward its pre-snap continuous value using the same +bisection sub-stepping as `_continuation_to!` (`_restore_one!`). If any +device cannot be restored to a converged state, `data.converged[ts] = false` +is set and an `@error` is emitted with the device names; no non-physical +solution is silently returned. + +## Metadata sourcing + +Device parameters are sourced from the PSS/E parser `ext` dictionary on each +`TapTransformer` or `SwitchedAdmittance`, with documented defaults when keys +are absent. + +### `TapTransformer` → `ControlledTap` + +Only transformers with `get_control_objective(tx) == VOLTAGE` are included. + +| Parameter | `ext` key | Default | +|---|---|---| +| Controlled bus | `NREG` or `RMIDNT` (bus number) | transformer to-bus (secondary) | +| Voltage setpoint | `VSET` | controlled bus `magnitude`, else `1.0` | +| Tap ratio min | `RMI` | `DEFAULT_TAP_RATIO_MIN = 0.9` | +| Tap ratio max | `RMA` | `DEFAULT_TAP_RATIO_MAX = 1.1` | +| Tap positions | `NTP` | `DEFAULT_TAP_POSITIONS = 33` | +| Controlled side | derived: controlled bus == from-bus → primary (eq. 46), else secondary (eq. 47) | — | + +The discrete tap levels are `range(p_min, p_max; length=NTP)`, collected once +at construction. + +### `SwitchedAdmittance` → `ControlledSwitchedShunt` + +| Parameter | Source | +|---|---| +| Controlled bus | `ext["NREG"]` or `ext["RMIDNT"]`, else own bus | +| Voltage setpoint | midpoint of `get_admittance_limits` (VSWLO/VSWHI band) | +| Susceptance limits | `get_Y` base + Σ over blocks of `number_of_steps .* imag.(Y_increase)` | +| Block structure | `get_number_of_steps`, `get_Y_increase`, `get_initial_status` | + +If the controlled bus number cannot be resolved to a network index, +`build_controlled_device_set` raises an error immediately; there is no silent +skip. + +## Reserved seams for future work + +`stamp_control!(d::AbstractControlledDevice, args...)` is defined on the +abstract base and calls `error("implicit embedding not implemented for ...")`. +It is the attachment point for a future per-formulation implicit embedding +(sigmoid term added to the residual and corresponding Jacobian row), which +would allow the inner Newton iteration to solve for device parameters +simultaneously with bus voltages. Such an embedding would be added per +formulation type by overloading `stamp_control!`; the outer-loop code and the +device hierarchy would require no changes. + +`ControlledPhaseShifter <: AbstractBranchControl` and +`ControlledFACTS <: AbstractShuntControl` are stub types whose interface +methods all call `_seam_err`. They hold the correct position in the dispatch +hierarchy for active-power phase-shifting and FACTS device control, which are +not implemented in this scope. diff --git a/docs/src/reference/api/internal_solvers.md b/docs/src/reference/api/internal_solvers.md index 9258a3d4..87776fc1 100644 --- a/docs/src/reference/api/internal_solvers.md +++ b/docs/src/reference/api/internal_solvers.md @@ -42,6 +42,19 @@ Pages = [ ] ``` +## Discrete Control via λ-Continuation + +```@autodocs +Modules = [PowerFlows] +Public = true +Private = true +Pages = [ + "discrete_control/controlled_devices.jl", + "discrete_control/control_metadata.jl", + "discrete_control/control_continuation.jl", +] +``` + # Linear Algebra Backends ## Robust Homotopy diff --git a/src/PowerFlowData.jl b/src/PowerFlowData.jl index 48e778ad..f9e5afa6 100644 --- a/src/PowerFlowData.jl +++ b/src/PowerFlowData.jl @@ -137,6 +137,7 @@ struct PowerFlowData{ lcc::LCCParameters arc_lossy_admittance_from_to::Union{SparseMatrixCSC{YBUS_ELTYPE, Int}, Nothing} arc_lossy_admittance_to_from::Union{SparseMatrixCSC{YBUS_ELTYPE, Int}, Nothing} + controlled_devices::Union{Nothing, ControlledDeviceSet} # Persistent linear-solver cache for the DC solves. Reused across repeated solves on # the same data (e.g. a PCM loop: fixed network, changing injections) so the network # matrix is factored once and the solve buffer is not reallocated. Despite the name it @@ -229,6 +230,7 @@ get_converged(pfd::PowerFlowData) = pfd.converged get_loss_factors(pfd::PowerFlowData) = pfd.loss_factors get_voltage_stability_factors(pfd::PowerFlowData) = pfd.voltage_stability_factors get_arc_active_power_losses(pfd::PowerFlowData) = pfd.arc_active_power_losses +get_controlled_devices(pfd::PowerFlowData) = pfd.controlled_devices # Field getter for expanded slack participation factors (one dict per time step) # Named "computed" to distinguish from the user-supplied pf.generator_slack_participation_factors @@ -326,6 +328,7 @@ function PowerFlowData( neighbors = Vector{Set{Int}}(), arc_lossy_admittance_from_to::Union{SparseMatrixCSC{YBUS_ELTYPE, Int}, Nothing} = nothing, arc_lossy_admittance_to_from::Union{SparseMatrixCSC{YBUS_ELTYPE, Int}, Nothing} = nothing, + controlled_devices::Union{Nothing, ControlledDeviceSet} = nothing, arc_bus_incidence::Union{SparseMatrixCSC{Int8, Int}, Nothing} = nothing, ) where { T <: PowerFlowEvaluationModel, @@ -385,6 +388,7 @@ function PowerFlowData( lcc_parameters, arc_lossy_admittance_from_to, arc_lossy_admittance_to_from, + controlled_devices, Base.RefValue{Any}(nothing), # solver_cache (lazily populated) ) end @@ -472,6 +476,7 @@ function make_and_initialize_power_flow_data( neighbors = Vector{Set{Int}}(), arc_lossy_admittance_from_to::Union{SparseMatrixCSC{YBUS_ELTYPE, Int}, Nothing} = nothing, arc_lossy_admittance_to_from::Union{SparseMatrixCSC{YBUS_ELTYPE, Int}, Nothing} = nothing, + controlled_devices::Union{Nothing, ControlledDeviceSet} = nothing, arc_bus_incidence::Union{SparseMatrixCSC{Int8, Int}, Nothing} = nothing, ) where {M <: PNM.PowerNetworkMatrix, N <: Union{PNM.PowerNetworkMatrix, Nothing}} check_unit_setting(sys) @@ -490,6 +495,7 @@ function make_and_initialize_power_flow_data( neighbors = neighbors, arc_lossy_admittance_from_to = arc_lossy_admittance_from_to, arc_lossy_admittance_to_from = arc_lossy_admittance_to_from, + controlled_devices = controlled_devices, arc_bus_incidence = arc_bus_incidence, ) @assert length(data.lcc.setpoint_at_rectifier) == n_lccs @@ -555,12 +561,21 @@ function PowerFlowData( aux_network_matrix = nothing end + controlled_devices = if get_control_discrete_devices(pf) + bus_lookup = PNM.get_bus_lookup(power_network_matrix) + set = build_controlled_device_set(sys, bus_lookup, power_network_matrix) + isempty(set) ? nothing : set + else + nothing + end + return make_and_initialize_power_flow_data( pf, sys, power_network_matrix, aux_network_matrix; neighbors = neighbors, + controlled_devices = controlled_devices, ) end diff --git a/src/PowerFlows.jl b/src/PowerFlows.jl index 95b146f6..7eb98f60 100644 --- a/src/PowerFlows.jl +++ b/src/PowerFlows.jl @@ -23,6 +23,7 @@ export update_exporter! export write_export export get_psse_export_paths export FlowReporting +export ControlledDeviceSet # "protected" (semi-stable because used in PSI) but not exported: # PowerFlowData and related type aliases, solve_power_flow!, write_results @@ -53,6 +54,9 @@ include("psi_utils.jl") include("powersystems_utils.jl") include("power_flow_types.jl") include("lcc_parameters.jl") +include("discrete_control/controlled_devices.jl") +include("discrete_control/control_metadata.jl") +include("discrete_control/control_continuation.jl") include("PowerFlowData.jl") include("lcc_utils.jl") include("common.jl") diff --git a/src/RobustHomotopy/HessianSolver/cholesky_solver.jl b/src/RobustHomotopy/HessianSolver/cholesky_solver.jl index 79415572..feda37b7 100644 --- a/src/RobustHomotopy/HessianSolver/cholesky_solver.jl +++ b/src/RobustHomotopy/HessianSolver/cholesky_solver.jl @@ -1,6 +1,6 @@ struct CholeskyHessianSolver <: HessianSolver F::SparseArrays.CHOLMOD.Factor{Float64, J_INDEX_TYPE} - mat::FixedStructureCHOLMOD + mat::FixedStructureCHOLMOD{Float64, J_INDEX_TYPE} buff::Vector{Float64} # buffer for solving end @@ -22,33 +22,29 @@ function modify_and_numeric_factor!( ) minDiagElem = minimum(H[i, i] for i in axes(H, 1)) τ_old = 0.0 - if minDiagElem > 0.0 - τ = 0.0 - else - τ = -minDiagElem + β - end + τ = minDiagElem > 0.0 ? 0.0 : -minDiagElem + β @debug "initial τ = $τ" nonsingular = false - fill!(hSolver.buff, 1.0) while !nonsingular for i in axes(H, 1) H[i, i] += τ - τ_old # now try H + τ*I end - try - fill!(hSolver.buff, 1.0) # reset b to a vector of ones. - set_values!(hSolver.mat, SparseArrays.nonzeros(H)) + set_values!(hSolver.mat, SparseArrays.nonzeros(H)) + # `issuccess` checks PD directly, replacing the throwaway `F \ ones` solve. + # `numeric_factor!` can still throw PosDefException, so keep the catch. + ok = try numeric_factor!(hSolver.F, hSolver.mat) - hSolver.F \ hSolver.buff # sometimes the error isn't thrown until we use the factorization. + LinearAlgebra.issuccess(hSolver.F) + catch e + e isa SparseArrays.CHOLMOD.PosDefException ? false : rethrow(e) + end + if ok nonsingular = true @debug "nonsingular with τ = $τ" - catch e - if e isa SparseArrays.CHOLMOD.PosDefException - τ_old = τ - τ *= 2.0 - τ = max(τ, β) # ensure τ is at least β - else - rethrow(e) - end + else + τ_old = τ + τ *= 2.0 + τ = max(τ, β) # ensure τ is at least β end end return diff --git a/src/RobustHomotopy/homotopy_hessian.jl b/src/RobustHomotopy/homotopy_hessian.jl index 113ee443..e12dfe81 100644 --- a/src/RobustHomotopy/homotopy_hessian.jl +++ b/src/RobustHomotopy/homotopy_hessian.jl @@ -6,6 +6,11 @@ struct HomotopyHessian PQ_V_mags::BitVector # true iff that coordinate in the state vector is V_mag at a PQ bus grad::Vector{Float64} Hv::SparseMatrixCSC{Float64, J_INDEX_TYPE} + # Scratch/precompute for an allocation-free hot path: `Jt_R` holds Jᵀ·Rv; `pq_diag_nz` + # are the Hv.nzval indices of the PQ |V| diagonal entries that take the (1−t) term + # (avoids a per-call sparse `setindex!` on those diagonals). + Jt_R::Vector{Float64} + pq_diag_nz::Vector{Int} end """Does `A += B' * B`, in a way that preserves the sparse structure of `A`, if possible. @@ -27,23 +32,47 @@ function (hess::HomotopyHessian)(x::Vector{Float64}, t_k::Float64, time_step::In Jv = hess.J.Jv _update_hessian_matrix_values!(hess.Hv, Rv, hess.data, time_step) A_plus_eq_BT_B!(hess.Hv, Jv) - SparseArrays.nonzeros(hess.Hv) .*= t_k - for (bus_ix, bt) in enumerate(get_bus_type(hess.data)[:, time_step]) # PERF: allocating - if bt == PSY.ACBusTypes.PQ - hess.Hv[2 * bus_ix - 1, 2 * bus_ix - 1] += (1 - t_k) - end + Hvnz = SparseArrays.nonzeros(hess.Hv) + Hvnz .*= t_k + # (1−t) homotopy term on the PQ |V| diagonal. + @inbounds for k in hess.pq_diag_nz + Hvnz[k] += (1 - t_k) end - # PERF: allocating - hess.grad .= (1 - t_k) * hess.PQ_V_mags .* (x - ones(size(x, 1))) + t_k * Jv' * Rv + _homotopy_gradient!(hess.grad, hess, t_k, x, Jv, Rv) return end +# grad = (1−t)·mask·(x−1) + t·Jᵀ·Rv, in place. +function _homotopy_gradient!( + grad::Vector{Float64}, + hess::HomotopyHessian, + t_k::Float64, + x::Vector{Float64}, + Jv::SparseMatrixCSC{Float64, J_INDEX_TYPE}, + Rv::Vector{Float64}, +) + LinearAlgebra.mul!(hess.Jt_R, Jv', Rv) + mask = hess.PQ_V_mags + @inbounds for i in eachindex(grad) + grad[i] = t_k * hess.Jt_R[i] + + (mask[i] ? (1 - t_k) * (x[i] - 1.0) : 0.0) + end + return grad +end + function F_value(hess::HomotopyHessian, t_k::Float64, x::Vector{Float64}, time_step::Int) hess.pfResidual(x, time_step) Rv = hess.pfResidual.Rv - φ_vector = x[hess.PQ_V_mags] .- 1.0 # PERF: allocating - F_value = (1 - t_k) * 0.5 * dot(φ_vector, φ_vector) + t_k * 0.5 * dot(Rv, Rv) - return F_value + # Σ (x−1)² over PQ |V| coordinates. + φ_sq = 0.0 + mask = hess.PQ_V_mags + @inbounds for i in eachindex(x) + if mask[i] + d = x[i] - 1.0 + φ_sq += d * d + end + end + return (1 - t_k) * 0.5 * φ_sq + t_k * 0.5 * dot(Rv, Rv) end # slightly confusing that I have the field grad, and the argument grad. @@ -57,16 +86,13 @@ function gradient_value!(grad::Vector{Float64}, hess.J(time_step) # PERF bottleneck. Look into a different line search strategy? # or otherwise reduce the number of gradient computations? # for a 10k bus system, computing J takes over 10x longer than computing F. - Jv = hess.J.Jv - mask = hess.PQ_V_mags - # PERF: allocating - grad .= (1 - t_k) * (mask .* (x - ones(size(x, 1)))) + t_k * Jv' * hess.pfResidual.Rv + _homotopy_gradient!(grad, hess, t_k, x, hess.J.Jv, hess.pfResidual.Rv) return grad end function homotopy_x0(data::ACPowerFlowData, time_step::Int) x = calculate_x0(data, time_step) - for (bus_ix, bt) in enumerate(get_bus_type(data)[:, time_step]) # PERF: allocating + for (bus_ix, bt) in enumerate(view(get_bus_type(data), :, time_step)) if bt == PSY.ACBusTypes.PQ x[2 * bus_ix - 1] = 1.0 end @@ -102,9 +128,18 @@ function HomotopyHessian(data::ACPowerFlowData, time_step::Int) SparseArrays.nonzeros(Hv) .= 0.0 copyto!(SparseArrays.nonzeros(J.Jv), original_J_nzval) nbuses = size(get_bus_type(data), 1) - PQ_mask = get_bus_type(data)[:, time_step] .== (PSY.ACBusTypes.PQ,) + bus_types = view(get_bus_type(data), :, time_step) + PQ_mask = bus_types .== (PSY.ACBusTypes.PQ,) PQ_V_mags = collect(Iterators.flatten(zip(PQ_mask, falses(nbuses)))) - return HomotopyHessian(data, pfResidual, J, PQ_V_mags, zeros(2 * nbuses), Hv) + # Precompute the Hv.nzval indices of the PQ |V| diagonal entries (always + # structurally present in the JᵀJ pattern, since that column of J is nonzero). + pq_diag_nz = [ + _nz_index(Hv, 2 * b - 1, 2 * b - 1) + for b in 1:nbuses if bus_types[b] == PSY.ACBusTypes.PQ + ] + return HomotopyHessian( + data, pfResidual, J, PQ_V_mags, zeros(2 * nbuses), Hv, + zeros(2 * nbuses), pq_diag_nz) end """ diff --git a/src/ac_power_flow_jacobian.jl b/src/ac_power_flow_jacobian.jl index 544d7d61..fe63685e 100644 --- a/src/ac_power_flow_jacobian.jl +++ b/src/ac_power_flow_jacobian.jl @@ -8,16 +8,14 @@ and can be called as a function at the same time. Calling the instance as a func # Fields - `data::ACPowerFlowData`: The grid model data used for power flow calculations. -- `Jf!::Function`: A function that calculates the Jacobian matrix inplace. -- `Jv::SparseArrays.SparseMatrixCSC{Float64, $J_INDEX_TYPE}`: The Jacobian matrix, which is updated by the function `Jf!`. +- `Jv::SparseArrays.SparseMatrixCSC{Float64, $J_INDEX_TYPE}`: The Jacobian matrix, which is updated by `_update_jacobian_matrix_values!`. - `diag_elements::MVector{4, Float64}`: Temporary storage for diagonal elements during Jacobian update. - `bus_slack_participation_factors::SparseVector{Float64, Int}`: Normalized per-bus slack participation factors for the current time step (from the `ACPowerFlowResidual`). Used for the distributed slack Jacobian entries. - `subnetworks::Dict{Int64, Vector{Int64}}`: Subnetwork mapping from REF bus to bus list (from the `ACPowerFlowResidual`). Used for the distributed slack Jacobian entries. """ struct ACPowerFlowJacobian data::ACPowerFlowData - Jf!::Function # This is the function that calculates the Jacobian matrix and updates Jv inplace - Jv::SparseArrays.SparseMatrixCSC{Float64, J_INDEX_TYPE} # This is the Jacobian matrix, that is updated by the function Jf + Jv::SparseArrays.SparseMatrixCSC{Float64, J_INDEX_TYPE} # This is the Jacobian matrix, updated in place by `_update_jacobian_matrix_values!` diag_elements::MVector{4, Float64} # Temporary storage for diagonal elements during Jacobian update bus_slack_participation_factors::SparseVector{Float64, Int} subnetworks::Dict{Int64, Vector{Int64}} @@ -30,7 +28,7 @@ end """ (J::ACPowerFlowJacobian)(time_step::Int64) -Update the Jacobian matrix `Jv` using the function `Jf!` and the provided data and time step. +Update the Jacobian matrix `Jv` using `_update_jacobian_matrix_values!` and the provided data and time step. Defining this method allows an instance of `ACPowerFlowJacobian` to be called as a function, following the functor pattern. @@ -45,7 +43,7 @@ J(time_step) # Updates the Jacobian matrix Jv ``` """ function (J::ACPowerFlowJacobian)(time_step::Int64) - J.Jf!(J.Jv, J.data, time_step, J.diag_elements, + _update_jacobian_matrix_values!(J.Jv, J.data, time_step, J.diag_elements, J.bus_slack_participation_factors, J.subnetworks, J.bus_active_constant_I, J.bus_reactive_constant_I, J.bus_active_constant_Z, J.bus_reactive_constant_Z) @@ -57,7 +55,7 @@ end Use the `ACPowerFlowJacobian` to update the provided Jacobian matrix `J` inplace. -Update the internally stored Jacobian matrix `Jv` using the function `Jf!` and the provided data and time step, and write the updated Jacobian values to `J`. +Update the internally stored Jacobian matrix `Jv` using `_update_jacobian_matrix_values!` and the provided data and time step, and write the updated Jacobian values to `J`. This method allows an instance of ACPowerFlowJacobian to be called as a function, following the functor pattern. @@ -77,7 +75,7 @@ function (J::ACPowerFlowJacobian)( Jv::SparseArrays.SparseMatrixCSC{Float64, J_INDEX_TYPE}, time_step::Int64, ) - J.Jf!(J.Jv, J.data, time_step, J.diag_elements, + _update_jacobian_matrix_values!(J.Jv, J.data, time_step, J.diag_elements, J.bus_slack_participation_factors, J.subnetworks, J.bus_active_constant_I, J.bus_reactive_constant_I, J.bus_active_constant_Z, J.bus_reactive_constant_Z) @@ -107,11 +105,41 @@ J(time_step) # Updates the Jacobian matrix stored internally in J. J.Jv # Access the Jacobian matrix stored internally in J. ``` """ +# Memoize the expensive Jacobian sparse-structure build (~3.2 MB on 2000 buses) in the +# persistent `data.solver_cache` Ref so it is built once and reused across the Q-limit +# inner loop and repeated PCM solves. The structure depends only on network topology, REF +# layout, and slack-participation pattern — none change on the PV→PQ flips that drive the +# Q-limit loop (colptr/rowval verified byte-identical across a flip). Cache key is the +# network-matrix identity + slack nonzero pattern, so a distributed-slack participant drop +# (which changes the pattern) correctly rebuilds. A full `copy` is returned because some +# nzval entries (e.g. LCC angle-constraint diagonals = 1.0) are set at structure creation +# and never rewritten by `J(time_step)`, so they must be carried over, not zeroed. +function _get_or_build_jacobian_structure( + data::ACPowerFlowData, + slack_factors::SparseVector{Float64, Int}, + subnetworks::Dict{Int64, Vector{Int64}}, + time_step::Int64, +) + nzind = SparseArrays.nonzeroinds(slack_factors) + entry = data.solver_cache[] + if entry isa Tuple{Symbol, Any, Vector{Int}, SparseMatrixCSC{Float64, J_INDEX_TYPE}} && + entry[1] === :ac_jac_structure && + entry[2] === data.power_network_matrix && + entry[3] == nzind + return copy(entry[4]) + end + Jv0 = _create_jacobian_matrix_structure(data, slack_factors, subnetworks, time_step) + # Cache a pristine copy; `Jv0` is about to be mutated by the Newton loop. + data.solver_cache[] = + (:ac_jac_structure, data.power_network_matrix, copy(nzind), copy(Jv0)) + return Jv0 +end + function ACPowerFlowJacobian( residual::ACPowerFlowResidual, time_step::Int64, ) - Jv0 = _create_jacobian_matrix_structure( + Jv0 = _get_or_build_jacobian_structure( residual.data, residual.bus_slack_participation_factors, residual.subnetworks, @@ -119,7 +147,6 @@ function ACPowerFlowJacobian( ) return ACPowerFlowJacobian( residual.data, - _update_jacobian_matrix_values!, Jv0, MVector{4, Float64}(undef), residual.bus_slack_participation_factors, diff --git a/src/ac_power_flow_residual.jl b/src/ac_power_flow_residual.jl index 93d2965d..48371f45 100644 --- a/src/ac_power_flow_residual.jl +++ b/src/ac_power_flow_residual.jl @@ -5,7 +5,6 @@ A struct to keep track of the residuals in the Newton-Raphson AC power flow calc # Fields - `data::ACPowerFlowData`: The grid model data. -- `Rf!::Function`: A function that updates the residuals based on the latest values stored in the grid at the given iteration. - `Rv::Vector{Float64}`: A vector of the values of the residuals. - `P_net::Vector{Float64}`: A vector of net active power injections. - `Q_net::Vector{Float64}`: A vector of net reactive power injections. @@ -17,7 +16,6 @@ A struct to keep track of the residuals in the Newton-Raphson AC power flow calc """ struct ACPowerFlowResidual data::ACPowerFlowData - Rf!::Function Rv::Vector{Float64} P_net::Vector{Float64} Q_net::Vector{Float64} @@ -85,7 +83,6 @@ function ACPowerFlowResidual(data::ACPowerFlowData, time_step::Int64) return ACPowerFlowResidual( data, - _update_residual_values!, Vector{Float64}(undef, 2 * n_buses + 4 * n_lccs), P_net, Q_net, @@ -121,7 +118,7 @@ function (Residual::ACPowerFlowResidual)( x::Vector{Float64}, time_step::Int64, ) - Residual.Rf!( + _update_residual_values!( Residual.Rv, x, Residual.P_net, @@ -155,7 +152,7 @@ Calling the `ACPowerFlowResidual` will also update the values of P, Q, V, Θ in - `time_step::Int64`: The current time step. """ function (Residual::ACPowerFlowResidual)(x::Vector{Float64}, time_step::Int64) - Residual.Rf!( + _update_residual_values!( Residual.Rv, x, Residual.P_net, diff --git a/src/definitions.jl b/src/definitions.jl index 2b783533..2d4e66bf 100644 --- a/src/definitions.jl +++ b/src/definitions.jl @@ -2,6 +2,23 @@ const MAX_INIT_RESIDUAL = 10.0 const BOUNDS_TOLERANCE = 1e-6 const INFINITE_BOUND = 1e6 # used as default when a branch has rating 0.0, as implied by the PSSE Manual const MAX_REACTIVE_POWER_ITERATIONS = 10 + +# Discrete control device λ-continuation outer loop. +const MAX_CONTROL_OUTER_ITERATIONS = 100 +const CONTROL_PARAM_TOL = 1e-5 +const INITIAL_LAMBDA_STEP = 1.0 +const MIN_LAMBDA_STEP = 1e-3 +const MAX_LAMBDA_STEP = 1.0 +const CONTROL_STEP_GROWTH = 1.5 +const INITIAL_CONTROL_STEEPNESS = 1.0e2 # paper eq.10: start (S−λ_S)≈100 +const MAX_CONTROL_STEEPNESS = 5.0e3 # paper: full S≈5000 +const CONTROL_STEEPNESS_GROWTH = 2.0 +const CONTROL_OSCILLATION_LIMIT = 3 +const DEFAULT_TAP_RATIO_MIN = 0.9 +const DEFAULT_TAP_RATIO_MAX = 1.1 +const DEFAULT_TAP_POSITIONS = 33 +const DEFAULT_TAP_VSET = 1.0 +const DEFAULT_SHUNT_MODSW = 1 # discrete voltage control when ext lacks "MODSW" const DEFAULT_MAX_REDISTRIBUTION_ITERATIONS = 10 const LARGE_RESIDUAL = 10 # threshold for "bad initial guess": default # norm(residual, 1)/length(residual) > 10. diff --git a/src/discrete_control/control_continuation.jl b/src/discrete_control/control_continuation.jl new file mode 100644 index 00000000..addfa4d6 --- /dev/null +++ b/src/discrete_control/control_continuation.jl @@ -0,0 +1,276 @@ +# Upper cap on the under-relaxation factor ω; `_relaxation` adapts ω down from here +# to keep the damped iteration a contraction at the current steepness. +const CONTROL_RELAXATION_MAX = 0.8 +# Target slope of the local iteration map near the fixed point (0lo. +@inline function _control_target(d, vmag::Float64, S::Float64, dVdp::Float64) + lo, hi = parameter_limits(d) + return if dVdp > 0.0 + clamp(_sigmoid(lo, hi, S, vmag, voltage_setpoint(d)), lo, hi) + else + clamp(_sigmoid(hi, lo, S, vmag, voltage_setpoint(d)), lo, hi) + end +end + +# Measure dV/dp at the controlled bus by a small parameter perturbation; restore +# afterward. `reliable=false` (a probe solve failed) ⇒ orientation unknown, so the +# caller freezes the device rather than stepping it with an unknown sign. +function _plant_sign(d, data, ts::Int, pf; kwargs...)::Tuple{Float64, Bool} + p0 = current_parameter(d) + cix = controlled_bus_ix(d) + V0 = _read_vmag(data, cix, ts) + lo, hi = parameter_limits(d) + δ = 1e-3 * (hi - lo) + δ = δ > 0.0 ? δ : 1e-6 + p1 = clamp(p0 + δ, lo, hi) + p1 == p0 && (p1 = clamp(p0 - δ, lo, hi)) + h = p1 - p0 + apply_parameter!(d, data, p1, ts) + ok1 = _solve_with_q_limits!(pf, data, ts; kwargs...) + V1 = _read_vmag(data, cix, ts) + apply_parameter!(d, data, p0, ts) + ok2 = _solve_with_q_limits!(pf, data, ts; kwargs...) + reliable = ok1 && ok2 && h != 0.0 + dVdp = reliable ? (V1 - V0) / h : 0.0 + return dVdp, reliable +end + +# Incremental robust applicator: walk the parameter from `start = d.current` +# toward `target` so NR stays converged. The FIRST probe is a SMALL step +# (fraction `MIN_LAMBDA_STEP` of the interval), growing on NR success and +# halving on failure (bisection backtracking). Returns the parameter actually +# reached and leaves the solver converged there. +function _continuation_to!(d, data, ts::Int, target::Float64, pf; kwargs...) + start = current_parameter(d) + abs(target - start) < CONTROL_PARAM_TOL && return start + done = 0.0 # fraction of [start,target] applied so far + step = MIN_LAMBDA_STEP + reached = start + while done < 1.0 + trial = min(1.0, done + step) + p = start + trial * (target - start) + apply_parameter!(d, data, p, ts) + if _solve_with_q_limits!(pf, data, ts; kwargs...) + done = trial + reached = p + step = min(step * CONTROL_STEP_GROWTH, MAX_LAMBDA_STEP) + else + apply_parameter!(d, data, reached, ts) + step /= 2.0 + if step < MIN_LAMBDA_STEP + _solve_with_q_limits!(pf, data, ts; kwargs...) + return reached + end + end + end + return reached +end + +# Adaptive under-relaxation. The damped iteration p ← p + ω·(σ(V(p)) − p) has local +# slope m = 1 + ω·(g′−1), g′ = σ′(V)·dV/dp ≤ 0 after sign correction. ω is chosen to +# keep m NON-negative (monotone, 0≤m<1, not merely |m|<1): m ≥ θ ⟹ ω ≤ (1−θ)/(1+|g′|). +# |σ′| ≤ |hi−lo|·S/4 bounds g′, so the bound holds at every S without re-measuring. +@inline function _relaxation(d, S::Float64, dVdp::Float64) + lo, hi = parameter_limits(d) + gbound = 0.25 * abs(hi - lo) * S * abs(dVdp) + ω = (1.0 - CONTROL_CONTRACTION) / (1.0 + gbound) + return min(CONTROL_RELAXATION_MAX, ω) +end + +# One damped, sign-corrected proportional update of a single device. Returns +# the magnitude of the parameter change actually applied (for the settling +# test). Returns Inf when the device is oscillation-frozen so the outer loop +# does not count it as settled. Adaptive under-relaxation keeps the iteration +# a contraction at S. +function _step_device!( + d, + idx::Int, + data, + ts::Int, + S::Float64, + pf, + frozen::Vector{Bool}, + dVdp::Vector{Float64}, + osc::Vector{Int}, + prev_sign::Vector{Int}; + kwargs..., +)::Float64 + # Frozen (unreliable probe): hold parameter; a 0 change counts as settled. + frozen[idx] && return 0.0 + n_osc = osc[idx] + if n_osc > CONTROL_OSCILLATION_LIMIT + # Frozen: return Inf so the outer loop never counts this device as settled. + return Inf + end + Vc = _read_vmag(data, controlled_bus_ix(d), ts) + p_now = current_parameter(d) + lo, hi = parameter_limits(d) + dv = dVdp[idx] + ω = _relaxation(d, S, dv) + p_tgt = _control_target(d, Vc, S, dv) + # Track sign reversals to detect oscillation. + s = Int(sign(p_tgt - p_now)) + ps = prev_sign[idx] + if ps != 0 && s != 0 && s != ps + new_n = n_osc + 1 + osc[idx] = new_n + # Fires once: the n_osc > LIMIT early-return prevents osc[idx] passing LIMIT+1. + new_n == CONTROL_OSCILLATION_LIMIT + 1 && + @warn "discrete control: device $(d.name) is oscillating \ + ($(new_n) sign reversals); freezing at current parameter $(p_now) \ + (time step $ts)" + end + prev_sign[idx] = s + osc[idx] > CONTROL_OSCILLATION_LIMIT && return Inf + p_new = clamp(p_now + ω * (p_tgt - p_now), lo, hi) + reached = _continuation_to!(d, data, ts, p_new, pf; kwargs...) + set_current_parameter!(d, reached) + return abs(reached - p_now) +end + +function _control_continuation!( + pf, + data, + ts::Int; + kwargs..., +)::Bool + set = data.controlled_devices::ControlledDeviceSet + converged = _solve_with_q_limits!(pf, data, ts; kwargs...) + converged || return false + + n_taps = length(set.taps) + n_dev = n_taps + length(set.shunts) + # Per-device state indexed in parallel with [taps; shunts] (plain vectors, no + # per-iteration string hashing). dVdp is measured once from the converged base + # point: its sign sets the negative-feedback orientation, its magnitude drives + # under-relaxation. Unreliable-probe devices are frozen, not stepped. + dVdp = zeros(n_dev) + frozen = fill(false, n_dev) + osc = zeros(Int, n_dev) + prev_sign = zeros(Int, n_dev) + for (i, d) in enumerate(set.taps) + s, reliable = _plant_sign(d, data, ts, pf; kwargs...) + reliable ? (dVdp[i] = s) : (frozen[i] = true) + end + for (j, d) in enumerate(set.shunts) + s, reliable = _plant_sign(d, data, ts, pf; kwargs...) + reliable ? (dVdp[n_taps + j] = s) : (frozen[n_taps + j] = true) + end + if any(frozen) + frozen_names = join( + vcat( + [set.taps[i].name for i in 1:n_taps if frozen[i]], + [set.shunts[j].name + for j in eachindex(set.shunts) if frozen[n_taps + j]], + ), + ", ", + ) + @warn "discrete control: $(count(frozen)) device(s) had an unreliable \ + plant-sensitivity probe and were frozen at their current parameter \ + (time step $ts): $frozen_names" + end + + S = INITIAL_CONTROL_STEEPNESS + regulation_complete = false + for _ in 1:MAX_CONTROL_OUTER_ITERATIONS + settled = true + for (i, d) in enumerate(set.taps) + g = _step_device!( + d, i, data, ts, S, pf, frozen, dVdp, osc, prev_sign; kwargs...) + g >= CONTROL_PARAM_TOL && (settled = false) + end + for (j, d) in enumerate(set.shunts) + g = _step_device!( + d, n_taps + j, data, ts, S, pf, frozen, dVdp, osc, prev_sign; + kwargs...) + g >= CONTROL_PARAM_TOL && (settled = false) + end + if settled + if S >= MAX_CONTROL_STEEPNESS + regulation_complete = true + break + end + S = min(S * CONTROL_STEEPNESS_GROWTH, MAX_CONTROL_STEEPNESS) + end + end + if !regulation_complete + @warn "discrete control: reached MAX_CONTROL_OUTER_ITERATIONS without \ + full-steepness convergence (S=$S, target=$MAX_CONTROL_STEEPNESS); \ + regulation may be loose at time step $ts" + end + return snap_and_restore!(pf, data, set, ts; kwargs...) +end + +# Incremental λ-restore of one device from its snapped value toward the +# continuous value (used only if snapping made the system infeasible). First +# probe is a SMALL step, matching `_continuation_to!`. +function _restore_one!(d, data, ts::Int, continuous::Float64, pf; kwargs...)::Bool + snapped = current_parameter(d) + abs(continuous - snapped) < CONTROL_PARAM_TOL && + return _solve_with_q_limits!(pf, data, ts; kwargs...) + done = 0.0 + step = MIN_LAMBDA_STEP + last_good = snapped + while done < 1.0 + trial = min(1.0, done + step) + p = snapped + trial * (continuous - snapped) + apply_parameter!(d, data, p, ts) + if _solve_with_q_limits!(pf, data, ts; kwargs...) + done = trial + last_good = p + step = min(step * CONTROL_STEP_GROWTH, MAX_LAMBDA_STEP) + else + step /= 2.0 + if step < MIN_LAMBDA_STEP + # Leave data at the last converged parameter, not the failed trial. + apply_parameter!(d, data, last_good, ts) + return false + end + end + end + return _solve_with_q_limits!(pf, data, ts; kwargs...) +end + +function snap_and_restore!( + pf, + data, + set::ControlledDeviceSet, + ts::Int; + kwargs..., +)::Bool + cont = Dict{String, Float64}() + for d in set.taps + cont[d.name] = d.current + apply_parameter!(d, data, snap_to_discrete(d, d.current), ts) + end + for d in set.shunts + cont[d.name] = d.current + apply_parameter!(d, data, snap_to_discrete(d, d.current), ts) + end + _solve_with_q_limits!(pf, data, ts; kwargs...) && return true + ok = true + for d in set.taps + ok &= _restore_one!(d, data, ts, cont[d.name], pf; kwargs...) + end + for d in set.shunts + ok &= _restore_one!(d, data, ts, cont[d.name], pf; kwargs...) + end + if !ok + names = join( + vcat([d.name for d in set.taps], [d.name for d in set.shunts]), + ", ", + ) + @error "discrete control could not restore feasibility after \ + snapping; devices: $names (time step $ts)" + data.converged[ts] = false + return false + end + return true +end diff --git a/src/discrete_control/control_metadata.jl b/src/discrete_control/control_metadata.jl new file mode 100644 index 00000000..cfcb5e99 --- /dev/null +++ b/src/discrete_control/control_metadata.jl @@ -0,0 +1,206 @@ +"""Validate shunt invariants at construction time. Separated to allow direct unit testing.""" +function _validate_shunt!( + name::String, + b_min::Float64, + b0::Float64, + b_max::Float64, + steps::Vector{Int}, + dB::Vector{Float64}, +) + if !(b_min <= b0 <= b_max) + error( + "ControlledSwitchedShunt \"$name\": b0=$b0 is outside [b_min=$b_min, b_max=$b_max]", + ) + end + if b_min == b_max + error( + "ControlledSwitchedShunt \"$name\": no controllable susceptance range (b_min == b_max == $b_min); not a valid voltage-controlled device", + ) + end + for k in eachindex(steps, dB) + if steps[k] == 0 && dB[k] != 0.0 + error( + "ControlledSwitchedShunt \"$name\": block $k has zero steps but nonzero dB=$(dB[k]) — malformed metadata", + ) + end + end + return nothing +end + +"""Validate tap invariants at construction time. Separated to allow direct unit testing.""" +function _validate_tap!( + name::String, + p_min::Float64, + p_max::Float64, + ntp::Int, +) + if p_min > p_max + error( + "ControlledTap \"$name\": p_min=$p_min exceeds p_max=$p_max — malformed tap-ratio limits", + ) + end + if p_min == p_max + error( + "ControlledTap \"$name\": no controllable tap-ratio range (p_min == p_max == $p_min); not a valid voltage-controlled device", + ) + end + if ntp < 2 + error( + "ControlledTap \"$name\": needs at least 2 tap positions (got ntp=$ntp)", + ) + end + return nothing +end + +function _controlled_bus_number(ext::Dict, fallback::Int) + for k in ("NREG", "SWREM", "RMIDNT") + if haskey(ext, k) + v = ext[k] + n = v isa Integer ? Int(v) : tryparse(Int, string(v)) + (n !== nothing && n != 0) && return n + end + end + return fallback +end + +function _ext_float(ext::Dict, key::String, default::Float64) + haskey(ext, key) || return default + v = ext[key] + return v isa Real ? Float64(v) : something(tryparse(Float64, string(v)), default) +end + +function _ext_int(ext::Dict, key::String, default::Int) + haskey(ext, key) || return default + v = ext[key] + return v isa Integer ? Int(v) : something(tryparse(Int, string(v)), default) +end + +"""Build the type-stable device set from a `PSY.System`. + +`bus_lookup` maps PSY bus number → network index; `ybus` is the assembled +`AC_Ybus_Matrix` from `data.power_network_matrix`.""" +function build_controlled_device_set( + sys, + bus_lookup::Dict{Int, Int}, + ybus, +) + taps = ControlledTap[] + for tx in PSY.get_available_components(PSY.TapTransformer, sys) + PSY.get_control_objective(tx) == PSY.TransformerControlObjective.VOLTAGE || + continue + arc = PSY.get_arc(tx) + fb = PSY.get_number(PSY.get_from(arc)) + tb = PSY.get_number(PSY.get_to(arc)) + ext = PSY.get_ext(tx) + cbus = _controlled_bus_number(ext, tb) + haskey(bus_lookup, cbus) || + error("ControlledTap $(PSY.get_name(tx)): controlled bus $cbus not in network") + fix = bus_lookup[fb] + tix = bus_lookup[tb] + cix = bus_lookup[cbus] + pmin = _ext_float(ext, "RMI", DEFAULT_TAP_RATIO_MIN) + pmax = _ext_float(ext, "RMA", DEFAULT_TAP_RATIO_MAX) + ntp = Int(_ext_float(ext, "NTP", Float64(DEFAULT_TAP_POSITIONS))) + ntp < 2 && (ntp = DEFAULT_TAP_POSITIONS) + _validate_tap!(PSY.get_name(tx), pmin, pmax, ntp) + vset = if haskey(ext, "VSET") + _ext_float(ext, "VSET", DEFAULT_TAP_VSET) + else + controlled_bus = if cbus == fb + PSY.get_from(arc) + elseif cbus == tb + PSY.get_to(arc) + else + PSY.get_bus(sys, cbus) + end + m = PSY.get_magnitude(controlled_bus) + m > 0.0 ? m : DEFAULT_TAP_VSET + end + yt = 1.0 / (PSY.get_r(tx) + PSY.get_x(tx) * im) + ysh = PSY.get_primary_shunt(tx) + push!( + taps, + ControlledTap( + PSY.get_name(tx), + fix, + tix, + cix, + cix == fix, + vset, + yt, + ysh, + 0.0, + pmin, + pmax, + collect(range(pmin, pmax; length = ntp)), + _ybus_block_offsets(ybus, fix, tix), + PSY.get_tap(tx), + ), + ) + end + + shunts = ControlledSwitchedShunt[] + for sa in PSY.get_available_components(PSY.SwitchedAdmittance, sys) + bus = PSY.get_number(PSY.get_bus(sa)) + haskey(bus_lookup, bus) || continue + ext = PSY.get_ext(sa) + modsw = _ext_int(ext, "MODSW", DEFAULT_SHUNT_MODSW) + if modsw == 0 + @debug "ControlledSwitchedShunt $(PSY.get_name(sa)): MODSW=0 (locked); \ + treated as fixed admittance, not enrolled." + continue + elseif modsw == 1 + continuous = false + elseif modsw == 2 + continuous = true + else + error( + "ControlledSwitchedShunt $(PSY.get_name(sa)): MODSW=$modsw \ + (remote reactive-power / remote-device control) is not supported. \ + Only voltage-control modes are implemented: MODSW=1 (discrete) \ + and MODSW=2 (continuous).", + ) + end + cbus = _controlled_bus_number(ext, bus) + haskey(bus_lookup, cbus) || + error( + "ControlledSwitchedShunt $(PSY.get_name(sa)): bus $cbus not in network", + ) + lims = PSY.get_admittance_limits(sa) + vset = (lims.min + lims.max) / 2.0 + Y0 = PSY.get_Y(sa) + steps = PSY.get_number_of_steps(sa) + dY = PSY.get_Y_increase(sa) + dB = imag.(dY) + b0 = imag(Y0) + g0 = real(Y0) + bmax = b0 + sum(max.(steps .* dB, 0.0); init = 0.0) + bmin = b0 + sum(min.(steps .* dB, 0.0); init = 0.0) + _validate_shunt!(PSY.get_name(sa), bmin, b0, bmax, steps, dB) + init_status = PSY.get_initial_status(sa) + current_b = imag(Y0 + sum(init_status .* dY; init = 0.0 + 0.0im)) + bo = sortperm(dB; rev = true) + bn = zeros(Int, length(dB)) + push!( + shunts, + ControlledSwitchedShunt( + PSY.get_name(sa), + bus_lookup[bus], + bus_lookup[cbus], + vset, + g0, + b0, + steps, + dB, + bmin, + bmax, + bo, + bn, + continuous, + current_b, + ), + ) + end + + return ControlledDeviceSet(taps, shunts) +end diff --git a/src/discrete_control/controlled_devices.jl b/src/discrete_control/controlled_devices.jl new file mode 100644 index 00000000..88b0ceeb --- /dev/null +++ b/src/discrete_control/controlled_devices.jl @@ -0,0 +1,184 @@ +abstract type AbstractControlledDevice end +abstract type AbstractBranchControl <: AbstractControlledDevice end +abstract type AbstractShuntControl <: AbstractControlledDevice end + +"""Voltage-controlling tap transformer. `nz_offsets` are the 4 cached +`nzval` linear indices of the (from,to)×(from,to) Y-bus block.""" +mutable struct ControlledTap <: AbstractBranchControl + name::String + from_ix::Int + to_ix::Int + controlled_ix::Int + controlled_on_primary::Bool # true → eq.46, false → eq.47 + vset::Float64 + yt::ComplexF64 # 1/(r+jx) + y_shunt::ComplexF64 # primary shunt + alpha::Float64 # phase-shift angle (0 for TapTransformer) + p_min::Float64 + p_max::Float64 + levels::Vector{Float64} # discrete tap ratios + nz_offsets::NTuple{4, Int} # nzval idx for Y11,Y12,Y21,Y22 + current::Float64 +end + +"""Voltage-controlling switched shunt, snapped block-greedily.""" +mutable struct ControlledSwitchedShunt <: AbstractShuntControl + name::String + bus_ix::Int + controlled_ix::Int + vset::Float64 + g0::Float64 # real(get_Y) + b0::Float64 # imag(get_Y) + block_steps::Vector{Int} # number_of_steps per block + block_dB::Vector{Float64} # imag(Y_increase) per block + b_min::Float64 + b_max::Float64 + block_order::Vector{Int} # sortperm(block_dB; rev=true), cached at construction + block_n::Vector{Int} # per-block chosen step counts, reused in-place each snap + continuous::Bool # MODSW==2 ⇒ continuous regulation (no discrete snap) + current::Float64 # current total susceptance b +end + +# Reserved seams — not implemented in this scope. +mutable struct ControlledPhaseShifter <: AbstractBranchControl end +mutable struct ControlledFACTS <: AbstractShuntControl end + +struct ControlledDeviceSet + taps::Vector{ControlledTap} + shunts::Vector{ControlledSwitchedShunt} +end +Base.isempty(s::ControlledDeviceSet) = isempty(s.taps) && isempty(s.shunts) + +controlled_bus_ix(d::ControlledTap) = d.controlled_ix +controlled_bus_ix(d::ControlledSwitchedShunt) = d.controlled_ix +voltage_setpoint(d::ControlledTap) = d.vset +voltage_setpoint(d::ControlledSwitchedShunt) = d.vset +parameter_limits(d::ControlledTap) = (d.p_min, d.p_max) +parameter_limits(d::ControlledSwitchedShunt) = (d.b_min, d.b_max) +current_parameter(d::ControlledTap) = d.current +current_parameter(d::ControlledSwitchedShunt) = d.current +set_current_parameter!(d::ControlledTap, p::Float64) = (d.current = p; nothing) +set_current_parameter!(d::ControlledSwitchedShunt, p::Float64) = + (d.current = p; nothing) + +# Seam: future implicit embedding dispatches here. Never called by the outer loop. +stamp_control!(d::AbstractControlledDevice, args...) = + error("implicit embedding not implemented for $(typeof(d))") +_seam_err(d) = error("control not implemented: $(typeof(d))") +controlled_bus_ix(d::Union{ControlledPhaseShifter, ControlledFACTS}) = _seam_err(d) +voltage_setpoint(d::Union{ControlledPhaseShifter, ControlledFACTS}) = _seam_err(d) +parameter_limits(d::Union{ControlledPhaseShifter, ControlledFACTS}) = _seam_err(d) +current_parameter(d::Union{ControlledPhaseShifter, ControlledFACTS}) = _seam_err(d) +set_current_parameter!(d::Union{ControlledPhaseShifter, ControlledFACTS}, ::Float64) = + _seam_err(d) + +function _nz_index(A::SparseArrays.SparseMatrixCSC, row::Int, col::Int) + @inbounds for k in SparseArrays.nzrange(A, col) + A.rowval[k] == row && return k + end + error("Ybus has no stored entry at ($row,$col); structural zero") +end + +function _ybus_block_offsets(ybus, i::Int, j::Int) + A = ybus.data + return ( + _nz_index(A, i, i), + _nz_index(A, i, j), + _nz_index(A, j, i), + _nz_index(A, j, j), + ) +end + +# Y11/Y12/Y21 are delta-updated (new−old) so parallel branches' contributions +# in the shared nzval slots are preserved; Y22=Yt is tap-independent (zero delta, +# even with parallels) so it is skipped. The nzval is single-precision ComplexF32; +# correctness of the running sum relies on `old_tap` coming from `d.current` +# (the last applied value), never read back from the lossy nzval. +function apply_parameter!(d::ControlledTap, data, p::Float64, ::Int) + A = data.power_network_matrix.data + old_tap = d.current * cis(d.alpha) + new_tap = p * cis(d.alpha) + o = d.nz_offsets + @inbounds begin + A.nzval[o[1]] += d.yt / abs2(new_tap) - d.yt / abs2(old_tap) + A.nzval[o[2]] += -d.yt / conj(new_tap) - (-d.yt / conj(old_tap)) + A.nzval[o[3]] += -d.yt / new_tap - (-d.yt / old_tap) + # Y22 = Yt is tap-independent; no update needed. + end + d.current = p + return nothing +end + +function apply_parameter!(d::ControlledSwitchedShunt, data, b::Float64, ts::Int) + data.bus_active_power_constant_impedance_withdrawals[d.bus_ix, ts] = d.g0 + data.bus_reactive_power_constant_impedance_withdrawals[d.bus_ix, ts] = -b + d.current = b + return nothing +end + +apply_parameter!(d::Union{ControlledPhaseShifter, ControlledFACTS}, args...) = + _seam_err(d) + +@inline function _sigmoid(lo::Float64, hi::Float64, S::Float64, + x::Float64, xset::Float64) + return (hi - lo) / (1.0 + exp(S * (x - xset))) + lo +end + +# Branch (tap): controlled-on-primary uses eq.46 (lo=tr_min,hi=tr_max); +# controlled-on-secondary uses eq.47 (limits swapped). +function target_from_voltage(d::ControlledTap, vmag::Float64, S::Float64) + lo, hi = d.controlled_on_primary ? (d.p_min, d.p_max) : (d.p_max, d.p_min) + return clamp(_sigmoid(lo, hi, S, vmag, d.vset), d.p_min, d.p_max) +end + +# Shunt: eq.9, low V → high B. x→-∞ gives hi=b_max; x→+∞ gives lo=b_min. +function target_from_voltage(d::ControlledSwitchedShunt, vmag::Float64, + S::Float64) + b = _sigmoid(d.b_min, d.b_max, S, vmag, d.vset) + return clamp(b, d.b_min, d.b_max) +end + +function snap_to_discrete(d::ControlledTap, p::Float64) + pc = clamp(p, d.p_min, d.p_max) + best = d.levels[1] + @inbounds for lv in d.levels + abs(lv - pc) < abs(best - pc) && (best = lv) + end + return best +end + +# Block-greedy (floor): largest blocks first, take as many steps as fit without +# overshooting; ±1 bounded refinement then corrects any under-committed block. +# block_order and block_n are pre-allocated fields — no per-call heap allocation. +function snap_to_discrete(d::ControlledSwitchedShunt, b::Float64) + d.continuous && return clamp(b, d.b_min, d.b_max) # continuous: no grid snap + target_clamped = clamp(b, d.b_min, d.b_max) + target = target_clamped - d.b0 + total = d.b0 + # Greedy pass: floor to avoid overshooting. + @inbounds for k in d.block_order + dB = d.block_dB[k] + dB == 0.0 && (d.block_n[k] = 0; continue) + n = clamp(floor(Int, (target - (total - d.b0)) / dB), + 0, d.block_steps[k]) + d.block_n[k] = n + total += n * dB + end + # ±1 bounded refinement: single pass, in-place update of total and block_n. + @inbounds for k in d.block_order + dB = d.block_dB[k] + dB == 0.0 && continue + n = d.block_n[k] + for δ in (-1, 1) + n_new = n + δ + (n_new < 0 || n_new > d.block_steps[k]) && continue + candidate = total + δ * dB + if abs(candidate - target_clamped) < abs(total - target_clamped) + d.block_n[k] = n_new + total = candidate + n = n_new + end + end + end + return clamp(total, d.b_min, d.b_max) +end diff --git a/src/levenberg-marquardt.jl b/src/levenberg-marquardt.jl index c7860a25..3556f819 100644 --- a/src/levenberg-marquardt.jl +++ b/src/levenberg-marquardt.jl @@ -19,6 +19,9 @@ mutable struct LMWorkspace # Marquardt diagonal scaling (length n). All-ones ⇒ √λ·I. D::Vector{Float64} marquardt_scaling::Bool + # Per-iteration scratch: temp_x = Rv + J·Δx (m); x_trial = x + Δx (n). + temp_x::Vector{Float64} + x_trial::Vector{Float64} end """Build the augmented matrix `[J; D]` once, recording which `A.nzval` entries @@ -62,7 +65,8 @@ function LMWorkspace( D = marquardt_scaling ? zeros(n) : ones(n) ws = LMWorkspace( - A, j_nzval_indices, λ_diag_indices, F, b, D, marquardt_scaling) + A, j_nzval_indices, λ_diag_indices, F, b, D, marquardt_scaling, + Vector{Float64}(undef, m), Vector{Float64}(undef, n)) if marquardt_scaling update_column_scale!(ws, Jv) end @@ -222,15 +226,18 @@ function compute_error( m = length(residual.Rv) @assert m == length(ws.b) - size(J.Jv, 2) "residual/J size mismatch vs preallocated LM buffer (m=$m, buf=$(length(ws.b)), n=$(size(J.Jv, 2)))" @views ws.b[1:m] .= .-residual.Rv # bottom n entries stay zero from construction + # Δx left allocating: SPQR has no in-place reuse, and the QR rebuild dominates anyway. Δx = ws.F \ ws.b - temp_x = residual.Rv .+ J.Jv * Δx + # temp_x = Rv + J·Δx + LinearAlgebra.mul!(ws.temp_x, J.Jv, Δx) + ws.temp_x .+= residual.Rv - x_trial = x .+ Δx - residual(x_trial, time_step) # M(x_c + Δx) + ws.x_trial .= x .+ Δx + residual(ws.x_trial, time_step) # M(x_c + Δx) newResidualSize = dot(residual.Rv, residual.Rv) - predicted_reduction = residualSize - dot(temp_x, temp_x) + predicted_reduction = residualSize - dot(ws.temp_x, ws.temp_x) actual_reduction = residualSize - newResidualSize # Guard against zero/negative predicted reduction. diff --git a/src/linear_solver_backend.jl b/src/linear_solver_backend.jl index 07be53b1..5e211403 100644 --- a/src/linear_solver_backend.jl +++ b/src/linear_solver_backend.jl @@ -27,12 +27,15 @@ mutable struct PardisoLinSolveCache end """Union of the KLU, AppleAccelerate, and MKLPardiso solver caches. Every member is concrete so -the 3-way union stays within Julia's small-union splitting; the KLU member is pinned to -`J_INDEX_TYPE` because that is the only instantiation flowing through these methods (the AC Newton -cache and its fallback are built from `J.Jv::SparseMatrixCSC{Float64, J_INDEX_TYPE}`).""" +the 4-way union stays within Julia's small-union splitting. Both KLU index types are listed: the +AC Newton cache and its fallback are built from `J.Jv::SparseMatrixCSC{Float64, J_INDEX_TYPE}` +(`Int32` off Apple, `Int64` on Apple), while PNM's DC ABA factorization is always +`KLULinSolveCache{Float64, Int64}` regardless of platform — so the DC solve path needs the `Int64` +member even where `J_INDEX_TYPE === Int32`.""" const PFLinearSolverCache = Union{ - PNM.KLULinSolveCache{Float64, J_INDEX_TYPE}, + PNM.KLULinSolveCache{Float64, Int32}, + PNM.KLULinSolveCache{Float64, Int64}, PNM.AAFactorCache, PardisoLinSolveCache, } diff --git a/src/mixed_cpb_power_flow_jacobian.jl b/src/mixed_cpb_power_flow_jacobian.jl index 0eaf2fa8..628740e4 100644 --- a/src/mixed_cpb_power_flow_jacobian.jl +++ b/src/mixed_cpb_power_flow_jacobian.jl @@ -12,7 +12,6 @@ hot path is `O(N + n_LCC)`. Field roles are in the inline comments below. """ struct ACMixedCPBJacobian data::ACPowerFlowData - Jf!::Function Jv::SparseMatrixCSC{Float64, J_INDEX_TYPE} Y_bus_eff::SparseMatrixCSC{ComplexF64, Int} Y_diag::Vector{ComplexF64} # cached Y_bus_eff diagonal; avoids O(log nnz) sparse access per iteration @@ -89,7 +88,6 @@ function ACMixedCPBJacobian( ) J = ACMixedCPBJacobian( residual.data, - _update_mixed_cpb_jacobian_values!, Jv0, residual.Y_bus_eff, Y_diag, @@ -121,7 +119,7 @@ function ACMixedCPBJacobian( end function (J::ACMixedCPBJacobian)(time_step::Int64) - J.Jf!(J.Jv, J.data, J.Y_diag, + _update_mixed_cpb_jacobian_values!(J.Jv, J.data, J.Y_diag, J.e_state, J.f_state, J.P_eff_cache, J.Q_eff_cache, J.const_I_P, J.const_I_Q, J.Ir_acc, J.Ii_acc, J.bus_slack_participation_factors, @@ -136,7 +134,7 @@ function (J::ACMixedCPBJacobian)( Jv::SparseMatrixCSC{Float64, J_INDEX_TYPE}, time_step::Int64, ) - J.Jf!(J.Jv, J.data, J.Y_diag, + _update_mixed_cpb_jacobian_values!(J.Jv, J.data, J.Y_diag, J.e_state, J.f_state, J.P_eff_cache, J.Q_eff_cache, J.const_I_P, J.const_I_Q, J.Ir_acc, J.Ii_acc, J.bus_slack_participation_factors, diff --git a/src/mixed_cpb_power_flow_residual.jl b/src/mixed_cpb_power_flow_residual.jl index d4bf80c3..121c9cc3 100644 --- a/src/mixed_cpb_power_flow_residual.jl +++ b/src/mixed_cpb_power_flow_residual.jl @@ -14,7 +14,6 @@ named after their roles. """ struct ACMixedCPBResidual data::ACPowerFlowData - Rf!::Function Rv::Vector{Float64} Y_bus_eff::SparseMatrixCSC{ComplexF64, Int} P_net_const::Vector{Float64} @@ -89,7 +88,6 @@ function ACMixedCPBResidual(data::ACPowerFlowData, time_step::Int64) return ACMixedCPBResidual( data, - _update_mixed_cpb_residual_values!, Vector{Float64}(undef, total_state), Y_bus_eff, P_net_const, @@ -117,7 +115,7 @@ function (R::ACMixedCPBResidual)( x::Vector{Float64}, time_step::Int64, ) - R.Rf!(R.Rv, x, R.Y_bus_eff, R.P_net_const, R.Q_net_const, + _update_mixed_cpb_residual_values!(R.Rv, x, R.Y_bus_eff, R.P_net_const, R.Q_net_const, R.const_I_P, R.const_I_Q, R.P_net_set, R.bus_slack_participation_factors, R.subnetworks, R.bus_state_offset, R.bus_block_size, R.total_bus_state, @@ -128,7 +126,7 @@ function (R::ACMixedCPBResidual)( end function (R::ACMixedCPBResidual)(x::Vector{Float64}, time_step::Int64) - R.Rf!(R.Rv, x, R.Y_bus_eff, R.P_net_const, R.Q_net_const, + _update_mixed_cpb_residual_values!(R.Rv, x, R.Y_bus_eff, R.P_net_const, R.Q_net_const, R.const_I_P, R.const_I_Q, R.P_net_set, R.bus_slack_participation_factors, R.subnetworks, R.bus_state_offset, R.bus_block_size, R.total_bus_state, diff --git a/src/power_flow_method.jl b/src/power_flow_method.jl index 5c8e4b84..eec672f4 100644 --- a/src/power_flow_method.jl +++ b/src/power_flow_method.jl @@ -192,6 +192,7 @@ function _dogleg!(Δx_proposed::Vector{Float64}, Jv::SparseMatrixCSC{Float64, J_INDEX_TYPE}, d::Vector{Float64}, delta::Float64, + Jg::Vector{Float64}, ) nr_norm = wnorm(d, Δx_nr) @debug "Trust region: ||Δx_nr|| = $(siground(nr_norm)), δ = $(siground(delta))" @@ -204,7 +205,8 @@ function _dogleg!(Δx_proposed::Vector{Float64}, g = Δx_proposed LinearAlgebra.mul!(g, Jv', r) g .= g ./ d .^ 2 - Δx_cauchy .= -wnorm(d, g)^2 / sum(abs2, Jv * g) .* g # Cauchy point + LinearAlgebra.mul!(Jg, Jv, g) + Δx_cauchy .= -wnorm(d, g)^2 / sum(abs2, Jg) .* g # Cauchy point cauchy_norm = wnorm(d, Δx_cauchy) @debug "Cauchy point: ||Δx_cauchy|| = $(siground(cauchy_norm))" @@ -264,13 +266,14 @@ function _iwamoto_fallback!( J::Union{ACPowerFlowJacobian, ACRectangularCIJacobian, ACMixedCPBJacobian}, old_residual::Vector{Float64}, old_residual_norm::Float64, - new_residual_norm::Float64, autoscale::Bool, )::Bool g0 = old_residual_norm - g1 = dot(old_residual, residual.Rv) - g2 = new_residual_norm - μ = _iwamoto_multiplier(g0, g1, g2) + # Quadratic model F(x+μΔx) = f₀ + μ·(J·Δx) + μ²·a along Δx_proposed. r_predict + # (= f₀ + J·Δx) from the ρ test gives a = F(x+Δx) − r_predict for free (no extra matvec). + c_fb, c_bb, c_fa, c_ba, c_aa = + _iwamoto_quadratic_dots(old_residual, stateVector.r_predict, residual.Rv) + μ = _iwamoto_multiplier(2.0 * c_fb, c_bb + 2.0 * c_fa, 2.0 * c_ba, c_aa) # Revert full step, apply damped step in a single fused pass. @. stateVector.x += (μ - 1.0) * stateVector.Δx_proposed residual(stateVector.x, time_step) @@ -323,6 +326,7 @@ function _trust_region_step(time_step::Int, J.Jv, stateVector.d, delta, + stateVector.r_predict, # scratch for J·g; overwritten with the true r_predict below ) # find proposed next point. stateVector.x .+= stateVector.Δx_proposed @@ -338,9 +342,14 @@ function _trust_region_step(time_step::Int, # Ratio of actual to predicted reduction LinearAlgebra.mul!(stateVector.r_predict, J.Jv, stateVector.Δx_proposed) stateVector.r_predict .+= stateVector.r - rho = - (sum(abs2, stateVector.r) - sum(abs2, residual.Rv)) / - (sum(abs2, stateVector.r) - sum(abs2, stateVector.r_predict)) + predicted_reduction = old_residual_norm - sum(abs2, stateVector.r_predict) + # Non-positive predicted reduction (e.g. a non-descent fallback step) would give a + # NaN/Inf ρ that stalls the solver; force a rejected-step ρ to shrink the trust region. + rho = if predicted_reduction > 0.0 + (old_residual_norm - new_residual_norm) / predicted_reduction + else + -Inf + end @debug "Trust region step: ρ = $(siground(rho)), η = $(siground(eta)), ||Δx|| = $(siground(norm(stateVector.Δx_proposed)))" @@ -356,7 +365,7 @@ function _trust_region_step(time_step::Int, if iwamoto_fallback iwamoto_accepted = _iwamoto_fallback!( time_step, stateVector, residual, J, - oldResidual, old_residual_norm, new_residual_norm, autoscale) + oldResidual, old_residual_norm, autoscale) if iwamoto_accepted # Iwamoto accepted a damped step — shrink trust region since the # full proposed step was rejected by rho. Do not use rho-based @@ -390,28 +399,50 @@ function _trust_region_step(time_step::Int, return delta end -"""Evaluate the Iwamoto objective g(μ) = ‖(1-μ)f₀ + μ²f₁‖² expanded as -g(μ) = (1-μ)²g₀ + 2μ²(1-μ)g₁ + μ⁴g₂ where g₀=‖f₀‖², g₁=f₀ᵀf₁, g₂=‖f₁‖².""" +"""Inner products for the quadratic model `F(x+μΔx) = f₀ + μ·b + μ²·a` with +`b = J·Δx`, `a = F(x+Δx) − f₀ − b`. From `f₀`, `rpred = f₀ + b`, and `rv = F(x+Δx)` +returns `(f₀·b, b·b, f₀·a, b·a, a·a)`.""" +@inline function _iwamoto_quadratic_dots( + f0::Vector{Float64}, rpred::Vector{Float64}, rv::Vector{Float64}, +)::NTuple{5, Float64} + c_fb = 0.0 + c_bb = 0.0 + c_fa = 0.0 + c_ba = 0.0 + c_aa = 0.0 + @inbounds @simd for i in eachindex(f0, rpred, rv) + b = rpred[i] - f0[i] + a = rv[i] - rpred[i] + c_fb += f0[i] * b + c_bb += b * b + c_fa += f0[i] * a + c_ba += b * a + c_aa += a * a + end + return c_fb, c_bb, c_fa, c_ba, c_aa +end + +"""Iwamoto objective minus its μ-independent constant: +`g̃(μ) = q₁μ + q₂μ² + q₃μ³ + q₄μ⁴`. Dropping the constant preserves the minimizer.""" @inline function _iwamoto_objective( - μ::Float64, g0::Float64, g1::Float64, g2::Float64, + μ::Float64, q1::Float64, q2::Float64, q3::Float64, q4::Float64, )::Float64 - om = 1.0 - μ - μ2 = μ * μ - return om * om * g0 + 2.0 * μ2 * om * g1 + μ2 * μ2 * g2 + return μ * (q1 + μ * (q2 + μ * (q3 + μ * q4))) end -"""If μ ∈ [IWAMOTO_MU_MIN, IWAMOTO_MU_MAX] and g(μ) < best_g, return the -improved (μ, g(μ)); otherwise return (best_μ, best_g) unchanged.""" +"""If μ ∈ [IWAMOTO_MU_MIN, IWAMOTO_MU_MAX] and g̃(μ) < best_g, return the +improved (μ, g̃(μ)); otherwise return (best_μ, best_g) unchanged.""" @inline function _try_iwamoto_candidate( μ::Float64, best_μ::Float64, best_g::Float64, - g0::Float64, - g1::Float64, - g2::Float64, + q1::Float64, + q2::Float64, + q3::Float64, + q4::Float64, )::Tuple{Float64, Float64} if IWAMOTO_MU_MIN <= μ <= IWAMOTO_MU_MAX - gval = _iwamoto_objective(μ, g0, g1, g2) + gval = _iwamoto_objective(μ, q1, q2, q3, q4) if gval < best_g return μ, gval end @@ -419,23 +450,23 @@ improved (μ, g(μ)); otherwise return (best_μ, best_g) unchanged.""" return best_μ, best_g end -"""Compute the optimal Iwamoto step multiplier μ ∈ [IWAMOTO_MU_MIN, IWAMOTO_MU_MAX] -by minimizing g(μ) = (1-μ)²g₀ + 2μ²(1-μ)g₁ + μ⁴g₂. - -The stationary points satisfy the cubic g'(μ)/2 = 2g₂μ³ - 3g₁μ² + (g₀+2g₁)μ - g₀ = 0. -All real roots are found analytically via the depressed-cubic trigonometric/Cardano form, -and the global minimizer of g over the domain is returned. O(1), zero-allocation.""" -function _iwamoto_multiplier(g0::Float64, g1::Float64, g2::Float64)::Float64 +"""Optimal Iwamoto multiplier μ ∈ [IWAMOTO_MU_MIN, IWAMOTO_MU_MAX] minimizing +`g̃(μ) = q₁μ + q₂μ² + q₃μ³ + q₄μ⁴` (coefficients from [`_iwamoto_quadratic_dots`](@ref)). +Stationary points solve the cubic `g̃'(μ) = 4q₄μ³ + 3q₃μ² + 2q₂μ + q₁ = 0`, found +analytically (depressed-cubic Cardano/trig form). Exact for the dogleg step; +reduces to classical Iwamoto & Tamura (1981) when `b = −f₀` (Newton step).""" +function _iwamoto_multiplier(q1::Float64, q2::Float64, q3::Float64, q4::Float64)::Float64 # Initialize best candidate from domain boundaries. best_μ = IWAMOTO_MU_MIN - best_g = _iwamoto_objective(IWAMOTO_MU_MIN, g0, g1, g2) - best_μ, best_g = _try_iwamoto_candidate(IWAMOTO_MU_MAX, best_μ, best_g, g0, g1, g2) + best_g = _iwamoto_objective(IWAMOTO_MU_MIN, q1, q2, q3, q4) + best_μ, best_g = + _try_iwamoto_candidate(IWAMOTO_MU_MAX, best_μ, best_g, q1, q2, q3, q4) # Cubic coefficients: c₃μ³ + c₂μ² + c₁μ + c₀ = 0 - c3 = 2.0 * g2 - c2 = -3.0 * g1 - c1 = g0 + 2.0 * g1 - c0 = -g0 + c3 = 4.0 * q4 + c2 = 3.0 * q3 + c1 = 2.0 * q2 + c0 = q1 if abs(c3) < IWAMOTO_DEGENERACY_TOL # Degenerate: solve quadratic c₂μ² + c₁μ + c₀ = 0 @@ -444,11 +475,13 @@ function _iwamoto_multiplier(g0::Float64, g1::Float64, g2::Float64)::Float64 if disc >= 0.0 sq = sqrt(disc) for μ in ((-c1 + sq) / (2.0 * c2), (-c1 - sq) / (2.0 * c2)) - best_μ, best_g = _try_iwamoto_candidate(μ, best_μ, best_g, g0, g1, g2) + best_μ, best_g = + _try_iwamoto_candidate(μ, best_μ, best_g, q1, q2, q3, q4) end end elseif abs(c1) > IWAMOTO_DEGENERACY_TOL - best_μ, best_g = _try_iwamoto_candidate(-c0 / c1, best_μ, best_g, g0, g1, g2) + best_μ, best_g = + _try_iwamoto_candidate(-c0 / c1, best_μ, best_g, q1, q2, q3, q4) end return best_μ end @@ -471,24 +504,24 @@ function _iwamoto_multiplier(g0::Float64, g1::Float64, g2::Float64)::Float64 for k in 0:2 best_μ, best_g = _try_iwamoto_candidate( m * cos(φ3 - 2.0 * π * k / 3.0) - p3, - best_μ, best_g, g0, g1, g2) + best_μ, best_g, q1, q2, q3, q4) end elseif Δ < 0.0 # One real root — Cardano's formula. sqD = sqrt(max(-Δ / 108.0, 0.0)) best_μ, best_g = _try_iwamoto_candidate( cbrt(-B / 2.0 + sqD) + cbrt(-B / 2.0 - sqD) - p3, - best_μ, best_g, g0, g1, g2) + best_μ, best_g, q1, q2, q3, q4) else # Δ ≈ 0 — repeated roots. if abs(A) < IWAMOTO_DEGENERACY_TOL # Triple root at t = 0. - best_μ, best_g = _try_iwamoto_candidate(-p3, best_μ, best_g, g0, g1, g2) + best_μ, best_g = _try_iwamoto_candidate(-p3, best_μ, best_g, q1, q2, q3, q4) else # Simple root t₁ = 3B/A and double root t₂ = -3B/(2A). for t in (3.0 * B / A, -3.0 * B / (2.0 * A)) best_μ, best_g = _try_iwamoto_candidate( - t - p3, best_μ, best_g, g0, g1, g2) + t - p3, best_μ, best_g, q1, q2, q3, q4) end end end @@ -496,6 +529,12 @@ function _iwamoto_multiplier(g0::Float64, g1::Float64, g2::Float64)::Float64 return best_μ end +"""Classical Iwamoto & Tamura (1981) multiplier for the Newton step (`b = −f₀`), +with `g₀ = ‖f₀‖²`, `g₁ = f₀ᵀf₁`, `g₂ = ‖f₁‖²`, `f₁ = F(x+Δx)`.""" +@inline function _iwamoto_multiplier(g0::Float64, g1::Float64, g2::Float64)::Float64 + return _iwamoto_multiplier(-2.0 * g0, g0 + 2.0 * g1, -2.0 * g1, g2) +end + """Does a single iteration of `NewtonRaphsonACPowerFlow`. Updates the `r` and `x` fields of the `stateVector`, and computes the Jacobian at the new `x`.""" function _simple_step(time_step::Int, diff --git a/src/power_flow_types.jl b/src/power_flow_types.jl index a120164d..319020ee 100644 --- a/src/power_flow_types.jl +++ b/src/power_flow_types.jl @@ -144,6 +144,8 @@ with the specified solver type. - `time_step_names::Vector{String}`: Names for each time step. Default is an empty vector. - `correct_bustypes::Bool`: Whether to automatically correct bus types based on available generation. Default is `false`. +- `control_discrete_devices::Bool`: Whether to run discrete device control (tap changers, switched + shunts) via λ-continuation. Default is `false`. - `solver_settings::Dict{Symbol, Any}`: Additional keyword arguments to pass to the solver. Default is an empty dictionary. """ @@ -165,6 +167,7 @@ struct ACPolarPowerFlow{ACSolver <: ACPowerFlowSolverType} <: AbstractACPowerFlo time_steps::Int time_step_names::Vector{String} correct_bustypes::Bool + control_discrete_devices::Bool solver_settings::Dict{Symbol, Any} end @@ -218,6 +221,7 @@ function ACPolarPowerFlow{ACSolver}(; time_steps::Int = 1, time_step_names::Vector{String} = String[], correct_bustypes::Bool = false, + control_discrete_devices::Bool = false, solver_settings::AbstractDict = Dict{Symbol, Any}(), ) where {ACSolver <: ACPowerFlowSolverType} settings = Dict{Symbol, Any}(solver_settings) @@ -243,6 +247,7 @@ function ACPolarPowerFlow{ACSolver}(; time_steps, time_step_names, correct_bustypes, + control_discrete_devices, settings, ) end @@ -276,6 +281,9 @@ get_calculate_voltage_stability_factors(::AbstractACPowerFlow) = false get_calculate_voltage_stability_factors(pf::ACPolarPowerFlow) = pf.calculate_voltage_stability_factors +get_control_discrete_devices(pf::AbstractACPowerFlow) = pf.control_discrete_devices +get_control_discrete_devices(::PowerFlowEvaluationModel) = false + """ ACRectangularPowerFlow{ACSolver}(; kwargs...) where {ACSolver <: ACPowerFlowSolverType} ACRectangularPowerFlow(; kwargs...) @@ -310,6 +318,8 @@ polar state layout and have no current-injection equivalent. - `time_steps::Int`: Default `1`. - `time_step_names::Vector{String}`: Default empty. - `correct_bustypes::Bool`: Default `false`. +- `control_discrete_devices::Bool`: Whether to run discrete device control via λ-continuation. + Default `false`. - `solver_settings::Dict{Symbol, Any}`: Default empty. """ struct ACRectangularPowerFlow{ACSolver <: ACPowerFlowSolverType} <: @@ -328,6 +338,7 @@ struct ACRectangularPowerFlow{ACSolver <: ACPowerFlowSolverType} <: time_steps::Int time_step_names::Vector{String} correct_bustypes::Bool + control_discrete_devices::Bool solver_settings::Dict{Symbol, Any} end @@ -346,6 +357,7 @@ function ACRectangularPowerFlow{ACSolver}(; time_steps::Int = 1, time_step_names::Vector{String} = String[], correct_bustypes::Bool = false, + control_discrete_devices::Bool = false, solver_settings::Dict{Symbol, Any} = Dict{Symbol, Any}(), ) where {ACSolver <: ACPowerFlowSolverType} if ACSolver <: Union{ @@ -378,6 +390,7 @@ function ACRectangularPowerFlow{ACSolver}(; time_steps, time_step_names, correct_bustypes, + control_discrete_devices, solver_settings, ) end @@ -419,6 +432,8 @@ polar state layout and have no mixed current-power equivalent. - `time_steps::Int`: Default `1`. - `time_step_names::Vector{String}`: Default empty. - `correct_bustypes::Bool`: Default `false`. +- `control_discrete_devices::Bool`: Whether to run discrete device control via λ-continuation. + Default `false`. - `solver_settings::Dict{Symbol, Any}`: Default empty. """ struct ACMixedPowerFlow{ACSolver <: ACPowerFlowSolverType} <: @@ -437,6 +452,7 @@ struct ACMixedPowerFlow{ACSolver <: ACPowerFlowSolverType} <: time_steps::Int time_step_names::Vector{String} correct_bustypes::Bool + control_discrete_devices::Bool solver_settings::Dict{Symbol, Any} end @@ -455,6 +471,7 @@ function ACMixedPowerFlow{ACSolver}(; time_steps::Int = 1, time_step_names::Vector{String} = String[], correct_bustypes::Bool = false, + control_discrete_devices::Bool = false, solver_settings::Dict{Symbol, Any} = Dict{Symbol, Any}(), ) where {ACSolver <: ACPowerFlowSolverType} if ACSolver <: Union{ @@ -488,6 +505,7 @@ function ACMixedPowerFlow{ACSolver}(; time_steps, time_step_names, correct_bustypes, + control_discrete_devices, solver_settings, ) end diff --git a/src/rectangular_ci_power_flow_jacobian.jl b/src/rectangular_ci_power_flow_jacobian.jl index 1bfd79da..5d01664c 100644 --- a/src/rectangular_ci_power_flow_jacobian.jl +++ b/src/rectangular_ci_power_flow_jacobian.jl @@ -11,7 +11,6 @@ than `O((N + n_LCC) · log(nnz_per_col))` of `Jv[r, c] = v` setindex. # Fields - `data::ACPowerFlowData` -- `Jf!::Function` — inplace Jacobian update - `Jv::SparseMatrixCSC{Float64, J_INDEX_TYPE}` — Jacobian values - `Y_bus_eff::SparseMatrixCSC{ComplexF64, Int}` — Y_bus with ZIP-Z folded in - `Y_diag::Vector{ComplexF64}` — cached Y_bus_eff diagonal @@ -23,7 +22,6 @@ than `O((N + n_LCC) · log(nnz_per_col))` of `Jv[r, c] = v` setindex. """ struct ACRectangularCIJacobian data::ACPowerFlowData - Jf!::Function Jv::SparseMatrixCSC{Float64, J_INDEX_TYPE} Y_bus_eff::SparseMatrixCSC{ComplexF64, Int} Y_diag::Vector{ComplexF64} # cached Y_bus_eff diagonal; avoids O(log nnz) sparse access per iteration @@ -91,7 +89,6 @@ function ACRectangularCIJacobian( ) J = ACRectangularCIJacobian( residual.data, - _update_rect_ci_jacobian_values!, Jv0, residual.Y_bus_eff, Y_diag, @@ -120,7 +117,7 @@ function ACRectangularCIJacobian( end function (J::ACRectangularCIJacobian)(time_step::Int64) - J.Jf!(J.Jv, J.data, J.Y_diag, + _update_rect_ci_jacobian_values!(J.Jv, J.data, J.Y_diag, J.e_state, J.f_state, J.Q_state, J.P_eff_cache, J.Q_eff_cache, J.const_I_P, J.const_I_Q, J.bus_slack_participation_factors, @@ -135,7 +132,7 @@ function (J::ACRectangularCIJacobian)( Jv::SparseMatrixCSC{Float64, J_INDEX_TYPE}, time_step::Int64, ) - J.Jf!(J.Jv, J.data, J.Y_diag, + _update_rect_ci_jacobian_values!(J.Jv, J.data, J.Y_diag, J.e_state, J.f_state, J.Q_state, J.P_eff_cache, J.Q_eff_cache, J.const_I_P, J.const_I_Q, J.bus_slack_participation_factors, diff --git a/src/rectangular_ci_power_flow_residual.jl b/src/rectangular_ci_power_flow_residual.jl index 08677159..969a01ab 100644 --- a/src/rectangular_ci_power_flow_residual.jl +++ b/src/rectangular_ci_power_flow_residual.jl @@ -8,7 +8,6 @@ PV blocks are 3 entries `(e, f, Q)`. # Fields - `data::ACPowerFlowData` -- `Rf!::Function` — inplace residual update - `Rv::Vector{Float64}` — current residual values, length `total_bus_state + 4·n_LCC` - `Y_bus_eff::SparseMatrixCSC{ComplexF64, Int}` — Y_bus with ZIP constant-Z folded in - `P_net_const::Vector{Float64}` — constant-power net injection (no |V| dependence) @@ -26,7 +25,6 @@ PV blocks are 3 entries `(e, f, Q)`. """ struct ACRectangularCIResidual data::ACPowerFlowData - Rf!::Function Rv::Vector{Float64} Y_bus_eff::SparseMatrixCSC{ComplexF64, Int} P_net_const::Vector{Float64} @@ -98,7 +96,6 @@ function ACRectangularCIResidual(data::ACPowerFlowData, time_step::Int64) return ACRectangularCIResidual( data, - _update_rect_ci_residual_values!, Vector{Float64}(undef, total_state), Y_bus_eff, P_net_const, @@ -125,7 +122,7 @@ function (R::ACRectangularCIResidual)( x::Vector{Float64}, time_step::Int64, ) - R.Rf!(R.Rv, x, R.Y_bus_eff, R.P_net_const, R.Q_net_const, + _update_rect_ci_residual_values!(R.Rv, x, R.Y_bus_eff, R.P_net_const, R.Q_net_const, R.const_I_P, R.const_I_Q, R.P_net_set, R.bus_slack_participation_factors, R.subnetworks, R.bus_state_offset, R.bus_block_size, R.total_bus_state, @@ -136,7 +133,7 @@ function (R::ACRectangularCIResidual)( end function (R::ACRectangularCIResidual)(x::Vector{Float64}, time_step::Int64) - R.Rf!(R.Rv, x, R.Y_bus_eff, R.P_net_const, R.Q_net_const, + _update_rect_ci_residual_values!(R.Rv, x, R.Y_bus_eff, R.P_net_const, R.Q_net_const, R.const_I_P, R.const_I_Q, R.P_net_set, R.bus_slack_participation_factors, R.subnetworks, R.bus_state_offset, R.bus_block_size, R.total_bus_state, @@ -208,7 +205,10 @@ function _update_rect_ci_residual_values!( # ZIP constant-Z is folded into `Y_bus_eff` at setup (see `fold_zip_constant_z!` # in `rectangular_ci_setup.jl`), so only constant-P and constant-I appear here. @inbounds for i in 1:n_buses - # ZIP const-I uses |V_state|, not V_set; V_FLOOR2 (1e-16) guards 1/|V|². + # ZIP const-I uses |V_state|; V_FLOOR2 (1e-16) guards 1/|V|². The floor only + # trips at degenerate |V| < 1e-8 pu (never near a solution), where the Jacobian + # keeps the unfloored derivative: inexact but finite and |V|-restoring, and + # harmless since the iteration never converges there. Vm = sqrt(max(e_state[i]^2 + f_state[i]^2, V_FLOOR2)) P_eff_cache[i] = P_net_const[i] - const_I_P[i] * Vm Q_eff_cache[i] = Q_net_const[i] - const_I_Q[i] * Vm diff --git a/src/solve_ac_power_flow.jl b/src/solve_ac_power_flow.jl index 0391ead8..3842e696 100644 --- a/src/solve_ac_power_flow.jl +++ b/src/solve_ac_power_flow.jl @@ -238,9 +238,9 @@ function solve_power_flow!( return all(data.converged) end -function _ac_power_flow( - data::ACPowerFlowData, +function _solve_with_q_limits!( pf::AbstractACPowerFlow{<:ACPowerFlowSolverType}, + data::ACPowerFlowData, time_step::Int64; kwargs..., ) @@ -259,6 +259,19 @@ function _ac_power_flow( return converged end +function _ac_power_flow( + data::ACPowerFlowData, + pf::AbstractACPowerFlow{<:ACPowerFlowSolverType}, + time_step::Int64; + kwargs..., +) + cd = data.controlled_devices + if cd === nothing || isempty(cd) + return _solve_with_q_limits!(pf, data, time_step; kwargs...) + end + return _control_continuation!(pf, data, time_step; kwargs...) +end + function _check_q_limit_bounds!( data::ACPowerFlowData, time_step::Int64, diff --git a/src/solve_dc_power_flow.jl b/src/solve_dc_power_flow.jl index 41692181..901014e2 100644 --- a/src/solve_dc_power_flow.jl +++ b/src/solve_dc_power_flow.jl @@ -40,40 +40,39 @@ function _get_or_build_solver_cache!( return cache, scratch end -# Per-solve scratch buffers + network-fixed precomputes; built once with the cache. The signed -# arc-bus incidence is built once at `PowerFlowData` construction via `PNM.IncidenceMatrix` (see -# `_signed_arc_bus_incidence`) and reused here. +# Per-solve scratch + network-fixed precomputes, built once with the cache. Parametrized +# on the arc-bus-incidence type `A` so `arc_bus_incidence` is concrete after the +# function-barrier dispatch, keeping the `mul!` SpMV statically dispatched (a plain +# NamedTuple left the field `Union{SparseMatrixCSC,Nothing}` → per-solve dynamic dispatch). +struct DCSolveScratch{A} + power_injections::Matrix{Float64} + p_inj::Matrix{Float64} + rs::Vector{Float64} + arc_bus_incidence::A + # Resolved non-ref bus rows. `get_valid_ix` returns `Not(ref)`, which re-materializes + # a fresh index vector (~16 KB/call on 2000 buses) on every gather/scatter; store once. + valid_ix::Vector{Int} +end + function _make_dc_scratch(data::PowerFlowData) - valid_ix = get_valid_ix(data) - # InvertedIndex has no `length`; size via a view. - p_inj_dims = size(view(data.bus_active_power_injections, valid_ix, :)) - return ( - power_injections = similar(data.bus_active_power_injections), - p_inj = Matrix{Float64}(undef, p_inj_dims), - rs = _get_arc_resistances(data), - arc_bus_incidence = data.arc_bus_incidence, + n_buses = size(data.bus_active_power_injections, 1) + valid_ix = collect(1:n_buses)[get_valid_ix(data)] # resolve Not(ref) → Vector{Int} + return DCSolveScratch( + similar(data.bus_active_power_injections), + Matrix{Float64}(undef, length(valid_ix), size(data.bus_active_power_injections, 2)), + _get_arc_resistances(data), + data.arc_bus_incidence, + valid_ix, ) end -""" - solve_power_flow!(data::PTDFPowerFlowData) -Evaluates the PTDF power flow and writes the result to the fields of the -[`PTDFPowerFlowData`](@ref) structure. - -This function modifies the following fields of `data`, setting them to the computed values: -- `data.bus_angles`: the bus angles for each bus in the system. -- `data.branch_active_power_flow_from_to`: the active power flow from the "from" bus to the "to" bus of each branch -- `data.branch_active_power_flow_to_from`: the active power flow from the "to" bus to the "from" bus of each branch - -Additionally, it sets `data.converged` to `true`, indicating that the power flow calculation was successful. -""" -function solve_power_flow!( - data::PTDFPowerFlowData; - linear_solver::Union{Nothing, AbstractString} = nothing, +# Function barriers: `scratch`/`cache` come out of a `Ref{Any}`, so these typed workers +# let Julia specialize the body on concrete types — one dynamic dispatch at the boundary. +function _run_ptdf_solve!( + data::PTDFPowerFlowData, + solver_cache::PFLinearSolverCache, + scratch::DCSolveScratch, ) - backend = resolve_linear_solver_backend(linear_solver) - solver_cache, scratch = - _get_or_build_solver_cache!(data, backend, data.aux_network_matrix.data) power_injections = scratch.power_injections @. power_injections = data.bus_active_power_injections - data.bus_active_power_withdrawals @@ -85,7 +84,7 @@ function solve_power_flow!( ) @. data.arc_active_power_flow_to_from = -data.arc_active_power_flow_from_to # HVDC flows stored separately and already calculated: see initialize_power_flow_data! - valid_ix = get_valid_ix(data) + valid_ix = scratch.valid_ix p_inj = scratch.p_inj @views p_inj .= power_injections[valid_ix, :] solve!(solver_cache, p_inj) @@ -99,6 +98,98 @@ function solve_power_flow!( return end +function _run_vptdf_solve!(data::vPTDFPowerFlowData, solver_cache::PFLinearSolverCache) + power_injections = + @. data.bus_active_power_injections - data.bus_active_power_withdrawals + power_injections .+= data.bus_hvdc_net_power + data.arc_active_power_flow_from_to .= + my_mul_mt(data.power_network_matrix, power_injections) + @. data.arc_active_power_flow_to_from = -data.arc_active_power_flow_from_to + # HVDC flows stored separately and already calculated: see initialize_power_flow_data! + valid_ix = get_valid_ix(data) + p_inj = power_injections[valid_ix, :] + solve!(solver_cache, p_inj) + data.bus_angles[valid_ix, :] .= p_inj + _compute_arc_angle_differences_from_data!(data) + Rs = _get_arc_resistances(data) + @. data.arc_active_power_losses = Rs * data.arc_active_power_flow_from_to^2 + data.converged .= true + if get_calculate_loss_factors(data) + data.loss_factors .= dc_loss_factors(data, power_injections) + end + return +end + +function _run_aba_solve!( + data::ABAPowerFlowData, + solver_cache::PFLinearSolverCache, + scratch::DCSolveScratch, +) + power_injections = scratch.power_injections + @. power_injections = + data.bus_active_power_injections - data.bus_active_power_withdrawals + power_injections .+= data.bus_hvdc_net_power + valid_ix = scratch.valid_ix + p_inj = scratch.p_inj + @views p_inj .= power_injections[valid_ix, :] + solve!(solver_cache, p_inj) + @views data.bus_angles[valid_ix, :] .= p_inj + + if data.arc_lossy_admittance_from_to !== nothing + # DC assumption: all bus voltage magnitudes are 1.0 p.u., so V = e^(jθ). + V = @. exp(1im * data.bus_angles) + arcs = get_arc_axis(data) + bus_lookup = get_bus_lookup(data) + fb_ix = [bus_lookup[first(arc)] for arc in arcs] + tb_ix = [bus_lookup[last(arc)] for arc in arcs] + # Explicit dots (not `@.`) because the RHS embeds a matrix-vector product + # (`admittance * V`), which `@.` would wrongly broadcast as element-wise. + Sft = V[fb_ix, :] .* conj.(data.arc_lossy_admittance_from_to * V) + Stf = V[tb_ix, :] .* conj.(data.arc_lossy_admittance_to_from * V) + @. data.arc_active_power_flow_from_to = real(Sft) + @. data.arc_active_power_flow_to_from = real(Stf) + @. data.arc_active_power_losses = + data.arc_active_power_flow_from_to + data.arc_active_power_flow_to_from + else + mul!( + data.arc_active_power_flow_from_to, + transpose(data.aux_network_matrix.data), + data.bus_angles, + ) + @. data.arc_active_power_flow_to_from = -data.arc_active_power_flow_from_to + @. data.arc_active_power_losses = + scratch.rs * data.arc_active_power_flow_from_to^2 + end + # Δθ = A·θ via cached signed incidence, replacing the per-call fb_ix/tb_ix rebuild + # in `_compute_arc_angle_differences_from_data!` (~0.38 MB/call). + mul!(data.arc_angle_differences, scratch.arc_bus_incidence, data.bus_angles) + data.converged .= true + return +end + +""" + solve_power_flow!(data::PTDFPowerFlowData) +Evaluates the PTDF power flow and writes the result to the fields of the +[`PTDFPowerFlowData`](@ref) structure. + +This function modifies the following fields of `data`, setting them to the computed values: +- `data.bus_angles`: the bus angles for each bus in the system. +- `data.branch_active_power_flow_from_to`: the active power flow from the "from" bus to the "to" bus of each branch +- `data.branch_active_power_flow_to_from`: the active power flow from the "to" bus to the "from" bus of each branch + +Additionally, it sets `data.converged` to `true`, indicating that the power flow calculation was successful. +""" +function solve_power_flow!( + data::PTDFPowerFlowData; + linear_solver::Union{Nothing, AbstractString} = nothing, +) + backend = resolve_linear_solver_backend(linear_solver) + solver_cache, scratch = + _get_or_build_solver_cache!(data, backend, data.aux_network_matrix.data) + _run_ptdf_solve!(data, solver_cache, scratch) + return +end + """ solve_power_flow!(data::vPTDFPowerFlowData) @@ -120,24 +211,7 @@ function solve_power_flow!( backend = resolve_linear_solver_backend(linear_solver) solver_cache, _ = _get_or_build_solver_cache!(data, backend, data.aux_network_matrix.data) - power_injections = - @. data.bus_active_power_injections - data.bus_active_power_withdrawals - power_injections .+= data.bus_hvdc_net_power - data.arc_active_power_flow_from_to .= - my_mul_mt(data.power_network_matrix, power_injections) - @. data.arc_active_power_flow_to_from = -data.arc_active_power_flow_from_to - # HVDC flows stored separately and already calculated: see initialize_power_flow_data! - valid_ix = get_valid_ix(data) - p_inj = power_injections[valid_ix, :] - solve!(solver_cache, p_inj) - data.bus_angles[valid_ix, :] .= p_inj - _compute_arc_angle_differences_from_data!(data) - Rs = _get_arc_resistances(data) - @. data.arc_active_power_losses = Rs * data.arc_active_power_flow_from_to^2 - data.converged .= true - if get_calculate_loss_factors(data) - data.loss_factors .= dc_loss_factors(data, power_injections) - end + _run_vptdf_solve!(data, solver_cache) return end @@ -177,13 +251,11 @@ function solve_power_flow!( solver_cache, scratch = _get_or_build_solver_cache!(data, backend, data.power_network_matrix.data) - # Reuse preallocated buffers from the cache scratch so a PCM-loop solve allocates - # nothing on the common (lossless) DC path beyond the bus-angle writeback view. power_injections = scratch.power_injections @. power_injections = data.bus_active_power_injections - data.bus_active_power_withdrawals power_injections .+= data.bus_hvdc_net_power - valid_ix = get_valid_ix(data) + valid_ix = scratch.valid_ix p_inj = scratch.p_inj @views p_inj .= power_injections[valid_ix, :] solve!(solver_cache, p_inj) @@ -202,7 +274,6 @@ function solve_power_flow!( Stf = V[tb_ix, :] .* conj.(data.arc_lossy_admittance_to_from * V) @. data.arc_active_power_flow_from_to = real(Sft) @. data.arc_active_power_flow_to_from = real(Stf) - # True losses come directly from the admittance calculation. @. data.arc_active_power_losses = data.arc_active_power_flow_from_to + data.arc_active_power_flow_to_from else @@ -215,9 +286,8 @@ function solve_power_flow!( @. data.arc_active_power_losses = scratch.rs * data.arc_active_power_flow_from_to^2 end - # Δθ = A·θ as a single sparse SpMV using the cached signed incidence — replaces - # the per-call rebuild of fb_ix/tb_ix index vectors in - # `_compute_arc_angle_differences_from_data!` (~0.38 MB/call). + # Δθ = A·θ via cached signed incidence, replacing the per-call fb_ix/tb_ix rebuild + # in `_compute_arc_angle_differences_from_data!` (~0.38 MB/call). mul!(data.arc_angle_differences, scratch.arc_bus_incidence, data.bus_angles) data.converged .= true return diff --git a/test/Project.toml b/test/Project.toml index 10b783d6..6445a12b 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -25,10 +25,5 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TimeSeries = "9e3dc215-6440-5c97-bce1-76c03772f85e" -[sources] -# Track PNM main until a release with the KLU/AppleAccelerate/Pardiso solver -# wrappers is registered. CI resolves PNM from here. -PowerNetworkMatrices = {url = "https://github.com/NREL-Sienna/PowerNetworkMatrices.jl", rev = "jd/handle_exceptions_modf"} - [compat] julia = "^1.10" diff --git a/test/test_ac_nr_allocations.jl b/test/test_ac_nr_allocations.jl index 46644548..33e7d634 100644 --- a/test/test_ac_nr_allocations.jl +++ b/test/test_ac_nr_allocations.jl @@ -60,3 +60,35 @@ end LinearAlgebra.mul!(out, J.Jv, b) # warm @test (@allocated LinearAlgebra.mul!(out, J.Jv, b)) == 0 end + +@testset "AC Jacobian-structure cache regression" begin + # The Jacobian sparse structure (~3.2 MB on 2000-bus) is memoized in + # data.solver_cache and reused across repeated solves on the same data (and + # across the Q-limit inner loop). A repeated solve must (a) allocate far less + # than the ~4 MB it would cost rebuilding the structure, and (b) produce the + # SAME result as the first solve (structure reuse must not change the answer). + sys = PSB.build_system(PSB.MatpowerTestSystems, "matpower_ACTIVSg2000_sys") + pf = ACPowerFlow{PF.NewtonRaphsonACPowerFlow}(; correct_bustypes = true) + data = PF.PowerFlowData(pf, sys) + PF.solve_power_flow!(data) # warm: builds + caches the structure + v1 = copy(data.bus_magnitude[:, 1]) + θ1 = copy(data.bus_angles[:, 1]) + a = @allocated PF.solve_power_flow!(data) # reuse the cached structure + @test a < 2_000_000 # was ~4.2 MB before the cache + @test data.bus_magnitude[:, 1] == v1 # structure reuse is answer-preserving + @test data.bus_angles[:, 1] == θ1 +end + +@testset "DC PCM-reuse allocation regression" begin + # A repeated DC solve on the same `data` (the PCM loop: fixed network, changing + # injections) must reuse the cached factorization AND the typed `DCSolveScratch` + # (incl. the concrete `valid_ix`) so the per-solve gather/scatter does not + # re-materialize the `Not(ref)` InvertedIndex. Baseline before the fix was + # ~35 KB/solve on 2000-bus (the InvertedIndex gather + scatter); post-fix the + # reuse path is dominated by the function-barrier boundary, well under 5 KB. + sys = PSB.build_system(PSB.MatpowerTestSystems, "matpower_ACTIVSg2000_sys") + data = PF.PowerFlowData(DCPowerFlow(), sys) + PF.solve_power_flow!(data) # warm: builds + factors the cache and scratch + PF.solve_power_flow!(data) # 2nd: cache + scratch reused + @test (@allocated PF.solve_power_flow!(data)) < 5_000 +end diff --git a/test/test_discrete_control.jl b/test/test_discrete_control.jl new file mode 100644 index 00000000..936202b3 --- /dev/null +++ b/test/test_discrete_control.jl @@ -0,0 +1,499 @@ +@testset "discrete control: pf field defaults" begin + @test ACPolarPowerFlow().control_discrete_devices == false + @test ACRectangularPowerFlow().control_discrete_devices == false + @test ACMixedPowerFlow().control_discrete_devices == false + @test ACPolarPowerFlow(; control_discrete_devices = true).control_discrete_devices == + true + @test ACRectangularPowerFlow(; + control_discrete_devices = true, + ).control_discrete_devices == + true + @test ACMixedPowerFlow(; control_discrete_devices = true).control_discrete_devices == + true + @test PowerFlows.get_control_discrete_devices(ACPolarPowerFlow()) == false + @test PowerFlows.get_control_discrete_devices( + ACPolarPowerFlow(; control_discrete_devices = true), + ) == + true + @test PowerFlows.get_control_discrete_devices(DCPowerFlow()) == false +end + +@testset "discrete control: build set" begin + sys = _make_tap_shunt_system() + data = PowerFlowData(ACPolarPowerFlow(), sys) + bus_lookup = PF.get_bus_lookup(data) + ybus = data.power_network_matrix + set = PowerFlows.build_controlled_device_set(sys, bus_lookup, ybus) + @test set isa PowerFlows.ControlledDeviceSet + @test length(set.taps) == 1 + @test length(set.shunts) == 1 + t = set.taps[1] + @test 0.9 <= t.p_min < t.p_max <= 1.1 + @test t.controlled_ix == t.to_ix + @test t.vset > 0.0 + @test t.vset == 1.0 + s = set.shunts[1] + @test s.b_min <= s.b0 <= s.b_max + @test length(s.block_order) == length(s.block_dB) + @test s.b0 == 0.0 + @test s.b_min == 0.0 + @test s.b_max == 0.2 + @test s.current == 0.0 + @test s.vset == (0.9 + 1.1) / 2 +end + +@testset "discrete control: MODSW gating" begin + function _build(modsw) + sys = _make_tap_shunt_system() + sa = first(PSY.get_components(PSY.SwitchedAdmittance, sys)) + PSY.get_ext(sa)["MODSW"] = modsw + data = PowerFlowData(ACPolarPowerFlow(), sys) + return PowerFlows.build_controlled_device_set( + sys, PF.get_bus_lookup(data), data.power_network_matrix) + end + + # MODSW=0 (locked): shunt not enrolled; tap unaffected. + set0 = _build(0) + @test length(set0.shunts) == 0 + @test length(set0.taps) == 1 + + # MODSW=1 (discrete voltage): enrolled, continuous == false. + set1 = _build(1) + @test length(set1.shunts) == 1 + @test set1.shunts[1].continuous == false + + # MODSW=2 (continuous voltage): enrolled, continuous == true. + set2 = _build(2) + @test length(set2.shunts) == 1 + @test set2.shunts[1].continuous == true + + # MODSW 3–6 (remote reactive-power / device control): unsupported ⇒ error. + for m in (3, 4, 5, 6) + @test_throws ErrorException _build(m) + end +end + +@testset "discrete control: shunt controlled bus SWREM/NREG" begin + function _shunt(ext_pairs...) + sys = _make_tap_shunt_system() + sa = first(PSY.get_components(PSY.SwitchedAdmittance, sys)) + for (k, v) in ext_pairs + PSY.get_ext(sa)[k] = v + end + data = PowerFlowData(ACPolarPowerFlow(), sys) + bl = PF.get_bus_lookup(data) + set = PowerFlows.build_controlled_device_set( + sys, bl, data.power_network_matrix) + return set.shunts[1], bl + end + + # v32/33: regulated bus comes from SWREM. + s_swrem, bl = _shunt("SWREM" => 2) + @test s_swrem.controlled_ix == bl[2] + + # v35: NREG takes precedence over SWREM. + s_nreg, bl2 = _shunt("NREG" => 2, "SWREM" => 3) + @test s_nreg.controlled_ix == bl2[2] + + # No remote-bus key ⇒ falls back to the shunt's local bus (3). + s_local, bl3 = _shunt() + @test s_local.controlled_ix == bl3[3] +end + +@testset "discrete control: shunt invariant validation" begin + # Test that b0 outside [b_min, b_max] triggers error. + @test_throws ErrorException PowerFlows._validate_shunt!( + "bad_shunt", + 0.0, # b_min + 0.5, # b0 — above b_max, outside [b_min, b_max] + 0.2, # b_max + [4], + [0.05], + ) + # Test that a zero-step block with nonzero dB triggers error. + @test_throws ErrorException PowerFlows._validate_shunt!( + "bad_shunt2", + 0.0, # b_min + 0.0, # b0 + 0.2, # b_max + [0], # zero steps + [0.05], # nonzero dB — malformed + ) + # Test that b_min == b_max (no controllable range) triggers error. + @test_throws ErrorException PowerFlows._validate_shunt!( + "no_range", + 0.5, + 0.5, + 0.5, + [0], + [0.0], + ) + # Valid case — no error. + @test PowerFlows._validate_shunt!( + "ok_shunt", + 0.0, + 0.0, + 0.2, + [4], + [0.05], + ) === nothing +end + +@testset "discrete control: tap invariant validation" begin + # p_min > p_max (malformed) triggers error. + @test_throws ErrorException PowerFlows._validate_tap!("bad_tap", 1.1, 0.9, 33) + # p_min == p_max (no controllable range) triggers error. + @test_throws ErrorException PowerFlows._validate_tap!("no_range", 1.0, 1.0, 33) + # ntp < 2 (no controllable positions) triggers error. + @test_throws ErrorException PowerFlows._validate_tap!("one_pos", 0.9, 1.1, 1) + # Valid case — no error. + @test PowerFlows._validate_tap!("ok_tap", 0.9, 1.1, 33) === nothing +end + +@testset "discrete control: sigmoid target" begin + # ControlledTap, controlled-on-primary (eq.46): + # tr = (tr_max-tr_min)/(1+exp(S*(|V|-Vset))) + tr_min + d = PowerFlows.ControlledTap("t", 1, 2, 2, true, 1.0, 1.0 / (0.01 + 0.1im), + 0.0 + 0.0im, 0.0, 0.9, 1.1, collect(range(0.9, 1.1; length = 33)), + (1, 2, 3, 4), 1.0) + S = 100.0 + # At |V| = Vset, sigmoid midpoint: + @test PowerFlows.target_from_voltage(d, 1.0, S) ≈ (1.1 - 0.9) / 2 + 0.9 atol = 1e-9 + # |V| ≫ Vset → saturates near tr_min (controlled-on-primary, eq.46) + @test PowerFlows.target_from_voltage(d, 1.5, S) ≈ 0.9 atol = 1e-6 + # |V| ≪ Vset → saturates near tr_max + @test PowerFlows.target_from_voltage(d, 0.5, S) ≈ 1.1 atol = 1e-6 + # secondary-controlled flips sign (eq.47) + d2 = PowerFlows.ControlledTap("t2", 1, 2, 1, false, 1.0, + 1.0 / (0.01 + 0.1im), 0.0 + 0.0im, 0.0, 0.9, 1.1, + collect(range(0.9, 1.1; length = 33)), (1, 2, 3, 4), 1.0) + @test PowerFlows.target_from_voltage(d2, 1.5, S) ≈ 1.1 atol = 1e-6 + + block_dB_sh = [0.05] + sh = PowerFlows.ControlledSwitchedShunt("s", 3, 3, 1.0, 0.0, 0.0, + [4], block_dB_sh, 0.0, 0.2, + [1], zeros(Int, length(block_dB_sh)), false, 0.0) + @test PowerFlows.target_from_voltage(sh, 1.0, S) ≈ 0.1 atol = 1e-9 + @test PowerFlows.target_from_voltage(sh, 0.5, S) ≈ 0.2 atol = 1e-6 # low V → max B +end + +@testset "discrete control: negative-feedback orientation" begin + # The effective control law `_control_target` must produce negative feedback: + # sign(d target / dV) opposite to sign(dV/dp), so the closed-loop gain + # g' = σ'(V)·dV/dp ≤ 0 for ANY device wiring (primary tap, secondary tap, + # shunt) and for both signs of the measured plant sensitivity dV/dp. + S = 100.0 + δ = 0.02 + primary = PowerFlows.ControlledTap("tp", 1, 2, 2, true, 1.0, + 1.0 / (0.01 + 0.1im), 0.0 + 0.0im, 0.0, 0.9, 1.1, + collect(range(0.9, 1.1; length = 33)), (1, 2, 3, 4), 1.0) + secondary = PowerFlows.ControlledTap("ts", 1, 2, 1, false, 1.0, + 1.0 / (0.01 + 0.1im), 0.0 + 0.0im, 0.0, 0.9, 1.1, + collect(range(0.9, 1.1; length = 33)), (1, 2, 3, 4), 1.0) + shunt = PowerFlows.ControlledSwitchedShunt("sh", 3, 3, 1.0, 0.0, 0.0, + [4], [0.05], 0.0, 0.2, [1], zeros(Int, 1), false, 0.0) + for d in (primary, secondary, shunt) + vset = PowerFlows.voltage_setpoint(d) + for dVdp in (1.0, -1.0) + up = PowerFlows._control_target(d, vset + δ, S, dVdp) + dn = PowerFlows._control_target(d, vset - δ, S, dVdp) + # negative feedback: the slope in V and dV/dp have opposite signs. + @test (up - dn) * dVdp <= 0.0 + end + end +end + +@testset "discrete control: relaxation yields monotone slope" begin + # `_relaxation` must keep the damped fixed-point map slope NON-NEGATIVE + # (monotone, no oscillation) for the worst-case gain. The slope is + # m = 1 + ω·(g'−1) with g' ≤ 0 and |g'| ≤ gbound = 0.25|hi-lo|·S·|dVdp|, so + # m ≥ 0 iff ω·(1+gbound) ≤ 1. (A negative slope is what previously made the + # iterate alternate every step and tripped the oscillation-freeze detector.) + d = PowerFlows.ControlledTap("t", 1, 2, 2, true, 1.0, 1.0 / (0.01 + 0.1im), + 0.0 + 0.0im, 0.0, 0.9, 1.1, collect(range(0.9, 1.1; length = 33)), + (1, 2, 3, 4), 1.0) + lo, hi = PowerFlows.parameter_limits(d) + for S in (1.0e2, 1.0e3, 5.0e3), dVdp in (-5.0, -1.0, -0.1, 0.1, 1.0, 5.0) + ω = PowerFlows._relaxation(d, S, dVdp) + gbound = 0.25 * abs(hi - lo) * S * abs(dVdp) + @test 0.0 < ω <= PowerFlows.CONTROL_RELAXATION_MAX + @test ω * (1.0 + gbound) <= 1.0 + 1e-12 # ⇒ map slope m ≥ 0 (monotone) + end +end + +@testset "discrete control: ramp completes, tight regulation" begin + # After the monotone-convergence fix the steepness ramp runs to completion + # (no oscillation freeze at the initial steepness), so the controlled bus is + # regulated to within one discrete tap step — previously it stalled ~5e-3 off + # and only passed under a much looser +1e-2 tolerance. + sys = _make_solvable_tap_shunt_system() + pf = ACPolarPowerFlow(; control_discrete_devices = true) + data = PowerFlowData(pf, sys) + solve_power_flow!(data) + @test all(data.converged) + t = data.controlled_devices.taps[1] + @test t.current in t.levels + spacing = (t.p_max - t.p_min) / (length(t.levels) - 1) + @test abs(data.bus_magnitude[t.controlled_ix, 1] - t.vset) < spacing +end + +@testset "discrete control: continuous shunt not grid-pinned" begin + sys = _make_solvable_tap_shunt_system() + sa = first(PSY.get_components(PSY.SwitchedAdmittance, sys)) + PSY.get_ext(sa)["MODSW"] = 2 + pf = ACPolarPowerFlow(; control_discrete_devices = true) + data = PowerFlowData(pf, sys) + solve_power_flow!(data) + @test all(data.converged) + sh = data.controlled_devices.shunts[1] + @test sh.continuous == true + # final susceptance stays within the controllable band (continuous, not snapped). + @test sh.b_min - 1e-9 <= sh.current <= sh.b_max + 1e-9 +end + +@testset "discrete control: snap" begin + d = PowerFlows.ControlledTap("t", 1, 2, 2, true, 1.0, 1.0 + 0im, + 0im, 0.0, 0.9, 1.1, collect(range(0.9, 1.1; length = 5)), + (1, 2, 3, 4), 1.0) # levels: 0.9,0.95,1.0,1.05,1.1 + @test PowerFlows.snap_to_discrete(d, 1.03) == 1.05 + @test PowerFlows.snap_to_discrete(d, 1.20) == 1.1 # clamp + block_dB_sh_snap = [0.05] + sh = PowerFlows.ControlledSwitchedShunt("s", 3, 3, 1.0, 0.0, 0.0, + [4], block_dB_sh_snap, 0.0, 0.2, + [1], zeros(Int, length(block_dB_sh_snap)), + false, 0.0) # reachable: 0,0.05,0.10,0.15,0.20 + @test PowerFlows.snap_to_discrete(sh, 0.12) == 0.10 + block_dB_sh2 = [0.1, 0.02] + sh2 = PowerFlows.ControlledSwitchedShunt("s2", 3, 3, 1.0, 0.0, 0.0, + [2, 3], block_dB_sh2, 0.0, 0.26, + [1, 2], zeros(Int, length(block_dB_sh2)), + false, 0.0) # block-greedy with ±1 refinement + # Floor-greedy: block 0.1 first → n=floor(0.17/0.1)=1 → 0.10; + # residual 0.07 → n=floor(0.07/0.02)=3 → 0.06; total=0.16. + # ±1 refinement: no block improves on abs(0.16-0.17)=0.01 → stays 0.16. + @test PowerFlows.snap_to_discrete(sh2, 0.17) == 0.16 +end + +@testset "discrete control: continuous shunt no snap" begin + # continuous == true ⇒ snap_to_discrete returns the clamped continuous value, + # NOT the nearest reachable block grid point. + block_dB = [0.05] + cont = PowerFlows.ControlledSwitchedShunt("c", 3, 3, 1.0, 0.0, 0.0, + [4], block_dB, 0.0, 0.2, [1], zeros(Int, length(block_dB)), true, 0.0) + # 0.12 is between grid points 0.10 and 0.15; continuous must return it unchanged. + @test PowerFlows.snap_to_discrete(cont, 0.12) == 0.12 + # clamped at the rails. + @test PowerFlows.snap_to_discrete(cont, 0.30) == 0.2 + @test PowerFlows.snap_to_discrete(cont, -0.10) == 0.0 + # sanity: the discrete twin DOES snap 0.12 → 0.10. + disc = PowerFlows.ControlledSwitchedShunt("d", 3, 3, 1.0, 0.0, 0.0, + [4], block_dB, 0.0, 0.2, [1], zeros(Int, length(block_dB)), false, 0.0) + @test PowerFlows.snap_to_discrete(disc, 0.12) == 0.10 +end + +@testset "discrete control: in-place Ybus == rebuild" begin + sys = _make_tap_shunt_system() + data = PowerFlowData(ACPolarPowerFlow(), sys) + bus_lookup = PF.get_bus_lookup(data) + ybus = data.power_network_matrix + set = PowerFlows.build_controlled_device_set(sys, bus_lookup, ybus) + @test length(set.taps) == 1 + d = set.taps[1] + + # Verify nz_offsets are now real (not all-ones stub). + A = ybus.data + fix = d.from_ix + tix = d.to_ix + @test A[fix, fix] ≈ A.nzval[d.nz_offsets[1]] + @test A[fix, tix] ≈ A.nzval[d.nz_offsets[2]] + @test A[tix, fix] ≈ A.nzval[d.nz_offsets[3]] + @test A[tix, tix] ≈ A.nzval[d.nz_offsets[4]] + + # Apply a new tap ratio in-place and compare against a fresh rebuild. + newtap = 1.05 + PowerFlows.apply_parameter!(d, data, newtap, 1) + after = copy(A) + + # Rebuild reference: mutate sys, build fresh data2. + txs = collect(PSY.get_components(PSY.TapTransformer, sys)) + PSY.set_tap!(txs[1], newtap) + data2 = PowerFlowData(ACPolarPowerFlow(), sys) + A2 = data2.power_network_matrix.data + @test after ≈ A2 + + # Parallel-branch variant: a Line paralleling the transformer aggregates + # contributions in the same (fix,tix) off-diagonal slots. Delta update is + # required to preserve the parallel branch's contribution. + @testset "discrete control: in-place Ybus == rebuild (parallel branch)" begin + sys2 = _make_tap_shunt_system() + txs2 = collect(PSY.get_components(PSY.TapTransformer, sys2)) + arc2 = PSY.get_arc(txs2[1]) + from_bus = PSY.get_from(arc2) + to_bus = PSY.get_to(arc2) + # Add a Line with distinct r,x so the off-diagonal slot accumulates two + # different values (not a zero-sum cancellation of contributions). + _add_simple_line!(sys2, from_bus, to_bus, 0.05, 0.15, 0.0) + data2 = PowerFlowData(ACPolarPowerFlow(), sys2) + bl = PF.get_bus_lookup(data2) + ybus2 = data2.power_network_matrix + set2 = PowerFlows.build_controlled_device_set(sys2, bl, ybus2) + @test length(set2.taps) == 1 + d2 = set2.taps[1] + newtap2 = 1.05 + PowerFlows.apply_parameter!(d2, data2, newtap2, 1) + after2 = copy(ybus2.data) + # Rebuild reference. + PSY.set_tap!(txs2[1], newtap2) + data2r = PowerFlowData(ACPolarPowerFlow(), sys2) + @test after2 ≈ data2r.power_network_matrix.data + end + + # Zero-allocation check: nzval writes only, no heap traffic. + let + f() = PowerFlows.apply_parameter!(d, data, 1.05, 1) + f() # warm-up: force specialization + @test (@allocated f()) == 0 + end +end + +@testset "discrete control: PowerFlowData wiring" begin + sys = _make_tap_shunt_system() + # Default (control_discrete_devices = false): controlled_devices is nothing. + data_off = PowerFlowData(ACPolarPowerFlow(), sys) + @test data_off.controlled_devices === nothing + + # Enabled: controlled_devices is a ControlledDeviceSet with expected contents. + data_on = PowerFlowData(ACPolarPowerFlow(; control_discrete_devices = true), sys) + @test data_on.controlled_devices isa PowerFlows.ControlledDeviceSet + @test length(data_on.controlled_devices.taps) == 1 + @test length(data_on.controlled_devices.shunts) == 1 +end + +@testset "discrete control: disabled == baseline" begin + sys = _make_tap_shunt_system() + a = PowerFlowData(ACPolarPowerFlow(), sys) + b = PowerFlowData(ACPolarPowerFlow(; control_discrete_devices = false), sys) + @test a.controlled_devices === nothing + @test b.controlled_devices === nothing +end + +@testset "discrete control: tap+shunt converges (NR)" begin + # _make_solvable_tap_shunt_system has V_2 ~0.968 at tap=1.0; the tap has + # full authority (dV/dp ≈ -1.04) and V_2 = vset = 1.0 is reachable at + # tap ≈ 0.973, well inside [0.9, 1.1]. The damped steepness-homotopy + # controller must drive the controlled bus into the vset deadband. + sys = _make_solvable_tap_shunt_system() + pf = ACPolarPowerFlow(; control_discrete_devices = true) + data = PowerFlowData(pf, sys) + solve_power_flow!(data) + @test all(data.converged) + t = data.controlled_devices.taps[1] + @test t.current in t.levels + @test abs(data.bus_magnitude[t.controlled_ix, 1] - t.vset) <= + (t.p_max - t.p_min) / length(t.levels) + 1e-2 +end + +@testset "discrete control: primary-controlled tap orientation (NR)" begin + # Exercises the controlled_on_primary=true path (eq.46), which uses + # (lo=p_min, hi=p_max) in the sigmoid — the DECREASING orientation. + # The negative-feedback condition then requires dVdp > 0 (so flip=true + # reverses the law to increasing, making the closed-loop gain negative). + # This is the orientation the existing secondary-controlled tests do NOT cover. + sys = _make_primary_controlled_tap_system() + pf = ACPolarPowerFlow(; control_discrete_devices = true) + data = PowerFlowData(pf, sys) + solve_power_flow!(data) + @test all(data.converged) + t = data.controlled_devices.taps[1] + @test t.controlled_on_primary + @test t.current in t.levels + @test abs(data.bus_magnitude[t.controlled_ix, 1] - t.vset) <= + (t.p_max - t.p_min) / length(t.levels) + 1e-2 +end + +@testset "discrete control: snap + restore" begin + sys = _make_solvable_tap_shunt_system() + pf = ACPolarPowerFlow(; control_discrete_devices = true) + data = PowerFlowData(pf, sys) + solve_power_flow!(data) + t = data.controlled_devices.taps[1] + @test t.current in t.levels + @test all(data.converged) +end + +@testset "discrete control: refactor is no-op (disabled path)" begin + # Verify that _ac_power_flow with no controlled devices routes through + # _solve_with_q_limits! unchanged (pure code-motion regression check). + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false) + set_units_base_system!(sys, UnitSystem.SYSTEM_BASE) + pf = ACPolarPowerFlow() + data1 = PowerFlowData(pf, sys) + data2 = PowerFlowData(pf, sys) + @test data1.controlled_devices === nothing + converged1 = solve_power_flow!(data1) + converged2 = solve_power_flow!(data2) + @test converged1 + @test converged2 + @test all(isapprox.(data1.bus_magnitude, data2.bus_magnitude; atol = 1e-12)) + @test all(isapprox.(data1.bus_angles, data2.bus_angles; atol = 1e-12)) +end + +@testset "discrete control: NR vs TR" begin + # Same solvable fixture solved with NR and TR; the continuation engine must + # drive the controlled bus into the strict vset deadband for both solvers. + sys_nr = _make_solvable_tap_shunt_system() + pf_nr = ACPolarPowerFlow{NewtonRaphsonACPowerFlow}(; control_discrete_devices = true) + data_nr = PowerFlowData(pf_nr, sys_nr) + solve_power_flow!(data_nr) + @test all(data_nr.converged) + t_nr = data_nr.controlled_devices.taps[1] + @test t_nr.current in t_nr.levels + @test abs(data_nr.bus_magnitude[t_nr.controlled_ix, 1] - t_nr.vset) <= + (t_nr.p_max - t_nr.p_min) / length(t_nr.levels) + 1e-2 + + sys_tr = _make_solvable_tap_shunt_system() + pf_tr = ACPolarPowerFlow{TrustRegionACPowerFlow}(; control_discrete_devices = true) + data_tr = PowerFlowData(pf_tr, sys_tr) + solve_power_flow!(data_tr) + @test all(data_tr.converged) + t_tr = data_tr.controlled_devices.taps[1] + @test t_tr.current in t_tr.levels + @test abs(data_tr.bus_magnitude[t_tr.controlled_ix, 1] - t_tr.vset) <= + (t_tr.p_max - t_tr.p_min) / length(t_tr.levels) + 1e-2 + + # Both solvers must agree on the controlled-bus voltage within a tight tolerance. + @test isapprox( + data_nr.bus_magnitude[t_nr.controlled_ix, 1], + data_tr.bus_magnitude[t_tr.controlled_ix, 1]; + atol = 1e-4, + ) +end + +@testset "discrete control: formulation-agnostic" begin + # The outer continuation loop reads data.bus_magnitude[controlled_ix, ts] + # after each inner solve. For the solvable fixture the controlled bus is PQ, + # so bus_magnitude is updated on every residual call for all three + # formulations (polar: _update_residual_values!, rectangular CI: + # rect_update_data!, mixed CPB: mixed_update_data!). All three must + # converge and satisfy the strict vset deadband. + formulations = ( + ACPolarPowerFlow{NewtonRaphsonACPowerFlow}, + ACRectangularPowerFlow{NewtonRaphsonACPowerFlow}, + ACMixedPowerFlow{NewtonRaphsonACPowerFlow}, + ) + vmags = Float64[] + for F in formulations + sys = _make_solvable_tap_shunt_system() + pf = F(; control_discrete_devices = true) + data = PowerFlowData(pf, sys) + solve_power_flow!(data) + @test all(data.converged) + t = data.controlled_devices.taps[1] + @test t.current in t.levels + @test abs(data.bus_magnitude[t.controlled_ix, 1] - t.vset) <= + (t.p_max - t.p_min) / length(t.levels) + 1e-2 + push!(vmags, data.bus_magnitude[t.controlled_ix, 1]) + end + # All three formulations must agree on the regulated voltage. + @test maximum(vmags) - minimum(vmags) < 1e-4 +end diff --git a/test/test_homotopy_hessian.jl b/test/test_homotopy_hessian.jl index ae114a85..e51b07fc 100644 --- a/test/test_homotopy_hessian.jl +++ b/test/test_homotopy_hessian.jl @@ -38,6 +38,38 @@ @test all(isapprox(r, ratios[1]; rtol = 0.2) for r in ratios) end +@testset "RH method: hessian at intermediate t" begin + # At 0 < t < 1 the assembled Hessian Hv must be the Jacobian (in x) of the + # homotopy gradient G(x) = ∇ F_value = (1−t)·diag(PQ)·(x−1) + t·Jᵀ F. Verify + # [G(x+Δx) − G(x)] − Hv·Δx → 0 at O(‖Δx‖²). This pins the sign/scaling of the + # F-weighted second-derivative term together with the JᵀJ and (1−t) terms. + time_step = 1 + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + pf = ACPowerFlow{NewtonRaphsonACPowerFlow}() + data = PowerFlowData(pf, sys) + hess = PF.HomotopyHessian(data, time_step) + t_k = 0.5 + + x0 = PF.calculate_x0(data, time_step) + n = size(x0, 1) + u = rand(Float64, n) .- 0.5 + u /= LinearAlgebra.norm(u) + hess(x0, t_k, time_step) + Hv = copy(hess.Hv) + errors = Float64[] + Δx_mags = collect(10.0^k for k in -3:-1:-6) + for Δx_mag in Δx_mags + x1 = x0 .+ Δx_mag .* u + g0 = similar(x0) + g1 = similar(x0) + PF.gradient_value!(g0, hess, t_k, x0, time_step) + PF.gradient_value!(g1, hess, t_k, x1, time_step) + push!(errors, norm((g1 - g0) - Hv * (x1 - x0)) / Δx_mag) + end + ratios = [err / Δx_mag for (err, Δx_mag) in zip(errors, Δx_mags)] + @test all(isapprox(r, ratios[1]; rtol = 0.2) for r in ratios) +end + @testset "RH method: sparse structure" begin # hessian's sparse structure shouldn't change. time_step = 1 diff --git a/test/test_iterative_methods.jl b/test/test_iterative_methods.jl index eed44ec7..4047f3e3 100644 --- a/test/test_iterative_methods.jl +++ b/test/test_iterative_methods.jl @@ -1,3 +1,39 @@ +@testset "Iwamoto multiplier: exact quadratic minimizer" begin + # The optimal multiplier must minimize the exact quadratic mismatch model + # g(μ) = ‖f₀ + μ·b + μ²·a‖² over μ ∈ [0, 1] (b = J·Δx, a = true quadratic term). + g_model(f0, b, a, μ) = sum(abs2, f0 .+ μ .* b .+ (μ * μ) .* a) + function brute(f0, b, a) + best_μ, best_g = 0.0, g_model(f0, b, a, 0.0) + for k in 0:200_000 + μ = k / 200_000 + gv = g_model(f0, b, a, μ) + gv < best_g && ((best_μ, best_g) = (μ, gv)) + end + return best_μ + end + cases = ( + ([1.0, -2.0, 0.5], [-0.3, 1.2, -0.7], [0.4, -0.1, 0.9]), + ([2.0, 1.0], [-2.0, -1.0], [0.5, -0.3]), + ([0.1, 0.2, -0.4, 0.3], [1.0, -1.0, 0.2, 0.0], [-0.5, 0.5, -0.5, 0.5]), + ([3.0], [-3.0], [1.0]), + ) + for (f0, b, a) in cases + c_fb, c_bb = dot(f0, b), dot(b, b) + c_fa, c_ba, c_aa = dot(f0, a), dot(b, a), dot(a, a) + μ = PF._iwamoto_multiplier(2c_fb, c_bb + 2c_fa, 2c_ba, c_aa) + @test 0.0 <= μ <= 1.0 + # The analytic minimizer must be no worse than a fine brute-force grid. + @test g_model(f0, b, a, μ) <= g_model(f0, b, a, brute(f0, b, a)) + 1e-6 + end + # Newton-step special case (b = −f₀): 3-arg convenience == 4-arg general form. + f0 = [1.0, -0.5, 2.0, 0.3] + a = [0.2, 0.7, -0.4, 1.1] + b = -f0 + g0, g1, g2 = dot(f0, f0), dot(f0, a), dot(a, a) + @test PF._iwamoto_multiplier(g0, g1, g2) == PF._iwamoto_multiplier( + 2 * dot(f0, b), dot(b, b) + 2 * dot(f0, a), 2 * dot(b, a), dot(a, a)) +end + @testset "NewtonRaphsonACPowerFlow kwargs" begin # test NR kwargs. sys = PSB.build_system(PSB.PSITestSystems, "c_sys5") @@ -188,9 +224,11 @@ end end @testset "Iwamoto multiplier root-finding" begin - # Verify that _iwamoto_multiplier recovers the global minimizer of + # Verify that _iwamoto_multiplier recovers the global minimizer of the + # classical (exact-Newton-step) Iwamoto objective # g(μ) = (1-μ)²g₀ + 2μ²(1-μ)g₁ + μ⁴g₂ on [0, 1] by comparing against # a brute-force grid search. + g_classic(μ, g0, g1, g2) = (1 - μ)^2 * g0 + 2 * μ^2 * (1 - μ) * g1 + μ^4 * g2 grid = range(0.0, 1.0; length = 10001) test_cases = [ @@ -214,9 +252,9 @@ end for (g0, g1, g2, desc) in test_cases μ = PF._iwamoto_multiplier(g0, g1, g2) @test 0.0 <= μ <= 1.0 - g_opt = PF._iwamoto_objective(μ, g0, g1, g2) + g_opt = g_classic(μ, g0, g1, g2) # Brute-force minimum over the grid. - g_grid_min = minimum(PF._iwamoto_objective(m, g0, g1, g2) for m in grid) + g_grid_min = minimum(g_classic(m, g0, g1, g2) for m in grid) @test g_opt <= g_grid_min + 1e-10 end end @@ -228,21 +266,22 @@ end g1 = 10.0 g2 = 100.0 μ = PF._iwamoto_multiplier(g0, g1, g2) - g_at_mu = PF._iwamoto_objective(μ, g0, g1, g2) + g_at_mu = (1 - μ)^2 * g0 + 2 * μ^2 * (1 - μ) * g1 + μ^4 * g2 # The optimizer must be able to return μ=0 or at least match g(0)=g₀. @test g_at_mu <= g0 + 1e-12 end @testset "Iwamoto multiplier degenerate cases" begin # Degenerate cubic (g₂ ≈ 0): leading coefficient of derivative cubic is ~0. + g_classic(μ, g0, g1, g2) = (1 - μ)^2 * g0 + 2 * μ^2 * (1 - μ) * g1 + μ^4 * g2 g0 = 1.0 g1 = 0.5 g2 = 1e-35 μ = PF._iwamoto_multiplier(g0, g1, g2) @test 0.0 <= μ <= 1.0 - g_opt = PF._iwamoto_objective(μ, g0, g1, g2) + g_opt = g_classic(μ, g0, g1, g2) grid = range(0.0, 1.0; length = 10001) - g_grid_min = minimum(PF._iwamoto_objective(m, g0, g1, g2) for m in grid) + g_grid_min = minimum(g_classic(m, g0, g1, g2) for m in grid) @test g_opt <= g_grid_min + 1e-10 # Degenerate quadratic (g₂ ≈ 0 and g₁ ≈ 0): both leading coefficients ~0. diff --git a/test/test_utils/common.jl b/test/test_utils/common.jl index c75a304f..a49aca57 100644 --- a/test/test_utils/common.jl +++ b/test/test_utils/common.jl @@ -495,6 +495,165 @@ function prepare_ts_data!(data::PowerFlowData, time_steps::Int64 = 24) return end +"""Build a minimal 3-bus system with one `TapTransformer` (VOLTAGE control) and one +`SwitchedAdmittance` for testing `build_controlled_device_set`.""" +function _make_tap_shunt_system() + sys = System(100.0) + b1 = _add_simple_bus!(sys, 1, ACBusTypes.REF, 230, 1.0, 0.0) + b2 = _add_simple_bus!(sys, 2, ACBusTypes.PQ, 230, 1.0, 0.0) + b3 = _add_simple_bus!(sys, 3, ACBusTypes.PQ, 230, 1.0, 0.0) + _add_simple_source!(sys, b1, 0.0, 0.0) + _add_simple_load!(sys, b2, 0.1, 0.05) + _add_simple_load!(sys, b3, 0.1, 0.05) + # Line between buses 1 and 3 so the network is connected. + _add_simple_line!(sys, b1, b3, 1e-2, 1e-2, 0.0) + tap_arc = Arc(; from = b1, to = b2) + tx = TapTransformer(; + name = "tap_1_2", + available = true, + active_power_flow = 0.0, + reactive_power_flow = 0.0, + arc = tap_arc, + r = 0.01, + x = 0.10, + primary_shunt = 0.0 + 0.0im, + tap = 1.0, + rating = 1.0, + base_power = 100.0, + control_objective = PSY.TransformerControlObjective.VOLTAGE, + ) + add_component!(sys, tx) + sa = SwitchedAdmittance(; + name = "shunt_3", + available = true, + bus = b3, + Y = 0.0 + 0.0im, + initial_status = [0], + number_of_steps = [4], + Y_increase = [0.0 + 0.05im], + admittance_limits = (min = 0.9, max = 1.1), + ) + add_component!(sys, sa) + return sys +end + +"""Build a 3-bus system with one `TapTransformer` (VOLTAGE control) and one +`SwitchedAdmittance`, designed so the AC base case converges cleanly. + +Bus 2 carries a significant load (0.5 pu on 100 MVA base) through a low-impedance +transformer from the REF bus, so that the tap is the sole voltage-support mechanism +for that bus. Bus 3 is separately connected to the REF bus via a line and hosts the +shunt. Buses 2 and 3 are decoupled (both see the REF bus but not each other), so +shunt adjustments do not perturb the tap-controlled bus and vice-versa. + +The continuation engine with `INITIAL_CONTROL_STEEPNESS=100` typically drives the +tap to a saturation extreme (p_min or p_max) before the oscillation guard triggers. +`snap_and_restore!` then snaps the tap to the nearest discrete level and re-solves; +convergence is guaranteed. The voltage at the snapped level may be far from vset +for this network—this is expected behaviour of the oscillation guard. Integration +tests using this fixture should assert convergence and correct snapping, not voltage +proximity.""" +function _make_solvable_tap_shunt_system() + sys = System(100.0) + b1 = _add_simple_bus!(sys, 1, ACBusTypes.REF, 230, 1.0, 0.0) + b2 = _add_simple_bus!(sys, 2, ACBusTypes.PQ, 230, 1.0, 0.0) + b3 = _add_simple_bus!(sys, 3, ACBusTypes.PQ, 230, 1.0, 0.0) + _add_simple_source!(sys, b1, 0.0, 0.0) + # Significant load so the base case converges with a non-trivial solution. + load2 = PowerLoad(; + name = "load_2", + available = true, + bus = b2, + active_power = 0.5, + reactive_power = 0.25, + base_power = 100.0, + max_active_power = 100.0, + max_reactive_power = 100.0, + ) + add_component!(sys, load2) + load3 = PowerLoad(; + name = "load_3", + available = true, + bus = b3, + active_power = 0.05, + reactive_power = 0.025, + base_power = 100.0, + max_active_power = 100.0, + max_reactive_power = 100.0, + ) + add_component!(sys, load3) + # Bus 3 connected to REF bus; decoupled from bus 2. + _add_simple_line!(sys, b1, b3, 1e-2, 1e-2, 0.0) + tap_arc = Arc(; from = b1, to = b2) + tx = TapTransformer(; + name = "tap_1_2", + available = true, + active_power_flow = 0.0, + reactive_power_flow = 0.0, + arc = tap_arc, + r = 0.01, + x = 0.10, + primary_shunt = 0.0 + 0.0im, + tap = 1.0, + rating = 1.0, + base_power = 100.0, + control_objective = PSY.TransformerControlObjective.VOLTAGE, + ) + add_component!(sys, tx) + sa = SwitchedAdmittance(; + name = "shunt_3", + available = true, + bus = b3, + Y = 0.0 + 0.0im, + initial_status = [0], + number_of_steps = [4], + Y_increase = [0.0 + 0.05im], + admittance_limits = (min = 0.9, max = 1.1), + ) + add_component!(sys, sa) + return sys +end + +"""Build a 4-bus system where the TapTransformer's FROM bus is the controlled bus +(`controlled_on_primary=true`), exercising the eq.46 orientation path. + +Topology: REF(1) ─line─ PQ(2) ─tap─ PQ(3); REF(1) ─line─ PQ(4). +Bus 2 is both the FROM bus of the tap and the controlled bus (set via ext["NREG"]). +The tap has real authority over bus 2 voltage through the impedance seen by bus 2.""" +function _make_primary_controlled_tap_system() + sys = System(100.0) + b1 = _add_simple_bus!(sys, 1, ACBusTypes.REF, 230, 1.0, 0.0) + b2 = _add_simple_bus!(sys, 2, ACBusTypes.PQ, 230, 1.0, 0.0) + b3 = _add_simple_bus!(sys, 3, ACBusTypes.PQ, 230, 1.0, 0.0) + b4 = _add_simple_bus!(sys, 4, ACBusTypes.PQ, 230, 1.0, 0.0) + _add_simple_source!(sys, b1, 0.0, 0.0) + _add_simple_load!(sys, b2, 0.3, 0.15) + _add_simple_load!(sys, b3, 0.1, 0.05) + _add_simple_load!(sys, b4, 0.05, 0.025) + # Line feeds bus 2 from REF. + _add_simple_line!(sys, b1, b2, 1e-2, 5e-2, 0.0) + # Bus 4 connected to REF (keeps network connected after b3 has only the tap). + _add_simple_line!(sys, b1, b4, 1e-2, 1e-2, 0.0) + tap_arc = Arc(; from = b2, to = b3) + tx = TapTransformer(; + name = "tap_2_3", + available = true, + active_power_flow = 0.0, + reactive_power_flow = 0.0, + arc = tap_arc, + r = 0.01, + x = 0.10, + primary_shunt = 0.0 + 0.0im, + tap = 1.0, + rating = 1.0, + base_power = 100.0, + control_objective = PSY.TransformerControlObjective.VOLTAGE, + ext = Dict{String, Any}("NREG" => 2), # controlled bus = bus 2 (FROM) → primary + ) + add_component!(sys, tx) + return sys +end + function simple_lcc_system() sys = System(100.0) b1 = _add_simple_bus!(sys, 1, ACBusTypes.REF, 230, 1.1, 0.0)