From 1414ab1e032220641e39bb77e47f5351667d62e8 Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara <16385323+jd-lara@users.noreply.github.com> Date: Tue, 11 Aug 2026 12:08:31 -0600 Subject: [PATCH 01/11] Disable reserve code referencing PSY types removed on psy6 (#232) PowerOperationsModels does not load against the psy6 line of PowerSystems: UndefVarError: `ReserveNonSpinning` not defined in `PowerSystems` The psy6 reserve rework replaced ConstantReserve, ConstantReserveGroup, ReserveNonSpinning and VariableReserveNonSpinning with OnlineReserve, OfflineReserve and GroupReserve. POM still dispatches on the old names, so the module fails at precompile and nothing downstream can be built or solved. Comment out every path that references the four removed types so the package loads again. Each site is marked `PSY6-PORT-DISABLED:` with the symbol that went away, and the original expression is quoted where a clause was dropped from a condition or a Union, so the port back is mechanical. This is deliberately the minimum needed to compile, not a port. Reserve modelling on psy6 is unimplemented until the new reserve semantics are wired up, which is the follow-up PR. One behaviour change beyond commenting out: the empty-contributing-devices guard in `_populate_contributing_devices!` previously excepted ConstantReserveGroup, whose empty map was by design. With that type gone the guard now applies to every service type. No currently existing type is affected, but a future group-reserve type will need the exception restored. Not addressed here: the test suite and `ext/PowerFlowsExt` still reference these types, so `Pkg.test()` does not run. Only `using PowerOperationsModels` is fixed. Verified against PowerSystems jd/schema_matching: an RTS-GMLC unit commitment parsed by PowerTableDataParser builds and solves on CopperPlate, and on PTDF with both StaticBranch and StaticBranchBounds. --- src/core/problem_template.jl | 12 +++++++--- src/services_models/reserve_group.jl | 6 +++++ src/services_models/reserves.jl | 25 +++++++++++++++++++-- src/services_models/service_slacks.jl | 4 +++- src/services_models/services_constructor.jl | 9 ++++++++ 5 files changed, 50 insertions(+), 6 deletions(-) diff --git a/src/core/problem_template.jl b/src/core/problem_template.jl index 80c72a9..697ed9a 100644 --- a/src/core/problem_template.jl +++ b/src/core/problem_template.jl @@ -251,8 +251,9 @@ function _populate_contributing_devices!( # A reserve or interface with no available provider can never meet its requirement, # so error rather than let it silently force slacks or go infeasible. # ConstantReserveGroup aggregates other services, so its empty map is by design. - if !(service_type <: PSY.ConstantReserveGroup) && - isempty(get_contributing_devices_map(service_model, service_name)) + # PSY6-PORT-DISABLED: PSY.ConstantReserveGroup removed on jd/schema_matching + # (dropped from this condition; original: `!(service_type <: PSY.ConstantReserveGroup) &&`) + if isempty(get_contributing_devices_map(service_model, service_name)) error( "Service \"$(service_name)\" of type $(typeof(service)) has no available contributing devices/branches. Assign available contributing devices/branches to it in the system data, or remove its service model from the template.", ) @@ -283,6 +284,8 @@ function _modify_device_model!( return end +# PSY6-PORT-DISABLED: PSY.ReserveNonSpinning removed on jd/schema_matching +#= function _modify_device_model!( ::Dict{Symbol, DeviceModel}, ::ServiceModel{<:PSY.ReserveNonSpinning, <:AbstractReservesFormulation}, @@ -290,6 +293,7 @@ function _modify_device_model!( ) return end +=# function _modify_device_model!( ::Dict{Symbol, DeviceModel}, @@ -312,7 +316,9 @@ function _add_services_to_device_model!(template::PowerOperationsProblemTemplate devices_template = get_device_models(template) for (service_key, service_model) in service_models S = get_component_type(service_model) - (S <: PSY.AGC || S <: PSY.ConstantReserveGroup) && continue + # PSY6-PORT-DISABLED: PSY.ConstantReserveGroup removed on jd/schema_matching + # (dropped from this disjunction; original: `S <: PSY.AGC || S <: PSY.ConstantReserveGroup`) + S <: PSY.AGC && continue contributing_devices = get_contributing_devices(service_model) isempty(contributing_devices) && continue _modify_device_model!(devices_template, service_model, contributing_devices) diff --git a/src/services_models/reserve_group.jl b/src/services_models/reserve_group.jl index c7b4c7f..ad8b0aa 100644 --- a/src/services_models/reserve_group.jl +++ b/src/services_models/reserve_group.jl @@ -1,3 +1,5 @@ +# PSY6-PORT-DISABLED: PSY.ConstantReserveGroup removed on jd/schema_matching +#= function get_default_time_series_names( ::Type{PSY.ConstantReserveGroup{T}}, ::Type{GroupReserve}) where {T <: PSY.ReserveDirection} @@ -9,6 +11,7 @@ function get_default_attributes( ::Type{GroupReserve}) where {T <: PSY.ReserveDirection} return Dict{String, Any}() end +=# ############################### Reserve Variables` ######################################### """ @@ -32,6 +35,8 @@ function check_activeservice_variables( end ################################## Reserve Requirement Constraint ########################## +# PSY6-PORT-DISABLED: PSY.ConstantReserveGroup removed on jd/schema_matching +#= """ This function creates the requirement constraint that will be attained by the appropriate services """ @@ -65,6 +70,7 @@ function add_constraints!( return end +=# # Collect the group's contributing reserve variables into one bucket per time step, so the # constraint loop above indexes straight in rather than re-scanning per `(group, t)`. Services diff --git a/src/services_models/reserves.jl b/src/services_models/reserves.jl index 0d35a15..e3e69d9 100644 --- a/src/services_models/reserves.jl +++ b/src/services_models/reserves.jl @@ -11,11 +11,14 @@ get_variable_upper_bound(::Type{ActivePowerReserveVariable}, r::Union{PSY.Reserv get_variable_lower_bound(::Type{ActivePowerReserveVariable}, ::PSY.Reserve, ::PSY.Device, ::Type) = 0.0 ############################### ActivePowerReserveVariable, ReserveNonSpinning ######################################### +# PSY6-PORT-DISABLED: PSY.ReserveNonSpinning removed on jd/schema_matching +#= get_variable_binary(::Type{ActivePowerReserveVariable}, ::Type{<:PSY.ReserveNonSpinning}, ::Type{<:AbstractReservesFormulation}) = false function get_variable_upper_bound(::Type{ActivePowerReserveVariable}, r::PSY.ReserveNonSpinning, d::PSY.Device, ::Type{<:AbstractReservesFormulation}) return PSY.get_max_output_fraction(r) * PSY.get_max_active_power(d, PSY.SU) end get_variable_lower_bound(::Type{ActivePowerReserveVariable}, ::PSY.ReserveNonSpinning, ::PSY.Device, ::Type) = 0.0 +=# ############################### ServiceRequirementVariable, ReserveDemandCurve ################################ @@ -29,7 +32,8 @@ get_variable_lower_bound(::Type{ServiceRequirementVariable}, ::Union{PSY.Reserve _get_requirement(service) = PSY.get_requirement(service, PSY.SU) get_multiplier_value(::Type{RequirementTimeSeriesParameter}, d::PSY.Reserve, ::Type{<:AbstractReservesFormulation}) = _get_requirement(d) -get_multiplier_value(::Type{RequirementTimeSeriesParameter}, d::PSY.ReserveNonSpinning, ::Type{<:AbstractReservesFormulation}) = _get_requirement(d) +# PSY6-PORT-DISABLED: PSY.ReserveNonSpinning removed on jd/schema_matching +# get_multiplier_value(::Type{RequirementTimeSeriesParameter}, d::PSY.ReserveNonSpinning, ::Type{<:AbstractReservesFormulation}) = _get_requirement(d) get_parameter_multiplier(::Type{<:VariableValueParameter}, d::Type{<:PSY.AbstractReserve}, ::Type{<:AbstractReservesFormulation}) = 1.0 get_initial_parameter_value(::Type{<:VariableValueParameter}, d::Type{<:PSY.AbstractReserve}, ::Type{<:AbstractReservesFormulation}) = 0.0 @@ -53,12 +57,15 @@ function get_initial_conditions_service_model( return ServiceModel(T, D) end +# PSY6-PORT-DISABLED: PSY.VariableReserveNonSpinning removed on jd/schema_matching +#= function get_initial_conditions_service_model( ::IOM.AbstractOptimizationModel, ::ServiceModel{T, D}, ) where {T <: PSY.VariableReserveNonSpinning, D <: AbstractReservesFormulation} return ServiceModel(T, D) end +=# function get_default_time_series_names( ::Type{<:PSY.Reserve}, @@ -69,6 +76,8 @@ function get_default_time_series_names( ) end +# PSY6-PORT-DISABLED: PSY.ReserveNonSpinning removed on jd/schema_matching +#= function get_default_time_series_names( ::Type{<:PSY.ReserveNonSpinning}, ::Type{NonSpinningReserve}, @@ -77,6 +86,7 @@ function get_default_time_series_names( RequirementTimeSeriesParameter => "requirement", ) end +=# function get_default_time_series_names( ::Type{T}, @@ -92,12 +102,15 @@ function get_default_attributes( return Dict{String, Any}() end +# PSY6-PORT-DISABLED: PSY.ReserveNonSpinning removed on jd/schema_matching +#= function get_default_attributes( ::Type{<:PSY.ReserveNonSpinning}, ::Type{<:AbstractReservesFormulation}, ) return Dict{String, Any}() end +=# """ Add variables for ServiceRequirementVariable for StepWiseCostReserve @@ -268,6 +281,8 @@ function add_constraints!( return end +# PSY6-PORT-DISABLED: PSY.ConstantReserve removed on jd/schema_matching +#= function add_constraints!( container::OptimizationContainer, T::Type{RequirementConstraint}, @@ -302,6 +317,7 @@ function add_constraints!( return end +=# function add_to_objective_function!( container::OptimizationContainer, @@ -461,6 +477,8 @@ function add_constraints!( return end +# PSY6-PORT-DISABLED: PSY.VariableReserveNonSpinning removed on jd/schema_matching +#= function add_constraints!( container::OptimizationContainer, T::Type{ReservePowerConstraint}, @@ -510,6 +528,7 @@ function add_constraints!( end return end +=# function _add_reserve_power_constraint_device!( cons, @@ -658,7 +677,9 @@ function add_reserves_proportional_cost!( contributing_names::Vector{String}; skip_devices = Set{String}(), ) where { - T <: Union{PSY.Reserve, PSY.ReserveNonSpinning}, + # PSY6-PORT-DISABLED: PSY.ReserveNonSpinning removed on jd/schema_matching + # (dropped from this Union; original: `T <: Union{PSY.Reserve, PSY.ReserveNonSpinning}`) + T <: PSY.Reserve, U <: ActivePowerReserveVariable, V <: AbstractReservesFormulation, } diff --git a/src/services_models/service_slacks.jl b/src/services_models/service_slacks.jl index 2463696..21db429 100644 --- a/src/services_models/service_slacks.jl +++ b/src/services_models/service_slacks.jl @@ -1,8 +1,10 @@ +# PSY6-PORT-DISABLED: PSY.ReserveNonSpinning removed on jd/schema_matching +# (dropped from this Union; original: `T <: Union{PSY.Reserve, PSY.ReserveNonSpinning}`) function add_reserve_slacks!( container::OptimizationContainer, ::Type{T}, service_names::Vector{String}, -) where {T <: Union{PSY.Reserve, PSY.ReserveNonSpinning}} +) where {T <: PSY.Reserve} time_steps = get_time_steps(container) # Dense 2D container keyed `[service_name, time]`, built once per service type over all # the type's services (`use_slacks` is per type). Lower bound 0, penalty in objective. diff --git a/src/services_models/services_constructor.jl b/src/services_models/services_constructor.jl index ae39b2a..3f45c3a 100644 --- a/src/services_models/services_constructor.jl +++ b/src/services_models/services_constructor.jl @@ -132,6 +132,8 @@ function construct_service!( return end +# PSY6-PORT-DISABLED: PSY.ConstantReserve removed on jd/schema_matching +#= # ConstantReserve has no requirement time series, so its argument stage skips the # parameter add. The model stage is the shared `SR <: PSY.AbstractReserve` method below. function construct_service!( @@ -165,6 +167,7 @@ function construct_service!( end return end +=# # Shared RangeReserve model stage for both `PSY.Reserve` and `PSY.ConstantReserve`; # the inner `add_constraints!` calls resolve per service type. @@ -390,6 +393,8 @@ function construct_service!( end =# +# PSY6-PORT-DISABLED: PSY.ConstantReserveGroup removed on jd/schema_matching +#= """ Constructs a service for ConstantReserveGroup. """ @@ -441,6 +446,7 @@ function construct_service!( add_constraint_dual!(container, sys, model) return end +=# function construct_service!( container::OptimizationContainer, @@ -520,6 +526,8 @@ function construct_service!( return end +# PSY6-PORT-DISABLED: PSY.ReserveNonSpinning removed on jd/schema_matching +#= function construct_service!( container::OptimizationContainer, sys::PSY.System, @@ -596,6 +604,7 @@ function construct_service!( add_constraint_dual!(container, sys, model) return end +=# function construct_service!( container::OptimizationContainer, From ea8e374d8ed22b253ed2831a78f7ad016940bbc5 Mon Sep 17 00:00:00 2001 From: kdayday Date: Wed, 12 Aug 2026 09:57:00 -0600 Subject: [PATCH 02/11] Add PSY-style docs workflow concurrency guard. (#229) PR previews get their own cancelable group; main and tag deploys share docs-deploy and queue instead of canceling. Co-authored-by: Cursor --- .github/workflows/docs.yml | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index a94b6e5..5ca0717 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -8,9 +8,12 @@ on: tags: '*' pull_request: +# Prevent simultaneous gh-pages deployments (tag + branch race condition). +# PR previews each get their own group and cancel on new pushes. +# All other deploys (main, tags) share one group and queue instead of cancel. concurrency: - group: docs-deploy - cancel-in-progress: false + group: ${{ github.event_name == 'pull_request' && format('docs-pr-{0}', github.event.pull_request.number) || 'docs-deploy' }} + cancel-in-progress: ${{ github.event_name == 'pull_request' }} jobs: build: From a4c068e7f3a5945c64e197e5d963f1ff2722f70a Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Thu, 13 Aug 2026 14:14:06 -0400 Subject: [PATCH 03/11] refactor voltage-only ac network flow code --- src/ac_transmission_models/AC_branches.jl | 396 +++------ .../voltage_control_tap_models.jl | 788 ------------------ test/test_native_transformer_tap.jl | 79 +- 3 files changed, 139 insertions(+), 1124 deletions(-) delete mode 100644 src/ac_transmission_models/voltage_control_tap_models.jl diff --git a/src/ac_transmission_models/AC_branches.jl b/src/ac_transmission_models/AC_branches.jl index 76a7a8d..c71f293 100644 --- a/src/ac_transmission_models/AC_branches.jl +++ b/src/ac_transmission_models/AC_branches.jl @@ -1411,30 +1411,6 @@ function _add_flow_constraint_containers!( return cons_pft, cons_qft, cons_ptf, cons_qtf end -# Pure, tap-free π-model coefficients shared by the polar (ACP) and rectangular (ACR) -# Ohm's law, for both the fixed-tap StaticBranch path and the variable-tap VoltageControlTap -# path. `cs`/`sn` are the phase-shift trig; `a_cos`/`a_sin`/`c_cos`/`d_sin` are the tm-free -# coupling coefficients (each divided by the live tap at the constraint site — `tm` for -# fixed tap, `TapRatioVariable[name, t]` for variable tap). ACR uses `e_sin = -d_sin`. -function _tap_flow_coefficients(g, b, g_fr, b_fr, g_to, b_to, shift) - cs = cos(shift) - sn = sin(shift) - return ( - cs = cs, - sn = sn, - g = g, - b = b, - g_fr = g_fr, - b_fr = b_fr, - gg_to = g + g_to, - bb_to = b + b_to, - a_cos = -g * cs + b * sn, - a_sin = -b * cs - g * sn, - c_cos = -g * cs - b * sn, - d_sin = b * cs - g * sn, - ) -end - # Slack holders for the equality/limit rows. `_SlackPair` carries a metaed upper/lower pair # (equality relaxation, term `up - lo`); `_UpperSlack` carries a one-sided upper slack # (quadratic-limit relaxation, term `up`). The no-slack twins contribute a constant 0.0 so @@ -1576,107 +1552,105 @@ function _current_magnitude_slacks( return _UpperSlack(get_variable(container, FlowActivePowerSlackUpperBound, T, meta)) end -# Polar (ACP) π-model Ohm's law for one branch, one time step. `coef` from -# `_tap_flow_coefficients`; `tap` is the constant `tm` (fixed tap) or the -# `TapRatioVariable` (variable tap). The constraints reduce term-for-term to the -# fixed-tap StaticBranch form when `tap == tm`. -function _add_tap_acp_flow!( - jump_model, cons_pft, cons_qft, cons_ptf, cons_qtf, pft, qft, ptf, qtf, - name, t, vmf, vmt, θ, coef, tap, p_ft_slack, q_ft_slack, p_tf_slack, q_tf_slack, +function _voltage_products( + container::OptimizationContainer, + ::NetworkModel{ACPNetworkModel}, + ::Type{<:PSY.ACTransmission}, + ::String, + from_bus::String, + to_bus::String, + t::Int, ) - cons_pft[name, t] = JuMP.@constraint( - jump_model, - pft[name, t] == - (coef.g / tap^2 + coef.g_fr) * vmf^2 + - coef.a_cos / tap * vmf * vmt * cos(θ) + - coef.a_sin / tap * vmf * vmt * sin(θ) + p_ft_slack, - ) - cons_qft[name, t] = JuMP.@constraint( - jump_model, - qft[name, t] == - -(coef.b / tap^2 + coef.b_fr) * vmf^2 + - (-coef.a_sin) / tap * vmf * vmt * cos(θ) + - coef.a_cos / tap * vmf * vmt * sin(θ) + q_ft_slack, - ) - cons_ptf[name, t] = JuMP.@constraint( - jump_model, - ptf[name, t] == - coef.gg_to * vmt^2 + - coef.c_cos / tap * vmt * vmf * cos(θ) + - coef.d_sin / tap * vmt * vmf * sin(θ) + p_tf_slack, - ) - cons_qtf[name, t] = JuMP.@constraint( - jump_model, - qtf[name, t] == - -coef.bb_to * vmt^2 + - coef.d_sin / tap * vmt * vmf * cos(θ) + - (-coef.c_cos) / tap * vmt * vmf * sin(θ) + q_tf_slack, + vm = get_variable(container, VoltageMagnitude, PSY.ACBus) + va = get_variable(container, VoltageAngle, PSY.ACBus) + vmf = vm[from_bus, t] + vmt = vm[to_bus, t] + θ = va[from_bus, t] - va[to_bus, t] + return ( + v2_fr = vmf^2, + v2_to = vmt^2, + vv_cos = vmf * vmt * cos(θ), + vv_sin = vmf * vmt * sin(θ), ) - return end -# Rectangular (ACR) π-model Ohm's law for one branch, one time step. Same coefficients as -# ACP; the rectangular substitution replaces vmf²/vmf·vmt·cos/vmf·vmt·sin with the -# pre-built bilinears `vsq_fr`/`vv_cos`/`vv_sin`. `e_sin = -d_sin` (rectangular sin sign). -function _add_tap_acr_flow!( - jump_model, cons_pft, cons_qft, cons_ptf, cons_qtf, pft, qft, ptf, qtf, - name, t, vsq_fr, vsq_to, vv_cos, vv_sin, coef, tap, - p_ft_slack, q_ft_slack, p_tf_slack, q_tf_slack, +function _voltage_products( + container::OptimizationContainer, + ::NetworkModel{ACRNetworkModel}, + ::Type{<:PSY.ACTransmission}, + ::String, + from_bus::String, + to_bus::String, + t::Int, ) - e_sin = -coef.d_sin - cons_pft[name, t] = JuMP.@constraint( - jump_model, - pft[name, t] == - (coef.g / tap^2 + coef.g_fr) * vsq_fr + - coef.a_cos / tap * vv_cos + - coef.a_sin / tap * vv_sin + p_ft_slack, - ) - cons_qft[name, t] = JuMP.@constraint( - jump_model, - qft[name, t] == - -(coef.b / tap^2 + coef.b_fr) * vsq_fr + - (-coef.a_sin) / tap * vv_cos + - coef.a_cos / tap * vv_sin + q_ft_slack, - ) - cons_ptf[name, t] = JuMP.@constraint( - jump_model, - ptf[name, t] == - coef.gg_to * vsq_to + - coef.c_cos / tap * vv_cos - - e_sin / tap * vv_sin + p_tf_slack, - ) - cons_qtf[name, t] = JuMP.@constraint( - jump_model, - qtf[name, t] == - -coef.bb_to * vsq_to - - e_sin / tap * vv_cos - - coef.c_cos / tap * vv_sin + q_tf_slack, + vr = get_variable(container, VoltageReal, PSY.ACBus) + vi = get_variable(container, VoltageImaginary, PSY.ACBus) + vr_fr = vr[from_bus, t] + vr_to = vr[to_bus, t] + vi_fr = vi[from_bus, t] + vi_to = vi[to_bus, t] + return ( + v2_fr = vr_fr^2 + vi_fr^2, + v2_to = vr_to^2 + vi_to^2, + vv_cos = vr_fr * vr_to + vi_fr * vi_to, + vv_sin = vi_fr * vr_to - vr_fr * vi_to, ) - return end -""" -Add full π-model rectangular AC Ohm's law constraints for ACBranch under ACRNetworkModel. - -Four constraints per branch per time step (p_ft, q_ft, p_tf, q_tf) relate the four -directional flow variables to rectangular voltage components (vr, vi) via the -π-equivalent circuit. Rectangular identity applied to the ACP polar expressions: - vmf^2 → vr_fr^2 + vi_fr^2 - vmf*vmt*cos(θ) → vr_fr*vr_to + vi_fr*vi_to - vmf*vmt*sin(θ) → vi_fr*vr_to - vr_fr*vi_to -""" +function _voltage_products( + container::OptimizationContainer, + ::NetworkModel{LPACCNetworkModel}, + ::Type{T}, + name::String, + from_bus::String, + to_bus::String, + t::Int, +) where {T <: PSY.ACTransmission} + va = get_variable(container, VoltageAngle, PSY.ACBus) + phi = get_variable(container, VoltageDeviation, PSY.ACBus) + cs = get_variable(container, CosineApproximation, T) + phi_fr = phi[from_bus, t] + phi_to = phi[to_bus, t] + return ( + v2_fr = 1.0 + 2.0 * phi_fr, + v2_to = 1.0 + 2.0 * phi_to, + vv_cos = cs[name, t] + phi_fr + phi_to, + vv_sin = va[from_bus, t] - va[to_bus, t], + ) +end + +# Ybus terms, supporting Float64 and VariableRef taps. PNM's ybus functions +# use imaginary numbers which VariableRef doesn't support. +function _tapped_admittance(adm, tap) + cs = cos(adm.shift) + sn = sin(adm.shift) + return ( + g11 = adm.g / tap^2 + adm.g_fr, + b11 = adm.b / tap^2 + adm.b_fr, + g12 = (-adm.g * cs + adm.b * sn) / tap, + b12 = (-adm.b * cs - adm.g * sn) / tap, + g21 = (-adm.g * cs - adm.b * sn) / tap, + b21 = (adm.g * sn - adm.b * cs) / tap, + g22 = adm.g + adm.g_to, + b22 = adm.b + adm.b_to, + ) +end + +# Voltage-only AC networks. function add_constraints!( container::OptimizationContainer, sys::PSY.System, ::Type{NetworkFlowConstraint}, devices::IS.FlattenIteratorWrapper{T}, device_model::DeviceModel{T, U}, - network_model::NetworkModel{ACRNetworkModel}, -) where {T <: PSY.ACTransmission, U <: AbstractBranchFormulation} + network_model::NetworkModel{N}, +) where { + T <: PSY.ACTransmission, + U <: AbstractBranchFormulation, + N <: Union{ACPNetworkModel, ACRNetworkModel, LPACCNetworkModel}, +} time_steps = get_time_steps(container) - vr = get_variable(container, VoltageReal, PSY.ACBus) - vi = get_variable(container, VoltageImaginary, PSY.ACBus) pft = get_variable(container, FlowActivePowerFromToVariable, T) ptf = get_variable(container, FlowActivePowerToFromVariable, T) qft = get_variable(container, FlowReactivePowerFromToVariable, T) @@ -1689,34 +1663,42 @@ function add_constraints!( cons_pft, cons_qft, cons_ptf, cons_qtf = _add_flow_constraint_containers!(container, T, branch_names) jump_model = get_jump_model(container) + slacks = _flow_equality_slacks(container, device_model, T) + p_ft_slack = _slack_term(slacks.p_ft, name, t) + q_ft_slack = _slack_term(slacks.q_ft, name, t) + p_tf_slack = _slack_term(slacks.p_tf, name, t) + q_tf_slack = _slack_term(slacks.q_tf, name, t) for g_geom in geoms name = g_geom.name adm = g_geom.adm - tm = adm.tap from_bus = g_geom.from_name to_bus = g_geom.to_name - coef = _tap_flow_coefficients( - adm.g, adm.b, adm.g_fr, adm.b_fr, adm.g_to, adm.b_to, adm.shift, - ) for t in time_steps - vr_fr = vr[from_bus, t] - vr_to = vr[to_bus, t] - vi_fr = vi[from_bus, t] - vi_to = vi[to_bus, t] - vsq_fr = vr_fr^2 + vi_fr^2 - vsq_to = vr_to^2 + vi_to^2 - vv_cos = vr_fr * vr_to + vi_fr * vi_to - vv_sin = vi_fr * vr_to - vr_fr * vi_to - _add_tap_acr_flow!( - jump_model, cons_pft, cons_qft, cons_ptf, cons_qtf, pft, qft, ptf, qtf, - name, t, vsq_fr, vsq_to, vv_cos, vv_sin, coef, tm, - _slack_term(slacks.p_ft, name, t), - _slack_term(slacks.q_ft, name, t), - _slack_term(slacks.p_tf, name, t), - _slack_term(slacks.q_tf, name, t), + vp = _voltage_products(container, network_model, T, name, from_bus, to_bus, t) + y = _tapped_admittance(adm, adm.tap) + + cons_pft[name, t] = JuMP.@constraint( + jump_model, + pft[name, t] == + y.g11 * vp.v2_fr + y.g12 * vp.vv_cos + y.b12 * vp.vv_sin + p_ft_slack, + ) + cons_qft[name, t] = JuMP.@constraint( + jump_model, + qft[name, t] == + -y.b11 * vp.v2_fr - y.b12 * vp.vv_cos + y.g12 * vp.vv_sin + q_ft_slack, + ) + cons_ptf[name, t] = JuMP.@constraint( + jump_model, + ptf[name, t] == + y.g22 * vp.v2_to + y.g21 * vp.vv_cos - y.b21 * vp.vv_sin + p_tf_slack, + ) + cons_qtf[name, t] = JuMP.@constraint( + jump_model, + qtf[name, t] == + -y.b22 * vp.v2_to - y.b21 * vp.vv_cos - y.g21 * vp.vv_sin + q_tf_slack, ) end end @@ -1869,105 +1851,6 @@ function _entry_angle_limits(geometry, device_by_name::Dict{String, <:PSY.ACTran return (min = -π / 2, max = π / 2) end -""" -Add the LPAC-linearized π-model AC Ohm's law constraints for ACBranch under -LPACCNetworkModel. - -Four constraints per branch per time step (p_ft, q_ft, p_tf, q_tf) relate the directional -flow variables to the voltage-magnitude deviations (phi), the bus-pair cosine variable (cs), -and the voltage-angle difference (va_fr - va_to). Transcribed from PowerModels `lpac.jl` -`constraint_ohms_yt_from/to` for `AbstractLPACCNetworkModel`, with `tr = tm·cos(shift)`, -`ti = tm·sin(shift)`. -""" -function add_constraints!( - container::OptimizationContainer, - sys::PSY.System, - ::Type{NetworkFlowConstraint}, - devices::IS.FlattenIteratorWrapper{T}, - device_model::DeviceModel{T, U}, - network_model::NetworkModel{LPACCNetworkModel}, -) where {T <: PSY.ACTransmission, U <: AbstractBranchFormulation} - time_steps = get_time_steps(container) - - va = get_variable(container, VoltageAngle, PSY.ACBus) - phi = get_variable(container, VoltageDeviation, PSY.ACBus) - cs = get_variable(container, CosineApproximation, T) - pft = get_variable(container, FlowActivePowerFromToVariable, T) - ptf = get_variable(container, FlowActivePowerToFromVariable, T) - qft = get_variable(container, FlowReactivePowerFromToVariable, T) - qtf = get_variable(container, FlowReactivePowerToFromVariable, T) - - number_to_name = _retained_number_to_name(sys, network_model) - geoms = - _branch_geometries(number_to_name, network_model, devices, T, NetworkFlowConstraint) - branch_names = [g.name for g in geoms] - cons_pft, cons_qft, cons_ptf, cons_qtf = - _add_flow_constraint_containers!(container, T, branch_names) - - jump_model = get_jump_model(container) - slacks = _flow_equality_slacks(container, device_model, T) - for g_geom in geoms - name = g_geom.name - adm = g_geom.adm - g = adm.g - b = adm.b - g_fr = adm.g_fr - b_fr = adm.b_fr - g_to = adm.g_to - b_to = adm.b_to - tm = adm.tap - nominal_shift = adm.shift - from_bus = g_geom.from_name - to_bus = g_geom.to_name - tr = tm * cos(nominal_shift) - ti = tm * sin(nominal_shift) - # Coupling coefficients (identical to ACP / PowerModels lpac.jl). - c_cos_fr = (-g * tr + b * ti) / tm^2 - c_sin_fr = (-b * tr - g * ti) / tm^2 - c_cos_to = (-g * tr - b * ti) / tm^2 - c_sin_to = (-b * tr + g * ti) / tm^2 - - for t in time_steps - phi_fr = phi[from_bus, t] - phi_to = phi[to_bus, t] - vad = va[from_bus, t] - va[to_bus, t] - cs_b = cs[name, t] - - # Shared affine terms reused across the four flow constraints: - # cs_sum = cs + phi_fr + phi_to, dev_* = 1 + 2·phi_* - cs_sum = cs_b + phi_fr + phi_to - dev_fr = 1.0 + 2.0 * phi_fr - dev_to = 1.0 + 2.0 * phi_to - - cons_pft[name, t] = JuMP.@constraint( - jump_model, - pft[name, t] == - (g / tm^2 + g_fr) * dev_fr + c_cos_fr * cs_sum + c_sin_fr * vad + - _slack_term(slacks.p_ft, name, t), - ) - cons_qft[name, t] = JuMP.@constraint( - jump_model, - qft[name, t] == - -(b / tm^2 + b_fr) * dev_fr - c_sin_fr * cs_sum + c_cos_fr * vad + - _slack_term(slacks.q_ft, name, t), - ) - cons_ptf[name, t] = JuMP.@constraint( - jump_model, - ptf[name, t] == - (g + g_to) * dev_to + c_cos_to * cs_sum + c_sin_to * (-vad) + - _slack_term(slacks.p_tf, name, t), - ) - cons_qtf[name, t] = JuMP.@constraint( - jump_model, - qtf[name, t] == - -(b + b_to) * dev_to - c_sin_to * cs_sum + c_cos_to * (-vad) + - _slack_term(slacks.q_tf, name, t), - ) - end - end - return -end - ################################## IVRNetworkModel branch constraints ################## # Compute the per-unit current rating bound for an IVR branch variable. @@ -2769,65 +2652,6 @@ function add_constraints!( return end -""" -Add full π-model AC Ohm's law constraints for ACBranch under ACPNetworkModel. - -Four constraints per branch per time step (p_ft, q_ft, p_tf, q_tf) relate the four -directional flow variables to voltage magnitudes and angles via the π-equivalent circuit. -""" -function add_constraints!( - container::OptimizationContainer, - sys::PSY.System, - ::Type{NetworkFlowConstraint}, - devices::IS.FlattenIteratorWrapper{T}, - device_model::DeviceModel{T, U}, - network_model::NetworkModel{ACPNetworkModel}, -) where {T <: PSY.ACTransmission, U <: AbstractBranchFormulation} - time_steps = get_time_steps(container) - - va = get_variable(container, VoltageAngle, PSY.ACBus) - vm = get_variable(container, VoltageMagnitude, PSY.ACBus) - pft = get_variable(container, FlowActivePowerFromToVariable, T) - ptf = get_variable(container, FlowActivePowerToFromVariable, T) - qft = get_variable(container, FlowReactivePowerFromToVariable, T) - qtf = get_variable(container, FlowReactivePowerToFromVariable, T) - - number_to_name = _retained_number_to_name(sys, network_model) - geoms = - _branch_geometries(number_to_name, network_model, devices, T, NetworkFlowConstraint) - branch_names = [g.name for g in geoms] - cons_pft, cons_qft, cons_ptf, cons_qtf = - _add_flow_constraint_containers!(container, T, branch_names) - slacks = _flow_equality_slacks(container, device_model, T) - - for g_geom in geoms - name = g_geom.name - adm = g_geom.adm - tm = adm.tap - from_bus = g_geom.from_name - to_bus = g_geom.to_name - coef = _tap_flow_coefficients( - adm.g, adm.b, adm.g_fr, adm.b_fr, adm.g_to, adm.b_to, adm.shift, - ) - - for t in time_steps - θ = va[from_bus, t] - va[to_bus, t] - vmf = vm[from_bus, t] - vmt = vm[to_bus, t] - jump_model = get_jump_model(container) - _add_tap_acp_flow!( - jump_model, cons_pft, cons_qft, cons_ptf, cons_qtf, pft, qft, ptf, qtf, - name, t, vmf, vmt, θ, coef, tm, - _slack_term(slacks.p_ft, name, t), - _slack_term(slacks.q_ft, name, t), - _slack_term(slacks.p_tf, name, t), - _slack_term(slacks.q_tf, name, t), - ) - end - end - return -end - ################################## DCPLLNetworkModel branch constraints ################# # Tighten a flow variable to ±rate without loosening any bound it already carries (a diff --git a/src/ac_transmission_models/voltage_control_tap_models.jl b/src/ac_transmission_models/voltage_control_tap_models.jl deleted file mode 100644 index c694b8d..0000000 --- a/src/ac_transmission_models/voltage_control_tap_models.jl +++ /dev/null @@ -1,788 +0,0 @@ -################################################################################# -# Voltage-controlling tap transformer (Family B). -# -# `VoltageControlTap` models the off-nominal tap ratio of a `PSY.TwoWindingTransformer` -# as a bounded continuous decision variable `t ∈ [t_min, t_max]` -# (`TapRatioVariable`) that enters the AC π-model Ohm's law nonlinearly (the fixed -# tap `tm` of the StaticBranch law is replaced by the variable `t`, so the self -# terms scale as `1/t²` and the coupling terms as `1/t`). The control objective is -# applied count-invariantly with a single `JuMP.fix` on an already-created variable: -# VOLTAGE → fix the regulated-bus VoltageMagnitude to voltage_setpoint -# REACTIVE_POWER_FLOW → fix the from-to reactive flow to reactive_power_flow -# ACTIVE_POWER_FLOW → fix the from-to active flow to active_power_flow -# No per-mode constraint is ever added — the variable/constraint containers are the -# same in every mode (a `FixRef` lives at the variable level). -# -# Voltage-objective regulation: under ACP, the scalar VoltageMagnitude is pinned -# directly; under ACR/IVR, a per-device RegulatedVoltageMagnitude aux variable is -# tied to the rectangular components via RegulatedVoltageMagnitudeConstraint and then -# fixed. Reactive/active-flow objectives share a common path across ACP and ACR. -# The formulation is dropped from DC templates via `models_reactive_power`. -################################################################################# - -# Finite tap-ratio bounds (pu turns ratio) for the control variable `t`. A -# non-finite limit is a data error (Principle 0 / IPOPT). -function _tap_ratio_limits(d::PSY.TwoWindingTransformer) - lims = PSY.get_tap_limits(d) - lo = lims.min - hi = lims.max - if !(isfinite(lo) && isfinite(hi)) - error( - "TwoWindingTransformer $(PSY.get_name(d)) has non-finite tap_limits ", - "($(lo), $(hi)); cannot bound TapRatioVariable", - ) - end - if lo <= 0.0 - error( - "TwoWindingTransformer $(PSY.get_name(d)) has a non-positive tap lower limit ", - "($(lo)); the variable-tap Ohm's law divides by t and requires t > 0", - ) - end - if hi < lo - error( - "TwoWindingTransformer $(PSY.get_name(d)) has tap_limits.max < tap_limits.min ", - "($(hi) < $(lo))", - ) - end - return (min = lo, max = hi) -end - -################################################################################# -# TapRatioVariable traits -################################################################################# - -get_variable_binary( - ::Type{TapRatioVariable}, - ::Type{<:PSY.TwoWindingTransformer}, - ::Type{VoltageControlTap}, -) = false - -get_variable_multiplier( - ::Type{TapRatioVariable}, - ::Type{<:PSY.TwoWindingTransformer}, - ::Type{VoltageControlTap}, -) = 1.0 - -function get_variable_lower_bound( - ::Type{TapRatioVariable}, - d::PSY.TwoWindingTransformer, - ::Type{VoltageControlTap}, -) - return _tap_ratio_limits(d).min -end - -function get_variable_upper_bound( - ::Type{TapRatioVariable}, - d::PSY.TwoWindingTransformer, - ::Type{VoltageControlTap}, -) - return _tap_ratio_limits(d).max -end - -# Warm-start the tap at its current position so IPOPT begins inside the bounds. -function get_variable_warm_start_value( - ::Type{TapRatioVariable}, - d::PSY.TwoWindingTransformer, - ::Type{VoltageControlTap}, -) - return PSY.get_tap(d) -end - -requires_initialization(::VoltageControlTap) = false - -function get_default_attributes( - ::Type{<:PSY.TwoWindingTransformer}, - ::Type{VoltageControlTap}, -) - return Dict{String, Any}( - PARALLEL_BRANCH_MAX_RATING_KEY => "single_element_contingency", - ) -end - -function get_default_time_series_names( - ::Type{<:PSY.TwoWindingTransformer}, - ::Type{VoltageControlTap}, -) - return Dict{Type{<:TimeSeriesParameter}, String}() -end - -################################################################################# -# Variable-tap AC π-model Ohm's law constraints. -# -# The ACP/ACR variable-tap Ohm's law shares its π-model coefficients and constraint -# builders with the fixed-tap StaticBranch path in AC_branches.jl -# (`_tap_flow_coefficients`, `_add_tap_acp_flow!`, `_add_tap_acr_flow!`), passing -# `TapRatioVariable[name, t]` as the tap in place of the constant `tm`, so each constraint -# reduces to its StaticBranch counterpart when `t == tm`. IVR keeps its own form below. -################################################################################# - -# ACP (polar) variable-tap Ohm's law. -function add_constraints!( - container::OptimizationContainer, - sys::PSY.System, - ::Type{NetworkFlowConstraint}, - devices::IS.FlattenIteratorWrapper{T}, - device_model::DeviceModel{T, VoltageControlTap}, - network_model::NetworkModel{ACPNetworkModel}, -) where {T <: PSY.TwoWindingTransformer} - time_steps = get_time_steps(container) - - va = get_variable(container, VoltageAngle, PSY.ACBus) - vm = get_variable(container, VoltageMagnitude, PSY.ACBus) - tap = get_variable(container, TapRatioVariable, T) - pft = get_variable(container, FlowActivePowerFromToVariable, T) - ptf = get_variable(container, FlowActivePowerToFromVariable, T) - qft = get_variable(container, FlowReactivePowerFromToVariable, T) - qtf = get_variable(container, FlowReactivePowerToFromVariable, T) - - number_to_name = _retained_number_to_name(sys, network_model) - geoms = - _branch_geometries(number_to_name, network_model, devices, T, NetworkFlowConstraint) - branch_names = [g.name for g in geoms] - cons_pft, cons_qft, cons_ptf, cons_qtf = - _add_flow_constraint_containers!(container, T, branch_names) - jump_model = get_jump_model(container) - slacks = _flow_equality_slacks(container, device_model, T) - - for g_geom in geoms - name = g_geom.name - adm = g_geom.adm - from_bus = g_geom.from_name - to_bus = g_geom.to_name - coef = _tap_flow_coefficients( - adm.g, adm.b, adm.g_fr, adm.b_fr, adm.g_to, adm.b_to, adm.shift, - ) - for t in time_steps - θ = va[from_bus, t] - va[to_bus, t] - vmf = vm[from_bus, t] - vmt = vm[to_bus, t] - _add_tap_acp_flow!( - jump_model, cons_pft, cons_qft, cons_ptf, cons_qtf, pft, qft, ptf, qtf, - name, t, vmf, vmt, θ, coef, tap[name, t], - _slack_term(slacks.p_ft, name, t), - _slack_term(slacks.q_ft, name, t), - _slack_term(slacks.p_tf, name, t), - _slack_term(slacks.q_tf, name, t), - ) - end - end - return -end - -# ACR (rectangular) variable-tap Ohm's law. -function add_constraints!( - container::OptimizationContainer, - sys::PSY.System, - ::Type{NetworkFlowConstraint}, - devices::IS.FlattenIteratorWrapper{T}, - device_model::DeviceModel{T, VoltageControlTap}, - network_model::NetworkModel{ACRNetworkModel}, -) where {T <: PSY.TwoWindingTransformer} - time_steps = get_time_steps(container) - - vr = get_variable(container, VoltageReal, PSY.ACBus) - vi = get_variable(container, VoltageImaginary, PSY.ACBus) - tap = get_variable(container, TapRatioVariable, T) - pft = get_variable(container, FlowActivePowerFromToVariable, T) - ptf = get_variable(container, FlowActivePowerToFromVariable, T) - qft = get_variable(container, FlowReactivePowerFromToVariable, T) - qtf = get_variable(container, FlowReactivePowerToFromVariable, T) - - number_to_name = _retained_number_to_name(sys, network_model) - geoms = - _branch_geometries(number_to_name, network_model, devices, T, NetworkFlowConstraint) - branch_names = [g.name for g in geoms] - cons_pft, cons_qft, cons_ptf, cons_qtf = - _add_flow_constraint_containers!(container, T, branch_names) - jump_model = get_jump_model(container) - slacks = _flow_equality_slacks(container, device_model, T) - - for g_geom in geoms - name = g_geom.name - adm = g_geom.adm - from_bus = g_geom.from_name - to_bus = g_geom.to_name - coef = _tap_flow_coefficients( - adm.g, adm.b, adm.g_fr, adm.b_fr, adm.g_to, adm.b_to, adm.shift, - ) - for t in time_steps - vr_fr = vr[from_bus, t] - vr_to = vr[to_bus, t] - vi_fr = vi[from_bus, t] - vi_to = vi[to_bus, t] - vv_fr = vr_fr^2 + vi_fr^2 - vv_to = vr_to^2 + vi_to^2 - cosprod = vr_fr * vr_to + vi_fr * vi_to - sinprod = vi_fr * vr_to - vr_fr * vi_to - _add_tap_acr_flow!( - jump_model, cons_pft, cons_qft, cons_ptf, cons_qtf, pft, qft, ptf, qtf, - name, t, vv_fr, vv_to, cosprod, sinprod, coef, tap[name, t], - _slack_term(slacks.p_ft, name, t), - _slack_term(slacks.q_ft, name, t), - _slack_term(slacks.p_tf, name, t), - _slack_term(slacks.q_tf, name, t), - ) - end - end - return -end - -# IVR (current-injection, rectangular) variable-tap Ohm's law. -# -# Mirrors the fixed-tap IVR branch constraints in AC_branches.jl term-by-term, with -# the constant tap `tm` (and the derived `tr = tm·cos(shift)`, `ti = tm·sin(shift)`, -# `tm² = tm^2`) replaced by the variable tap `t = TapRatioVariable[name, ts]`: -# tr → t·cs, ti → t·sn, tm² → t² (cs = cos(shift), sn = sin(shift)). -# The series impedance Z = r + jx is tap-independent (unchanged). Ten constraints -# per branch per time step (the same ten as the fixed-tap IVR branch). Because -# every `tm`-bearing term carries the live `t` symbol, each constraint reduces -# EXACTLY to its fixed-tap counterpart when `t == tap_nominal` (= PSY.get_tap(d) = -# adm.tap). The multiplied-through form (LHS·t²) keeps the equations polynomial -# (no division) and makes the reduction term-identical; t > 0 (TapRatioVariable -# bounds) guarantees equivalence with the divided form. -function add_constraints!( - container::OptimizationContainer, - sys::PSY.System, - ::Type{NetworkFlowConstraint}, - devices::IS.FlattenIteratorWrapper{T}, - device_model::DeviceModel{T, VoltageControlTap}, - network_model::NetworkModel{IVRNetworkModel}, -) where {T <: PSY.TwoWindingTransformer} - time_steps = get_time_steps(container) - - vr = get_variable(container, VoltageReal, PSY.ACBus) - vi = get_variable(container, VoltageImaginary, PSY.ACBus) - tap = get_variable(container, TapRatioVariable, T) - pft = get_variable(container, FlowActivePowerFromToVariable, T) - ptf = get_variable(container, FlowActivePowerToFromVariable, T) - qft = get_variable(container, FlowReactivePowerFromToVariable, T) - qtf = get_variable(container, FlowReactivePowerToFromVariable, T) - cr_fr = get_variable(container, BranchCurrentFromToReal, T) - ci_fr = get_variable(container, BranchCurrentFromToImaginary, T) - cr_to = get_variable(container, BranchCurrentToFromReal, T) - ci_to = get_variable(container, BranchCurrentToFromImaginary, T) - csr = get_variable(container, BranchSeriesCurrentReal, T) - csi = get_variable(container, BranchSeriesCurrentImaginary, T) - - number_to_name = _retained_number_to_name(sys, network_model) - geoms = - _branch_geometries(number_to_name, network_model, devices, T, NetworkFlowConstraint) - branch_names = [g.name for g in geoms] - - cons_pft = add_constraints_container!( - container, NetworkFlowConstraint, T, branch_names, time_steps; meta = "p_ft", - ) - cons_qft = add_constraints_container!( - container, NetworkFlowConstraint, T, branch_names, time_steps; meta = "q_ft", - ) - cons_ptf = add_constraints_container!( - container, NetworkFlowConstraint, T, branch_names, time_steps; meta = "p_tf", - ) - cons_qtf = add_constraints_container!( - container, NetworkFlowConstraint, T, branch_names, time_steps; meta = "q_tf", - ) - cons_cr_fr = add_constraints_container!( - container, NetworkFlowConstraint, T, branch_names, time_steps; meta = "cr_fr", - ) - cons_ci_fr = add_constraints_container!( - container, NetworkFlowConstraint, T, branch_names, time_steps; meta = "ci_fr", - ) - cons_cr_to = add_constraints_container!( - container, NetworkFlowConstraint, T, branch_names, time_steps; meta = "cr_to", - ) - cons_ci_to = add_constraints_container!( - container, NetworkFlowConstraint, T, branch_names, time_steps; meta = "ci_to", - ) - cons_vr_to = add_constraints_container!( - container, NetworkFlowConstraint, T, branch_names, time_steps; meta = "vr_to", - ) - cons_vi_to = add_constraints_container!( - container, NetworkFlowConstraint, T, branch_names, time_steps; meta = "vi_to", - ) - - jump_model = get_jump_model(container) - for g_geom in geoms - name = g_geom.name - adm = g_geom.adm - g = adm.g - b = adm.b - g_fr = adm.g_fr - b_fr = adm.b_fr - g_to = adm.g_to - b_to = adm.b_to - from_bus = g_geom.from_name - to_bus = g_geom.to_name - cs = cos(adm.shift) - sn = sin(adm.shift) - - # Series impedance Z = r + jx = conj(y)/|y|² (tap-independent). - ymag2 = g^2 + b^2 - r = g / ymag2 - x = -b / ymag2 - - for t in time_steps - vr_f = vr[from_bus, t] - vi_f = vi[from_bus, t] - vr_t = vr[to_bus, t] - vi_t = vi[to_bus, t] - tt = tap[name, t] - tt2 = tt^2 - tr = tt * cs - ti = tt * sn - csr_b = csr[name, t] - csi_b = csi[name, t] - cr_f = cr_fr[name, t] - ci_f = ci_fr[name, t] - cr_t = cr_to[name, t] - ci_t = ci_to[name, t] - - # Bilinear power-current linking (tap-independent) - cons_pft[name, t] = JuMP.@constraint( - jump_model, pft[name, t] == vr_f * cr_f + vi_f * ci_f, - ) - cons_qft[name, t] = JuMP.@constraint( - jump_model, qft[name, t] == vi_f * cr_f - vr_f * ci_f, - ) - cons_ptf[name, t] = JuMP.@constraint( - jump_model, ptf[name, t] == vr_t * cr_t + vi_t * ci_t, - ) - cons_qtf[name, t] = JuMP.@constraint( - jump_model, qtf[name, t] == vi_t * cr_t - vr_t * ci_t, - ) - - # KCL at from terminal (tm → t) - cons_cr_fr[name, t] = JuMP.@constraint( - jump_model, - cr_f * tt2 == - tr * csr_b - ti * csi_b + (g_fr * vr_f - b_fr * vi_f) * tt2, - ) - cons_ci_fr[name, t] = JuMP.@constraint( - jump_model, - ci_f * tt2 == - tr * csi_b + ti * csr_b + (g_fr * vi_f + b_fr * vr_f) * tt2, - ) - - # KCL at to terminal (no tap) - cons_cr_to[name, t] = JuMP.@constraint( - jump_model, cr_t == -csr_b + g_to * vr_t - b_to * vi_t, - ) - cons_ci_to[name, t] = JuMP.@constraint( - jump_model, ci_t == -csi_b + g_to * vi_t + b_to * vr_t, - ) - - # Ohm's law across series impedance (tm → t) - cons_vr_to[name, t] = JuMP.@constraint( - jump_model, - vr_t * tt2 == - vr_f * tr + vi_f * ti - r * csr_b * tt2 + x * csi_b * tt2, - ) - cons_vi_to[name, t] = JuMP.@constraint( - jump_model, - vi_t * tt2 == - vi_f * tr - vr_f * ti - r * csi_b * tt2 - x * csr_b * tt2, - ) - end - end - return -end - -################################################################################# -# Control-objective application — count-invariant JuMP.fix on existing variables. -# Branch on the enum value (data, not type). -################################################################################# - -# Shared handler for REACTIVE_POWER_FLOW and ACTIVE_POWER_FLOW objectives — identical -# between ACP and ACR. VOLTAGE regulation is handled per-network (ACP: direct vm fix; -# ACR/IVR: fix via RegulatedVoltageMagnitude aux variable). -function _fix_tap_flow_objective!( - d::PSY.TwoWindingTransformer, - name::String, - qft, - pft, - objective, - time_steps, -) - if objective == PSY.TransformerControlObjective.REACTIVE_POWER_FLOW - target = PSY.get_reactive_power_flow(d, PSY.SU) - for t in time_steps - JuMP.fix(qft[name, t], target; force = true) - end - elseif objective == PSY.TransformerControlObjective.ACTIVE_POWER_FLOW - target = PSY.get_active_power_flow(d, PSY.SU) - for t in time_steps - JuMP.fix(pft[name, t], target; force = true) - end - end - return -end - -# Resolve the regulated-bus name for a transformer: `regulated_bus_number == 0` -# means the arc's to-bus (local control). -function _tap_regulated_bus_name(d::PSY.TwoWindingTransformer, geom, number_to_name) - reg = PSY.get_regulated_bus_number(d) - if iszero(reg) - return geom.to_name - end - if !haskey(number_to_name, reg) - error( - "TwoWindingTransformer $(PSY.get_name(d)) regulates bus number $(reg), which is \ - not a retained bus — it does not exist or was absorbed by a network \ - reduction. Fix the regulated_bus_number or exclude the bus from the \ - reduction with a PNM reduction filter.", - ) - end - return number_to_name[reg] -end - -# Resolve the regulated ACBus for a transformer (used to bound and tie the ACR/IVR -# RegulatedVoltageMagnitude aux variable). `regulated_bus_number == 0` means the -# arc's to-bus (local control); otherwise the ACBus carrying that number. -function _tap_regulated_bus(d::PSY.TwoWindingTransformer, bus_by_number) - reg = PSY.get_regulated_bus_number(d) - if iszero(reg) - return PSY.get_to(PSY.get_arc(d)) - end - if !haskey(bus_by_number, reg) - error( - "TwoWindingTransformer $(PSY.get_name(d)) regulates bus number $(reg), which does \ - not exist in the system. Fix the regulated_bus_number.", - ) - end - return bus_by_number[reg] -end - -_regulated_buses(d::PSY.TwoWindingTransformer, bus_by_number) = - [("1", _tap_regulated_bus(d, bus_by_number))] - -# Dispatch entry: the VOLTAGE objective is pinned differently depending on how the -# network expresses a regulated bus voltage magnitude (polar scalar vs rectangular aux). -function _apply_tap_control_objective!( - container::OptimizationContainer, - sys::PSY.System, - devices::IS.FlattenIteratorWrapper{T}, - network_model::NetworkModel{N}, -) where {T <: PSY.TwoWindingTransformer, N} - return _apply_tap_control_objective!( - regulated_voltage_form(N), - container, - sys, - devices, - network_model, - ) -end - -# Polar (ACP): VOLTAGE pins the regulated-bus VoltageMagnitude directly; -# REACTIVE/ACTIVE_POWER_FLOW pin the from-to terminal flow. Other objectives -# (UNDEFINED / disabled) free-float. -function _apply_tap_control_objective!( - ::PolarRegulatedVoltage, - container::OptimizationContainer, - sys::PSY.System, - devices::IS.FlattenIteratorWrapper{T}, - network_model::NetworkModel, -) where {T <: PSY.TwoWindingTransformer} - time_steps = get_time_steps(container) - vm = get_variable(container, VoltageMagnitude, PSY.ACBus) - qft = get_variable(container, FlowReactivePowerFromToVariable, T) - pft = get_variable(container, FlowActivePowerFromToVariable, T) - number_to_name = _retained_number_to_name(sys, network_model) - # Control objectives act on the device's own terminals; the reduction guard in the - # ArgumentConstructStage ensures every device here is a direct (un-aggregated) entry. - for d in devices - geom = _branch_geometry(d) - name = geom.name - objective = PSY.get_control_objective(d) - if objective == PSY.TransformerControlObjective.VOLTAGE - reg_name = _tap_regulated_bus_name(d, geom, number_to_name) - setpoint = PSY.get_voltage_setpoint(d) - for t in time_steps - JuMP.fix(vm[reg_name, t], setpoint; force = true) - end - else - _fix_tap_flow_objective!(d, name, qft, pft, objective, time_steps) - end - end - return -end - -# Rectangular (ACR/IVR): VOLTAGE pins the regulated-bus magnitude via the component-owned -# (component, "1") RegulatedVoltageMagnitude aux variable (see fix_regulated_voltage!); -# reactive/active-flow objectives pin the from-to terminal flow. The aux variable/ -# constraint are added unconditionally in the construction stages, so only the fix is -# objective-conditional (count-invariance). Under IVR the from-to power variables -# (pft/qft) are bilinear-linked to the branch currents in the IVR Ohm's law, so the -# flow objectives are well-defined in current space — identical control logic to ACR. -function _apply_tap_control_objective!( - ::RectangularRegulatedVoltage, - container::OptimizationContainer, - sys::PSY.System, - devices::IS.FlattenIteratorWrapper{T}, - network_model::NetworkModel, -) where {T <: PSY.TwoWindingTransformer} - time_steps = get_time_steps(container) - qft = get_variable(container, FlowReactivePowerFromToVariable, T) - pft = get_variable(container, FlowActivePowerFromToVariable, T) - bus_by_number = _bus_by_number(sys) - for d in devices - name = PSY.get_name(d) - objective = PSY.get_control_objective(d) - if objective == PSY.TransformerControlObjective.VOLTAGE - reg_bus = _tap_regulated_bus(d, bus_by_number) - fix_regulated_voltage!( - container, d, "1", reg_bus, PSY.get_voltage_setpoint(d), network_model, - ) - else - _fix_tap_flow_objective!(d, name, qft, pft, objective, time_steps) - end - end - return -end - -################################################################################# -# construct_device! — two-stage. ACP/ACR build the branch in power only; -# IVR adds explicit branch current variables and a CurrentLimitConstraint. The -# `tap_branch_current_form` trait selects between the two construction paths. -################################################################################# - -function construct_device!( - container::OptimizationContainer, - sys::PSY.System, - stage::ArgumentConstructStage, - device_model::DeviceModel{T, VoltageControlTap}, - network_model::NetworkModel{N}, -) where { - T <: PSY.TwoWindingTransformer, - N <: Union{ACPNetworkModel, ACRNetworkModel, IVRNetworkModel}, -} - return construct_device!( - tap_branch_current_form(N), - container, - sys, - stage, - device_model, - network_model, - ) -end - -function construct_device!( - container::OptimizationContainer, - sys::PSY.System, - stage::ModelConstructStage, - device_model::DeviceModel{T, VoltageControlTap}, - network_model::NetworkModel{N}, -) where { - T <: PSY.TwoWindingTransformer, - N <: Union{ACPNetworkModel, ACRNetworkModel, IVRNetworkModel}, -} - return construct_device!( - tap_branch_current_form(N), - container, - sys, - stage, - device_model, - network_model, - ) -end - -# Power-only branch construction (ACP/ACR), mirrors StaticBranch. -function construct_device!( - ::PowerOnlyTapBranch, - container::OptimizationContainer, - sys::PSY.System, - ::ArgumentConstructStage, - device_model::DeviceModel{T, VoltageControlTap}, - network_model::NetworkModel, -) where {T <: PSY.TwoWindingTransformer} - @debug "construct_device VoltageControlTap (ArgumentConstructStage)" _group = - LOG_GROUP_BRANCH_CONSTRUCTIONS - devices = get_available_components(device_model, sys) - _validate_controlled_branch_not_reduced(network_model, devices, "VoltageControlTap") - add_variables!(container, TapRatioVariable, devices, VoltageControlTap) - add_variables!(container, FlowActivePowerFromToVariable, devices, VoltageControlTap) - add_variables!(container, FlowActivePowerToFromVariable, devices, VoltageControlTap) - add_variables!(container, FlowReactivePowerFromToVariable, devices, VoltageControlTap) - add_variables!(container, FlowReactivePowerToFromVariable, devices, VoltageControlTap) - add_regulated_voltage_magnitude!( - container, devices, sys, network_model, - ) - add_to_expression!( - container, ActivePowerBalance, FlowActivePowerFromToVariable, - devices, device_model, network_model, - ) - add_to_expression!( - container, ActivePowerBalance, FlowActivePowerToFromVariable, - devices, device_model, network_model, - ) - add_to_expression!( - container, ReactivePowerBalance, FlowReactivePowerFromToVariable, - devices, device_model, network_model, - ) - add_to_expression!( - container, ReactivePowerBalance, FlowReactivePowerToFromVariable, - devices, device_model, network_model, - ) - add_feedforward_arguments!(container, device_model, devices) - return -end - -function construct_device!( - ::PowerOnlyTapBranch, - container::OptimizationContainer, - sys::PSY.System, - ::ModelConstructStage, - device_model::DeviceModel{T, VoltageControlTap}, - network_model::NetworkModel{N}, -) where {T <: PSY.TwoWindingTransformer, N} - @debug "construct_device VoltageControlTap (ModelConstructStage)" _group = - LOG_GROUP_BRANCH_CONSTRUCTIONS - devices = get_available_components(device_model, sys) - add_constraints!( - container, FlowRateConstraintFromTo, devices, device_model, network_model, - ) - add_constraints!( - container, FlowRateConstraintToFrom, devices, device_model, network_model, - ) - add_constraints!( - container, sys, NetworkFlowConstraint, devices, device_model, network_model, - ) - add_constraints!( - container, sys, AngleDifferenceConstraint, devices, device_model, network_model, - ) - add_regulated_voltage_magnitude_constraints!( - container, devices, sys, network_model, - ) - _apply_tap_control_objective!(container, sys, devices, network_model) - add_feedforward_constraints!(container, device_model, devices) - add_to_objective_function!(container, devices, device_model, N) - add_constraint_dual!(container, sys, device_model) - return -end - -################################################################################# -# construct_device! — IVR (current-injection) variable-tap branch. -# Mirrors StaticBranch under IVRNetworkModel (branch_constructor.jl) plus the -# TapRatioVariable and the RegulatedVoltageMagnitude aux variable / constraint. -################################################################################# - -function construct_device!( - ::CurrentInjectionTapBranch, - container::OptimizationContainer, - sys::PSY.System, - ::ArgumentConstructStage, - device_model::DeviceModel{T, VoltageControlTap}, - network_model::NetworkModel, -) where {T <: PSY.TwoWindingTransformer} - @debug "construct_device IVR VoltageControlTap (ArgumentConstructStage)" _group = - LOG_GROUP_BRANCH_CONSTRUCTIONS - devices = get_available_components(device_model, sys) - _validate_controlled_branch_not_reduced(network_model, devices, "VoltageControlTap") - add_variables!(container, TapRatioVariable, devices, VoltageControlTap) - add_variables!(container, FlowActivePowerFromToVariable, devices, VoltageControlTap) - add_variables!(container, FlowActivePowerToFromVariable, devices, VoltageControlTap) - add_variables!(container, FlowReactivePowerFromToVariable, devices, VoltageControlTap) - add_variables!(container, FlowReactivePowerToFromVariable, devices, VoltageControlTap) - add_variables!(container, BranchCurrentFromToReal, devices, device_model, network_model) - add_variables!( - container, - BranchCurrentFromToImaginary, - devices, - device_model, - network_model, - ) - add_variables!(container, BranchCurrentToFromReal, devices, device_model, network_model) - add_variables!( - container, - BranchCurrentToFromImaginary, - devices, - device_model, - network_model, - ) - add_variables!(container, BranchSeriesCurrentReal, devices, device_model, network_model) - add_variables!( - container, - BranchSeriesCurrentImaginary, - devices, - device_model, - network_model, - ) - add_regulated_voltage_magnitude!( - container, devices, sys, network_model, - ) - add_to_expression!( - container, ActivePowerBalance, FlowActivePowerFromToVariable, - devices, device_model, network_model, - ) - add_to_expression!( - container, ActivePowerBalance, FlowActivePowerToFromVariable, - devices, device_model, network_model, - ) - add_to_expression!( - container, ReactivePowerBalance, FlowReactivePowerFromToVariable, - devices, device_model, network_model, - ) - add_to_expression!( - container, ReactivePowerBalance, FlowReactivePowerToFromVariable, - devices, device_model, network_model, - ) - add_feedforward_arguments!(container, device_model, devices) - return -end - -function construct_device!( - ::CurrentInjectionTapBranch, - container::OptimizationContainer, - sys::PSY.System, - ::ModelConstructStage, - device_model::DeviceModel{T, VoltageControlTap}, - network_model::NetworkModel{N}, -) where {T <: PSY.TwoWindingTransformer, N} - @debug "construct_device IVR VoltageControlTap (ModelConstructStage)" _group = - LOG_GROUP_BRANCH_CONSTRUCTIONS - devices = get_available_components(device_model, sys) - add_constraints!( - container, FlowRateConstraintFromTo, devices, device_model, network_model, - ) - add_constraints!( - container, FlowRateConstraintToFrom, devices, device_model, network_model, - ) - add_constraints!( - container, sys, NetworkFlowConstraint, devices, device_model, network_model, - ) - add_constraints!( - container, sys, CurrentLimitConstraint, devices, device_model, network_model, - ) - add_constraints!( - container, sys, AngleDifferenceConstraint, devices, device_model, network_model, - ) - add_regulated_voltage_magnitude_constraints!( - container, devices, sys, network_model, - ) - _apply_tap_control_objective!(container, sys, devices, network_model) - add_feedforward_constraints!(container, device_model, devices) - add_to_objective_function!(container, devices, device_model, N) - add_constraint_dual!(container, sys, device_model) - return -end - -# Defensive no-ops for active-power-only networks. template_validation drops the -# reactive VoltageControlTap formulation before construction, so these are only -# reached if template validation is bypassed. -function construct_device!( - ::OptimizationContainer, - ::PSY.System, - ::ArgumentConstructStage, - ::DeviceModel{T, VoltageControlTap}, - ::NetworkModel{<:AbstractActivePowerModel}, -) where {T <: PSY.TwoWindingTransformer} - return -end - -function construct_device!( - ::OptimizationContainer, - ::PSY.System, - ::ModelConstructStage, - ::DeviceModel{T, VoltageControlTap}, - ::NetworkModel{<:AbstractActivePowerModel}, -) where {T <: PSY.TwoWindingTransformer} - return -end diff --git a/test/test_native_transformer_tap.jl b/test/test_native_transformer_tap.jl index c7988ca..5f65143 100644 --- a/test/test_native_transformer_tap.jl +++ b/test/test_native_transformer_tap.jl @@ -2,8 +2,8 @@ # Off-nominal transformer tap under the native network models. # # These testsets cover only code that ships in the current module: the DC susceptance -# `b_dc = 1/(tap*x)` used by `BThetaBranchFlow`/`NetworkFlowConstraint`, and the tap-free -# π coefficients in `_tap_flow_coefficients` (both in `ac_transmission_models/ +# `b_dc = 1/(tap*x)` used by `BThetaBranchFlow`/`NetworkFlowConstraint`, and the Ybus +# two-port terms in `_pi_flow_terms` (both in `ac_transmission_models/ # AC_branches.jl`). They deliberately do NOT touch `TapControl` / `VoltageControlTap`, # whose formulation files are not yet included — those live in # `test_native_tapcontrol.jl` / `test_voltage_control_tap_models.jl` and stay disabled @@ -61,54 +61,33 @@ @test tested_a_real_tap end -@testset "_tap_flow_coefficients ground truth (hand-computed)" begin - # No shift: cs=1, sn=0. Hand-computed π terms and coupling coefficients. - c0 = POM._tap_flow_coefficients(1.0, -2.0, 0.1, 0.3, 0.2, 0.4, 0.0) - @test c0.cs == 1.0 - @test c0.sn == 0.0 - # From side is returned UNFOLDED (series and shunt separate) because the two get - # different tap treatment at the constraint site; see the composition asserts below. - @test c0.g == 1.0 - @test c0.b == -2.0 - @test c0.g_fr == 0.1 - @test c0.b_fr == 0.3 - # To side stays folded: neither half is tap-referred, matching PNM's `Y22 = Y_l + y_to`. - @test c0.gg_to == 1.2 - @test c0.bb_to == -1.6 - @test c0.a_cos == -1.0 # -g*cs + b*sn - @test c0.a_sin == 2.0 # -b*cs - g*sn - @test c0.c_cos == -1.0 # -g*cs - b*sn - @test c0.d_sin == -2.0 # b*cs - g*sn +@testset "_tapped_admittance round-trips PNM.ybus_branch_entries" begin + function check_terms(y, ybus) + Y11, Y12, Y21, Y22 = ybus + @test isapprox(complex(y.g11, y.b11), Y11; rtol = 1e-10, atol = 1e-12) + @test isapprox(complex(y.g12, y.b12), Y12; rtol = 1e-10, atol = 1e-12) + @test isapprox(complex(y.g21, y.b21), Y21; rtol = 1e-10, atol = 1e-12) + @test isapprox(complex(y.g22, y.b22), Y22; rtol = 1e-10, atol = 1e-12) + end - # From-side composition exactly as the ACP/ACR constraint bodies build it: - # g/tm^2 + g_fr (the series is tap-referred; the magnetizing shunt is NOT) - # This mirrors PNM's Ybus stamp `Y11 = Y_series/abs2(tap) + y_shunt_from`. The folded - # convention `(g + g_fr)/tm^2` would give 0.704 / -1.088 instead, which is what made - # POM's AC solutions disagree with PowerFlows for off-nominal taps. - tm = 1.25 # tm^2 == 1.5625, so 1/tm^2 == 0.64 exactly - @test c0.g / tm^2 + c0.g_fr ≈ 0.74 # 0.64 + 0.1 - @test c0.b / tm^2 + c0.b_fr ≈ -0.98 # -1.28 + 0.3 - # Continuity: at nominal tap the split form reproduces the folded value, so the - # convention only bites for off-nominal taps. - @test c0.g / 1.0^2 + c0.g_fr == 1.1 - @test c0.b / 1.0^2 + c0.b_fr == -1.7 + sys = PSB.build_system(PSITestSystems, "c_sys14") + for br in Iterators.flatten(( + PSY.get_components(PSY.Line, sys), + PSY.get_components(PSY.TwoWindingTransformer, sys), + )) + adm = PNM.branch_admittance(br) + check_terms(POM._pi_flow_terms(adm, adm.tap), PNM.ybus_branch_entries(br)) + end - # Nonzero shift = π/6: cs=√3/2, sn=1/2 exercises the trig. - cs = cos(pi / 6) - sn = sin(pi / 6) - cS = POM._tap_flow_coefficients(1.0, -2.0, 0.1, 0.3, 0.2, 0.4, pi / 6) - @test cS.cs ≈ cs - @test cS.sn ≈ sn - @test cS.g == 1.0 - @test cS.b == -2.0 - @test cS.g_fr == 0.1 - @test cS.b_fr == 0.3 - @test cS.gg_to == 1.2 - @test cS.bb_to == -1.6 - @test cS.a_cos ≈ -1.0 * cs + (-2.0) * sn - @test cS.a_sin ≈ -(-2.0) * cs - 1.0 * sn - @test cS.c_cos ≈ -1.0 * cs - (-2.0) * sn - @test cS.d_sin ≈ (-2.0) * cs - 1.0 * sn - # ACR's e_sin sign relationship the constraint body relies on. - @test -cS.d_sin ≈ -(-2.0) * cs + 1.0 * sn + tr = PSY.get_component(PSY.TwoWindingTransformer, sys, "Trans1") + circuit = PSY.get_circuit(tr) + for shift in (-pi / 5, 0.0, pi / 6) + PSY.set_α!(circuit, shift) + PSY.set_tap!(circuit, 1.0) + adm = PNM.branch_admittance(tr) + for tap in (0.9, 1.0, 1.1, 1.25) + PSY.set_tap!(circuit, tap) + check_terms(POM._pi_flow_terms(adm, tap), PNM.ybus_branch_entries(tr)) + end + end end From c4cfdaec18682110afc6134d995d0adc35ebb846 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Thu, 13 Aug 2026 18:03:31 -0400 Subject: [PATCH 04/11] refactor irreducible bus code; make controllable transformers irreducible --- .../instantiate_network_model.jl | 146 ++++++++++-------- test/test_native_transformer_tap.jl | 2 +- .../test_postcontingency_mixed_outage_axes.jl | 12 +- 3 files changed, 86 insertions(+), 74 deletions(-) diff --git a/src/network_models/instantiate_network_model.jl b/src/network_models/instantiate_network_model.jl index cec315a..0ac8002 100644 --- a/src/network_models/instantiate_network_model.jl +++ b/src/network_models/instantiate_network_model.jl @@ -74,25 +74,16 @@ _assign_subnetworks_to_buses( ::PSY.System, ) where {T <: AbstractNetworkModel} = nothing -function _push_component_buses!(buses::Set{Int64}, branch::PSY.Branch) +function _push_component_buses!( + buses::Set{Int64}, + branch::Union{PSY.Branch, PSY.TransformerCircuit}, +) arc = PSY.get_arc(branch) push!(buses, PSY.get_number(PSY.get_from(arc))) push!(buses, PSY.get_number(PSY.get_to(arc))) return end -function _push_component_buses!(buses::Set{Int64}, branch::PSY.ThreeWindingTransformer) - for arc in ( - PSY.get_primary_star_arc(branch), - PSY.get_secondary_star_arc(branch), - PSY.get_tertiary_star_arc(branch), - ) - push!(buses, PSY.get_number(PSY.get_from(arc))) - push!(buses, PSY.get_number(PSY.get_to(arc))) - end - return -end - function _push_component_buses!(buses::Set{Int64}, device::PSY.StaticInjection) push!(buses, PSY.get_number(PSY.get_bus(device))) return @@ -113,7 +104,7 @@ function _push_component_buses!(::Set{Int64}, ::PSY.AreaInterchange) end # Fallback for monitored/outaged component types with no bus-pinning rule. Reached -# from `_add_outage_monitored_irreducible_buses!`, which iterates the raw +# from `_outage_irreducible_buses`, which iterates the raw # `PSY.get_monitored_components(outage)` UUIDs (unfiltered — unlike the # template-validation path), so any non-{Branch, ThreeWindingTransformer, # StaticInjection, ACBus} monitored type lands here. Warn and skip rather than @@ -129,16 +120,28 @@ function _push_component_buses!(::Set{Int64}, c::PSY.Component) return end -# Outages registered on an outage-aware branch DeviceModel pin both their -# monitored and their outaged (associated) component buses so the network -# reduction can't collapse them: the MODF column for a contingency is keyed by -# the outaged arc's endpoints, and post-contingency flow constraints reference -# the monitored components' real bus numbers. -function _add_outage_monitored_irreducible_buses!( - irreducible_buses::Set{Int64}, +function _get_irreducible_buses( + sys::PSY.System, + network_model::NetworkModel, + branch_models::BranchModelContainer, +) + @debug "Identifying buses that are irreducible due to monitored components" + return collect( + union( + _dynamic_rating_irreducible_buses(branch_models), + _outage_irreducible_buses(sys, branch_models), + _monitored_lines_irreducible_buses(branch_models), + _controllable_transformers_irreducbile_buses(branch_models), + ), + ) +end + +# Pin both buses with outage-monitored and outaged lines +function _outage_irreducible_buses( sys::PSY.System, branch_models::BranchModelContainer, ) + irreducible_buses = Set{Int64}() outage_uuids = Set{Base.UUID}() for m in values(branch_models) IOM.supports_outages(get_formulation(m)) || continue @@ -162,78 +165,89 @@ function _add_outage_monitored_irreducible_buses!( _push_component_buses!(irreducible_buses, component) end end - return + return irreducible_buses end -# Buses that must survive PNM network reductions because something monitored is -# pinned to them: branch endpoints carrying a `BranchRatingTimeSeriesParameter` -# (dynamic line ratings), and the monitored/outaged endpoints of outages -# registered on outage-aware (security-constrained) branch DeviceModels. -function _get_irreducible_buses_due_to_monitored_components( - sys::PSY.System, - network_model::NetworkModel, +_is_ac_transmission(::DeviceModel{PSY.ACTransmission}) = true +_is_ac_transmission(_) = false + +_is_3w(::DeviceModel{PSY.ThreeWindingTransformer}) = true +_is_3w(_) = false + +# Pin buses with DLR branches +function _dynamic_rating_irreducible_buses( branch_models::BranchModelContainer, ) - @debug "Identifying buses that are irreducible due to monitored components" irreducible_buses = Set{Int64}() - for branch_type in network_model.modeled_branch_types - branch_type <: PSY.ACTransmission || continue - device_model = branch_models[nameof(branch_type)] + for model in values(branch_models) + _is_ac_transmission(model) || continue if !haskey( - get_time_series_names(device_model), + get_time_series_names(model), BranchRatingTimeSeriesParameter, ) continue end - - if branch_type == PSY.ThreeWindingTransformer - @warn "Dynamic branch ratings for ThreeWindingTransformers are not implemented yet. Skipping it." + if _is_3w(model) + @warn "Dynamic branch ratings for ThreeWindingTransformers are not implemented yet. Allowing these devices to be reduced." continue end ts_name = - get_time_series_names(device_model)[BranchRatingTimeSeriesParameter] + get_time_series_names(model)[BranchRatingTimeSeriesParameter] ts_type = PSY.Deterministic #TODO workaround since we dont have the container - branches = PSY.get_available_components(branch_type, sys) - for branch in branches + for branch in get_device_cache(model) if !PSY.has_time_series(branch, ts_type, ts_name) continue end _push_component_buses!(irreducible_buses, branch) end end - _add_outage_monitored_irreducible_buses!(irreducible_buses, sys, branch_models) - # `model_all_branches` MonitoredLine models pin their lines so zero-impedance - # ones survive the reduction instead of being merged away. - _add_model_all_branches_irreducible_buses!(irreducible_buses, branch_models) - return collect(irreducible_buses) + return irreducible_buses end -# Pin both endpoint buses of every branch a `model_all_branches` MonitoredLine model -# covers. Dispatch on the model type so it is a no-op for other branch types. -function _add_model_all_branches_irreducible_buses!( - irreducible_buses::Set{Int64}, - branch_models::BranchModelContainer, -) +_is_monitored(m::DeviceModel{PSY.MonitoredLine}) = + get_attribute(m, MODEL_ALL_BRANCHES_KEY) === true +_is_monitored(_) = false + +# Pin buses with MonitoredLines with MODEL_ALL_BRANCHES_KEY === true +function _monitored_lines_irreducible_buses(branch_models::BranchModelContainer) + irreducible_buses = Set{Int64}() for m in values(branch_models) - _pin_model_all_branches!(irreducible_buses, m) + _is_monitored(m) || continue + for branch in get_device_cache(m) + _push_component_buses!(irreducible_buses, branch) + end end - return + return irreducible_buses end -_pin_model_all_branches!(::Set{Int64}, ::DeviceModel) = nothing +_is_controllable( + m::Union{ + DeviceModel{PSY.TwoWindingTransformer}, + DeviceModel{PSY.ThreeWindingTransformer}, + }, +) = get_attribute(m, ENABLE_CONTROLS_KEY) === true +_is_controllable(c::PSY.TransformerCircuit) = + !( + get_control_objective(c) in + (TransformerControlObjective.UNDEFINED, TransformerControlObjective.FIXED) + ) +_is_controllable(_) = false -function _pin_model_all_branches!( - irreducible_buses::Set{Int64}, - m::DeviceModel{PSY.MonitoredLine}, -) - get_attribute(m, MODEL_ALL_BRANCHES_KEY) === true || return - # The device cache is the modeled set (available + filter_function). - for branch in get_device_cache(m) - _push_component_buses!(irreducible_buses, branch) +# Pin buses with transformers with non-default control +function _controllable_transformers_irreducbile_buses(branch_models::BranchModelContainer) + irreducible_buses = Set{Int64}() + for m in values(branch_models) + _is_controllable(m) || continue + for transformer in get_device_cache(m) + for circuit in PSY.get_circuits(transformer) + _is_controllable(m) || continue + _push_component_buses!(irreducible_buses, branch) + end + end end - return + return irreducible_buses end # Drop (and warn about) any branch type whose components were all merged away by the @@ -357,7 +371,7 @@ function IOM.instantiate_network_model!( sys::PSY.System, ) where {T <: AbstractNetworkModel} _validate_network_and_branches(model, branch_models, sys) - irreducible_buses = _get_irreducible_buses_due_to_monitored_components( + irreducible_buses = _get_irreducible_buses( sys, model, branch_models, @@ -378,7 +392,7 @@ function IOM.instantiate_network_model!( sys::PSY.System, ) _validate_network_and_branches(model, branch_models, sys) - irreducible_buses = _get_irreducible_buses_due_to_monitored_components( + irreducible_buses = _get_irreducible_buses( sys, model, branch_models, @@ -518,7 +532,7 @@ function IOM.instantiate_network_model!( number_of_steps::Int, sys::PSY.System, ) - irreducible_buses = _get_irreducible_buses_due_to_monitored_components( + irreducible_buses = _get_irreducible_buses( sys, model, branch_models, diff --git a/test/test_native_transformer_tap.jl b/test/test_native_transformer_tap.jl index 5f65143..0d7a198 100644 --- a/test/test_native_transformer_tap.jl +++ b/test/test_native_transformer_tap.jl @@ -3,7 +3,7 @@ # # These testsets cover only code that ships in the current module: the DC susceptance # `b_dc = 1/(tap*x)` used by `BThetaBranchFlow`/`NetworkFlowConstraint`, and the Ybus -# two-port terms in `_pi_flow_terms` (both in `ac_transmission_models/ +# two-port terms in `_tapped_admittance` (both in `ac_transmission_models/ # AC_branches.jl`). They deliberately do NOT touch `TapControl` / `VoltageControlTap`, # whose formulation files are not yet included — those live in # `test_native_tapcontrol.jl` / `test_voltage_control_tap_models.jl` and stay disabled diff --git a/test/test_postcontingency_mixed_outage_axes.jl b/test/test_postcontingency_mixed_outage_axes.jl index f2647d6..5af24fd 100644 --- a/test/test_postcontingency_mixed_outage_axes.jl +++ b/test/test_postcontingency_mixed_outage_axes.jl @@ -1,7 +1,7 @@ # Ported from PowerSimulations.jl PR #1579 (MODF SCUC migration), adapted to # POM / PS6 APIs. Exercises the Phase-4 `_build_device_model_outages!` planned # vs. unplanned outage-axis selection and the -# `_add_outage_monitored_irreducible_buses!` bus-pinning (N3) logic. +# `_outage_irreducible_buses` bus-pinning (N3) logic. # # Adaptation notes: # * internal symbols are namespaced `POM.` / `IOM.` / `PNM.`; expression and @@ -132,7 +132,7 @@ end @testset "Outage pinning includes outaged-component buses (N3)" begin - # Regression: `_add_outage_monitored_irreducible_buses!` must pin both the + # Regression: `_outage_irreducible_buses` must pin both the # MONITORED components' buses AND the OUTAGED (associated) components' # buses. If only the monitored set is pinned, a degree-two reduction # between the outaged arc's endpoints can collapse the contingency arc out @@ -163,8 +163,7 @@ end ) branch_models[nameof(PSY.Line)] = dm - irreducible_buses = Set{Int64}() - POM._add_outage_monitored_irreducible_buses!(irreducible_buses, sys, branch_models) + irreducible_buses = POM._outage_irreducible_buses(sys, branch_models) monitored_arc = PSY.get_arc(monitored_only_line) outaged_arc = PSY.get_arc(outaged_line) @@ -199,11 +198,10 @@ end branch_models = IOM.BranchModelContainer() dm = DeviceModel(PSY.Line, POM.StaticBranch) # Even if an outage dict were present, StaticBranch is not outage-aware so - # `_add_outage_monitored_irreducible_buses!` skips it. + # `_outage_irreducible_buses` skips it. branch_models[nameof(PSY.Line)] = dm - irreducible_buses = Set{Int64}() - POM._add_outage_monitored_irreducible_buses!(irreducible_buses, sys, branch_models) + irreducible_buses = POM._outage_irreducible_buses(sys, branch_models) @test isempty(irreducible_buses) end From d95b43837bec6b7c1b3d91a171531b087de3f64b Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Thu, 13 Aug 2026 23:57:18 -0400 Subject: [PATCH 05/11] rough tap control implementation; still need to run tests and add voltage and q bounds --- src/ac_transmission_models/AC_branches.jl | 112 ++++++++++++------ .../branch_constructor.jl | 1 + 2 files changed, 75 insertions(+), 38 deletions(-) diff --git a/src/ac_transmission_models/AC_branches.jl b/src/ac_transmission_models/AC_branches.jl index c71f293..c9288fb 100644 --- a/src/ac_transmission_models/AC_branches.jl +++ b/src/ac_transmission_models/AC_branches.jl @@ -69,6 +69,14 @@ function get_default_time_series_names( return Dict{Type{<:TimeSeriesParameter}, String}() end +ENABLE_CONTROLS_KEY = "enable_controls" + +_control_attribute( + ::Union{Type{PSY.TwoWindingTransformer}, Type{PSY.ThreeWindingTransformer}}, +) = (ENABLE_CONTROLS_KEY => false,) + +_control_attribute(_) = () + """ DeviceModel attribute key selecting which `PowerNetworkMatrices` function aggregates the individual circuit ratings of a `PNM.BranchesParallel` into a single maximum flow @@ -84,6 +92,7 @@ function get_default_attributes( ) where {U <: PSY.ACTransmission, V <: AbstractBranchFormulation} return Dict{String, Any}( PARALLEL_BRANCH_MAX_RATING_KEY => "single_element_contingency", + _control_attribute(U)... ) end @@ -94,6 +103,7 @@ function get_default_attributes( return Dict{String, Any}( PARALLEL_BRANCH_MAX_RATING_KEY => "single_element_contingency", "include_planned_outages" => false, + _control_attribute(U)... ) end @@ -271,12 +281,48 @@ function add_variables!( return end -# Non-negative flow-definition slack container carrying a container META. StaticBranchBounds -# distinguishes its per-direction slack pairs ("p_ft"/"p_tf"/"q_ft"/"q_tf") by meta on the -# shared FlowActivePowerSlack{Upper,Lower}Bound types; `add_variables!` threads no meta, so -# build the container directly. One slack per representative arc — the equality is written -# once per arc. Axes are precomputed by the caller (shared across all metas of one device -# model). +# Matches the names returned by _branch_geometries +_circuit_arc_name(d::PSY.TwoWindingTransformer, ::PSY.TransformerCircuit, ::Int) = + PSY.get_name(d) +_circuit_arc_name(d::PSY.ThreeWindingTransformer, c::PSY.TransformerCircuit, i::Int) = + PNM.get_name(PNM.ThreeWindingTransformerCircuit(d, c, i)) + +function _add_transformer_control_variables!( + container::OptimizationContainer, + model::DeviceModel{U, F}, + devices::IS.FlattenIteratorWrapper{U}, + network_model::NetworkModel, +) where { + U <: Union{PSY.TwoWindingTransformer, PSY.ThreeWindingTransformer}, + F <: AbstractBranchFormulation, +} + get_attribute(model, ENABLE_CONTROLS_KEY) === true || return + time_steps = get_time_steps(container) + jump_model = get_jump_model(container) + names = [ + _circuit_arc_name(d, c, i) + for d in devices + for (i, c) in enumerate(PSY.get_circuits(d)) + if PSY.get_available(c) && PSY.get_control_objective(c) in (PSY.TransformerControlObjective.VOLTAGE, PSY.TransformerControlObjective.REACTIVE_POWER_FLOW) + ] + _validate_controlled_branch_not_reduced(network_model, U, names) + variable = add_variable_container!(container, TapRatioVariable, U, names, time_steps) + for name in names, t in time_steps + variable[name, t] = JuMP.@variable( + jump_model, + base_name = "TapRatioVariable_$(U)_{$(name), $(t)}", + ) + end + return +end + +_add_transformer_control_variables!( + ::OptimizationContainer, + ::DeviceModel, + ::IS.FlattenIteratorWrapper, + ::NetworkModel, +) = nothing + function _add_meta_flow_slack!( container::OptimizationContainer, ::Type{T}, @@ -1113,27 +1159,19 @@ function _branch_rating_entries( ] end -# Formulations that model a per-device control decision variable (variable tap ratio, -# phase-shifter angle) cannot be expressed on a PNM series/parallel equivalent — the -# reduction folds a FIXED device setting into the merged π-parameters. A controlled -# branch absorbed by a network reduction is a modeling conflict the user must resolve, -# not something to silently approximate. function _validate_controlled_branch_not_reduced( network_model::NetworkModel, - devices::IS.FlattenIteratorWrapper{T}, - formulation_name::String, + ::Type{T}, + controlled_names, ) where {T <: PSY.ACTransmission} network_reduction = get_network_reduction(network_model) isempty(network_reduction) && return arc_map = get_name_to_arc_map_entries(network_reduction, T) - for d in devices - name = PSY.get_name(d) - if !haskey(arc_map, name) || arc_map[name][2] != "direct_branch_map" + for name in controlled_names + entry = get(arc_map, name, nothing) + if entry === nothing || entry[2] != "direct_branch_map" error( - "$(formulation_name) branch $(name) was absorbed by a network \ - reduction (radial, degree-two, or parallel aggregation). Exclude it \ - from the reduction with a PNM reduction filter or model it with a \ - static branch formulation.", + "Controlled transformer circuit $(name) was merged with a parallel branch. Either remove the parallel branch or disable control for this circuit.", ) end end @@ -1172,6 +1210,9 @@ _dc_phase_shift(branch::PNM.AbstractReductionAggregate, nr::PNM.NetworkReduction _dc_phase_shift(branch::PSY.ACTransmission, ::PNM.NetworkReductionData) = PNM.get_series_phase_shift(branch) +_control_objective(c::PSY.TransformerCircuit) = PSY.get_control_objective(c) +_control_objective(_) = PSY.TransformerControlObjective.UNDEFINED + function _branch_geometry( nr::PNM.NetworkReductionData, number_to_name::Dict{Int, String}, @@ -1192,6 +1233,7 @@ function _branch_geometry( shift_dc = _dc_phase_shift(branch, nr), r_dc = PNM.arc_dc_resistance(nr, arc_tuple), direct = !_is_aggregate(branch), + objective = _control_objective(branch) ) end @@ -1382,14 +1424,6 @@ function add_constraints!( return end -""" -Create the four directional `NetworkFlowConstraint` containers shared by every AC -branch-flow formulation, fixed- and variable-tap alike: active and reactive power in -the from→to and to→from directions, keyed by branch name and time step. Thin factory -over `add_constraints_container!`; returns them in (p_ft, q_ft, p_tf, q_tf) order so a -caller can write `cons_pft, cons_qft, cons_ptf, cons_qtf = ...`. Keeping this in one -place lets each formulation's method show only the Ohm's-law math that actually differs. -""" function _add_flow_constraint_containers!( container::OptimizationContainer, ::Type{T}, @@ -1656,6 +1690,8 @@ function add_constraints!( qft = get_variable(container, FlowReactivePowerFromToVariable, T) qtf = get_variable(container, FlowReactivePowerToFromVariable, T) + tap_var = get_attribute(device_model, ENABLE_CONTROLS_KEY) ? get_variable(container, TapRatioVariable, T) : nothing + number_to_name = _retained_number_to_name(sys, network_model) geoms = _branch_geometries(number_to_name, network_model, devices, T, NetworkFlowConstraint) @@ -1663,12 +1699,7 @@ function add_constraints!( cons_pft, cons_qft, cons_ptf, cons_qtf = _add_flow_constraint_containers!(container, T, branch_names) jump_model = get_jump_model(container) - slacks = _flow_equality_slacks(container, device_model, T) - p_ft_slack = _slack_term(slacks.p_ft, name, t) - q_ft_slack = _slack_term(slacks.q_ft, name, t) - p_tf_slack = _slack_term(slacks.p_tf, name, t) - q_tf_slack = _slack_term(slacks.q_tf, name, t) for g_geom in geoms name = g_geom.name @@ -1678,27 +1709,32 @@ function add_constraints!( for t in time_steps vp = _voltage_products(container, network_model, T, name, from_bus, to_bus, t) - y = _tapped_admittance(adm, adm.tap) + tap = g_geom.control in TAP_CONTROL_OBJECTIVES ? tap_var[name, t] : adm.tap + y = _tapped_admittance(adm, tap) cons_pft[name, t] = JuMP.@constraint( jump_model, pft[name, t] == - y.g11 * vp.v2_fr + y.g12 * vp.vv_cos + y.b12 * vp.vv_sin + p_ft_slack, + y.g11 * vp.v2_fr + y.g12 * vp.vv_cos + y.b12 * vp.vv_sin + + p_slack_term(slacks.p_ft, name, t) ) cons_qft[name, t] = JuMP.@constraint( jump_model, qft[name, t] == - -y.b11 * vp.v2_fr - y.b12 * vp.vv_cos + y.g12 * vp.vv_sin + q_ft_slack, + -y.b11 * vp.v2_fr - y.b12 * vp.vv_cos + y.g12 * vp.vv_sin + + _slack_term(slacks.q_ft, name, t), ) cons_ptf[name, t] = JuMP.@constraint( jump_model, ptf[name, t] == - y.g22 * vp.v2_to + y.g21 * vp.vv_cos - y.b21 * vp.vv_sin + p_tf_slack, + y.g22 * vp.v2_to + y.g21 * vp.vv_cos - y.b21 * vp.vv_sin + + _slack_term(slacks.p_tf, name, t), ) cons_qtf[name, t] = JuMP.@constraint( jump_model, qtf[name, t] == - -y.b22 * vp.v2_to - y.b21 * vp.vv_cos - y.g21 * vp.vv_sin + q_tf_slack, + -y.b22 * vp.v2_to - y.b21 * vp.vv_cos - y.g21 * vp.vv_sin + + _slack_term(slacks.q_tf, name, t), ) end end diff --git a/src/ac_transmission_models/branch_constructor.jl b/src/ac_transmission_models/branch_constructor.jl index 270db4e..018b49e 100644 --- a/src/ac_transmission_models/branch_constructor.jl +++ b/src/ac_transmission_models/branch_constructor.jl @@ -288,6 +288,7 @@ function construct_device!( devices = get_available_components(device_model, sys) _add_static_branch_flow_variables!(container, devices, network_model) _add_static_branch_balance_arguments!(container, device_model, devices, network_model) + _add_transformer_control_variables!(container, device_model, devices, network_model) return end From 90b3a4b0d30b3de6abdf7f27a04c05bbcbeca32b Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Fri, 14 Aug 2026 09:52:45 -0400 Subject: [PATCH 06/11] finished tap controls --- src/ac_transmission_models/AC_branches.jl | 62 ++++++++++++------- .../branch_constructor.jl | 1 + 2 files changed, 42 insertions(+), 21 deletions(-) diff --git a/src/ac_transmission_models/AC_branches.jl b/src/ac_transmission_models/AC_branches.jl index c9288fb..106f402 100644 --- a/src/ac_transmission_models/AC_branches.jl +++ b/src/ac_transmission_models/AC_branches.jl @@ -287,6 +287,8 @@ _circuit_arc_name(d::PSY.TwoWindingTransformer, ::PSY.TransformerCircuit, ::Int) _circuit_arc_name(d::PSY.ThreeWindingTransformer, c::PSY.TransformerCircuit, i::Int) = PNM.get_name(PNM.ThreeWindingTransformerCircuit(d, c, i)) +# TODO: add other controls + refactor +# TODO: Change TransformerCircuit <: DeviceParameter -> <: Device function _add_transformer_control_variables!( container::OptimizationContainer, model::DeviceModel{U, F}, @@ -297,32 +299,30 @@ function _add_transformer_control_variables!( F <: AbstractBranchFormulation, } get_attribute(model, ENABLE_CONTROLS_KEY) === true || return - time_steps = get_time_steps(container) - jump_model = get_jump_model(container) - names = [ - _circuit_arc_name(d, c, i) + # TODO: collectzip is ugly, refactor now or do Circuit <: Device first + names, circuits = collect(zip([ + (_circuit_arc_name(d, c, i), c) for d in devices for (i, c) in enumerate(PSY.get_circuits(d)) if PSY.get_available(c) && PSY.get_control_objective(c) in (PSY.TransformerControlObjective.VOLTAGE, PSY.TransformerControlObjective.REACTIVE_POWER_FLOW) - ] + ])) _validate_controlled_branch_not_reduced(network_model, U, names) + + time_steps = get_time_steps(container) + jump_model = get_jump_model(container) variable = add_variable_container!(container, TapRatioVariable, U, names, time_steps) - for name in names, t in time_steps + for (i, name) in enumerate(names), t in time_steps + bounds = get_control_limits(circuits[i]) variable[name, t] = JuMP.@variable( jump_model, base_name = "TapRatioVariable_$(U)_{$(name), $(t)}", + lower_bound = bounds.min, + upper_bound = bounds.max ) end return end -_add_transformer_control_variables!( - ::OptimizationContainer, - ::DeviceModel, - ::IS.FlattenIteratorWrapper, - ::NetworkModel, -) = nothing - function _add_meta_flow_slack!( container::OptimizationContainer, ::Type{T}, @@ -1213,6 +1213,12 @@ _dc_phase_shift(branch::PSY.ACTransmission, ::PNM.NetworkReductionData) = _control_objective(c::PSY.TransformerCircuit) = PSY.get_control_objective(c) _control_objective(_) = PSY.TransformerControlObjective.UNDEFINED +_quantity_limits(c::PSY.TransformerCircuit) = PSY.get_controlled_quantity_limits(c) +_quantity_limits(_) = (min = -Inf, max = Inf) + +_regulated_bus(c::PSY.TransformerCircuit) = PSY.get_regulated_bus_number(c) === PSY.get_number(PSY.get_from(PSY.get_arc(c))) ? :from : :to +_regulated_bus(_) = nothing + function _branch_geometry( nr::PNM.NetworkReductionData, number_to_name::Dict{Int, String}, @@ -1234,6 +1240,8 @@ function _branch_geometry( r_dc = PNM.arc_dc_resistance(nr, arc_tuple), direct = !_is_aggregate(branch), objective = _control_objective(branch) + quantity_limits = _quantity_limits(branch) + regulated_bus = _regulated_bus(branch) ) end @@ -1718,24 +1726,36 @@ function add_constraints!( y.g11 * vp.v2_fr + y.g12 * vp.vv_cos + y.b12 * vp.vv_sin + p_slack_term(slacks.p_ft, name, t) ) - cons_qft[name, t] = JuMP.@constraint( - jump_model, - qft[name, t] == - -y.b11 * vp.v2_fr - y.b12 * vp.vv_cos + y.g12 * vp.vv_sin + - _slack_term(slacks.q_ft, name, t), - ) cons_ptf[name, t] = JuMP.@constraint( jump_model, ptf[name, t] == y.g22 * vp.v2_to + y.g21 * vp.vv_cos - y.b21 * vp.vv_sin + _slack_term(slacks.p_tf, name, t), ) - cons_qtf[name, t] = JuMP.@constraint( + qft_expr = JuMP.@expression( + jump_model, + -y.b11 * vp.v2_fr - y.b12 * vp.vv_cos + y.g12 * vp.vv_sin + + _slack_term(slacks.q_ft, name, t), + ) + qtf_expr = JuMP.@expression( jump_model, - qtf[name, t] == -y.b22 * vp.v2_to - y.b21 * vp.vv_cos - y.g21 * vp.vv_sin + _slack_term(slacks.q_tf, name, t), ) + cons_qft[name, t] = JuMP.@constraint(jump_model, qft[name, t] == qft_expr) + cons_qtf[name, t] = JuMP.@constraint(jump_model, qtf[name, t] == qtf_expr) + + if g_geom.control === PSY.TransformerControlObjective.VOLTAGE + voltage = g_geom.regulated_bus == :from ? vp.v2_fr : vp.v2_to + JuMP.@constraint(jump_model, voltage >= g_geom.quantity_limits.min) + JuMP.@constraint(jump_model, voltage <= g_geom.quantity_limits.max) + end + if g_geom.control === PSY.TransformerControlObjective.REACTIVE_POWER_FLOW + JuMP.@constraint(jump_model, qft_expr >= g_geom.quantity_limits.min) + JuMP.@constraint(jump_model, qft_expr <= g_geom.quantity_limits.max) + JuMP.@constraint(jump_model, qtf_expr >= g_geom.quantity_limits.min) + JuMP.@constraint(jump_model, qtf_expr <= g_geom.quantity_limits.max) + end end end return diff --git a/src/ac_transmission_models/branch_constructor.jl b/src/ac_transmission_models/branch_constructor.jl index 018b49e..39326da 100644 --- a/src/ac_transmission_models/branch_constructor.jl +++ b/src/ac_transmission_models/branch_constructor.jl @@ -320,6 +320,7 @@ function construct_device!( add_constraints!( container, sys, AngleDifferenceConstraint, devices, device_model, network_model, ) + _add_transformer_control_constraints!(container, device_model, devices, network_model) add_feedforward_constraints!(container, device_model, devices) add_to_objective_function!(container, devices, device_model, ACPNetworkModel) add_constraint_dual!(container, sys, device_model) From d6d631b58d26942dfae22d1fa02f30d2b3f3f56d Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Fri, 14 Aug 2026 11:00:09 -0400 Subject: [PATCH 07/11] claude reorganize tests and update tap tests --- test/runtests.jl | 9 +- test/test_device_branch_constructors.jl | 57 -- test/test_native_dcp_acp_models.jl | 37 +- test/test_native_lpacc_model.jl | 55 +- test/test_native_network_reductions.jl | 59 -- test/test_native_tapcontrol.jl | 72 -- test/test_power_flow_in_the_loop.jl | 60 -- test/test_transformer_controls.jl | 707 ++++++++++++++++++ ...r_tap.jl => test_transformer_fixed_tap.jl} | 15 +- test/test_voltage_control_tap_models.jl | 405 ---------- 10 files changed, 752 insertions(+), 724 deletions(-) delete mode 100644 test/test_native_tapcontrol.jl create mode 100644 test/test_transformer_controls.jl rename test/{test_native_transformer_tap.jl => test_transformer_fixed_tap.jl} (86%) delete mode 100644 test/test_voltage_control_tap_models.jl diff --git a/test/runtests.jl b/test/runtests.jl index b092cde..25c59ad 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -20,14 +20,7 @@ const TEST_DIR = @__DIR__ # helpers, and `test_data/` are shared infrastructure, not standalone testsets — they # must not be run as tests (ParallelTestRunner's default discovery would pick them up). -# psy6: disabled pending the transformer refactor. Both remaining entries dispatch on a -# formulation that is still commented out of the module — `TapControl` in -# `transformer_models.jl` and `VoltageControlTap` in `voltage_control_tap_models.jl` — so -# their bodies cannot even be compiled yet. Re-enable them with those formulations. -const DISABLED_TESTS = Set([ - "test_native_tapcontrol", - "test_voltage_control_tap_models", -]) +const DISABLED_TESTS = Set(String[]) testsuite = Dict{String, Expr}( splitext(f)[1] => :(include($(joinpath(TEST_DIR, f)))) for diff --git a/test/test_device_branch_constructors.jl b/test/test_device_branch_constructors.jl index 565f149..8b7c875 100644 --- a/test/test_device_branch_constructors.jl +++ b/test/test_device_branch_constructors.jl @@ -526,63 +526,6 @@ end ) end -@testset "DC Power Flow Models for phase-shifting TwoWindingTransformer and Line" begin - # system = build_system(PSITestSystems, "c_sys5_uc") - # - # line = get_component(Line, system, "1") - # - # ps = TwoWindingTransformer(; - # name = get_name(line), - # available = true, - # active_power_flow = 0.0, - # reactive_power_flow = 0.0, - # r = get_r(line, PSY.SU), - # x = get_r(line, PSY.SU), - # primary_shunt = 0.0, - # tap = 1.0, - # α = 0.0, - # rating = get_rating(line, PSY.SU), - # arc = get_arc(line), - # base_power = get_base_power(system, PSY.NU), - # ) - # - # add_component!(system, ps) - # remove_component!(system, line) - # - # template = get_template_dispatch_with_network( - # NetworkModel(PTDFNetworkModel; network_matrix = PTDF(system)), - # ) - # set_device_model!(template, DeviceModel(TwoWindingTransformer, PhaseAngleControl)) - # model_m = DecisionModel(template, system; optimizer = HiGHS_optimizer) - # @test build!(model_m; output_dir = mktempdir(; cleanup = true)) == - # IOM.ModelBuildStatus.BUILT - # - # @test check_variable_unbounded( - # model_m, - # FlowActivePowerVariable, - # TwoWindingTransformer, - # ) - # - # @test solve!(model_m) == IOM.RunStatus.SUCCESSFULLY_FINALIZED - # - # @test check_flow_variable_values( - # model_m, - # FlowActivePowerVariable, - # TwoWindingTransformer, - # "1", - # get_rating(ps, PSY.SU), - # ) - # - # @test check_flow_variable_values( - # model_m, - # PhaseShifterAngle, - # TwoWindingTransformer, - # "1", - # -π / 2, - # π / 2, - # ) -end - @testset "AC Power Flow Models for TwoTerminalGenericHVDCLine Flow Constraints and TwoWindingTransformer Unbounded" begin ratelimit_constraint_keys = [ IOM.ConstraintKey(FlowRateConstraintFromTo, TwoWindingTransformer), diff --git a/test/test_native_dcp_acp_models.jl b/test/test_native_dcp_acp_models.jl index 4335e15..7dd5f41 100644 --- a/test/test_native_dcp_acp_models.jl +++ b/test/test_native_dcp_acp_models.jl @@ -479,29 +479,20 @@ end end @testset "use_slacks on a no-machinery formulation fails template validation" begin - # # slack_spec defaults to NoBranchSlacks, so every pair whose constructors build no - # # slack containers now rejects the request instead of silently ignoring it. - # # StaticBranchUnbounded builds nothing at all; VoltageControlTap never creates slacks. - # sys = PSB.build_system(PSITestSystems, "c_sys5") - # for network_formulation in (DCPNetworkModel, ACPNetworkModel) - # template = - # get_thermal_dispatch_template_network(NetworkModel(network_formulation)) - # set_device_model!( - # template, - # DeviceModel(PSY.Line, StaticBranchUnbounded; use_slacks = true), - # ) - # model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) - # @test_throws IS.ConflictingInputsError POM.validate_template(model) - # end - - # sys14 = PSB.build_system(PSITestSystems, "c_sys14") - # template = get_thermal_dispatch_template_network(NetworkModel(ACPNetworkModel)) - # set_device_model!( - # template, - # DeviceModel(PSY.TwoWindingTransformer, VoltageControlTap; use_slacks = true), - # ) - # model = DecisionModel(template, sys14; optimizer = ipopt_optimizer) - # @test_throws IS.ConflictingInputsError POM.validate_template(model) + # slack_spec defaults to NoBranchSlacks, so every pair whose constructors build no + # slack containers rejects the request instead of silently ignoring it. + # StaticBranchUnbounded builds nothing at all. + sys = PSB.build_system(PSITestSystems, "c_sys5") + for network_formulation in (DCPNetworkModel, ACPNetworkModel) + template = + get_thermal_dispatch_template_network(NetworkModel(network_formulation)) + set_device_model!( + template, + DeviceModel(PSY.Line, StaticBranchUnbounded; use_slacks = true), + ) + model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + @test_throws IS.ConflictingInputsError POM.validate_template(model) + end end @testset "CopperPlateNetworkModel accepts use_slacks as inert with a validation warning" begin diff --git a/test/test_native_lpacc_model.jl b/test/test_native_lpacc_model.jl index 0a5d3c7..f3c96f5 100644 --- a/test/test_native_lpacc_model.jl +++ b/test/test_native_lpacc_model.jl @@ -44,36 +44,27 @@ end end @testset "LPACCNetworkModel rejects reactive control devices at validation" begin - # # LPACC is reactive-capable at the network level (network_has_reactive_power is - # # true), but VoltageControlTap/ShuntSusceptanceDispatch have no LPACC construct - # # path. The validation gate must reject the pairing with a - # # ConflictingInputsError. (build! swallows build/validation exceptions into a - # # FAILED status, so assert against validate_template directly — same pattern as - # # test_network_constructors_with_branch_rating_time_series.jl.) - # sys = PSB.build_system(PSITestSystems, "c_sys14") - # template = get_thermal_dispatch_template_network(NetworkModel(LPACCNetworkModel)) - # set_device_model!(template, PSY.TwoWindingTransformer, VoltageControlTap) - - # model = DecisionModel(template, sys; optimizer = ipopt_optimizer) - # @test_throws IS.ConflictingInputsError POM.validate_template(model) - - # # Same gate for the shunt controller. - # sys_shunt = PSB.build_system(PSITestSystems, "c_sys5") - # bus = PSY.get_component(PSY.ACBus, sys_shunt, "nodeA") - # PSY.add_component!( - # sys_shunt, - # PSY.SwitchedAdmittance(; - # name = "shunt_lpacc_gate", - # available = true, - # bus = bus, - # Y = 0.0 + 0.1im, - # number_of_steps = [2], - # Y_increase = [0.0 + 0.1im], - # ), - # ) - # template_shunt = - # get_thermal_dispatch_template_network(NetworkModel(LPACCNetworkModel)) - # set_device_model!(template_shunt, PSY.SwitchedAdmittance, ShuntSusceptanceDispatch) - # model_shunt = DecisionModel(template_shunt, sys_shunt; optimizer = ipopt_optimizer) - # @test_throws IS.ConflictingInputsError POM.validate_template(model_shunt) + # LPACC is reactive-capable at the network level (network_has_reactive_power is true), + # but ShuntSusceptanceDispatch has no LPACC construct path. The validation gate must + # reject the pairing with a ConflictingInputsError. (build! swallows build/validation + # exceptions into a FAILED status, so assert against validate_template directly — same + # pattern as test_network_constructors_with_branch_rating_time_series.jl.) + sys_shunt = PSB.build_system(PSITestSystems, "c_sys5") + bus = PSY.get_component(PSY.ACBus, sys_shunt, "nodeA") + PSY.add_component!( + sys_shunt, + PSY.SwitchedAdmittance(; + name = "shunt_lpacc_gate", + available = true, + bus = bus, + Y = 0.0 + 0.1im, + number_of_steps = [2], + Y_increase = [0.0 + 0.1im], + ), + ) + template_shunt = + get_thermal_dispatch_template_network(NetworkModel(LPACCNetworkModel)) + set_device_model!(template_shunt, PSY.SwitchedAdmittance, ShuntSusceptanceDispatch) + model_shunt = DecisionModel(template_shunt, sys_shunt; optimizer = ipopt_optimizer) + @test_throws IS.ConflictingInputsError POM.validate_template(model_shunt) end diff --git a/test/test_native_network_reductions.jl b/test/test_native_network_reductions.jl index c186dfc..acef398 100644 --- a/test/test_native_network_reductions.jl +++ b/test/test_native_network_reductions.jl @@ -347,65 +347,6 @@ end @test occursin("absorbed by a network reduction", log) end -# TODO: reenable with Phase angle control -@testset "PhaseAngleControl branch absorbed by a network reduction fails with a clear error" begin - # # "1-6-i_1" is one segment of the (1,2) series chain, so under reduction it has no - # # direct-branch entry of its own — the same _validate_controlled_branch_not_reduced - # # gate exercised above for VoltageControlTap also covers PhaseAngleControl. - # sys = _case11_with_forecast() - # line = PSY.get_component(Line, sys, "1-6-i_1") - # arc = PSY.get_arc(line) - # - # # TODO: phase_angle_limits? - # ps = PSY.TwoWindingTransformer(; - # name = PSY.get_name(line), - # circuit = PSY.TransformerCircuit(; - # available = true, - # active_power_flow = 0.0, - # reactive_power_flow = 0.0, - # r = PSY.get_r(line, PSY.SU), - # x = PSY.get_x(line, PSY.SU), - # tap = 1.0, - # α = 0.0, - # rating = PSY.get_rating(line, PSY.SU), - # arc = arc, - # base_power = PSY.get_base_power(sys, PSY.NU) - # ), - # magnetizing_shunt = 0.0 + 0.0im, - # shunt_location = TwoWindingTransformerShuntLocation.PRIMARY - # ) - # PSY.add_component!(sys, ps) - # PSY.remove_component!(sys, line) - # - # net = NetworkModel( - # DCPNetworkModel; - # reduce_radial_branches = true, - # reduce_degree_two_branches = true, - # ) - # template = get_thermal_dispatch_template_network(net) - # set_device_model!( - # template, DeviceModel(PSY.TwoWindingTransformer, PhaseAngleControl), - # ) - # model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) - # out = mktempdir(; cleanup = true) - # @test build!(model; output_dir = out, console_level = Logging.Error) == - # IOM.ModelBuildStatus.FAILED - # log = read(joinpath(out, "operation_problem.log"), String) - # @test occursin("absorbed by a network reduction", log) -end - -# TODO: reenable with tap control -@testset "tap regulated-bus resolution errors for non-retained bus numbers" begin - # sys = PSB.build_system(PSITestSystems, "c_sys14") - # tr = PSY.get_component(PSY.TwoWindingTransformer, sys, "Trans1") - # PSY.set_regulated_bus_number!(PSY.get_circuit(tr), 999) - # geom = POM._branch_geometry(tr) - # number_to_name = Dict(1 => "Bus 1") - # @test_throws ErrorException POM._tap_regulated_bus_name(tr, geom, number_to_name) - # bus_by_number = Dict(1 => PSY.get_from(PSY.get_arc(tr))) - # @test_throws ErrorException POM._tap_regulated_bus(tr, bus_by_number) -end - @testset "ACP + StaticBranchBounds use_slacks wires flow-definition slacks per reduced arc" begin # The c_sys5 coefficient tests exercise only the identity reduction; here the arcs # genuinely merge, so the flow-definition slack machinery is checked against the diff --git a/test/test_native_tapcontrol.jl b/test/test_native_tapcontrol.jl deleted file mode 100644 index 63e6962..0000000 --- a/test/test_native_tapcontrol.jl +++ /dev/null @@ -1,72 +0,0 @@ -######################################################################################### -# `TapControl` coverage. Disabled in `runtests.jl` until -# `ac_transmission_models/transformer_models.jl` is included again and the `TapControl` -# formulation exists. -# -# The tap physics that DOES ship today — `StaticBranch` under DCP, whose susceptance is -# tap-divided (`b_dc = 1/(tap*x)`) — is covered by `test_native_transformer_tap.jl`, which -# runs. Do not duplicate it here. -######################################################################################### - -@testset "TapControl models transformer tap ratio under DCP (c_sys14)" begin - sys = PSB.build_system(PSITestSystems, "c_sys14") - template = get_thermal_dispatch_template_network(NetworkModel(DCPNetworkModel)) - set_device_model!(template, PSY.TwoWindingTransformer, TapControl) - model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) - @test build!(model; output_dir = mktempdir()) == IOM.ModelBuildStatus.BUILT - @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED - - res = IOM.OptimizationProblemOutputs(model) - base = IOM.get_model_base_power(res) - flow = read_variable( - res, "FlowActivePowerVariable__TwoWindingTransformer"; - table_format = TableFormat.WIDE, - ) - va = read_variable(res, "VoltageAngle__ACBus"; table_format = TableFormat.WIDE) - - tested_a_real_tap = false - for tr in PSY.get_components(PSY.TwoWindingTransformer, sys) - name = PSY.get_name(tr) - @test name in names(flow) - adm = PNM.branch_admittance(tr) - x = -adm.b / (adm.g^2 + adm.b^2) - fr = PSY.get_name(PSY.get_from(PSY.get_arc(tr))) - to = PSY.get_name(PSY.get_to(PSY.get_arc(tr))) - if !isapprox(adm.tap, 1.0; atol = 1e-6) - tested_a_real_tap = true - end - for r in 1:nrow(flow) - p_pu = flow[r, name] / base - expected = (va[r, fr] - va[r, to] - adm.shift) / (x * adm.tap) - @test isapprox(p_pu, expected; atol = 1e-5) - end - end - # Guard: the test system must actually have a non-unit tap, else it proves nothing. - @test tested_a_real_tap -end - -@testset "TapControl differs from StaticBranch for non-unit-tap transformers (c_sys14)" begin - sys = PSB.build_system(PSITestSystems, "c_sys14") - - function _solve_obj(transformer_formulation) - template = get_thermal_dispatch_template_network(NetworkModel(DCPNetworkModel)) - set_device_model!(template, PSY.TwoWindingTransformer, transformer_formulation) - model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) - @test build!(model; output_dir = mktempdir()) == IOM.ModelBuildStatus.BUILT - @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED - return JuMP.objective_value(IOM.get_jump_model(model)) - end - - static_obj = _solve_obj(StaticBranch) - tap_obj = _solve_obj(TapControl) - # `StaticBranch` under DCP now takes its susceptance from `PNM.get_series_susceptance`, - # which is already tap-divided (`1/(tap*x)`), so a FIXED tap is modelled identically by - # both formulations and the two optima must AGREE. This testset previously asserted the - # opposite, from when the DC Ohm's law used the tap-free π susceptance and StaticBranch - # ignored the tap entirely. - # - # Before re-enabling: decide whether `TapControl` still earns its place under DCP at - # all, given StaticBranch subsumes the fixed-tap case. If it survives only to carry a - # variable tap, this comparison should be replaced by a test that moves the tap. - @test isapprox(static_obj, tap_obj; rtol = 1e-8) -end diff --git a/test/test_power_flow_in_the_loop.jl b/test/test_power_flow_in_the_loop.jl index 1f71fef..a15a7ec 100644 --- a/test/test_power_flow_in_the_loop.jl +++ b/test/test_power_flow_in_the_loop.jl @@ -187,72 +187,12 @@ end # ----------------------------------------------------------------------------- # Baseline PFitL coverage (ported from PowerSimulations.jl test file lines 1-548). # These exercise the regular non-headroom paths through the migrated code: -# - PhaseShiftingTransformer in PFitL # - Parallel-line aggregation # - Breaker-switch (DiscreteControlledACBranch) # - HVDCs with DC PowerFlow # - Line active power loss aux variable # ----------------------------------------------------------------------------- -@testset "AC Power Flow in the loop for PhaseShiftingTransformer" begin - # system = buid_system(PSITestSystems, "c_sys5_uc") - # - # line = get_component(Line, system, "1") - # arc = get_arc(line) - # - # ps = PhaseShiftingTransformer(; - # name = get_name(line), - # available = true, - # active_power_flow = 0.0, - # reactive_power_flow = 0.0, - # r = get_r(line, PSY.SU), - # x = get_x(line, PSY.SU), - # primary_shunt = 0.0, - # tap = 1.0, - # α = 0.0, - # rating = get_rating(line, PSY.SU), - # arc = arc, - # base_power = get_base_power(system, PSY.NU), - # ) - # add_component!(system, ps) - # remove_component!(system, line) - # - # template = get_template_dispatch_with_network( - # NetworkModel( - # PTDFNetworkModel; - # network_matrix = PTDF(system), - # evaluations = power_flow_evaluations(ACPowerFlow()), - # ), - # ) - # set_device_model!(template, DeviceModel(PhaseShiftingTransformer, PhaseAngleControl)) - # model_m = DecisionModel(template, system; optimizer = HiGHS_optimizer) - # @test build!(model_m; output_dir = mktempdir(; cleanup = true)) == - # ModelBuildStatus.BUILT - # @test solve!(model_m) == RunStatus.SUCCESSFULLY_FINALIZED - # - # container = get_optimization_container(model_m) - # pf_e_data = only(values(get_evaluation_data(get_evaluations(container)))) - # data = get_inner_data(pf_e_data) - # bus_lookup = PFS.get_bus_lookup(data) - # - # flow_key = VariableKey(FlowActivePowerVariable, PhaseShiftingTransformer) - # flow_values = lookup_value(container, flow_key) - # line_name = get_name(line) - # line_flows = - # [JuMP.value(flow_values[line_name, t]) for t in 1:length(get_time_steps(container))] - # - # # The PhaseShiftingTransformer flow contributes to the "to"-bus active power injection. - # # Both sides are in per-unit; lookup_value returns raw JuMP values in the model unit - # # system rather than the natural-unit conversion that `read_variables(...; WIDE)` - # # performs in PSI. - # @test isapprox( - # data.bus_active_power_injections[bus_lookup[get_number(get_to(arc))], :], - # line_flows; - # atol = 1e-9, - # rtol = 0, - # ) -end - @testset "AC Power Flow in the loop with parallel lines" begin original_line_flow, parallel_line_flow = zero(ComplexF64), zero(ComplexF64) for replace_line in (true, false) diff --git a/test/test_transformer_controls.jl b/test/test_transformer_controls.jl new file mode 100644 index 0000000..0eaa2d9 --- /dev/null +++ b/test/test_transformer_controls.jl @@ -0,0 +1,707 @@ +######################################################################################### +# Transformer tap controls: every `TransformerControlObjective` other than FIXED / +# UNDEFINED. Control is opted into per `DeviceModel` with the `enable_controls` attribute; +# when it is on, each `TransformerCircuit`'s `control_objective` decides what is built. +# A controlled circuit gets a `TapRatioVariable` bounded by `control_limits`, and the +# controlled quantity is held inside `controlled_quantity_limits`. +# +# Fixed / off-nominal tap physics (tap as a constant component property) lives in +# `test_transformer_fixed_tap.jl` — do not duplicate it here. +######################################################################################### + +const VOLTAGE_CONTROL = PSY.TransformerControlObjective.VOLTAGE +const Q_FLOW_CONTROL = PSY.TransformerControlObjective.REACTIVE_POWER_FLOW +const P_FLOW_CONTROL = PSY.TransformerControlObjective.ACTIVE_POWER_FLOW + +_control_attributes(enable::Bool) = + Dict{String, Any}(POM.ENABLE_CONTROLS_KEY => enable) + +""" +`c_sys14` with one transformer circuit put under `objective`. `regulated` picks which end +of the circuit's arc is regulated (the bus number, not a sentinel — the API takes the +number of either the from or the to bus). Returns the system, the transformer, its +circuit, and the regulated bus name. +""" +function _controlled_sys14( + objective; + name = "Trans1", + regulated = :to, + quantity_limits = (min = 0.9, max = 1.1), + control_limits = (min = 0.9, max = 1.1), +) + sys = PSB.build_system(PSITestSystems, "c_sys14") + transformer = PSY.get_component(PSY.TwoWindingTransformer, sys, name) + circuit = PSY.get_circuit(transformer) + arc = PSY.get_arc(circuit) + bus = regulated == :from ? PSY.get_from(arc) : PSY.get_to(arc) + PSY.set_control_objective!(circuit, objective) + PSY.set_regulated_bus_number!(circuit, PSY.get_number(bus)) + PSY.set_controlled_quantity_limits!(circuit, quantity_limits) + PSY.set_control_limits!(circuit, control_limits) + return sys, transformer, circuit, PSY.get_name(bus) +end + +function _controlled_template(network_formulation; enable = true, kwargs...) + template = get_thermal_dispatch_template_network(NetworkModel(network_formulation; kwargs...)) + set_device_model!( + template, + DeviceModel( + PSY.TwoWindingTransformer, + StaticBranch; + attributes = _control_attributes(enable), + ), + ) + return template +end + +function _build_controlled(sys, network_formulation; enable = true, optimizer, kwargs...) + template = _controlled_template(network_formulation; enable = enable, kwargs...) + model = DecisionModel(template, sys; optimizer = optimizer) + status = build!(model; output_dir = mktempdir(; cleanup = true)) + return model, status +end + +_has_tap_variable(container) = + any(k -> occursin("TapRatioVariable", string(k)), keys(IOM.get_variables(container))) + +################################### attribute plumbing ################################# + +@testset "enable_controls is a transformer-only attribute defaulting to false" begin + for T in (PSY.TwoWindingTransformer, PSY.ThreeWindingTransformer) + attributes = POM.get_default_attributes(T, StaticBranch) + @test haskey(attributes, POM.ENABLE_CONTROLS_KEY) + @test attributes[POM.ENABLE_CONTROLS_KEY] === false + end + # Non-transformer branches carry no control switch at all. + @test !haskey( + POM.get_default_attributes(PSY.Line, StaticBranch), + POM.ENABLE_CONTROLS_KEY, + ) + + # The attribute survives onto the DeviceModel and merges with the other defaults. + device_model = DeviceModel( + PSY.TwoWindingTransformer, + StaticBranch; + attributes = _control_attributes(true), + ) + @test IOM.get_attribute(device_model, POM.ENABLE_CONTROLS_KEY) === true + @test IOM.get_attribute(device_model, POM.PARALLEL_BRANCH_MAX_RATING_KEY) == + "single_element_contingency" +end + +@testset "a controlled circuit builds no tap variable while enable_controls is off" begin + sys, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL) + model, status = + _build_controlled(sys, ACPNetworkModel; enable = false, optimizer = ipopt_optimizer) + @test status == IOM.ModelBuildStatus.BUILT + @test !_has_tap_variable(IOM.get_optimization_container(model)) +end + +@testset "TapRatioVariable is created only for controlled circuits, bounded by control_limits" begin + limits = (min = 0.95, max = 1.05) + sys, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL; control_limits = limits) + model, status = _build_controlled(sys, ACPNetworkModel; optimizer = ipopt_optimizer) + @test status == IOM.ModelBuildStatus.BUILT + + container = IOM.get_optimization_container(model) + tap = IOM.get_variable(container, TapRatioVariable, PSY.TwoWindingTransformer) + # Trans2 / Trans3 are left UNDEFINED, so only the controlled circuit gets a variable. + @test axes(tap)[1] == ["Trans1"] + @test check_variable_bounded(model, TapRatioVariable, PSY.TwoWindingTransformer) + for v in tap + @test JuMP.lower_bound(v) == limits.min + @test JuMP.upper_bound(v) == limits.max + end +end + +@testset "REACTIVE_POWER_FLOW control also creates a tap variable" begin + sys, _, _, _ = _controlled_sys14(Q_FLOW_CONTROL) + model, status = _build_controlled(sys, ACPNetworkModel; optimizer = ipopt_optimizer) + @test status == IOM.ModelBuildStatus.BUILT + container = IOM.get_optimization_container(model) + @test axes(IOM.get_variable(container, TapRatioVariable, PSY.TwoWindingTransformer))[1] == + ["Trans1"] +end + +@testset "ACTIVE_POWER_FLOW is a phase-shift objective, not a tap control" begin + # The tap controls cover the voltage / reactive-power objectives; an active-power + # (phase-shifting) circuit must not silently acquire a tap variable. + sys, _, _, _ = _controlled_sys14(P_FLOW_CONTROL) + model, status = _build_controlled(sys, ACPNetworkModel; optimizer = ipopt_optimizer) + @test status == IOM.ModelBuildStatus.BUILT + @test !_has_tap_variable(IOM.get_optimization_container(model)) +end + +@testset "a DISABLED objective builds no control" begin + sys, _, _, _ = + _controlled_sys14(PSY.TransformerControlObjective.VOLTAGE_DISABLED) + model, status = _build_controlled(sys, ACPNetworkModel; optimizer = ipopt_optimizer) + @test status == IOM.ModelBuildStatus.BUILT + @test !_has_tap_variable(IOM.get_optimization_container(model)) +end + +################################### VOLTAGE objective ################################## + +# Solve `c_sys14` once with the transformer uncontrolled and report the regulated bus +# voltage, so each control test can aim its band away from the free-running solution and +# prove the constraint actually bites. +function _uncontrolled_voltage(bus_name; network_formulation = ACPNetworkModel) + sys = PSB.build_system(PSITestSystems, "c_sys14") + model, status = + _build_controlled(sys, network_formulation; enable = false, optimizer = ipopt_optimizer) + @test status == IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + res = IOM.OptimizationProblemOutputs(model) + vm = read_variable(res, "VoltageMagnitude__ACBus"; table_format = TableFormat.WIDE) + return vm[1, bus_name] +end + +@testset "VOLTAGE control holds the regulated bus inside its band (ACP, to-side)" begin + _, _, _, bus_name = _controlled_sys14(VOLTAGE_CONTROL) + free_vm = _uncontrolled_voltage(bus_name) + # A band the free-running solution violates, so holding it requires the tap to move. + band = (min = free_vm + 0.01, max = free_vm + 0.02) + + sys, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL; quantity_limits = band) + model, status = _build_controlled(sys, ACPNetworkModel; optimizer = ipopt_optimizer) + @test status == IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + res = IOM.OptimizationProblemOutputs(model) + vm = read_variable(res, "VoltageMagnitude__ACBus"; table_format = TableFormat.WIDE) + @test bus_name in names(vm) + for r in 1:nrow(vm) + @test vm[r, bus_name] >= band.min - 1e-6 + @test vm[r, bus_name] <= band.max + 1e-6 + end + # The band is on the voltage itself, not on its square. + @test !(free_vm >= band.min - 1e-6 && free_vm <= band.max + 1e-6) + + tap = read_variable( + res, "TapRatioVariable__TwoWindingTransformer"; table_format = TableFormat.WIDE, + ) + for r in 1:nrow(tap) + @test tap[r, "Trans1"] >= 0.9 - 1e-6 + @test tap[r, "Trans1"] <= 1.1 + 1e-6 + end +end + +@testset "VOLTAGE control regulates the from-side bus when its number is given" begin + _, _, _, bus_name = _controlled_sys14(VOLTAGE_CONTROL; regulated = :from) + free_vm = _uncontrolled_voltage(bus_name) + band = (min = free_vm + 0.01, max = free_vm + 0.02) + + sys, _, _, _ = + _controlled_sys14(VOLTAGE_CONTROL; regulated = :from, quantity_limits = band) + model, status = _build_controlled(sys, ACPNetworkModel; optimizer = ipopt_optimizer) + @test status == IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + res = IOM.OptimizationProblemOutputs(model) + vm = read_variable(res, "VoltageMagnitude__ACBus"; table_format = TableFormat.WIDE) + for r in 1:nrow(vm) + @test vm[r, bus_name] >= band.min - 1e-6 + @test vm[r, bus_name] <= band.max + 1e-6 + end +end + +@testset "the VOLTAGE band is a voltage, not a squared voltage" begin + # A band well away from 1.0 pu separates the two readings: v ∈ [0.80, 0.82] is + # satisfied by v² ∈ [0.64, 0.67], so a builder comparing the raw band against v² + # would land the voltage near 0.9 pu instead. + band = (min = 0.80, max = 0.82) + sys, _, _, bus_name = _controlled_sys14(VOLTAGE_CONTROL; quantity_limits = band) + model, status = _build_controlled(sys, ACPNetworkModel; optimizer = ipopt_optimizer) + @test status == IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + res = IOM.OptimizationProblemOutputs(model) + vm = read_variable(res, "VoltageMagnitude__ACBus"; table_format = TableFormat.WIDE) + @test vm[1, bus_name] >= band.min - 1e-6 + @test vm[1, bus_name] <= band.max + 1e-6 +end + +############################ REACTIVE_POWER_FLOW objective ############################# + +@testset "REACTIVE_POWER_FLOW control holds the terminal flow inside its band (ACP)" begin + # `controlled_quantity_limits` reaches the constraint builder unconverted, so it is + # read as system-base pu here; the reported flow is MVAR and divided back down. + band = (min = -0.05, max = 0.05) + sys, _, _, _ = _controlled_sys14(Q_FLOW_CONTROL; quantity_limits = band) + model, status = _build_controlled(sys, ACPNetworkModel; optimizer = ipopt_optimizer) + @test status == IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + res = IOM.OptimizationProblemOutputs(model) + base = IOM.get_model_base_power(res) + for key in ( + "FlowReactivePowerFromToVariable__TwoWindingTransformer", + "FlowReactivePowerToFromVariable__TwoWindingTransformer", + ) + flow = read_variable(res, key; table_format = TableFormat.WIDE) + for r in 1:nrow(flow) + @test flow[r, "Trans1"] / base >= band.min - 1e-6 + @test flow[r, "Trans1"] / base <= band.max + 1e-6 + end + end +end + +################################### model invariants ################################### + +@testset "the model is count-invariant across control objectives (ACP)" begin + function _container_for(objective) + sys, _, _, _ = _controlled_sys14(objective) + model, status = _build_controlled(sys, ACPNetworkModel; optimizer = ipopt_optimizer) + @test status == IOM.ModelBuildStatus.BUILT + return IOM.get_optimization_container(model) + end + + cv = _container_for(VOLTAGE_CONTROL) + cq = _container_for(Q_FLOW_CONTROL) + + var_v = IOM.get_variables(cv) + var_q = IOM.get_variables(cq) + @test Set(keys(var_v)) == Set(keys(var_q)) + for k in keys(var_v) + @test size(var_v[k]) == size(var_q[k]) + end + + con_v = IOM.get_constraints(cv) + con_q = IOM.get_constraints(cq) + @test Set(keys(con_v)) == Set(keys(con_q)) + for k in keys(con_v) + @test size(con_v[k]) == size(con_q[k]) + end +end + +@testset "a tap pinned at nominal reproduces the uncontrolled model (ACP)" begin + # White-box reduction gate: with the tap variable fixed at the circuit's nominal + # ratio and a band too wide to bind, the controlled Ohm's law is term-by-term the + # fixed-tap one, so both models must reach the same optimum and terminal flows. + wide = (min = 0.5, max = 1.5) + + sys_fixed, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL; quantity_limits = wide) + model_fixed, status_fixed = _build_controlled( + sys_fixed, ACPNetworkModel; enable = false, optimizer = ipopt_optimizer, + ) + @test status_fixed == IOM.ModelBuildStatus.BUILT + @test solve!(model_fixed) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + sys_var, transformer, circuit, _ = _controlled_sys14( + VOLTAGE_CONTROL; quantity_limits = wide, control_limits = wide, + ) + model_var, status_var = + _build_controlled(sys_var, ACPNetworkModel; optimizer = ipopt_optimizer) + @test status_var == IOM.ModelBuildStatus.BUILT + + container = IOM.get_optimization_container(model_var) + tap = IOM.get_variable(container, TapRatioVariable, PSY.TwoWindingTransformer) + for t in axes(tap, 2) + JuMP.fix(tap[PSY.get_name(transformer), t], PSY.get_tap(circuit); force = true) + end + @test solve!(model_var) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + res_fixed = IOM.OptimizationProblemOutputs(model_fixed) + res_var = IOM.OptimizationProblemOutputs(model_var) + @test isapprox( + IOM.get_objective_value(res_var), + IOM.get_objective_value(res_fixed); + rtol = 1e-3, + ) + for key in ( + "FlowActivePowerFromToVariable__TwoWindingTransformer", + "FlowReactivePowerFromToVariable__TwoWindingTransformer", + ) + flow_fixed = read_variable(res_fixed, key; table_format = TableFormat.WIDE) + flow_var = read_variable(res_var, key; table_format = TableFormat.WIDE) + for d in PSY.get_components(PSY.TwoWindingTransformer, sys_var) + name = PSY.get_name(d) + @test isapprox(flow_var[1, name], flow_fixed[1, name]; atol = 1e-3) + end + end +end + +@testset "NetworkFlowConstraint carries the live tap variable (ACP coefficients)" begin + # Ground truth: the built from-to flow constraint must use exactly the + # `_tapped_admittance` terms evaluated at the tap VARIABLE. Evaluate + # `constraint_object(con).func` at an arbitrary point and compare against the + # hand-assembled right-hand side. + sys, transformer, _, _ = _controlled_sys14(VOLTAGE_CONTROL) + model, status = _build_controlled(sys, ACPNetworkModel; optimizer = ipopt_optimizer) + @test status == IOM.ModelBuildStatus.BUILT + + container = IOM.get_optimization_container(model) + pft = IOM.get_variable(container, FlowActivePowerFromToVariable, PSY.TwoWindingTransformer) + vm = IOM.get_variable(container, VoltageMagnitude, PSY.ACBus) + va = IOM.get_variable(container, VoltageAngle, PSY.ACBus) + tap = IOM.get_variable(container, TapRatioVariable, PSY.TwoWindingTransformer) + con_pft = IOM.get_constraint( + container, POM.NetworkFlowConstraint, PSY.TwoWindingTransformer, "p_ft", + ) + + t = 1 + name = PSY.get_name(transformer) + arc = PSY.get_arc(PSY.get_circuit(transformer)) + fr = PSY.get_name(PSY.get_from(arc)) + to = PSY.get_name(PSY.get_to(arc)) + adm = PNM.branch_admittance(transformer) + + vals = Dict{JuMP.VariableRef, Float64}( + vm[fr, t] => 1.02, vm[to, t] => 0.98, + va[fr, t] => 0.05, va[to, t] => -0.03, + tap[name, t] => 1.05, + pft[name, t] => 0.7, + ) + lookup = z -> vals[z] + + y = POM._tapped_admittance(adm, vals[tap[name, t]]) + vmf = vals[vm[fr, t]] + vmt = vals[vm[to, t]] + θ = vals[va[fr, t]] - vals[va[to, t]] + rhs = y.g11 * vmf^2 + y.g12 * vmf * vmt * cos(θ) + y.b12 * vmf * vmt * sin(θ) + # `func` is stored as (lhs - rhs). + @test isapprox( + JuMP.value(lookup, JuMP.constraint_object(con_pft[name, t]).func), + vals[pft[name, t]] - rhs; + atol = 1e-10, + ) +end + +################################### network coverage ################################### + +@testset "VOLTAGE control is wired on every voltage-carrying AC network" begin + for network_formulation in (ACRNetworkModel, IVRNetworkModel, LPACCNetworkModel) + @testset "$network_formulation" begin + _, _, _, bus_name = _controlled_sys14(VOLTAGE_CONTROL) + band = (min = 1.00, max = 1.02) + sys, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL; quantity_limits = band) + model, status = + _build_controlled(sys, network_formulation; optimizer = ipopt_optimizer) + @test status == IOM.ModelBuildStatus.BUILT + @test _has_tap_variable(IOM.get_optimization_container(model)) + @test check_variable_bounded(model, TapRatioVariable, PSY.TwoWindingTransformer) + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + res = IOM.OptimizationProblemOutputs(model) + if network_formulation == LPACCNetworkModel + phi = read_variable( + res, "VoltageDeviation__ACBus"; table_format = TableFormat.WIDE, + ) + magnitude = 1.0 + phi[1, bus_name] + else + vr = read_variable(res, "VoltageReal__ACBus"; table_format = TableFormat.WIDE) + vi = read_variable( + res, "VoltageImaginary__ACBus"; table_format = TableFormat.WIDE, + ) + magnitude = sqrt(vr[1, bus_name]^2 + vi[1, bus_name]^2) + end + @test magnitude >= band.min - 1e-4 + @test magnitude <= band.max + 1e-4 + end + end +end + +@testset "DC networks carry no tap control" begin + # The DC network has no voltage magnitude or reactive flow to regulate, so an enabled + # control switch has nothing to build; the tap stays the component property that the + # tap-divided DC susceptance already carries. + sys, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL) + model, status = _build_controlled(sys, DCPNetworkModel; optimizer = HiGHS_optimizer) + @test status == IOM.ModelBuildStatus.BUILT + @test !_has_tap_variable(IOM.get_optimization_container(model)) + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED +end + +################################ reductions and conflicts ############################## + +@testset "a controlled circuit survives the network reduction" begin + # Controlled transformers pin their endpoint buses irreducible, so the circuit keeps + # its own arc (and therefore its own tap variable) even with reductions requested. + sys, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL) + model, status = _build_controlled( + sys, + ACPNetworkModel; + optimizer = ipopt_optimizer, + reduce_radial_branches = true, + reduce_degree_two_branches = true, + ) + @test status == IOM.ModelBuildStatus.BUILT + container = IOM.get_optimization_container(model) + @test axes(IOM.get_variable(container, TapRatioVariable, PSY.TwoWindingTransformer))[1] == + ["Trans1"] +end + +@testset "a controlled circuit merged with a parallel branch fails with a clear error" begin + # PNM collapses parallel branches onto one equivalent arc before POM sees them, which + # would leave the control acting on a flow that is not the transformer's own. + sys, transformer, _, _ = _controlled_sys14(VOLTAGE_CONTROL) + arc = PSY.get_arc(PSY.get_circuit(transformer)) + PSY.add_component!( + sys, + PSY.Line(; + name = "parallel_to_Trans1", + available = true, + active_power_flow = 0.0, + reactive_power_flow = 0.0, + arc = arc, + r = 0.01, + x = 0.1, + b = (from = 0.0, to = 0.0), + rating = 2.0, + angle_limits = (min = -π / 2, max = π / 2), + ), + ) + template = _controlled_template(ACPNetworkModel) + model = DecisionModel(template, sys; optimizer = ipopt_optimizer) + out = mktempdir(; cleanup = true) + @test build!(model; output_dir = out, console_level = Logging.Error) == + IOM.ModelBuildStatus.FAILED + log = read(joinpath(out, "operation_problem.log"), String) + @test occursin("Controlled transformer circuit", log) + @test occursin(PSY.get_name(transformer), log) +end + +# `case11_network_reductions` is the purpose-built reducible system (c_sys14 reduces +# nothing); it carries no forecast, which a DecisionModel build requires. +function _case11_with_forecast() + sys = PSB.build_system(PSITestSystems, "case11_network_reductions") + dummy_data = Dict( + DateTime("2020-01-01T08:00:00") => [5.0, 6, 7, 7, 7], + DateTime("2020-01-01T08:30:00") => [9.0, 9, 9, 9, 8], + DateTime("2020-01-01T09:00:00") => [6.0, 6, 5, 5, 4], + ) + dummy_forecast = Deterministic("max_active_power", dummy_data, Dates.Minute(5)) + load = first(PSY.get_components(PSY.StandardLoad, sys)) + PSY.add_time_series!(sys, load, dummy_forecast) + return sys +end + +@testset "a controlled circuit absorbed by a series reduction fails with a clear error" begin + # "1-6-i_1" is one segment of the (1,2) series chain, so under reduction it has no + # direct-branch entry of its own and its tap would control an equivalent arc it only + # partly owns. + sys = _case11_with_forecast() + line = PSY.get_component(PSY.Line, sys, "1-6-i_1") + arc = PSY.get_arc(line) + transformer = PSY.TwoWindingTransformer(; + name = PSY.get_name(line), + circuit = PSY.TransformerCircuit(; + available = true, + arc = arc, + r = PSY.get_r(line, PSY.SU), + x = PSY.get_x(line, PSY.SU), + tap = 1.0, + α = 0.0, + rating = PSY.get_rating(line, PSY.SU), + control_objective = VOLTAGE_CONTROL, + regulated_bus_number = PSY.get_number(PSY.get_to(arc)), + control_limits = (min = 0.9, max = 1.1), + controlled_quantity_limits = (min = 0.95, max = 1.05), + base_power = PSY.get_base_power(sys, PSY.NU), + ), + magnetizing_shunt = 0.0 + 0.0im, + shunt_location = PSY.TwoWindingTransformerShuntLocation.PRIMARY, + ) + PSY.remove_component!(sys, line) + PSY.add_component!(sys, transformer) + + template = _controlled_template( + ACPNetworkModel; + reduce_radial_branches = true, + reduce_degree_two_branches = true, + ) + model = DecisionModel(template, sys; optimizer = ipopt_optimizer) + out = mktempdir(; cleanup = true) + @test build!(model; output_dir = out, console_level = Logging.Error) == + IOM.ModelBuildStatus.FAILED + @test occursin( + "Controlled transformer circuit", + read(joinpath(out, "operation_problem.log"), String), + ) +end + +@testset "ACP rejects two transformers regulating the same bus" begin + # Trans1 (4 → 9) and Trans3 (4 → 7) share bus 4. Under ACP both would drive the one + # shared VoltageMagnitude, so the conflicting bands must be caught at validation. + sys = PSB.build_system(PSITestSystems, "c_sys14") + for name in ("Trans1", "Trans3") + circuit = PSY.get_circuit(PSY.get_component(PSY.TwoWindingTransformer, sys, name)) + PSY.set_control_objective!(circuit, VOLTAGE_CONTROL) + PSY.set_regulated_bus_number!(circuit, PSY.get_number(PSY.get_from(PSY.get_arc(circuit)))) + PSY.set_controlled_quantity_limits!(circuit, (min = 0.98, max = 1.02)) + end + template = _controlled_template(ACPNetworkModel) + model = DecisionModel(template, sys; optimizer = ipopt_optimizer) + @test_throws IS.ConflictingInputsError POM.validate_template(model) +end + +######################################################################################### +# Phase control (the ACTIVE_POWER_FLOW / ASYMMETRIC_ACTIVE_POWER_FLOW objectives, where the +# phase shift α rather than the tap ratio is the decision variable) is NOT supported yet. +# The testsets below are the coverage that existed for the old `PhaseAngleControl` +# formulation, kept commented until the objective is implemented. They still name +# `PhaseAngleControl` / `PhaseShiftingTransformer`; port them onto the control-objective +# framework when phase control lands. +######################################################################################### + +# @testset "PhaseAngleControl branch absorbed by a network reduction fails with a clear error" begin +# # "1-6-i_1" is one segment of the (1,2) series chain, so under reduction it has no +# # direct-branch entry of its own — the same _validate_controlled_branch_not_reduced +# # gate exercised above for tap control also covers PhaseAngleControl. +# sys = _case11_with_forecast() +# line = PSY.get_component(Line, sys, "1-6-i_1") +# arc = PSY.get_arc(line) +# +# # TODO: phase_angle_limits? +# ps = PSY.TwoWindingTransformer(; +# name = PSY.get_name(line), +# circuit = PSY.TransformerCircuit(; +# available = true, +# active_power_flow = 0.0, +# reactive_power_flow = 0.0, +# r = PSY.get_r(line, PSY.SU), +# x = PSY.get_x(line, PSY.SU), +# tap = 1.0, +# α = 0.0, +# rating = PSY.get_rating(line, PSY.SU), +# arc = arc, +# base_power = PSY.get_base_power(sys, PSY.NU) +# ), +# magnetizing_shunt = 0.0 + 0.0im, +# shunt_location = TwoWindingTransformerShuntLocation.PRIMARY +# ) +# PSY.add_component!(sys, ps) +# PSY.remove_component!(sys, line) +# +# net = NetworkModel( +# DCPNetworkModel; +# reduce_radial_branches = true, +# reduce_degree_two_branches = true, +# ) +# template = get_thermal_dispatch_template_network(net) +# set_device_model!( +# template, DeviceModel(PSY.TwoWindingTransformer, PhaseAngleControl), +# ) +# model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) +# out = mktempdir(; cleanup = true) +# @test build!(model; output_dir = out, console_level = Logging.Error) == +# IOM.ModelBuildStatus.FAILED +# log = read(joinpath(out, "operation_problem.log"), String) +# @test occursin("absorbed by a network reduction", log) +# end + +# @testset "DC Power Flow Models for phase-shifting TwoWindingTransformer and Line" begin +# system = build_system(PSITestSystems, "c_sys5_uc") +# +# line = get_component(Line, system, "1") +# +# ps = TwoWindingTransformer(; +# name = get_name(line), +# available = true, +# active_power_flow = 0.0, +# reactive_power_flow = 0.0, +# r = get_r(line, PSY.SU), +# x = get_r(line, PSY.SU), +# primary_shunt = 0.0, +# tap = 1.0, +# α = 0.0, +# rating = get_rating(line, PSY.SU), +# arc = get_arc(line), +# base_power = get_base_power(system, PSY.NU), +# ) +# +# add_component!(system, ps) +# remove_component!(system, line) +# +# template = get_template_dispatch_with_network( +# NetworkModel(PTDFNetworkModel; network_matrix = PTDF(system)), +# ) +# set_device_model!(template, DeviceModel(TwoWindingTransformer, PhaseAngleControl)) +# model_m = DecisionModel(template, system; optimizer = HiGHS_optimizer) +# @test build!(model_m; output_dir = mktempdir(; cleanup = true)) == +# IOM.ModelBuildStatus.BUILT +# +# @test check_variable_unbounded( +# model_m, +# FlowActivePowerVariable, +# TwoWindingTransformer, +# ) +# +# @test solve!(model_m) == IOM.RunStatus.SUCCESSFULLY_FINALIZED +# +# @test check_flow_variable_values( +# model_m, +# FlowActivePowerVariable, +# TwoWindingTransformer, +# "1", +# get_rating(ps, PSY.SU), +# ) +# +# @test check_flow_variable_values( +# model_m, +# PhaseShifterAngle, +# TwoWindingTransformer, +# "1", +# -π / 2, +# π / 2, +# ) +# end + +# @testset "AC Power Flow in the loop for PhaseShiftingTransformer" begin +# system = buid_system(PSITestSystems, "c_sys5_uc") +# +# line = get_component(Line, system, "1") +# arc = get_arc(line) +# +# ps = PhaseShiftingTransformer(; +# name = get_name(line), +# available = true, +# active_power_flow = 0.0, +# reactive_power_flow = 0.0, +# r = get_r(line, PSY.SU), +# x = get_x(line, PSY.SU), +# primary_shunt = 0.0, +# tap = 1.0, +# α = 0.0, +# rating = get_rating(line, PSY.SU), +# arc = arc, +# base_power = get_base_power(system, PSY.NU), +# ) +# add_component!(system, ps) +# remove_component!(system, line) +# +# template = get_template_dispatch_with_network( +# NetworkModel( +# PTDFNetworkModel; +# network_matrix = PTDF(system), +# evaluations = power_flow_evaluations(ACPowerFlow()), +# ), +# ) +# set_device_model!(template, DeviceModel(PhaseShiftingTransformer, PhaseAngleControl)) +# model_m = DecisionModel(template, system; optimizer = HiGHS_optimizer) +# @test build!(model_m; output_dir = mktempdir(; cleanup = true)) == +# ModelBuildStatus.BUILT +# @test solve!(model_m) == RunStatus.SUCCESSFULLY_FINALIZED +# +# container = get_optimization_container(model_m) +# pf_e_data = only(values(get_evaluation_data(get_evaluations(container)))) +# data = get_inner_data(pf_e_data) +# bus_lookup = PFS.get_bus_lookup(data) +# +# flow_key = VariableKey(FlowActivePowerVariable, PhaseShiftingTransformer) +# flow_values = lookup_value(container, flow_key) +# line_name = get_name(line) +# line_flows = +# [JuMP.value(flow_values[line_name, t]) for t in 1:length(get_time_steps(container))] +# +# # The PhaseShiftingTransformer flow contributes to the "to"-bus active power injection. +# # Both sides are in per-unit; lookup_value returns raw JuMP values in the model unit +# # system rather than the natural-unit conversion that `read_variables(...; WIDE)` +# # performs in PSI. +# @test isapprox( +# data.bus_active_power_injections[bus_lookup[get_number(get_to(arc))], :], +# line_flows; +# atol = 1e-9, +# rtol = 0, +# ) +# end diff --git a/test/test_native_transformer_tap.jl b/test/test_transformer_fixed_tap.jl similarity index 86% rename from test/test_native_transformer_tap.jl rename to test/test_transformer_fixed_tap.jl index 0d7a198..a9bd882 100644 --- a/test/test_native_transformer_tap.jl +++ b/test/test_transformer_fixed_tap.jl @@ -1,13 +1,12 @@ ######################################################################################### -# Off-nominal transformer tap under the native network models. +# Fixed off-nominal transformer tap under the native network models: the tap as a constant +# component property, with no control block (FIXED / UNDEFINED control objectives). Covers +# the DC susceptance `b_dc = 1/(tap*x)` used by `BThetaBranchFlow`/`NetworkFlowConstraint` +# and the Ybus two-port terms in `_tapped_admittance` (both in +# `ac_transmission_models/AC_branches.jl`). # -# These testsets cover only code that ships in the current module: the DC susceptance -# `b_dc = 1/(tap*x)` used by `BThetaBranchFlow`/`NetworkFlowConstraint`, and the Ybus -# two-port terms in `_tapped_admittance` (both in `ac_transmission_models/ -# AC_branches.jl`). They deliberately do NOT touch `TapControl` / `VoltageControlTap`, -# whose formulation files are not yet included — those live in -# `test_native_tapcontrol.jl` / `test_voltage_control_tap_models.jl` and stay disabled -# until the formulations are re-enabled. +# Tap CONTROL — the tap as a decision variable under a `TransformerControlObjective` — is +# covered by `test_transformer_controls.jl`. ######################################################################################### @testset "StaticBranch models transformer off-nominal tap under DCP (c_sys14)" begin diff --git a/test/test_voltage_control_tap_models.jl b/test/test_voltage_control_tap_models.jl deleted file mode 100644 index b081e4c..0000000 --- a/test/test_voltage_control_tap_models.jl +++ /dev/null @@ -1,405 +0,0 @@ -@testset "VoltageControlTap tap bounds are finite (Principle 0)" begin - sys = PSB.build_system(PSITestSystems, "c_sys14") - for tr in PSY.get_components(PSY.TwoWindingTransformer, sys) - lim = POM._tap_ratio_limits(tr) - @test isfinite(lim.min) - @test isfinite(lim.max) - @test lim.min > 0.0 - @test lim.max >= lim.min - end -end - -@testset "VoltageControlTap VOLTAGE objective pins regulated bus voltage (c_sys14)" begin - sys = PSB.build_system(PSITestSystems, "c_sys14") - # Regulate the to-bus of Trans1 (Bus 9) to 1.0 pu via a local (regbus 0) tap. - tr = PSY.get_component(PSY.TwoWindingTransformer, sys, "Trans1") - PSY.set_control_objective!(tr, PSY.TransformerControlObjective.VOLTAGE) - PSY.set_regulated_bus_number!(tr, 0) - PSY.set_voltage_setpoint!(tr, 1.0) - regulated_bus = PSY.get_name(PSY.get_to(PSY.get_arc(tr))) - setpoint = PSY.get_voltage_setpoint(tr) - - template = get_thermal_dispatch_template_network(NetworkModel(ACPNetworkModel)) - set_device_model!(template, PSY.TwoWindingTransformer, VoltageControlTap) - - model = DecisionModel(template, sys; optimizer = ipopt_optimizer) - @test build!(model; output_dir = mktempdir(; cleanup = true)) == - IOM.ModelBuildStatus.BUILT - @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED - - res = IOM.OptimizationProblemOutputs(model) - vm = read_variable(res, "VoltageMagnitude__ACBus"; table_format = TableFormat.WIDE) - @test regulated_bus in names(vm) - for r in 1:nrow(vm) - @test isapprox(vm[r, regulated_bus], setpoint; atol = 1e-6) - end - - # The tap floats within its bounds to hold the setpoint. - @test check_variable_bounded(model, TapRatioVariable, PSY.TwoWindingTransformer) -end - -@testset "VoltageControlTap is count-invariant across control objectives (c_sys14)" begin - function _container_for_objective(objective) - sys = PSB.build_system(PSITestSystems, "c_sys14") - tr = PSY.get_component(PSY.TwoWindingTransformer, sys, "Trans1") - PSY.set_control_objective!(tr, objective) - PSY.set_regulated_bus_number!(tr, 0) - PSY.set_voltage_setpoint!(tr, 1.0) - template = get_thermal_dispatch_template_network(NetworkModel(ACPNetworkModel)) - set_device_model!(template, PSY.TwoWindingTransformer, VoltageControlTap) - model = DecisionModel(template, sys; optimizer = ipopt_optimizer) - @test build!(model; output_dir = mktempdir(; cleanup = true)) == - IOM.ModelBuildStatus.BUILT - return IOM.get_optimization_container(model) - end - - cv = _container_for_objective(PSY.TransformerControlObjective.VOLTAGE) - cq = _container_for_objective(PSY.TransformerControlObjective.REACTIVE_POWER_FLOW) - - var_v = IOM.get_variables(cv) - var_q = IOM.get_variables(cq) - @test Set(keys(var_v)) == Set(keys(var_q)) - for k in keys(var_v) - @test size(var_v[k]) == size(var_q[k]) - end - - con_v = IOM.get_constraints(cv) - con_q = IOM.get_constraints(cq) - @test Set(keys(con_v)) == Set(keys(con_q)) - for k in keys(con_v) - @test size(con_v[k]) == size(con_q[k]) - end -end - -@testset "VoltageControlTap VOLTAGE objective pins regulated bus (ACR, c_sys14)" begin - sys = PSB.build_system(PSITestSystems, "c_sys14") - tr = PSY.get_component(PSY.TwoWindingTransformer, sys, "Trans1") - PSY.set_control_objective!(tr, PSY.TransformerControlObjective.VOLTAGE) - PSY.set_regulated_bus_number!(tr, 0) - PSY.set_voltage_setpoint!(tr, 1.0) - regulated_bus = PSY.get_name(PSY.get_to(PSY.get_arc(tr))) - setpoint = PSY.get_voltage_setpoint(tr) - - template = get_thermal_dispatch_template_network(NetworkModel(ACRNetworkModel)) - set_device_model!(template, PSY.TwoWindingTransformer, VoltageControlTap) - - model = DecisionModel(template, sys; optimizer = ipopt_optimizer) - @test build!(model; output_dir = mktempdir(; cleanup = true)) == - IOM.ModelBuildStatus.BUILT - @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED - - res = IOM.OptimizationProblemOutputs(model) - vr = read_variable(res, "VoltageReal__ACBus"; table_format = TableFormat.WIDE) - vi = read_variable(res, "VoltageImaginary__ACBus"; table_format = TableFormat.WIDE) - @test regulated_bus in names(vr) - for r in 1:nrow(vr) - mag = sqrt(vr[r, regulated_bus]^2 + vi[r, regulated_bus]^2) - @test isapprox(mag, setpoint; atol = 1e-4) - end - - @test check_variable_bounded(model, TapRatioVariable, PSY.TwoWindingTransformer) -end - -@testset "VoltageControlTap is count-invariant across control objectives (ACR, c_sys14)" begin - function _acr_container_for_objective(objective) - sys = PSB.build_system(PSITestSystems, "c_sys14") - tr = PSY.get_component(PSY.TwoWindingTransformer, sys, "Trans1") - PSY.set_control_objective!(tr, objective) - PSY.set_regulated_bus_number!(tr, 0) - PSY.set_voltage_setpoint!(tr, 1.0) - template = get_thermal_dispatch_template_network(NetworkModel(ACRNetworkModel)) - set_device_model!(template, PSY.TwoWindingTransformer, VoltageControlTap) - model = DecisionModel(template, sys; optimizer = ipopt_optimizer) - @test build!(model; output_dir = mktempdir(; cleanup = true)) == - IOM.ModelBuildStatus.BUILT - return IOM.get_optimization_container(model) - end - - cv = _acr_container_for_objective(PSY.TransformerControlObjective.VOLTAGE) - cq = _acr_container_for_objective(PSY.TransformerControlObjective.REACTIVE_POWER_FLOW) - - var_v = IOM.get_variables(cv) - var_q = IOM.get_variables(cq) - @test any(k -> occursin("RegulatedVoltageMagnitude", string(k)), keys(var_v)) - @test Set(keys(var_v)) == Set(keys(var_q)) - for k in keys(var_v) - @test size(var_v[k]) == size(var_q[k]) - end - - con_v = IOM.get_constraints(cv) - con_q = IOM.get_constraints(cq) - @test any(k -> occursin("RegulatedVoltageMagnitudeConstraint", string(k)), keys(con_v)) - @test Set(keys(con_v)) == Set(keys(con_q)) - for k in keys(con_v) - @test size(con_v[k]) == size(con_q[k]) - end -end - -@testset "VoltageControlTap VOLTAGE objective pins regulated bus (IVR, c_sys14)" begin - sys = PSB.build_system(PSITestSystems, "c_sys14") - tr = PSY.get_component(PSY.TwoWindingTransformer, sys, "Trans1") - PSY.set_control_objective!(tr, PSY.TransformerControlObjective.VOLTAGE) - PSY.set_regulated_bus_number!(tr, 0) - PSY.set_voltage_setpoint!(tr, 1.0) - regulated_bus = PSY.get_name(PSY.get_to(PSY.get_arc(tr))) - setpoint = PSY.get_voltage_setpoint(tr) - - template = get_thermal_dispatch_template_network(NetworkModel(IVRNetworkModel)) - set_device_model!(template, PSY.TwoWindingTransformer, VoltageControlTap) - - model = DecisionModel(template, sys; optimizer = ipopt_optimizer) - @test build!(model; output_dir = mktempdir(; cleanup = true)) == - IOM.ModelBuildStatus.BUILT - @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED - - res = IOM.OptimizationProblemOutputs(model) - vr = read_variable(res, "VoltageReal__ACBus"; table_format = TableFormat.WIDE) - vi = read_variable(res, "VoltageImaginary__ACBus"; table_format = TableFormat.WIDE) - @test regulated_bus in names(vr) - for r in 1:nrow(vr) - mag = sqrt(vr[r, regulated_bus]^2 + vi[r, regulated_bus]^2) - @test isapprox(mag, setpoint; atol = 1e-4) - end - - @test check_variable_bounded(model, TapRatioVariable, PSY.TwoWindingTransformer) -end - -@testset "VoltageControlTap IVR currents reduce to fixed-tap at t==tap_nominal" begin - # White-box reduction gate: with the tap variable pinned at its nominal value - # (PSY.get_tap), the variable-tap IVR Ohm's law is term-by-term identical to the - # fixed-tap (StaticBranch) IVR branch, so the two models must converge to the same - # optimum and the same physical (gauge-invariant) terminal power flows. - sys = PSB.build_system(PSITestSystems, "c_sys14") - - template_fixed = get_thermal_dispatch_template_network(NetworkModel(IVRNetworkModel)) - model_fixed = DecisionModel(template_fixed, sys; optimizer = ipopt_optimizer) - @test build!(model_fixed; output_dir = mktempdir(; cleanup = true)) == - IOM.ModelBuildStatus.BUILT - @test solve!(model_fixed) == IOM.RunStatus.SUCCESSFULLY_FINALIZED - - template_var = get_thermal_dispatch_template_network(NetworkModel(IVRNetworkModel)) - set_device_model!(template_var, PSY.TwoWindingTransformer, VoltageControlTap) - model_var = DecisionModel(template_var, sys; optimizer = ipopt_optimizer) - @test build!(model_var; output_dir = mktempdir(; cleanup = true)) == - IOM.ModelBuildStatus.BUILT - - # Pin every tap variable at its nominal ratio before solving. - container = IOM.get_optimization_container(model_var) - tapvar = IOM.get_variable(container, TapRatioVariable, PSY.TwoWindingTransformer) - for d in PSY.get_components(PSY.TwoWindingTransformer, sys) - name = PSY.get_name(d) - for t in axes(tapvar, 2) - JuMP.fix(tapvar[name, t], PSY.get_tap(d); force = true) - end - end - @test solve!(model_var) == IOM.RunStatus.SUCCESSFULLY_FINALIZED - - obj_fixed = IOM.get_objective_value(IOM.OptimizationProblemOutputs(model_fixed)) - obj_var = IOM.get_objective_value(IOM.OptimizationProblemOutputs(model_var)) - @test isapprox(obj_var, obj_fixed; rtol = 1e-3) - - # Compare physical terminal flows on the TwoWindingTransformers (reference-invariant). - res_fixed = IOM.OptimizationProblemOutputs(model_fixed) - res_var = IOM.OptimizationProblemOutputs(model_var) - pft_fixed = read_variable( - res_fixed, "FlowActivePowerFromToVariable__TwoWindingTransformer"; - table_format = TableFormat.WIDE, - ) - pft_var = read_variable( - res_var, "FlowActivePowerFromToVariable__TwoWindingTransformer"; - table_format = TableFormat.WIDE, - ) - qft_fixed = read_variable( - res_fixed, "FlowReactivePowerFromToVariable__TwoWindingTransformer"; - table_format = TableFormat.WIDE, - ) - qft_var = read_variable( - res_var, "FlowReactivePowerFromToVariable__TwoWindingTransformer"; - table_format = TableFormat.WIDE, - ) - for d in PSY.get_components(PSY.TwoWindingTransformer, sys) - name = PSY.get_name(d) - @test isapprox(pft_var[1, name], pft_fixed[1, name]; atol = 1e-3) - @test isapprox(qft_var[1, name], qft_fixed[1, name]; atol = 1e-3) - end -end - -@testset "VoltageControlTap is count-invariant across control objectives (IVR, c_sys14)" begin - function _ivr_container_for_objective(objective) - sys = PSB.build_system(PSITestSystems, "c_sys14") - tr = PSY.get_component(PSY.TwoWindingTransformer, sys, "Trans1") - PSY.set_control_objective!(tr, objective) - PSY.set_regulated_bus_number!(tr, 0) - PSY.set_voltage_setpoint!(tr, 1.0) - template = get_thermal_dispatch_template_network(NetworkModel(IVRNetworkModel)) - set_device_model!(template, PSY.TwoWindingTransformer, VoltageControlTap) - model = DecisionModel(template, sys; optimizer = ipopt_optimizer) - @test build!(model; output_dir = mktempdir(; cleanup = true)) == - IOM.ModelBuildStatus.BUILT - return IOM.get_optimization_container(model) - end - - cv = _ivr_container_for_objective(PSY.TransformerControlObjective.VOLTAGE) - cq = _ivr_container_for_objective(PSY.TransformerControlObjective.REACTIVE_POWER_FLOW) - cp = _ivr_container_for_objective(PSY.TransformerControlObjective.ACTIVE_POWER_FLOW) - - var_v = IOM.get_variables(cv) - var_q = IOM.get_variables(cq) - var_p = IOM.get_variables(cp) - @test any(k -> occursin("RegulatedVoltageMagnitude", string(k)), keys(var_v)) - @test Set(keys(var_v)) == Set(keys(var_q)) == Set(keys(var_p)) - for k in keys(var_v) - @test size(var_v[k]) == size(var_q[k]) == size(var_p[k]) - end - - con_v = IOM.get_constraints(cv) - con_q = IOM.get_constraints(cq) - con_p = IOM.get_constraints(cp) - @test any(k -> occursin("RegulatedVoltageMagnitudeConstraint", string(k)), keys(con_v)) - @test Set(keys(con_v)) == Set(keys(con_q)) == Set(keys(con_p)) - for k in keys(con_v) - @test size(con_v[k]) == size(con_q[k]) == size(con_p[k]) - end -end - -@testset "VoltageControlTap @info-drop under DCPNetworkModel (c_sys14)" begin - sys = PSB.build_system(PSITestSystems, "c_sys14") - template = get_thermal_dispatch_template_network(NetworkModel(DCPNetworkModel)) - set_device_model!(template, PSY.TwoWindingTransformer, VoltageControlTap) - - model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) - # The voltage-controlling tap formulation is reactive-only, so it is dropped - # with an @info from the (active-power-only) DC template during validation. - # A TwoWindingTransformer is a branch the DC network still requires to be modeled, so - # the build then fails on the now-unmodeled branch (unlike a droppable shunt - # injection). Both facts are asserted: the drop happened, and the build failed. - @test build!(model; output_dir = mktempdir(; cleanup = true)) == - IOM.ModelBuildStatus.FAILED - @test !haskey(get_branch_models(get_template(model)), :TwoWindingTransformer) -end - -@testset "ACP rejects two voltage regulators on one bus" begin - # Two TwoWindingTransformers both set to VOLTAGE control regulating bus 9. Under ACP each - # pins the shared network VoltageMagnitude via JuMP.fix(force=true), so the second - # silently overrides the first. validate_template! must reject this. (build! - # swallows the throw into FAILED, so assert against validate_template directly.) - sys = PSB.build_system(PSITestSystems, "c_sys14") - for nm in ("Trans1", "Trans2") - tr = PSY.get_component(PSY.TwoWindingTransformer, sys, nm) - PSY.set_control_objective!(tr, PSY.TransformerControlObjective.VOLTAGE) - PSY.set_regulated_bus_number!(tr, 9) - PSY.set_voltage_setpoint!(tr, 1.0) - end - template = get_thermal_dispatch_template_network(NetworkModel(ACPNetworkModel)) - set_device_model!(template, PSY.TwoWindingTransformer, VoltageControlTap) - model = DecisionModel(template, sys; optimizer = ipopt_optimizer) - @test_throws IS.ConflictingInputsError POM.validate_template(model) -end - -@testset "ACR does not validation-reject two regulators on one bus" begin - # Under ACR each regulator owns a (component, tag) RegulatedVoltageMagnitude aux - # variable tied by vm_reg^2 == vr^2 + vi^2, so the conflict is solver-infeasibility, - # not a validation error. validate_template must NOT throw. - sys = PSB.build_system(PSITestSystems, "c_sys14") - for nm in ("Trans1", "Trans2") - tr = PSY.get_component(PSY.TwoWindingTransformer, sys, nm) - PSY.set_control_objective!(tr, PSY.TransformerControlObjective.VOLTAGE) - PSY.set_regulated_bus_number!(tr, 9) - PSY.set_voltage_setpoint!(tr, 1.0) - end - template = get_thermal_dispatch_template_network(NetworkModel(ACRNetworkModel)) - set_device_model!(template, PSY.TwoWindingTransformer, VoltageControlTap) - model = DecisionModel(template, sys; optimizer = ipopt_optimizer) - @test POM.validate_template(model) === nothing -end - -# The `_tap_flow_coefficients` hand-computed ground truth lives in -# `test_native_transformer_tap.jl`: it exercises only `AC_branches.jl`, which is included, -# so it must not sit behind this file's `DISABLED_TESTS` entry. - -@testset "ACR NetworkFlowConstraint coefficients equal _tap_flow_coefficients" begin - # Ground-truth: the built ACR to-from flow constraint (ptf/qtf) must use exactly the - # pure-function coefficients. Evaluate constraint_object(con).func at chosen variable - # values and compare to the hand RHS assembled from _tap_flow_coefficients. - sys = PSB.build_system(PSITestSystems, "c_sys14") - template = get_thermal_dispatch_template_network(NetworkModel(ACRNetworkModel)) - set_device_model!(template, PSY.TwoWindingTransformer, VoltageControlTap) - model = DecisionModel(template, sys; optimizer = ipopt_optimizer) - @test build!(model; output_dir = mktempdir(; cleanup = true)) == - IOM.ModelBuildStatus.BUILT - - container = IOM.get_optimization_container(model) - ptf = IOM.get_variable( - container, - FlowActivePowerToFromVariable, - PSY.TwoWindingTransformer, - ) - qft = IOM.get_variable( - container, - FlowReactivePowerFromToVariable, - PSY.TwoWindingTransformer, - ) - vr = IOM.get_variable(container, VoltageReal, PSY.ACBus) - vi = IOM.get_variable(container, VoltageImaginary, PSY.ACBus) - tap = IOM.get_variable(container, TapRatioVariable, PSY.TwoWindingTransformer) - con_ptf = - IOM.get_constraint( - container, - POM.NetworkFlowConstraint, - PSY.TwoWindingTransformer, - "p_tf", - ) - con_qft = - IOM.get_constraint( - container, - POM.NetworkFlowConstraint, - PSY.TwoWindingTransformer, - "q_ft", - ) - - t = 1 - for d in Iterators.take(PSY.get_components(PSY.TwoWindingTransformer, sys), 3) - name = PSY.get_name(d) - # Read the π-model and endpoints straight from PNM/PSY rather than through - # `_branch_geometry`, which is reduction-keyed and takes the reduction data. - # c_sys14 reduces nothing, so the device's own arc is the retained arc. - adm = PNM.branch_admittance(d) - coef = POM._tap_flow_coefficients( - adm.g, adm.b, adm.g_fr, adm.b_fr, adm.g_to, adm.b_to, adm.shift, - ) - e_sin = -coef.d_sin - arc = PSY.get_arc(PSY.get_circuit(d)) - fr = PSY.get_name(PSY.get_from(arc)) - to = PSY.get_name(PSY.get_to(arc)) - - # Arbitrary evaluation point for the (nonlinear) constraint functions. - vals = Dict{JuMP.VariableRef, Float64}( - vr[fr, t] => 1.02, vi[fr, t] => 0.05, - vr[to, t] => 0.98, vi[to, t] => -0.03, - tap[name, t] => 1.05, - ptf[name, t] => 0.7, qft[name, t] => -0.2, - ) - lookup = z -> vals[z] - - vv_to = vals[vr[to, t]]^2 + vals[vi[to, t]]^2 - cosprod = vals[vr[fr, t]] * vals[vr[to, t]] + vals[vi[fr, t]] * vals[vi[to, t]] - sinprod = vals[vi[fr, t]] * vals[vr[to, t]] - vals[vr[fr, t]] * vals[vi[to, t]] - tt = vals[tap[name, t]] - - # func is stored as (lhs - rhs); assert it matches the hand-assembled (lhs - rhs). - rhs_ptf = - coef.gg_to * vv_to + coef.c_cos / tt * cosprod + e_sin / tt * (-sinprod) - want_ptf = vals[ptf[name, t]] - rhs_ptf - got_ptf = JuMP.value(lookup, JuMP.constraint_object(con_ptf[name, t]).func) - @test isapprox(got_ptf, want_ptf; atol = 1e-10) - - vv_fr = vals[vr[fr, t]]^2 + vals[vi[fr, t]]^2 - # From side: only the series term is tap-referred, the magnetizing shunt is not. - rhs_qft = - -(coef.b / tt^2 + coef.b_fr) * vv_fr + (-coef.a_sin) / tt * cosprod + - coef.a_cos / tt * sinprod - want_qft = vals[qft[name, t]] - rhs_qft - got_qft = JuMP.value(lookup, JuMP.constraint_object(con_qft[name, t]).func) - @test isapprox(got_qft, want_qft; atol = 1e-10) - end -end From 7df5eb9daefc05bc46ba94359a4a38a588ad81e1 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Fri, 14 Aug 2026 17:08:11 -0400 Subject: [PATCH 08/11] claude ported tests; tests now pass --- src/ac_transmission_models/AC_branches.jl | 296 +++++++++++------- .../branch_constructor.jl | 4 +- .../instantiate_network_model.jl | 26 +- test/Project.toml | 2 - test/includes.jl | 4 +- test/test_transformer_controls.jl | 241 +++++++------- 6 files changed, 324 insertions(+), 249 deletions(-) diff --git a/src/ac_transmission_models/AC_branches.jl b/src/ac_transmission_models/AC_branches.jl index 106f402..a584c06 100644 --- a/src/ac_transmission_models/AC_branches.jl +++ b/src/ac_transmission_models/AC_branches.jl @@ -74,9 +74,24 @@ ENABLE_CONTROLS_KEY = "enable_controls" _control_attribute( ::Union{Type{PSY.TwoWindingTransformer}, Type{PSY.ThreeWindingTransformer}}, ) = (ENABLE_CONTROLS_KEY => false,) - _control_attribute(_) = () +_control_enabled(m::Union{DeviceModel{PSY.TwoWindingTransformer}, DeviceModel{PSY.ThreeWindingTransformer}}) = get_attribute(m, ENABLE_CONTROLS_KEY) === true +_control_enabled(c::PSY.TransformerCircuit) = PSY.get_available(c) && !(PSY.get_control_objective(c) in (PSY.TransformerControlObjective.UNDEFINED, PSY.TransformerControlObjective.FIXED)) +_control_enabled(_) = false + +_tap_controlled(c::PSY.TransformerControlObjective) = c in (PSY.TransformerControlObjective.VOLTAGE, PSY.TransformerControlObjective.REACTIVE_POWER_FLOW) +_tap_controlled(c::PSY.TransformerCircuit) = PSY.get_available(c) && _tap_controlled(PSY.get_control_objective(c)) + +_voltage_controlled(c::PSY.TransformerControlObjective) = c === PSY.TransformerControlObjective.VOLTAGE +_voltage_controlled(c::PSY.TransformerCircuit) = PSY.get_available(c) && _voltage_controlled(PSY.get_control_objective(c)) + +_reactive_controlled(c::PSY.TransformerControlObjective) = c === PSY.TransformerControlObjective.REACTIVE_POWER_FLOW + +_tap_controlled(m::DeviceModel, d) = _control_enabled(m) && _tap_controlled(d) +_voltage_controlled(m::DeviceModel, d) = _control_enabled(m) && _voltage_controlled(d) +_reactive_controlled(m::DeviceModel, d) = _control_enabled(m) && _reactive_controlled(d) + """ DeviceModel attribute key selecting which `PowerNetworkMatrices` function aggregates the individual circuit ratings of a `PNM.BranchesParallel` into a single maximum flow @@ -92,7 +107,7 @@ function get_default_attributes( ) where {U <: PSY.ACTransmission, V <: AbstractBranchFormulation} return Dict{String, Any}( PARALLEL_BRANCH_MAX_RATING_KEY => "single_element_contingency", - _control_attribute(U)... + _control_attribute(U)..., ) end @@ -103,7 +118,7 @@ function get_default_attributes( return Dict{String, Any}( PARALLEL_BRANCH_MAX_RATING_KEY => "single_element_contingency", "include_planned_outages" => false, - _control_attribute(U)... + _control_attribute(U)..., ) end @@ -287,9 +302,14 @@ _circuit_arc_name(d::PSY.TwoWindingTransformer, ::PSY.TransformerCircuit, ::Int) _circuit_arc_name(d::PSY.ThreeWindingTransformer, c::PSY.TransformerCircuit, i::Int) = PNM.get_name(PNM.ThreeWindingTransformerCircuit(d, c, i)) -# TODO: add other controls + refactor -# TODO: Change TransformerCircuit <: DeviceParameter -> <: Device -function _add_transformer_control_variables!( +_add_tap_control_variables!( + ::OptimizationContainer, + ::DeviceModel, + ::IS.FlattenIteratorWrapper, + ::NetworkModel, +) = nothing + +function _add_tap_control_variables!( container::OptimizationContainer, model::DeviceModel{U, F}, devices::IS.FlattenIteratorWrapper{U}, @@ -299,21 +319,22 @@ function _add_transformer_control_variables!( F <: AbstractBranchFormulation, } get_attribute(model, ENABLE_CONTROLS_KEY) === true || return - # TODO: collectzip is ugly, refactor now or do Circuit <: Device first - names, circuits = collect(zip([ - (_circuit_arc_name(d, c, i), c) - for d in devices - for (i, c) in enumerate(PSY.get_circuits(d)) - if PSY.get_available(c) && PSY.get_control_objective(c) in (PSY.TransformerControlObjective.VOLTAGE, PSY.TransformerControlObjective.REACTIVE_POWER_FLOW) - ])) + names = String[] + circuits = PSY.TransformerCircuit[] + for d in devices, (i, c) in enumerate(PSY.get_circuits(d)) + _tap_controlled(c) || continue + push!(names, _circuit_arc_name(d, c, i)) + push!(circuits, c) + end + isempty(names) && return _validate_controlled_branch_not_reduced(network_model, U, names) time_steps = get_time_steps(container) jump_model = get_jump_model(container) - variable = add_variable_container!(container, TapRatioVariable, U, names, time_steps) + tap_var = add_variable_container!(container, TapRatioVariable, U, names, time_steps) for (i, name) in enumerate(names), t in time_steps - bounds = get_control_limits(circuits[i]) - variable[name, t] = JuMP.@variable( + bounds = PSY.get_control_limits(circuits[i]) + tap_var[name, t] = JuMP.@variable( jump_model, base_name = "TapRatioVariable_$(U)_{$(name), $(t)}", lower_bound = bounds.min, @@ -1178,25 +1199,6 @@ function _validate_controlled_branch_not_reduced( return end -# Concrete element type for `_branch_geometries` so constraint builders stay type-stable -# and empty axes still yield `String` name comprehensions (an axis can be empty when the -# other branch type's constructor claimed every shared reduced arc first). -const BranchGeometry = @NamedTuple{ - name::String, - from_name::String, - to_name::String, - from_number::Int, - to_number::Int, - adm::NamedTuple{ - (:g, :b, :g_fr, :b_fr, :g_to, :b_to, :tap, :shift), - NTuple{8, Float64}, - }, - b_dc::Float64, - shift_dc::Float64, - r_dc::Float64, - direct::Bool, -} - _is_aggregate(::PNM.AbstractReductionAggregate) = true _is_aggregate(::PSY.ACTransmission) = false @@ -1210,16 +1212,46 @@ _dc_phase_shift(branch::PNM.AbstractReductionAggregate, nr::PNM.NetworkReduction _dc_phase_shift(branch::PSY.ACTransmission, ::PNM.NetworkReductionData) = PNM.get_series_phase_shift(branch) -_control_objective(c::PSY.TransformerCircuit) = PSY.get_control_objective(c) -_control_objective(_) = PSY.TransformerControlObjective.UNDEFINED +_get_circuit( + b::Union{PSY.TwoWindingTransformer, PNM.ThreeWindingTransformerCircuit}, +) = PSY.get_circuit(b) +_get_circuit(_) = nothing +_control_objective(branch) = _control_objective(_get_circuit(branch)) +_control_objective(::Nothing) = PSY.TransformerControlObjective.UNDEFINED +_control_objective(c::PSY.TransformerCircuit) = + PSY.get_available(c) ? PSY.get_control_objective(c) : + PSY.TransformerControlObjective.UNDEFINED + +_quantity_limits(branch) = _quantity_limits(_get_circuit(branch)) +_quantity_limits(::Nothing) = (min = -Inf, max = Inf) _quantity_limits(c::PSY.TransformerCircuit) = PSY.get_controlled_quantity_limits(c) -_quantity_limits(_) = (min = -Inf, max = Inf) -_regulated_bus(c::PSY.TransformerCircuit) = PSY.get_regulated_bus_number(c) === PSY.get_number(PSY.get_from(PSY.get_arc(c))) ? :from : :to -_regulated_bus(_) = nothing +_regulated_number(branch) = _regulated_number(_get_circuit(branch)) +_regulated_number(::Nothing) = -1 +_regulated_number(c::PSY.TransformerCircuit) = PSY.get_regulated_bus_number(c) + +Base.@kwdef struct BranchGeometry + name::String + from_name::String + to_name::String + from_number::Int + to_number::Int + adm::NamedTuple{ + (:g, :b, :g_fr, :b_fr, :g_to, :b_to, :tap, :shift), + NTuple{8, Float64}, + } + b_dc::Float64 + shift_dc::Float64 + r_dc::Float64 + direct::Bool + control::PSY.TransformerControlObjective + quantity_limits::MinMax + regulated_number::Int -function _branch_geometry( +end + +function BranchGeometry( nr::PNM.NetworkReductionData, number_to_name::Dict{Int, String}, name::String, @@ -1228,7 +1260,7 @@ function _branch_geometry( ) from_no = arc_tuple[1] to_no = arc_tuple[2] - return ( + return BranchGeometry(; name = name, from_name = number_to_name[from_no], to_name = number_to_name[to_no], @@ -1239,17 +1271,18 @@ function _branch_geometry( shift_dc = _dc_phase_shift(branch, nr), r_dc = PNM.arc_dc_resistance(nr, arc_tuple), direct = !_is_aggregate(branch), - objective = _control_objective(branch) - quantity_limits = _quantity_limits(branch) - regulated_bus = _regulated_bus(branch) + control = _control_objective(branch), + quantity_limits = _quantity_limits(branch), + regulated_number = _regulated_number(branch), ) end +_tap_controlled(g::BranchGeometry) = _tap_controlled(g.control) +_voltage_controlled(g::BranchGeometry) = _voltage_controlled(g.control) +_reactive_controlled(g::BranchGeometry) = _reactive_controlled(g.control) """ -Per-branch network geometry for the native nodal constraint builders. - -One geometry per arc of `T` not yet claimed for the constraint family `C` — the -representative axis from [`get_branch_argument_constraint_axis`](@ref) — with PNM's +One [`BranchGeometry`](@ref) per arc of `T` not yet claimed for the constraint family `C` — +the representative axis from [`get_branch_argument_constraint_axis`](@ref) — with PNM's reduction-aware equivalent admittance. Every member of a reduced arc (series segments, parallel groups, across branch types) @@ -1270,7 +1303,7 @@ function _branch_geometries( arc_map = get_name_to_arc_map_entries(nr, T) all_branch_maps_by_type = PNM.get_all_branch_maps_by_type(nr) geoms = BranchGeometry[ - _branch_geometry( + BranchGeometry( nr, number_to_name, name, @@ -1603,16 +1636,16 @@ function _voltage_products( to_bus::String, t::Int, ) + jump_model = get_jump_model(container) vm = get_variable(container, VoltageMagnitude, PSY.ACBus) va = get_variable(container, VoltageAngle, PSY.ACBus) - vmf = vm[from_bus, t] - vmt = vm[to_bus, t] - θ = va[from_bus, t] - va[to_bus, t] + vmf, vmt = vm[from_bus, t], vm[to_bus, t] + vaf, vat = va[from_bus, t], va[to_bus, t] return ( - v2_fr = vmf^2, - v2_to = vmt^2, - vv_cos = vmf * vmt * cos(θ), - vv_sin = vmf * vmt * sin(θ), + v2_fr = JuMP.@expression(jump_model, vmf^2), + v2_to = JuMP.@expression(jump_model, vmt^2), + vv_cos = JuMP.@expression(jump_model, vmf * vmt * cos(vaf - vat)), + vv_sin = JuMP.@expression(jump_model, vmf * vmt * sin(vaf - vat)), ) end @@ -1625,17 +1658,16 @@ function _voltage_products( to_bus::String, t::Int, ) + jump_model = get_jump_model(container) vr = get_variable(container, VoltageReal, PSY.ACBus) vi = get_variable(container, VoltageImaginary, PSY.ACBus) - vr_fr = vr[from_bus, t] - vr_to = vr[to_bus, t] - vi_fr = vi[from_bus, t] - vi_to = vi[to_bus, t] + vr_fr, vr_to = vr[from_bus, t], vr[to_bus, t] + vi_fr, vi_to = vi[from_bus, t], vi[to_bus, t] return ( - v2_fr = vr_fr^2 + vi_fr^2, - v2_to = vr_to^2 + vi_to^2, - vv_cos = vr_fr * vr_to + vi_fr * vi_to, - vv_sin = vi_fr * vr_to - vr_fr * vi_to, + v2_fr = JuMP.@expression(jump_model, vr_fr^2 + vi_fr^2), + v2_to = JuMP.@expression(jump_model, vr_to^2 + vi_to^2), + vv_cos = JuMP.@expression(jump_model, vr_fr * vr_to + vi_fr * vi_to), + vv_sin = JuMP.@expression(jump_model, vi_fr * vr_to - vr_fr * vi_to), ) end @@ -1648,31 +1680,31 @@ function _voltage_products( to_bus::String, t::Int, ) where {T <: PSY.ACTransmission} + jump_model = get_jump_model(container) va = get_variable(container, VoltageAngle, PSY.ACBus) phi = get_variable(container, VoltageDeviation, PSY.ACBus) cs = get_variable(container, CosineApproximation, T) - phi_fr = phi[from_bus, t] - phi_to = phi[to_bus, t] + phi_fr, phi_to = phi[from_bus, t], phi[to_bus, t] return ( - v2_fr = 1.0 + 2.0 * phi_fr, - v2_to = 1.0 + 2.0 * phi_to, - vv_cos = cs[name, t] + phi_fr + phi_to, - vv_sin = va[from_bus, t] - va[to_bus, t], + v2_fr = JuMP.@expression(jump_model, 1.0 + 2.0 * phi_fr), + v2_to = JuMP.@expression(jump_model, 1.0 + 2.0 * phi_to), + vv_cos = JuMP.@expression(jump_model, cs[name, t] + phi_fr + phi_to), + vv_sin = JuMP.@expression(jump_model, va[from_bus, t] - va[to_bus, t]), ) end # Ybus terms, supporting Float64 and VariableRef taps. PNM's ybus functions # use imaginary numbers which VariableRef doesn't support. -function _tapped_admittance(adm, tap) - cs = cos(adm.shift) - sn = sin(adm.shift) +function _tapped_admittance(jump_model, adm, tap) + g_cos, g_sin = adm.g * cos(adm.shift), adm.g * sin(adm.shift) + b_cos, b_sin = adm.b * cos(adm.shift), adm.b * sin(adm.shift) return ( - g11 = adm.g / tap^2 + adm.g_fr, - b11 = adm.b / tap^2 + adm.b_fr, - g12 = (-adm.g * cs + adm.b * sn) / tap, - b12 = (-adm.b * cs - adm.g * sn) / tap, - g21 = (-adm.g * cs - adm.b * sn) / tap, - b21 = (adm.g * sn - adm.b * cs) / tap, + g11 = JuMP.@expression(jump_model, adm.g / tap^2 + adm.g_fr), + b11 = JuMP.@expression(jump_model, adm.b / tap^2 + adm.b_fr), + g12 = JuMP.@expression(jump_model, (-g_cos + b_sin) / tap), + b12 = JuMP.@expression(jump_model, (-b_cos - g_sin) / tap), + g21 = JuMP.@expression(jump_model, (-g_cos - b_sin) / tap), + b21 = JuMP.@expression(jump_model, (g_sin - b_cos) / tap), g22 = adm.g + adm.g_to, b22 = adm.b + adm.b_to, ) @@ -1698,8 +1730,6 @@ function add_constraints!( qft = get_variable(container, FlowReactivePowerFromToVariable, T) qtf = get_variable(container, FlowReactivePowerToFromVariable, T) - tap_var = get_attribute(device_model, ENABLE_CONTROLS_KEY) ? get_variable(container, TapRatioVariable, T) : nothing - number_to_name = _retained_number_to_name(sys, network_model) geoms = _branch_geometries(number_to_name, network_model, devices, T, NetworkFlowConstraint) @@ -1717,14 +1747,14 @@ function add_constraints!( for t in time_steps vp = _voltage_products(container, network_model, T, name, from_bus, to_bus, t) - tap = g_geom.control in TAP_CONTROL_OBJECTIVES ? tap_var[name, t] : adm.tap - y = _tapped_admittance(adm, tap) + tap = _tap_controlled(device_model, g_geom) ? get_variable(container, TapRatioVariable, T)[name, t] : adm.tap + y = _tapped_admittance(jump_model, adm, tap) cons_pft[name, t] = JuMP.@constraint( jump_model, pft[name, t] == y.g11 * vp.v2_fr + y.g12 * vp.vv_cos + y.b12 * vp.vv_sin + - p_slack_term(slacks.p_ft, name, t) + _slack_term(slacks.p_ft, name, t) ) cons_ptf[name, t] = JuMP.@constraint( jump_model, @@ -1745,12 +1775,8 @@ function add_constraints!( cons_qft[name, t] = JuMP.@constraint(jump_model, qft[name, t] == qft_expr) cons_qtf[name, t] = JuMP.@constraint(jump_model, qtf[name, t] == qtf_expr) - if g_geom.control === PSY.TransformerControlObjective.VOLTAGE - voltage = g_geom.regulated_bus == :from ? vp.v2_fr : vp.v2_to - JuMP.@constraint(jump_model, voltage >= g_geom.quantity_limits.min) - JuMP.@constraint(jump_model, voltage <= g_geom.quantity_limits.max) - end - if g_geom.control === PSY.TransformerControlObjective.REACTIVE_POWER_FLOW + # TODO: register in containers + if _reactive_controlled(device_model, g_geom) JuMP.@constraint(jump_model, qft_expr >= g_geom.quantity_limits.min) JuMP.@constraint(jump_model, qft_expr <= g_geom.quantity_limits.max) JuMP.@constraint(jump_model, qtf_expr >= g_geom.quantity_limits.min) @@ -1761,6 +1787,62 @@ function add_constraints!( return end +_voltage_magnitude(::Type{<:Union{ACPNetworkModel, IVRNetworkModel}}) = VoltageMagnitude +_voltage_magnitude(_) = RegulatedVoltageMagnitude + +# TODO: It might benefit solvers to reset VariableRef bounds instead of using singleton constraints (or even to use fix()), but then we lose the dual. Worth it? +function _add_voltage_control_constraints!( + container::OptimizationContainer, + sys::PSY.System, + devices::IS.FlattenIteratorWrapper{T}, + device_model::DeviceModel{T}, + ::NetworkModel{N} +) where {T <: Union{PSY.TwoWindingTransformer, PSY.ThreeWindingTransformer}, N <: AbstractNetworkModel} + _control_enabled(device_model) || return + + cons = add_constraints_container!( + container, + VoltageMagnitudeConstraint, + T, + String[], + Int[], + Int[]; + sparse = true + ) + vm = get_variable(container, _voltage_magnitude(N), PSY.ACBus) + + time_steps = get_time_steps(container) + jump_model = get_jump_model(container) + for d in devices + for (i, circuit) in enumerate(PSY.get_circuits(d)) + _voltage_controlled(device_model, circuit) || continue + circuit_name = _circuit_arc_name(d, circuit, i) + + bus = PSY.get_bus(sys, PSY.get_regulated_bus_number(circuit)) + bus_name = PSY.get_name(bus) + bus_limits = PSY.get_voltage_limits(bus) + ctl_limits = PSY.get_controlled_quantity_limits(circuit) + + # TODO: temporary pending PSY#1755 + (bus_limits.min <= ctl_limits.min <= ctl_limits.max <= bus_limits.max) || error("Bus limits for $bus_name disagree with control limits for circuit $circuit_name.") + + for t in time_steps + cons[circuit_name, 1, t] = JuMP.@constraint(jump_model, vm[bus_name, t] >= ctl_limits.min) + cons[circuit_name, 2, t] = JuMP.@constraint(jump_model, vm[bus_name, t] <= ctl_limits.max) + end + end + end + return +end + +_add_voltage_control_constraints!( + ::OptimizationContainer, + ::PSY.System, + ::IS.FlattenIteratorWrapper{T}, + ::DeviceModel{T}, + ::NetworkModel, +) where {T} = nothing + ################################## LPACCNetworkModel branch constraints ############### # Branch voltage-angle-difference bounds (angmin, angmax). Only Line / MonitoredLine @@ -1909,6 +1991,22 @@ end ################################## IVRNetworkModel branch constraints ################## +_branch_arc(d::PSY.ACTransmission) = PSY.get_arc(d) +_branch_arc(d::PSY.TwoWindingTransformer) = PSY.get_arc(PSY.get_circuit(d)) + +function _min_endpoint_voltage_limit(branch::PSY.ACTransmission) + arc = _branch_arc(branch) + # bus voltage limits are already per-unit + vmin_fr = PSY.get_voltage_limits(PSY.get_from(arc)).min + vmin_to = PSY.get_voltage_limits(PSY.get_to(arc)).min + return min(vmin_fr, vmin_to) +end + +# Series segments may themselves be parallel groups; recursion bottoms out at devices. +function _min_endpoint_voltage_limit(entry::PNM.AbstractReductionAggregate) + return minimum(_min_endpoint_voltage_limit(member) for member in entry) +end + # Compute the per-unit current rating bound for an IVR branch variable. # c_rating_a = rate_a / vmin (system-base power / per-unit voltage → per-unit current). function _ivr_current_rating(branch::PSY.ACTransmission) @@ -1946,22 +2044,6 @@ function _ivr_current_rating( return rate_a / vmin end -_branch_arc(d::PSY.ACTransmission) = PSY.get_arc(d) -_branch_arc(d::PSY.TwoWindingTransformer) = PSY.get_arc(PSY.get_circuit(d)) - -function _min_endpoint_voltage_limit(branch::PSY.ACTransmission) - arc = _branch_arc(branch) - # bus voltage limits are already per-unit - vmin_fr = PSY.get_voltage_limits(PSY.get_from(arc)).min - vmin_to = PSY.get_voltage_limits(PSY.get_to(arc)).min - return min(vmin_fr, vmin_to) -end - -# Series segments may themselves be parallel groups; recursion bottoms out at devices. -function _min_endpoint_voltage_limit(entry::PNM.AbstractReductionAggregate) - return minimum(_min_endpoint_voltage_limit(member) for member in entry) -end - function add_variables!( container::OptimizationContainer, ::Type{V}, diff --git a/src/ac_transmission_models/branch_constructor.jl b/src/ac_transmission_models/branch_constructor.jl index 39326da..9a3c48b 100644 --- a/src/ac_transmission_models/branch_constructor.jl +++ b/src/ac_transmission_models/branch_constructor.jl @@ -288,7 +288,7 @@ function construct_device!( devices = get_available_components(device_model, sys) _add_static_branch_flow_variables!(container, devices, network_model) _add_static_branch_balance_arguments!(container, device_model, devices, network_model) - _add_transformer_control_variables!(container, device_model, devices, network_model) + _add_tap_control_variables!(container, device_model, devices, network_model) return end @@ -320,7 +320,7 @@ function construct_device!( add_constraints!( container, sys, AngleDifferenceConstraint, devices, device_model, network_model, ) - _add_transformer_control_constraints!(container, device_model, devices, network_model) + _add_voltage_control_constraints!(container, sys, devices, device_model, network_model) add_feedforward_constraints!(container, device_model, devices) add_to_objective_function!(container, devices, device_model, ACPNetworkModel) add_constraint_dual!(container, sys, device_model) diff --git a/src/network_models/instantiate_network_model.jl b/src/network_models/instantiate_network_model.jl index 0ac8002..a5627c2 100644 --- a/src/network_models/instantiate_network_model.jl +++ b/src/network_models/instantiate_network_model.jl @@ -222,28 +222,16 @@ function _monitored_lines_irreducible_buses(branch_models::BranchModelContainer) return irreducible_buses end -_is_controllable( - m::Union{ - DeviceModel{PSY.TwoWindingTransformer}, - DeviceModel{PSY.ThreeWindingTransformer}, - }, -) = get_attribute(m, ENABLE_CONTROLS_KEY) === true -_is_controllable(c::PSY.TransformerCircuit) = - !( - get_control_objective(c) in - (TransformerControlObjective.UNDEFINED, TransformerControlObjective.FIXED) - ) -_is_controllable(_) = false - -# Pin buses with transformers with non-default control +# Pin buses with transformers with non-default control and their regulated buses. function _controllable_transformers_irreducbile_buses(branch_models::BranchModelContainer) irreducible_buses = Set{Int64}() - for m in values(branch_models) - _is_controllable(m) || continue - for transformer in get_device_cache(m) + for model in values(branch_models) + _control_enabled(model) || continue + for transformer in get_device_cache(model) for circuit in PSY.get_circuits(transformer) - _is_controllable(m) || continue - _push_component_buses!(irreducible_buses, branch) + _control_enabled(circuit) || continue + _push_component_buses!(irreducible_buses, transformer) + push!(irreducible_buses, PSY.get_regulated_bus_number(circuit)) end end end diff --git a/test/Project.toml b/test/Project.toml index 89209fd..0b8e83a 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -17,7 +17,6 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -PowerFlows = "94fada2c-fd9a-4e89-8d82-81405f5cb4f6" PowerNetworkMatrices = "bed98974-b02a-5e2f-9fe0-a103f5c450dd" PowerOperationsModels = "bed98974-b02a-5e2f-9ee0-a103f5c450dd" PowerSystemCaseBuilder = "f00506e0-b84f-492a-93c2-c0a9afc4364e" @@ -38,7 +37,6 @@ InfrastructureSystems = {rev = "IS4", url = "https://github.com/Sienna-Platform/ PowerNetworkMatrices = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerNetworkMatrices.jl"} PowerSystemCaseBuilder = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerSystemCaseBuilder.jl"} PowerSystems = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerSystems.jl"} -PowerFlows = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerFlows.jl"} [compat] HiGHS = "1" diff --git a/test/includes.jl b/test/includes.jl index 0fa7e74..92a254f 100644 --- a/test/includes.jl +++ b/test/includes.jl @@ -7,7 +7,7 @@ using InfrastructureSystems import InfrastructureSystems: TableFormat using PowerNetworkMatrices import PowerSystemCaseBuilder: PSITestSystems -using PowerFlows +#using PowerFlows using DataFramesMeta # Test Packages @@ -34,7 +34,7 @@ import LinearAlgebra const PSY = PowerSystems const POM = PowerOperationsModels const IOM = InfrastructureOptimizationModels -const PFS = PowerFlows +#const PFS = PowerFlows const PSB = PowerSystemCaseBuilder const PNM = PowerNetworkMatrices const ISOPT = InfrastructureSystems.Optimization diff --git a/test/test_transformer_controls.jl b/test/test_transformer_controls.jl index 0eaa2d9..90ccc1b 100644 --- a/test/test_transformer_controls.jl +++ b/test/test_transformer_controls.jl @@ -26,7 +26,9 @@ function _controlled_sys14( objective; name = "Trans1", regulated = :to, - quantity_limits = (min = 0.9, max = 1.1), + # c_sys14 buses carry (0.94, 1.06) voltage limits; a VOLTAGE band has to sit inside + # the regulated bus's own limits. + quantity_limits = (min = 0.95, max = 1.05), control_limits = (min = 0.9, max = 1.1), ) sys = PSB.build_system(PSITestSystems, "c_sys14") @@ -42,7 +44,8 @@ function _controlled_sys14( end function _controlled_template(network_formulation; enable = true, kwargs...) - template = get_thermal_dispatch_template_network(NetworkModel(network_formulation; kwargs...)) + template = + get_thermal_dispatch_template_network(NetworkModel(network_formulation; kwargs...)) set_device_model!( template, DeviceModel( @@ -148,7 +151,12 @@ end function _uncontrolled_voltage(bus_name; network_formulation = ACPNetworkModel) sys = PSB.build_system(PSITestSystems, "c_sys14") model, status = - _build_controlled(sys, network_formulation; enable = false, optimizer = ipopt_optimizer) + _build_controlled( + sys, + network_formulation; + enable = false, + optimizer = ipopt_optimizer, + ) @test status == IOM.ModelBuildStatus.BUILT @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED res = IOM.OptimizationProblemOutputs(model) @@ -160,7 +168,9 @@ end _, _, _, bus_name = _controlled_sys14(VOLTAGE_CONTROL) free_vm = _uncontrolled_voltage(bus_name) # A band the free-running solution violates, so holding it requires the tap to move. - band = (min = free_vm + 0.01, max = free_vm + 0.02) + # It sits below `free_vm`: the free-running voltage rides near the bus's 1.06 upper + # limit, and the band may not reach outside the bus's own limits. + band = (min = free_vm - 0.02, max = free_vm - 0.01) sys, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL; quantity_limits = band) model, status = _build_controlled(sys, ACPNetworkModel; optimizer = ipopt_optimizer) @@ -189,7 +199,7 @@ end @testset "VOLTAGE control regulates the from-side bus when its number is given" begin _, _, _, bus_name = _controlled_sys14(VOLTAGE_CONTROL; regulated = :from) free_vm = _uncontrolled_voltage(bus_name) - band = (min = free_vm + 0.01, max = free_vm + 0.02) + band = (min = free_vm - 0.02, max = free_vm - 0.01) sys, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL; regulated = :from, quantity_limits = band) @@ -205,22 +215,6 @@ end end end -@testset "the VOLTAGE band is a voltage, not a squared voltage" begin - # A band well away from 1.0 pu separates the two readings: v ∈ [0.80, 0.82] is - # satisfied by v² ∈ [0.64, 0.67], so a builder comparing the raw band against v² - # would land the voltage near 0.9 pu instead. - band = (min = 0.80, max = 0.82) - sys, _, _, bus_name = _controlled_sys14(VOLTAGE_CONTROL; quantity_limits = band) - model, status = _build_controlled(sys, ACPNetworkModel; optimizer = ipopt_optimizer) - @test status == IOM.ModelBuildStatus.BUILT - @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED - - res = IOM.OptimizationProblemOutputs(model) - vm = read_variable(res, "VoltageMagnitude__ACBus"; table_format = TableFormat.WIDE) - @test vm[1, bus_name] >= band.min - 1e-6 - @test vm[1, bus_name] <= band.max + 1e-6 -end - ############################ REACTIVE_POWER_FLOW objective ############################# @testset "REACTIVE_POWER_FLOW control holds the terminal flow inside its band (ACP)" begin @@ -248,39 +242,16 @@ end ################################### model invariants ################################### -@testset "the model is count-invariant across control objectives (ACP)" begin - function _container_for(objective) - sys, _, _, _ = _controlled_sys14(objective) - model, status = _build_controlled(sys, ACPNetworkModel; optimizer = ipopt_optimizer) - @test status == IOM.ModelBuildStatus.BUILT - return IOM.get_optimization_container(model) - end - - cv = _container_for(VOLTAGE_CONTROL) - cq = _container_for(Q_FLOW_CONTROL) - - var_v = IOM.get_variables(cv) - var_q = IOM.get_variables(cq) - @test Set(keys(var_v)) == Set(keys(var_q)) - for k in keys(var_v) - @test size(var_v[k]) == size(var_q[k]) - end - - con_v = IOM.get_constraints(cv) - con_q = IOM.get_constraints(cq) - @test Set(keys(con_v)) == Set(keys(con_q)) - for k in keys(con_v) - @test size(con_v[k]) == size(con_q[k]) - end -end - @testset "a tap pinned at nominal reproduces the uncontrolled model (ACP)" begin # White-box reduction gate: with the tap variable fixed at the circuit's nominal # ratio and a band too wide to bind, the controlled Ohm's law is term-by-term the # fixed-tap one, so both models must reach the same optimum and terminal flows. - wide = (min = 0.5, max = 1.5) + # The band is the regulated bus's own (0.94, 1.06) limits — the widest a VOLTAGE band + # may be — so the control constrains nothing the bus does not already. + band = (min = 0.94, max = 1.06) + tap_range = (min = 0.5, max = 1.5) - sys_fixed, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL; quantity_limits = wide) + sys_fixed, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL; quantity_limits = band) model_fixed, status_fixed = _build_controlled( sys_fixed, ACPNetworkModel; enable = false, optimizer = ipopt_optimizer, ) @@ -288,7 +259,7 @@ end @test solve!(model_fixed) == IOM.RunStatus.SUCCESSFULLY_FINALIZED sys_var, transformer, circuit, _ = _controlled_sys14( - VOLTAGE_CONTROL; quantity_limits = wide, control_limits = wide, + VOLTAGE_CONTROL; quantity_limits = band, control_limits = tap_range, ) model_var, status_var = _build_controlled(sys_var, ACPNetworkModel; optimizer = ipopt_optimizer) @@ -331,7 +302,11 @@ end @test status == IOM.ModelBuildStatus.BUILT container = IOM.get_optimization_container(model) - pft = IOM.get_variable(container, FlowActivePowerFromToVariable, PSY.TwoWindingTransformer) + pft = IOM.get_variable( + container, + FlowActivePowerFromToVariable, + PSY.TwoWindingTransformer, + ) vm = IOM.get_variable(container, VoltageMagnitude, PSY.ACBus) va = IOM.get_variable(container, VoltageAngle, PSY.ACBus) tap = IOM.get_variable(container, TapRatioVariable, PSY.TwoWindingTransformer) @@ -354,7 +329,7 @@ end ) lookup = z -> vals[z] - y = POM._tapped_admittance(adm, vals[tap[name, t]]) + y = POM._tapped_admittance(get_jump_model(container), adm, vals[tap[name, t]]) vmf = vals[vm[fr, t]] vmt = vals[vm[to, t]] θ = vals[va[fr, t]] - vals[va[to, t]] @@ -389,7 +364,11 @@ end ) magnitude = 1.0 + phi[1, bus_name] else - vr = read_variable(res, "VoltageReal__ACBus"; table_format = TableFormat.WIDE) + vr = read_variable( + res, + "VoltageReal__ACBus"; + table_format = TableFormat.WIDE, + ) vi = read_variable( res, "VoltageImaginary__ACBus"; table_format = TableFormat.WIDE, ) @@ -434,21 +413,26 @@ end @testset "a controlled circuit merged with a parallel branch fails with a clear error" begin # PNM collapses parallel branches onto one equivalent arc before POM sees them, which # would leave the control acting on a flow that is not the transformer's own. - sys, transformer, _, _ = _controlled_sys14(VOLTAGE_CONTROL) - arc = PSY.get_arc(PSY.get_circuit(transformer)) + # Trans1's arc spans a voltage change, so the parallel branch has to be another + # transformer: PSY rejects a Line whose endpoints differ in base voltage. + sys, transformer, circuit, _ = _controlled_sys14(VOLTAGE_CONTROL) + arc = PSY.get_arc(circuit) PSY.add_component!( sys, - PSY.Line(; + PSY.TwoWindingTransformer(; name = "parallel_to_Trans1", - available = true, - active_power_flow = 0.0, - reactive_power_flow = 0.0, - arc = arc, - r = 0.01, - x = 0.1, - b = (from = 0.0, to = 0.0), - rating = 2.0, - angle_limits = (min = -π / 2, max = π / 2), + circuit = PSY.TransformerCircuit(; + available = true, + arc = arc, + r = PSY.get_r(circuit, PSY.SU), + x = PSY.get_x(circuit, PSY.SU), + tap = 1.0, + α = 0.0, + rating = PSY.get_rating(circuit, PSY.SU), + base_power = PSY.get_base_power(sys, PSY.NU), + ), + magnetizing_shunt = 0.0 + 0.0im, + shunt_location = PSY.TwoWindingTransformerShuntLocation.PRIMARY, ), ) template = _controlled_template(ACPNetworkModel) @@ -476,65 +460,88 @@ function _case11_with_forecast() return sys end -@testset "a controlled circuit absorbed by a series reduction fails with a clear error" begin - # "1-6-i_1" is one segment of the (1,2) series chain, so under reduction it has no - # direct-branch entry of its own and its tap would control an equivalent arc it only - # partly owns. - sys = _case11_with_forecast() - line = PSY.get_component(PSY.Line, sys, "1-6-i_1") - arc = PSY.get_arc(line) - transformer = PSY.TwoWindingTransformer(; - name = PSY.get_name(line), - circuit = PSY.TransformerCircuit(; - available = true, - arc = arc, - r = PSY.get_r(line, PSY.SU), - x = PSY.get_x(line, PSY.SU), - tap = 1.0, - α = 0.0, - rating = PSY.get_rating(line, PSY.SU), - control_objective = VOLTAGE_CONTROL, - regulated_bus_number = PSY.get_number(PSY.get_to(arc)), - control_limits = (min = 0.9, max = 1.1), - controlled_quantity_limits = (min = 0.95, max = 1.05), - base_power = PSY.get_base_power(sys, PSY.NU), - ), - magnetizing_shunt = 0.0 + 0.0im, - shunt_location = PSY.TwoWindingTransformerShuntLocation.PRIMARY, - ) - PSY.remove_component!(sys, line) - PSY.add_component!(sys, transformer) +@testset "StaticBranch models transformer off-nominal tap under DCP (c_sys14)" begin + sys = PSB.build_system(PSITestSystems, "c_sys14") + template = get_thermal_dispatch_template_network(NetworkModel(DCPNetworkModel)) + set_device_model!(template, PSY.TwoWindingTransformer, StaticBranch) + model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED - template = _controlled_template( - ACPNetworkModel; - reduce_radial_branches = true, - reduce_degree_two_branches = true, - ) - model = DecisionModel(template, sys; optimizer = ipopt_optimizer) - out = mktempdir(; cleanup = true) - @test build!(model; output_dir = out, console_level = Logging.Error) == - IOM.ModelBuildStatus.FAILED - @test occursin( - "Controlled transformer circuit", - read(joinpath(out, "operation_problem.log"), String), + res = IOM.OptimizationProblemOutputs(model) + base = IOM.get_model_base_power(res) + # StaticBranch under DCP has no FlowActivePowerVariable: the flow IS the + # BThetaBranchFlow expression, reported in natural units (MW). VoltageAngle is + # unitless (radians, no conversion), so compare in per-unit. + pflow = read_expression( + res, + "BThetaBranchFlow__TwoWindingTransformer"; + table_format = TableFormat.WIDE, ) + va = read_variable(res, "VoltageAngle__ACBus"; table_format = TableFormat.WIDE) + + tested_a_real_tap = false + for tr in PSY.get_components(PSY.TwoWindingTransformer, sys) + name = PSY.get_name(tr) + @test name in names(pflow) + + # Recover the series reactance independently, from the π-model admittance, so the + # oracle does not simply re-call the susceptance helper the source uses. + adm = PNM.branch_admittance(tr) + x = -adm.b / (adm.g^2 + adm.b^2) + # The DC susceptance is tap-divided: b_dc == 1/(tap*x). Pin the equivalence of the + # independent recovery and PNM's DC entry point. + @test 1 / (x * adm.tap) ≈ PNM.get_series_susceptance(tr, PSY.SU) + + arc = PSY.get_arc(PSY.get_circuit(tr)) + fr = PSY.get_name(PSY.get_from(arc)) + to = PSY.get_name(PSY.get_to(arc)) + shift = PNM.get_series_phase_shift(tr) + if !isapprox(adm.tap, 1.0; atol = 1e-6) + tested_a_real_tap = true + end + for r in 1:nrow(pflow) + p_pu = pflow[r, name] / base + expected = (va[r, fr] - va[r, to] - shift) / (x * adm.tap) + @test isapprox(p_pu, expected; atol = 1e-5) + end + end + # Guard: the test system must actually carry a non-unit tap, else this proves nothing. + @test tested_a_real_tap end -@testset "ACP rejects two transformers regulating the same bus" begin - # Trans1 (4 → 9) and Trans3 (4 → 7) share bus 4. Under ACP both would drive the one - # shared VoltageMagnitude, so the conflicting bands must be caught at validation. +@testset "_tapped_admittance round-trips PNM.ybus_branch_entries" begin + function check_terms(y, ybus) + Y11, Y12, Y21, Y22 = ybus + @test isapprox(complex(y.g11, y.b11), Y11; rtol = 1e-10, atol = 1e-12) + @test isapprox(complex(y.g12, y.b12), Y12; rtol = 1e-10, atol = 1e-12) + @test isapprox(complex(y.g21, y.b21), Y21; rtol = 1e-10, atol = 1e-12) + @test isapprox(complex(y.g22, y.b22), Y22; rtol = 1e-10, atol = 1e-12) + end + + model = JuMP.Model() sys = PSB.build_system(PSITestSystems, "c_sys14") - for name in ("Trans1", "Trans3") - circuit = PSY.get_circuit(PSY.get_component(PSY.TwoWindingTransformer, sys, name)) - PSY.set_control_objective!(circuit, VOLTAGE_CONTROL) - PSY.set_regulated_bus_number!(circuit, PSY.get_number(PSY.get_from(PSY.get_arc(circuit)))) - PSY.set_controlled_quantity_limits!(circuit, (min = 0.98, max = 1.02)) + for br in Iterators.flatten(( + PSY.get_components(PSY.Line, sys), + PSY.get_components(PSY.TwoWindingTransformer, sys), + )) + adm = PNM.branch_admittance(br) + check_terms(POM._tapped_admittance(model, adm, adm.tap), PNM.ybus_branch_entries(br)) end - template = _controlled_template(ACPNetworkModel) - model = DecisionModel(template, sys; optimizer = ipopt_optimizer) - @test_throws IS.ConflictingInputsError POM.validate_template(model) -end + tr = PSY.get_component(PSY.TwoWindingTransformer, sys, "Trans1") + circuit = PSY.get_circuit(tr) + for shift in (-pi / 5, 0.0, pi / 6) + PSY.set_α!(circuit, shift) + PSY.set_tap!(circuit, 1.0) + adm = PNM.branch_admittance(tr) + for tap in (0.9, 1.0, 1.1, 1.25) + PSY.set_tap!(circuit, tap) + check_terms(POM._tapped_admittance(model, adm, tap), PNM.ybus_branch_entries(tr)) + end + end +end ######################################################################################### # Phase control (the ACTIVE_POWER_FLOW / ASYMMETRIC_ACTIVE_POWER_FLOW objectives, where the # phase shift α rather than the tap ratio is the decision variable) is NOT supported yet. From 30c0bfbe2ba2584ea517b17615c28b5797267e3c Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Fri, 14 Aug 2026 23:08:41 -0400 Subject: [PATCH 09/11] claude finished the impl --- src/PowerOperationsModels.jl | 1 + src/ac_transmission_models/AC_branches.jl | 248 ++++++++---- .../branch_constructor.jl | 19 +- src/core/constraints.jl | 7 + src/core/network_formulations.jl | 23 -- .../instantiate_network_model.jl | 4 +- test/test_transformer_controls.jl | 375 ++++++++++-------- test/test_transformer_fixed_tap.jl | 11 +- 8 files changed, 428 insertions(+), 260 deletions(-) diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl index e29654a..a73e8f2 100644 --- a/src/PowerOperationsModels.jl +++ b/src/PowerOperationsModels.jl @@ -770,6 +770,7 @@ export ActiveRangeICConstraint export NodalBalanceActiveConstraint export ReferenceBusConstraint export VoltageMagnitudeConstraint +export ReactivePowerFlowControlConstraint export RegulatedVoltageMagnitudeConstraint export CurrentLimitConstraint export AngleDifferenceConstraint diff --git a/src/ac_transmission_models/AC_branches.jl b/src/ac_transmission_models/AC_branches.jl index a584c06..07d5900 100644 --- a/src/ac_transmission_models/AC_branches.jl +++ b/src/ac_transmission_models/AC_branches.jl @@ -69,24 +69,37 @@ function get_default_time_series_names( return Dict{Type{<:TimeSeriesParameter}, String}() end -ENABLE_CONTROLS_KEY = "enable_controls" +const ENABLE_CONTROLS_KEY = "enable_controls" _control_attribute( ::Union{Type{PSY.TwoWindingTransformer}, Type{PSY.ThreeWindingTransformer}}, ) = (ENABLE_CONTROLS_KEY => false,) _control_attribute(_) = () -_control_enabled(m::Union{DeviceModel{PSY.TwoWindingTransformer}, DeviceModel{PSY.ThreeWindingTransformer}}) = get_attribute(m, ENABLE_CONTROLS_KEY) === true -_control_enabled(c::PSY.TransformerCircuit) = PSY.get_available(c) && !(PSY.get_control_objective(c) in (PSY.TransformerControlObjective.UNDEFINED, PSY.TransformerControlObjective.FIXED)) +_TRANSFORMERS = Union{PSY.TwoWindingTransformer, PSY.ThreeWindingTransformer} + +_control_enabled(m::DeviceModel{<:_TRANSFORMERS}) = + get_attribute(m, ENABLE_CONTROLS_KEY) === true +_control_enabled(c::PSY.TransformerCircuit) = + PSY.get_available(c) && !(PSY.get_control_objective(c) in (PSY.TransformerControlObjective.UNDEFINED, PSY.TransformerControlObjective.FIXED)) _control_enabled(_) = false -_tap_controlled(c::PSY.TransformerControlObjective) = c in (PSY.TransformerControlObjective.VOLTAGE, PSY.TransformerControlObjective.REACTIVE_POWER_FLOW) -_tap_controlled(c::PSY.TransformerCircuit) = PSY.get_available(c) && _tap_controlled(PSY.get_control_objective(c)) +_tap_controlled(c::PSY.TransformerControlObjective) = c in ( + PSY.TransformerControlObjective.VOLTAGE, + PSY.TransformerControlObjective.REACTIVE_POWER_FLOW, +) +_tap_controlled(c::PSY.TransformerCircuit) = + PSY.get_available(c) && _tap_controlled(PSY.get_control_objective(c)) -_voltage_controlled(c::PSY.TransformerControlObjective) = c === PSY.TransformerControlObjective.VOLTAGE -_voltage_controlled(c::PSY.TransformerCircuit) = PSY.get_available(c) && _voltage_controlled(PSY.get_control_objective(c)) +_voltage_controlled(c::PSY.TransformerControlObjective) = + c === PSY.TransformerControlObjective.VOLTAGE +_voltage_controlled(c::PSY.TransformerCircuit) = + PSY.get_available(c) && _voltage_controlled(PSY.get_control_objective(c)) -_reactive_controlled(c::PSY.TransformerControlObjective) = c === PSY.TransformerControlObjective.REACTIVE_POWER_FLOW +_reactive_controlled(c::PSY.TransformerControlObjective) = + c === PSY.TransformerControlObjective.REACTIVE_POWER_FLOW +_reactive_controlled(c::PSY.TransformerCircuit) = + PSY.get_available(c) && _reactive_controlled(PSY.get_control_objective(c)) _tap_controlled(m::DeviceModel, d) = _control_enabled(m) && _tap_controlled(d) _voltage_controlled(m::DeviceModel, d) = _control_enabled(m) && _voltage_controlled(d) @@ -309,6 +322,10 @@ _add_tap_control_variables!( ::NetworkModel, ) = nothing +_warn_tap_control_nonconvexity(::NetworkModel{N}) where {N <: Union{LPACCNetworkModel, DCPNetworkModel, DCPLLNetworkModel}} = + @warn "Tap control makes $N network models non-convex. Use Ipopt or change circuit controls." +_warn_tap_control_nonconvexity(_) = nothing + function _add_tap_control_variables!( container::OptimizationContainer, model::DeviceModel{U, F}, @@ -319,6 +336,8 @@ function _add_tap_control_variables!( F <: AbstractBranchFormulation, } get_attribute(model, ENABLE_CONTROLS_KEY) === true || return + _warn_tap_control_nonconvexity(network_model) + names = String[] circuits = PSY.TransformerCircuit[] for d in devices, (i, c) in enumerate(PSY.get_circuits(d)) @@ -1220,8 +1239,11 @@ _get_circuit(_) = nothing _control_objective(branch) = _control_objective(_get_circuit(branch)) _control_objective(::Nothing) = PSY.TransformerControlObjective.UNDEFINED _control_objective(c::PSY.TransformerCircuit) = - PSY.get_available(c) ? PSY.get_control_objective(c) : - PSY.TransformerControlObjective.UNDEFINED + if PSY.get_available(c) + PSY.get_control_objective(c) + else + PSY.TransformerControlObjective.UNDEFINED + end _quantity_limits(branch) = _quantity_limits(_get_circuit(branch)) _quantity_limits(::Nothing) = (min = -Inf, max = Inf) @@ -1248,7 +1270,6 @@ Base.@kwdef struct BranchGeometry control::PSY.TransformerControlObjective quantity_limits::MinMax regulated_number::Int - end function BranchGeometry( @@ -1747,7 +1768,11 @@ function add_constraints!( for t in time_steps vp = _voltage_products(container, network_model, T, name, from_bus, to_bus, t) - tap = _tap_controlled(device_model, g_geom) ? get_variable(container, TapRatioVariable, T)[name, t] : adm.tap + tap = if _tap_controlled(device_model, g_geom) + get_variable(container, TapRatioVariable, T)[name, t] + else + adm.tap + end y = _tapped_admittance(jump_model, adm, tap) cons_pft[name, t] = JuMP.@constraint( @@ -1762,42 +1787,41 @@ function add_constraints!( y.g22 * vp.v2_to + y.g21 * vp.vv_cos - y.b21 * vp.vv_sin + _slack_term(slacks.p_tf, name, t), ) - qft_expr = JuMP.@expression( + cons_qft[name, t] = JuMP.@constraint( jump_model, + qft[name, t] == -y.b11 * vp.v2_fr - y.b12 * vp.vv_cos + y.g12 * vp.vv_sin + _slack_term(slacks.q_ft, name, t), ) - qtf_expr = JuMP.@expression( + cons_qtf[name, t] = JuMP.@constraint( jump_model, + qtf[name, t] == -y.b22 * vp.v2_to - y.b21 * vp.vv_cos - y.g21 * vp.vv_sin + _slack_term(slacks.q_tf, name, t), ) - cons_qft[name, t] = JuMP.@constraint(jump_model, qft[name, t] == qft_expr) - cons_qtf[name, t] = JuMP.@constraint(jump_model, qtf[name, t] == qtf_expr) - - # TODO: register in containers - if _reactive_controlled(device_model, g_geom) - JuMP.@constraint(jump_model, qft_expr >= g_geom.quantity_limits.min) - JuMP.@constraint(jump_model, qft_expr <= g_geom.quantity_limits.max) - JuMP.@constraint(jump_model, qtf_expr >= g_geom.quantity_limits.min) - JuMP.@constraint(jump_model, qtf_expr <= g_geom.quantity_limits.max) - end end end return end -_voltage_magnitude(::Type{<:Union{ACPNetworkModel, IVRNetworkModel}}) = VoltageMagnitude -_voltage_magnitude(_) = RegulatedVoltageMagnitude +_voltage_magnitude(container, name, ::NetworkModel{ACPNetworkModel}) = + get_variable(container, VoltageMagnitude, PSY.ACBus)[name, :] +_voltage_magnitude(container, name, ::NetworkModel{<:Union{ACRNetworkModel, IVRNetworkModel}}) = + JuMP.@expression(get_jump_model(container), [t in get_time_steps(container)], get_variable(container, VoltageReal, PSY.ACBus)[name, t]^2 + get_variable(container, VoltageImaginary, PSY.ACBus)[name, t]^2) +_voltage_magnitude(container, name, ::NetworkModel{LPACCNetworkModel}) = + get_variable(container, VoltageDeviation, PSY.ACBus)[name, :] + +_voltage_limits(limits, ::NetworkModel{ACPNetworkModel}) = limits +_voltage_limits(limits, ::NetworkModel{<:Union{ACRNetworkModel, IVRNetworkModel}}) = (min = limits.min^2, max = limits.max^2) +_voltage_limits(limits, ::NetworkModel{LPACCNetworkModel}) = (min = limits.min - 1, max = limits.max - 1) -# TODO: It might benefit solvers to reset VariableRef bounds instead of using singleton constraints (or even to use fix()), but then we lose the dual. Worth it? function _add_voltage_control_constraints!( container::OptimizationContainer, sys::PSY.System, devices::IS.FlattenIteratorWrapper{T}, device_model::DeviceModel{T}, - ::NetworkModel{N} -) where {T <: Union{PSY.TwoWindingTransformer, PSY.ThreeWindingTransformer}, N <: AbstractNetworkModel} + ::NetworkModel{<:NativeACNetworkModel} +) where {T <: _TRANSFORMER_CONTROL_TYPES} _control_enabled(device_model) || return cons = add_constraints_container!( @@ -1807,28 +1831,28 @@ function _add_voltage_control_constraints!( String[], Int[], Int[]; - sparse = true + sparse = true, ) - vm = get_variable(container, _voltage_magnitude(N), PSY.ACBus) time_steps = get_time_steps(container) jump_model = get_jump_model(container) for d in devices for (i, circuit) in enumerate(PSY.get_circuits(d)) _voltage_controlled(device_model, circuit) || continue - circuit_name = _circuit_arc_name(d, circuit, i) bus = PSY.get_bus(sys, PSY.get_regulated_bus_number(circuit)) bus_name = PSY.get_name(bus) bus_limits = PSY.get_voltage_limits(bus) ctl_limits = PSY.get_controlled_quantity_limits(circuit) - # TODO: temporary pending PSY#1755 + circuit_name = _circuit_arc_name(d, circuit, i) (bus_limits.min <= ctl_limits.min <= ctl_limits.max <= bus_limits.max) || error("Bus limits for $bus_name disagree with control limits for circuit $circuit_name.") + lims = _voltage_limits(ctl_limits) + vm = _voltage_magnitude(container, bus_name, network_model) for t in time_steps - cons[circuit_name, 1, t] = JuMP.@constraint(jump_model, vm[bus_name, t] >= ctl_limits.min) - cons[circuit_name, 2, t] = JuMP.@constraint(jump_model, vm[bus_name, t] <= ctl_limits.max) + cons[circuit_name, 1, t] = JuMP.@constraint(jump_model, vm[t] >= lims.min) + cons[circuit_name, 2, t] = JuMP.@constraint(jump_model, vm[t] <= lims.max) end end end @@ -1843,6 +1867,68 @@ _add_voltage_control_constraints!( ::NetworkModel, ) where {T} = nothing +function _add_reactive_control_constraints!( + container::OptimizationContainer, + devices::IS.FlattenIteratorWrapper{T}, + device_model::DeviceModel{T}, + ::NetworkModel{<:NativeACNetworkModel}, +) where {T <: _TRANSFORMER_CONTROL_TYPES} + _control_enabled(device_model) || return + + cons = add_constraints_container!( + container, + ReactivePowerFlowControlConstraint, + T, + String[], + Int[], + Int[]; + sparse = true, + ) + qft = get_variable(container, FlowReactivePowerFromToVariable, T) + qtf = get_variable(container, FlowReactivePowerToFromVariable, T) + + time_steps = get_time_steps(container) + jump_model = get_jump_model(container) + for d in devices + for (i, circuit) in enumerate(PSY.get_circuits(d)) + _reactive_controlled(device_model, circuit) || continue + name = _circuit_arc_name(d, circuit, i) + lims = PSY.get_controlled_quantity_limits(circuit) + + for t in time_steps + cons[name, 1, t] = + JuMP.@constraint(jump_model, qft[name, t] >= lims.min) + cons[name, 2, t] = + JuMP.@constraint(jump_model, qft[name, t] <= lims.max) + cons[name, 3, t] = + JuMP.@constraint(jump_model, qtf[name, t] >= lims.min) + cons[name, 4, t] = + JuMP.@constraint(jump_model, qtf[name, t] <= lims.max) + end + end + end + return +end + +_add_reactive_control_constraints!( + ::OptimizationContainer, + ::IS.FlattenIteratorWrapper{T}, + ::DeviceModel{T}, + ::NetworkModel, +) where {T} = nothing + +function _add_transformer_control_constraints!( + container::OptimizationContainer, + sys::PSY.System, + devices::IS.FlattenIteratorWrapper{T}, + device_model::DeviceModel{T}, + network_model::NetworkModel, +) where {T <: PSY.ACTransmission} + _add_voltage_control_constraints!(container, sys, devices, device_model, network_model) + _add_reactive_control_constraints!(container, devices, device_model, network_model) + return +end + ################################## LPACCNetworkModel branch constraints ############### # Branch voltage-angle-difference bounds (angmin, angmax). Only Line / MonitoredLine @@ -2118,6 +2204,7 @@ Ten constraints per branch per time step: (9-10) Ohm's law across series impedance Z = r + jx = 1/(g + jb) (linear): vr_to·tm² = vr_fr·tr + vi_fr·ti - r·csr·tm² + x·csi·tm² vi_to·tm² = vi_fr·tr - vr_fr·ti - r·csi·tm² - x·csr·tm² + """ function add_constraints!( container::OptimizationContainer, @@ -2181,6 +2268,12 @@ function add_constraints!( jump_model = get_jump_model(container) slacks = _flow_equality_slacks(container, device_model, T) cslacks = _current_equality_slacks(container, device_model, T) + tap_var = + if has_container_key(container, TapRatioVariable, T) + get_variable(container, TapRatioVariable, T) + else + nothing + end for g_geom in geoms name = g_geom.name adm = g_geom.adm @@ -2190,20 +2283,20 @@ function add_constraints!( b_fr = adm.b_fr g_to = adm.g_to b_to = adm.b_to - tm = adm.tap from_bus = g_geom.from_name to_bus = g_geom.to_name - tr = tm * cos(adm.shift) - ti = tm * sin(adm.shift) - tm2 = tm^2 - # Series impedance Z = r + jx = conj(y)/|y|² ymag2 = g^2 + b^2 r = g / ymag2 x = -b / ymag2 for t in time_steps + tm = _tap_controlled(device_model, g_geom) ? tap_var[name, t] : adm.tap + tr = tm * cos(adm.shift) + ti = tm * sin(adm.shift) + tm2 = tm^2 + vr_f = vr[from_bus, t] vi_f = vi[from_bus, t] vr_t = vr[to_bus, t] @@ -2488,15 +2581,6 @@ function add_constraints!( return end -""" -Add branch Ohm's law (DC power flow) constraint for ACBranch under DCPNetworkModel: - - p_fr == b * (va_fr - va_to - shift) - -where `b` is the DC series susceptance `1/(a·x)` and `shift` is the DC phase-shift angle -(0 for non-PST branches) — the same pair PNM's `BA_Matrix` and `arc_dc_shift_injection` -use, not the π-recovery `adm.b`/`adm.shift`. -""" function add_constraints!( container::OptimizationContainer, sys::PSY.System, @@ -2517,25 +2601,40 @@ function add_constraints!( container, NetworkFlowConstraint, T, branch_names, time_steps, ) - # StaticBranchBounds relaxes the rating by slacking this defining equality: the bounded - # decision flow `p` stays within rating while the physical angle-implied flow may deviate - # by the signed slack. StaticBranch never reaches this method (it carries flow as the - # BThetaBranchFlow expression); only SBB does, so the slacks exist iff use_slacks. use_slacks = get_use_slacks(device_model) if use_slacks slack_ub = get_variable(container, FlowActivePowerSlackUpperBound, T) slack_lb = get_variable(container, FlowActivePowerSlackLowerBound, T) end - for g in geoms - for t in time_steps - rhs = g.b_dc * (va[g.from_name, t] - va[g.to_name, t] - g.shift_dc) + jump_model = get_jump_model(container) + tap_var = + if has_container_key(container, TapRatioVariable, T) + get_variable(container, TapRatioVariable, T) + else + nothing + end + + for g in geoms, t in time_steps + angle = va[g.from_name, t] - va[g.to_name, t] - g.shift_dc + flow = if use_slacks - rhs += slack_ub[g.name, t] - slack_lb[g.name, t] + JuMP.@expression( + jump_model, + p[g.name, t] - slack_ub[g.name, t] + slack_lb[g.name, t] + ) + else + p[g.name, t] + end + cons[g.name, t] = + if _tap_controlled(device_model, g) + JuMP.@constraint( + jump_model, + flow * tap_var[g.name, t] == g.b_dc * g.adm.tap * angle + ) + else + JuMP.@constraint(jump_model, flow == g.b_dc * angle) end - cons[g.name, t] = - JuMP.@constraint(get_jump_model(container), p[g.name, t] == rhs) - end end return end @@ -2913,13 +3012,6 @@ function add_constraints!( return end -""" -Add the DC Ohm's law for the from-to directional flow under DCPLLNetworkModel: - - p_fr == b * (va_fr - va_to - shift) - -identical to the DCP law; the to-from flow is determined by the quadratic loss constraint. -""" function add_constraints!( container::OptimizationContainer, sys::PSY.System, @@ -2941,14 +3033,24 @@ function add_constraints!( ) jump_model = get_jump_model(container) - for g in geoms - for t in time_steps - cons[g.name, t] = JuMP.@constraint( + tap_var = + if has_container_key(container, TapRatioVariable, T) + get_variable(container, TapRatioVariable, T) + else + nothing + end + + for g in geoms, t in time_steps + angle = va[g.from_name, t] - va[g.to_name, t] - g.shift_dc + cons[g.name, t] = + if _tap_controlled(device_model, g) + JuMP.@constraint( jump_model, - pft[g.name, t] == - g.b_dc * (va[g.from_name, t] - va[g.to_name, t] - g.shift_dc), + pft[g.name, t] * tap_var[g.name, t] == g.b_dc * g.adm.tap * angle ) - end + else + JuMP.@constraint(jump_model, pft[g.name, t] == g.b_dc * angle) + end end return end diff --git a/src/ac_transmission_models/branch_constructor.jl b/src/ac_transmission_models/branch_constructor.jl index 9a3c48b..a2b625d 100644 --- a/src/ac_transmission_models/branch_constructor.jl +++ b/src/ac_transmission_models/branch_constructor.jl @@ -320,7 +320,9 @@ function construct_device!( add_constraints!( container, sys, AngleDifferenceConstraint, devices, device_model, network_model, ) - _add_voltage_control_constraints!(container, sys, devices, device_model, network_model) + _add_transformer_control_constraints!( + container, sys, devices, device_model, network_model, + ) add_feedforward_constraints!(container, device_model, devices) add_to_objective_function!(container, devices, device_model, ACPNetworkModel) add_constraint_dual!(container, sys, device_model) @@ -415,6 +417,7 @@ function construct_device!( devices = get_available_components(device_model, sys) _add_static_branch_flow_variables!(container, devices, network_model) _add_static_branch_balance_arguments!(container, device_model, devices, network_model) + _add_tap_control_variables!(container, device_model, devices, network_model) return end @@ -447,6 +450,9 @@ function construct_device!( add_constraints!( container, sys, AngleDifferenceConstraint, devices, device_model, network_model, ) + _add_transformer_control_constraints!( + container, sys, devices, device_model, network_model, + ) add_feedforward_constraints!(container, device_model, devices) add_to_objective_function!(container, devices, device_model, ACRNetworkModel) add_constraint_dual!(container, sys, device_model) @@ -509,6 +515,7 @@ ArgumentConstructStage for StaticBranch under LPACCNetworkModel. Creates the four directional flow variables, the bus-pair cosine variable (cs), optional slacks, and registers each flow's contribution to the per-bus ActivePowerBalance and ReactivePowerBalance expressions. + """ function construct_device!( container::OptimizationContainer, @@ -523,6 +530,7 @@ function construct_device!( _add_static_branch_flow_variables!(container, devices, network_model) add_variables!(container, CosineApproximation, devices, network_model) _add_static_branch_balance_arguments!(container, device_model, devices, network_model) + _add_tap_control_variables!(container, device_model, devices, network_model) return end @@ -557,6 +565,9 @@ function construct_device!( add_constraints!( container, sys, AngleDifferenceConstraint, devices, device_model, network_model, ) + _add_transformer_control_constraints!( + container, sys, devices, device_model, network_model, + ) add_feedforward_constraints!(container, device_model, devices) add_to_objective_function!(container, devices, device_model, LPACCNetworkModel) add_constraint_dual!(container, sys, device_model) @@ -697,6 +708,7 @@ function construct_device!( network_model, ) end + _add_tap_control_variables!(container, device_model, devices, network_model) add_feedforward_arguments!(container, device_model, devices) return end @@ -734,6 +746,9 @@ function construct_device!( add_constraints!( container, sys, AngleDifferenceConstraint, devices, device_model, network_model, ) + _add_transformer_control_constraints!( + container, sys, devices, device_model, network_model, + ) add_feedforward_constraints!(container, device_model, devices) add_to_objective_function!(container, devices, device_model, IVRNetworkModel) add_constraint_dual!(container, sys, device_model) @@ -1125,6 +1140,7 @@ function construct_device!( container, ActivePowerBalance, FlowActivePowerToFromVariable, devices, device_model, network_model, ) + _add_tap_control_variables!(container, device_model, devices, network_model) add_feedforward_arguments!(container, device_model, devices) return end @@ -1197,6 +1213,7 @@ function construct_device!( device_model, network_model, ) + _add_tap_control_variables!(container, device_model, devices, network_model) add_feedforward_arguments!(container, device_model, devices) return end diff --git a/src/core/constraints.jl b/src/core/constraints.jl index 47d7b95..88d3393 100644 --- a/src/core/constraints.jl +++ b/src/core/constraints.jl @@ -200,6 +200,13 @@ struct ReferenceBusConstraint <: ConstraintType end """Rectangular-coordinate voltage magnitude bounds: vmin² ≤ vr² + vi² ≤ vmax².""" struct VoltageMagnitudeConstraint <: ConstraintType end """ +Terminal reactive-flow band for a transformer circuit whose control objective is +`REACTIVE_POWER_FLOW`. Both directional flows are held inside the circuit's +`controlled_quantity_limits`. Sparse, indexed by (circuit name, side, time step) with +side ∈ 1:4 = (from-to lower, from-to upper, to-from lower, to-from upper). +""" +struct ReactivePowerFlowControlConstraint <: ConstraintType end +""" Ties a component-owned [`RegulatedVoltageMagnitude`](@ref) auxiliary variable to the rectangular voltage components at its regulated bus under ACR/IVR formulations. One entry per regulating device per time step: diff --git a/src/core/network_formulations.jl b/src/core/network_formulations.jl index 18542b7..cbea8d0 100644 --- a/src/core/network_formulations.jl +++ b/src/core/network_formulations.jl @@ -177,26 +177,3 @@ voltage_form(::Type{DCPNetworkModel}) = AngleBasedVoltage() voltage_form(::Type{DCPLLNetworkModel}) = AngleBasedVoltage() voltage_form(::Type{ACPNetworkModel}) = AngleBasedVoltage() voltage_form(::Type{LPACCNetworkModel}) = AngleBasedVoltage() - -# --- How a regulated bus voltage magnitude is pinned by a controlling device --- -# Polar networks carry a scalar VoltageMagnitude that is fixed directly; rectangular -# networks (vr, vi) have no magnitude primitive, so a per-device RegulatedVoltageMagnitude -# aux variable is tied to the components and fixed instead. Selects the objective- -# application path for the voltage-controlling tap (and any future voltage regulator). -abstract type RegulatedVoltageForm end -struct PolarRegulatedVoltage <: RegulatedVoltageForm end -struct RectangularRegulatedVoltage <: RegulatedVoltageForm end - -regulated_voltage_form(::Type{<:AbstractNetworkModel}) = RectangularRegulatedVoltage() -regulated_voltage_form(::Type{ACPNetworkModel}) = PolarRegulatedVoltage() - -# --- Whether a tap branch is built with explicit current variables --- -# IVR carries branch terminal/series current variables (and a CurrentLimitConstraint); -# ACP/ACR model the branch in power only. Selects the tap-branch construction path. -abstract type TapBranchCurrentForm end -struct PowerOnlyTapBranch <: TapBranchCurrentForm end -struct CurrentInjectionTapBranch <: TapBranchCurrentForm end - -tap_branch_current_form(::Type{ACPNetworkModel}) = PowerOnlyTapBranch() -tap_branch_current_form(::Type{ACRNetworkModel}) = PowerOnlyTapBranch() -tap_branch_current_form(::Type{IVRNetworkModel}) = CurrentInjectionTapBranch() diff --git a/src/network_models/instantiate_network_model.jl b/src/network_models/instantiate_network_model.jl index a5627c2..e519547 100644 --- a/src/network_models/instantiate_network_model.jl +++ b/src/network_models/instantiate_network_model.jl @@ -168,7 +168,7 @@ function _outage_irreducible_buses( return irreducible_buses end -_is_ac_transmission(::DeviceModel{PSY.ACTransmission}) = true +_is_ac_transmission(::DeviceModel{<:PSY.ACTransmission}) = true _is_ac_transmission(_) = false _is_3w(::DeviceModel{PSY.ThreeWindingTransformer}) = true @@ -230,7 +230,7 @@ function _controllable_transformers_irreducbile_buses(branch_models::BranchModel for transformer in get_device_cache(model) for circuit in PSY.get_circuits(transformer) _control_enabled(circuit) || continue - _push_component_buses!(irreducible_buses, transformer) + _push_component_buses!(irreducible_buses, circuit) push!(irreducible_buses, PSY.get_regulated_bus_number(circuit)) end end diff --git a/test/test_transformer_controls.jl b/test/test_transformer_controls.jl index 90ccc1b..86b6361 100644 --- a/test/test_transformer_controls.jl +++ b/test/test_transformer_controls.jl @@ -43,22 +43,36 @@ function _controlled_sys14( return sys, transformer, circuit, PSY.get_name(bus) end -function _controlled_template(network_formulation; enable = true, kwargs...) +function _controlled_template( + network_formulation; + enable = true, + formulation = StaticBranch, + kwargs..., +) template = get_thermal_dispatch_template_network(NetworkModel(network_formulation; kwargs...)) set_device_model!( template, DeviceModel( PSY.TwoWindingTransformer, - StaticBranch; + formulation; attributes = _control_attributes(enable), ), ) return template end -function _build_controlled(sys, network_formulation; enable = true, optimizer, kwargs...) - template = _controlled_template(network_formulation; enable = enable, kwargs...) +function _build_controlled( + sys, + network_formulation; + enable = true, + optimizer, + formulation = StaticBranch, + kwargs..., +) + template = _controlled_template( + network_formulation; enable = enable, formulation = formulation, kwargs..., + ) model = DecisionModel(template, sys; optimizer = optimizer) status = build!(model; output_dir = mktempdir(; cleanup = true)) return model, status @@ -242,52 +256,66 @@ end ################################### model invariants ################################### -@testset "a tap pinned at nominal reproduces the uncontrolled model (ACP)" begin +@testset "a tap pinned at nominal reproduces the uncontrolled model" begin # White-box reduction gate: with the tap variable fixed at the circuit's nominal # ratio and a band too wide to bind, the controlled Ohm's law is term-by-term the # fixed-tap one, so both models must reach the same optimum and terminal flows. # The band is the regulated bus's own (0.94, 1.06) limits — the widest a VOLTAGE band # may be — so the control constrains nothing the bus does not already. + # + # IVR is the sharpest case: its law is multiplied through by the ratio, so every + # tm-bearing term has to reduce exactly for the flows to match. band = (min = 0.94, max = 1.06) tap_range = (min = 0.5, max = 1.5) - sys_fixed, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL; quantity_limits = band) - model_fixed, status_fixed = _build_controlled( - sys_fixed, ACPNetworkModel; enable = false, optimizer = ipopt_optimizer, + for network_formulation in ( + ACPNetworkModel, + ACRNetworkModel, + LPACCNetworkModel, + IVRNetworkModel, ) - @test status_fixed == IOM.ModelBuildStatus.BUILT - @test solve!(model_fixed) == IOM.RunStatus.SUCCESSFULLY_FINALIZED - - sys_var, transformer, circuit, _ = _controlled_sys14( - VOLTAGE_CONTROL; quantity_limits = band, control_limits = tap_range, - ) - model_var, status_var = - _build_controlled(sys_var, ACPNetworkModel; optimizer = ipopt_optimizer) - @test status_var == IOM.ModelBuildStatus.BUILT - - container = IOM.get_optimization_container(model_var) - tap = IOM.get_variable(container, TapRatioVariable, PSY.TwoWindingTransformer) - for t in axes(tap, 2) - JuMP.fix(tap[PSY.get_name(transformer), t], PSY.get_tap(circuit); force = true) - end - @test solve!(model_var) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + sys_fixed, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL; quantity_limits = band) + model_fixed, status_fixed = _build_controlled( + sys_fixed, network_formulation; enable = false, + optimizer = ipopt_optimizer, + ) + @test status_fixed == IOM.ModelBuildStatus.BUILT + @test solve!(model_fixed) == IOM.RunStatus.SUCCESSFULLY_FINALIZED - res_fixed = IOM.OptimizationProblemOutputs(model_fixed) - res_var = IOM.OptimizationProblemOutputs(model_var) - @test isapprox( - IOM.get_objective_value(res_var), - IOM.get_objective_value(res_fixed); - rtol = 1e-3, - ) - for key in ( - "FlowActivePowerFromToVariable__TwoWindingTransformer", - "FlowReactivePowerFromToVariable__TwoWindingTransformer", - ) - flow_fixed = read_variable(res_fixed, key; table_format = TableFormat.WIDE) - flow_var = read_variable(res_var, key; table_format = TableFormat.WIDE) - for d in PSY.get_components(PSY.TwoWindingTransformer, sys_var) - name = PSY.get_name(d) - @test isapprox(flow_var[1, name], flow_fixed[1, name]; atol = 1e-3) + sys_var, transformer, circuit, _ = _controlled_sys14( + VOLTAGE_CONTROL; quantity_limits = band, control_limits = tap_range, + ) + model_var, status_var = _build_controlled( + sys_var, network_formulation; optimizer = ipopt_optimizer, + ) + @test status_var == IOM.ModelBuildStatus.BUILT + + container = IOM.get_optimization_container(model_var) + tap = IOM.get_variable(container, TapRatioVariable, PSY.TwoWindingTransformer) + for t in axes(tap, 2) + JuMP.fix( + tap[PSY.get_name(transformer), t], PSY.get_tap(circuit); force = true, + ) + end + @test solve!(model_var) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + res_fixed = IOM.OptimizationProblemOutputs(model_fixed) + res_var = IOM.OptimizationProblemOutputs(model_var) + @test isapprox( + IOM.get_objective_value(res_var), + IOM.get_objective_value(res_fixed); + rtol = 1e-3, + ) + for key in ( + "FlowActivePowerFromToVariable__TwoWindingTransformer", + "FlowReactivePowerFromToVariable__TwoWindingTransformer", + ) + flow_fixed = read_variable(res_fixed, key; table_format = TableFormat.WIDE) + flow_var = read_variable(res_var, key; table_format = TableFormat.WIDE) + for d in PSY.get_components(PSY.TwoWindingTransformer, sys_var) + name = PSY.get_name(d) + @test isapprox(flow_var[1, name], flow_fixed[1, name]; atol = 1e-3) + end end end end @@ -346,49 +374,160 @@ end @testset "VOLTAGE control is wired on every voltage-carrying AC network" begin for network_formulation in (ACRNetworkModel, IVRNetworkModel, LPACCNetworkModel) - @testset "$network_formulation" begin - _, _, _, bus_name = _controlled_sys14(VOLTAGE_CONTROL) - band = (min = 1.00, max = 1.02) - sys, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL; quantity_limits = band) - model, status = - _build_controlled(sys, network_formulation; optimizer = ipopt_optimizer) - @test status == IOM.ModelBuildStatus.BUILT - @test _has_tap_variable(IOM.get_optimization_container(model)) - @test check_variable_bounded(model, TapRatioVariable, PSY.TwoWindingTransformer) - @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED - - res = IOM.OptimizationProblemOutputs(model) - if network_formulation == LPACCNetworkModel - phi = read_variable( - res, "VoltageDeviation__ACBus"; table_format = TableFormat.WIDE, - ) - magnitude = 1.0 + phi[1, bus_name] - else - vr = read_variable( - res, - "VoltageReal__ACBus"; - table_format = TableFormat.WIDE, - ) - vi = read_variable( - res, "VoltageImaginary__ACBus"; table_format = TableFormat.WIDE, - ) - magnitude = sqrt(vr[1, bus_name]^2 + vi[1, bus_name]^2) - end - @test magnitude >= band.min - 1e-4 - @test magnitude <= band.max + 1e-4 + _, _, _, bus_name = _controlled_sys14(VOLTAGE_CONTROL) + band = (min = 1.00, max = 1.02) + sys, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL; quantity_limits = band) + model, status = + _build_controlled(sys, network_formulation; optimizer = ipopt_optimizer) + @test status == IOM.ModelBuildStatus.BUILT + @test _has_tap_variable(IOM.get_optimization_container(model)) + @test check_variable_bounded(model, TapRatioVariable, PSY.TwoWindingTransformer) + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + res = IOM.OptimizationProblemOutputs(model) + if network_formulation == LPACCNetworkModel + phi = read_variable( + res, "VoltageDeviation__ACBus"; table_format = TableFormat.WIDE, + ) + magnitude = 1.0 + phi[1, bus_name] + else + vr = read_variable( + res, + "VoltageReal__ACBus"; + table_format = TableFormat.WIDE, + ) + vi = read_variable( + res, "VoltageImaginary__ACBus"; table_format = TableFormat.WIDE, + ) + magnitude = sqrt(vr[1, bus_name]^2 + vi[1, bus_name]^2) end + @test magnitude >= band.min - 1e-4 + @test magnitude <= band.max + 1e-4 + + # The band is written on the network's own voltage variables, so no per-device + # RegulatedVoltageMagnitude aux is introduced on any of these networks. + container = IOM.get_optimization_container(model) + @test !IOM.has_container_key( + container, RegulatedVoltageMagnitude, PSY.TwoWindingTransformer, + ) + # One controlled circuit, two rows (lower/upper) per time step. + @test length( + IOM.get_constraint( + container, VoltageMagnitudeConstraint, PSY.TwoWindingTransformer, + ), + ) == 2 * length(IOM.get_time_steps(container)) end end -@testset "DC networks carry no tap control" begin - # The DC network has no voltage magnitude or reactive flow to regulate, so an enabled - # control switch has nothing to build; the tap stays the component property that the - # tap-divided DC susceptance already carries. - sys, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL) - model, status = _build_controlled(sys, DCPNetworkModel; optimizer = HiGHS_optimizer) +@testset "REACTIVE_POWER_FLOW control holds the terminal flow inside its band (IVR)" begin + # IVR builds its own current-based flow constraints, so the band has to be applied + # there too and not only on the shared pi-model path. + band = (min = -0.05, max = 0.05) + sys, _, _, _ = _controlled_sys14(Q_FLOW_CONTROL; quantity_limits = band) + model, status = _build_controlled(sys, IVRNetworkModel; optimizer = ipopt_optimizer) @test status == IOM.ModelBuildStatus.BUILT - @test !_has_tap_variable(IOM.get_optimization_container(model)) @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + res = IOM.OptimizationProblemOutputs(model) + base = IOM.get_model_base_power(res) + for key in ( + "FlowReactivePowerFromToVariable__TwoWindingTransformer", + "FlowReactivePowerToFromVariable__TwoWindingTransformer", + ) + flow = read_variable(res, key; table_format = TableFormat.WIDE) + for r in 1:nrow(flow) + @test flow[r, "Trans1"] / base >= band.min - 1e-6 + @test flow[r, "Trans1"] / base <= band.max + 1e-6 + end + end +end + +################################### DC networks ######################################## + +@testset "DC tap control needs StaticBranchBounds; StaticBranch is rejected" begin + # StaticBranch under DCP carries its flow as the BThetaBranchFlow expression, which the + # affine nodal balance cannot hold once the ratio is a variable. The pair is refused up + # front rather than letting the control go silently inert. + sys, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL) + template = _controlled_template(DCPNetworkModel) + model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + out = mktempdir(; cleanup = true) + @test build!(model; output_dir = out, console_level = Logging.Error) == + IOM.ModelBuildStatus.FAILED + @test occursin( + "builds no tap ratio variable", + read(joinpath(out, "operation_problem.log"), String), + ) +end + +@testset "DC networks build a variable tap under StaticBranchBounds" begin + tap_range = (min = 0.9, max = 1.1) + for network_formulation in (DCPNetworkModel, DCPLLNetworkModel) + sys, _, _, _ = + _controlled_sys14(VOLTAGE_CONTROL; control_limits = tap_range) + model, status = _build_controlled( + sys, + network_formulation; + formulation = StaticBranchBounds, + optimizer = ipopt_optimizer, + ) + @test status == IOM.ModelBuildStatus.BUILT + container = IOM.get_optimization_container(model) + tap = IOM.get_variable(container, TapRatioVariable, PSY.TwoWindingTransformer) + @test axes(tap)[1] == ["Trans1"] + @test check_variable_bounded(model, TapRatioVariable, PSY.TwoWindingTransformer) + # Neither quantity exists on a DC network, so no band is built either way. + @test !IOM.has_container_key( + container, VoltageMagnitudeConstraint, PSY.TwoWindingTransformer, + ) + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + end +end + +@testset "the DC tap law reduces to the fixed-tap law at the nominal ratio" begin + # The bilinear DC law is multiplied through by the ratio, so evaluating the built row at + # the nominal ratio must reproduce `nominal * (p - b_dc * (va_fr - va_to - shift))`. + # Checked on the constraint itself rather than by comparing two solves: c_sys14's + # transformer ratings make a fixed-tap StaticBranchBounds DCP model infeasible (that + # formulation enforces the rating as hard variable bounds), so there is no fixed-tap + # reference solution to compare against on this system. + sys, transformer, circuit, _ = _controlled_sys14(VOLTAGE_CONTROL) + model, status = _build_controlled( + sys, DCPNetworkModel; + formulation = StaticBranchBounds, optimizer = ipopt_optimizer, + ) + @test status == IOM.ModelBuildStatus.BUILT + + container = IOM.get_optimization_container(model) + p = IOM.get_variable(container, FlowActivePowerVariable, PSY.TwoWindingTransformer) + va = IOM.get_variable(container, VoltageAngle, PSY.ACBus) + tap = IOM.get_variable(container, TapRatioVariable, PSY.TwoWindingTransformer) + cons = + IOM.get_constraint(container, POM.NetworkFlowConstraint, PSY.TwoWindingTransformer) + + t = 1 + name = PSY.get_name(transformer) + arc = PSY.get_arc(circuit) + fr = PSY.get_name(PSY.get_from(arc)) + to = PSY.get_name(PSY.get_to(arc)) + nominal = PNM.branch_admittance(transformer).tap + b_dc = PNM.get_series_susceptance(transformer, PSY.SU) + shift = PNM.get_series_phase_shift(transformer) + @test !isapprox(nominal, 1.0; atol = 1e-6) + + vals = Dict{JuMP.VariableRef, Float64}( + va[fr, t] => 0.05, va[to, t] => -0.03, + tap[name, t] => nominal, + p[name, t] => 0.7, + ) + lookup = z -> vals[z] + + angle = vals[va[fr, t]] - vals[va[to, t]] - shift + @test isapprox( + JuMP.value(lookup, JuMP.constraint_object(cons[name, t]).func), + nominal * (vals[p[name, t]] - b_dc * angle); + atol = 1e-10, + ) end ################################ reductions and conflicts ############################## @@ -460,88 +599,6 @@ function _case11_with_forecast() return sys end -@testset "StaticBranch models transformer off-nominal tap under DCP (c_sys14)" begin - sys = PSB.build_system(PSITestSystems, "c_sys14") - template = get_thermal_dispatch_template_network(NetworkModel(DCPNetworkModel)) - set_device_model!(template, PSY.TwoWindingTransformer, StaticBranch) - model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) - @test build!(model; output_dir = mktempdir(; cleanup = true)) == - IOM.ModelBuildStatus.BUILT - @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED - - res = IOM.OptimizationProblemOutputs(model) - base = IOM.get_model_base_power(res) - # StaticBranch under DCP has no FlowActivePowerVariable: the flow IS the - # BThetaBranchFlow expression, reported in natural units (MW). VoltageAngle is - # unitless (radians, no conversion), so compare in per-unit. - pflow = read_expression( - res, - "BThetaBranchFlow__TwoWindingTransformer"; - table_format = TableFormat.WIDE, - ) - va = read_variable(res, "VoltageAngle__ACBus"; table_format = TableFormat.WIDE) - - tested_a_real_tap = false - for tr in PSY.get_components(PSY.TwoWindingTransformer, sys) - name = PSY.get_name(tr) - @test name in names(pflow) - - # Recover the series reactance independently, from the π-model admittance, so the - # oracle does not simply re-call the susceptance helper the source uses. - adm = PNM.branch_admittance(tr) - x = -adm.b / (adm.g^2 + adm.b^2) - # The DC susceptance is tap-divided: b_dc == 1/(tap*x). Pin the equivalence of the - # independent recovery and PNM's DC entry point. - @test 1 / (x * adm.tap) ≈ PNM.get_series_susceptance(tr, PSY.SU) - - arc = PSY.get_arc(PSY.get_circuit(tr)) - fr = PSY.get_name(PSY.get_from(arc)) - to = PSY.get_name(PSY.get_to(arc)) - shift = PNM.get_series_phase_shift(tr) - if !isapprox(adm.tap, 1.0; atol = 1e-6) - tested_a_real_tap = true - end - for r in 1:nrow(pflow) - p_pu = pflow[r, name] / base - expected = (va[r, fr] - va[r, to] - shift) / (x * adm.tap) - @test isapprox(p_pu, expected; atol = 1e-5) - end - end - # Guard: the test system must actually carry a non-unit tap, else this proves nothing. - @test tested_a_real_tap -end - -@testset "_tapped_admittance round-trips PNM.ybus_branch_entries" begin - function check_terms(y, ybus) - Y11, Y12, Y21, Y22 = ybus - @test isapprox(complex(y.g11, y.b11), Y11; rtol = 1e-10, atol = 1e-12) - @test isapprox(complex(y.g12, y.b12), Y12; rtol = 1e-10, atol = 1e-12) - @test isapprox(complex(y.g21, y.b21), Y21; rtol = 1e-10, atol = 1e-12) - @test isapprox(complex(y.g22, y.b22), Y22; rtol = 1e-10, atol = 1e-12) - end - - model = JuMP.Model() - sys = PSB.build_system(PSITestSystems, "c_sys14") - for br in Iterators.flatten(( - PSY.get_components(PSY.Line, sys), - PSY.get_components(PSY.TwoWindingTransformer, sys), - )) - adm = PNM.branch_admittance(br) - check_terms(POM._tapped_admittance(model, adm, adm.tap), PNM.ybus_branch_entries(br)) - end - - tr = PSY.get_component(PSY.TwoWindingTransformer, sys, "Trans1") - circuit = PSY.get_circuit(tr) - for shift in (-pi / 5, 0.0, pi / 6) - PSY.set_α!(circuit, shift) - PSY.set_tap!(circuit, 1.0) - adm = PNM.branch_admittance(tr) - for tap in (0.9, 1.0, 1.1, 1.25) - PSY.set_tap!(circuit, tap) - check_terms(POM._tapped_admittance(model, adm, tap), PNM.ybus_branch_entries(tr)) - end - end -end ######################################################################################### # Phase control (the ACTIVE_POWER_FLOW / ASYMMETRIC_ACTIVE_POWER_FLOW objectives, where the # phase shift α rather than the tap ratio is the decision variable) is NOT supported yet. diff --git a/test/test_transformer_fixed_tap.jl b/test/test_transformer_fixed_tap.jl index a9bd882..a363c80 100644 --- a/test/test_transformer_fixed_tap.jl +++ b/test/test_transformer_fixed_tap.jl @@ -69,13 +69,17 @@ end @test isapprox(complex(y.g22, y.b22), Y22; rtol = 1e-10, atol = 1e-12) end + model = JuMP.Model() sys = PSB.build_system(PSITestSystems, "c_sys14") for br in Iterators.flatten(( PSY.get_components(PSY.Line, sys), PSY.get_components(PSY.TwoWindingTransformer, sys), )) adm = PNM.branch_admittance(br) - check_terms(POM._pi_flow_terms(adm, adm.tap), PNM.ybus_branch_entries(br)) + check_terms( + POM._tapped_admittance(model, adm, adm.tap), + PNM.ybus_branch_entries(br), + ) end tr = PSY.get_component(PSY.TwoWindingTransformer, sys, "Trans1") @@ -86,7 +90,10 @@ end adm = PNM.branch_admittance(tr) for tap in (0.9, 1.0, 1.1, 1.25) PSY.set_tap!(circuit, tap) - check_terms(POM._pi_flow_terms(adm, tap), PNM.ybus_branch_entries(tr)) + check_terms( + POM._tapped_admittance(model, adm, tap), + PNM.ybus_branch_entries(tr), + ) end end end From 6259895f500f813e1dc7a4cacca9e5f47c62da88 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Fri, 14 Aug 2026 23:18:21 -0400 Subject: [PATCH 10/11] remove unnecessary test --- src/ac_transmission_models/AC_branches.jl | 14 ++++++-------- test/test_transformer_controls.jl | 16 ---------------- 2 files changed, 6 insertions(+), 24 deletions(-) diff --git a/src/ac_transmission_models/AC_branches.jl b/src/ac_transmission_models/AC_branches.jl index 07d5900..6127457 100644 --- a/src/ac_transmission_models/AC_branches.jl +++ b/src/ac_transmission_models/AC_branches.jl @@ -332,7 +332,7 @@ function _add_tap_control_variables!( devices::IS.FlattenIteratorWrapper{U}, network_model::NetworkModel, ) where { - U <: Union{PSY.TwoWindingTransformer, PSY.ThreeWindingTransformer}, + U <: _TRANSFORMERS, F <: AbstractBranchFormulation, } get_attribute(model, ENABLE_CONTROLS_KEY) === true || return @@ -1231,9 +1231,7 @@ _dc_phase_shift(branch::PNM.AbstractReductionAggregate, nr::PNM.NetworkReduction _dc_phase_shift(branch::PSY.ACTransmission, ::PNM.NetworkReductionData) = PNM.get_series_phase_shift(branch) -_get_circuit( - b::Union{PSY.TwoWindingTransformer, PNM.ThreeWindingTransformerCircuit}, -) = PSY.get_circuit(b) +_get_circuit(b::_TRANSFORMERS) = PSY.get_circuit(b) _get_circuit(_) = nothing _control_objective(branch) = _control_objective(_get_circuit(branch)) @@ -1820,8 +1818,8 @@ function _add_voltage_control_constraints!( sys::PSY.System, devices::IS.FlattenIteratorWrapper{T}, device_model::DeviceModel{T}, - ::NetworkModel{<:NativeACNetworkModel} -) where {T <: _TRANSFORMER_CONTROL_TYPES} + network_model::NetworkModel{<:NativeACNetworkModel} +) where {T <: _TRANSFORMERS} _control_enabled(device_model) || return cons = add_constraints_container!( @@ -1848,7 +1846,7 @@ function _add_voltage_control_constraints!( circuit_name = _circuit_arc_name(d, circuit, i) (bus_limits.min <= ctl_limits.min <= ctl_limits.max <= bus_limits.max) || error("Bus limits for $bus_name disagree with control limits for circuit $circuit_name.") - lims = _voltage_limits(ctl_limits) + lims = _voltage_limits(ctl_limits, network_model) vm = _voltage_magnitude(container, bus_name, network_model) for t in time_steps cons[circuit_name, 1, t] = JuMP.@constraint(jump_model, vm[t] >= lims.min) @@ -1872,7 +1870,7 @@ function _add_reactive_control_constraints!( devices::IS.FlattenIteratorWrapper{T}, device_model::DeviceModel{T}, ::NetworkModel{<:NativeACNetworkModel}, -) where {T <: _TRANSFORMER_CONTROL_TYPES} +) where {T <: _TRANSFORMERS} _control_enabled(device_model) || return cons = add_constraints_container!( diff --git a/test/test_transformer_controls.jl b/test/test_transformer_controls.jl index 86b6361..93b6986 100644 --- a/test/test_transformer_controls.jl +++ b/test/test_transformer_controls.jl @@ -444,22 +444,6 @@ end ################################### DC networks ######################################## -@testset "DC tap control needs StaticBranchBounds; StaticBranch is rejected" begin - # StaticBranch under DCP carries its flow as the BThetaBranchFlow expression, which the - # affine nodal balance cannot hold once the ratio is a variable. The pair is refused up - # front rather than letting the control go silently inert. - sys, _, _, _ = _controlled_sys14(VOLTAGE_CONTROL) - template = _controlled_template(DCPNetworkModel) - model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) - out = mktempdir(; cleanup = true) - @test build!(model; output_dir = out, console_level = Logging.Error) == - IOM.ModelBuildStatus.FAILED - @test occursin( - "builds no tap ratio variable", - read(joinpath(out, "operation_problem.log"), String), - ) -end - @testset "DC networks build a variable tap under StaticBranchBounds" begin tap_range = (min = 0.9, max = 1.1) for network_formulation in (DCPNetworkModel, DCPLLNetworkModel) From dd4ad35e2a19f8ea547ea44ed476e0d036c7663c Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Fri, 14 Aug 2026 23:18:52 -0400 Subject: [PATCH 11/11] formatting --- src/ac_transmission_models/AC_branches.jl | 52 +++++++++++++++-------- 1 file changed, 35 insertions(+), 17 deletions(-) diff --git a/src/ac_transmission_models/AC_branches.jl b/src/ac_transmission_models/AC_branches.jl index 6127457..859dc48 100644 --- a/src/ac_transmission_models/AC_branches.jl +++ b/src/ac_transmission_models/AC_branches.jl @@ -81,7 +81,10 @@ _TRANSFORMERS = Union{PSY.TwoWindingTransformer, PSY.ThreeWindingTransformer} _control_enabled(m::DeviceModel{<:_TRANSFORMERS}) = get_attribute(m, ENABLE_CONTROLS_KEY) === true _control_enabled(c::PSY.TransformerCircuit) = - PSY.get_available(c) && !(PSY.get_control_objective(c) in (PSY.TransformerControlObjective.UNDEFINED, PSY.TransformerControlObjective.FIXED)) + PSY.get_available(c) && !( + PSY.get_control_objective(c) in + (PSY.TransformerControlObjective.UNDEFINED, PSY.TransformerControlObjective.FIXED) + ) _control_enabled(_) = false _tap_controlled(c::PSY.TransformerControlObjective) = c in ( @@ -322,7 +325,9 @@ _add_tap_control_variables!( ::NetworkModel, ) = nothing -_warn_tap_control_nonconvexity(::NetworkModel{N}) where {N <: Union{LPACCNetworkModel, DCPNetworkModel, DCPLLNetworkModel}} = +_warn_tap_control_nonconvexity( + ::NetworkModel{N}, +) where {N <: Union{LPACCNetworkModel, DCPNetworkModel, DCPLLNetworkModel}} = @warn "Tap control makes $N network models non-convex. Use Ipopt or change circuit controls." _warn_tap_control_nonconvexity(_) = nothing @@ -1804,21 +1809,32 @@ end _voltage_magnitude(container, name, ::NetworkModel{ACPNetworkModel}) = get_variable(container, VoltageMagnitude, PSY.ACBus)[name, :] -_voltage_magnitude(container, name, ::NetworkModel{<:Union{ACRNetworkModel, IVRNetworkModel}}) = - JuMP.@expression(get_jump_model(container), [t in get_time_steps(container)], get_variable(container, VoltageReal, PSY.ACBus)[name, t]^2 + get_variable(container, VoltageImaginary, PSY.ACBus)[name, t]^2) +_voltage_magnitude( + container, + name, + ::NetworkModel{<:Union{ACRNetworkModel, IVRNetworkModel}}, +) = + JuMP.@expression( + get_jump_model(container), + [t in get_time_steps(container)], + get_variable(container, VoltageReal, PSY.ACBus)[name, t]^2 + + get_variable(container, VoltageImaginary, PSY.ACBus)[name, t]^2 + ) _voltage_magnitude(container, name, ::NetworkModel{LPACCNetworkModel}) = get_variable(container, VoltageDeviation, PSY.ACBus)[name, :] _voltage_limits(limits, ::NetworkModel{ACPNetworkModel}) = limits -_voltage_limits(limits, ::NetworkModel{<:Union{ACRNetworkModel, IVRNetworkModel}}) = (min = limits.min^2, max = limits.max^2) -_voltage_limits(limits, ::NetworkModel{LPACCNetworkModel}) = (min = limits.min - 1, max = limits.max - 1) +_voltage_limits(limits, ::NetworkModel{<:Union{ACRNetworkModel, IVRNetworkModel}}) = + (min = limits.min^2, max = limits.max^2) +_voltage_limits(limits, ::NetworkModel{LPACCNetworkModel}) = + (min = limits.min - 1, max = limits.max - 1) function _add_voltage_control_constraints!( container::OptimizationContainer, sys::PSY.System, devices::IS.FlattenIteratorWrapper{T}, device_model::DeviceModel{T}, - network_model::NetworkModel{<:NativeACNetworkModel} + network_model::NetworkModel{<:NativeACNetworkModel}, ) where {T <: _TRANSFORMERS} _control_enabled(device_model) || return @@ -1844,7 +1860,9 @@ function _add_voltage_control_constraints!( ctl_limits = PSY.get_controlled_quantity_limits(circuit) # TODO: temporary pending PSY#1755 circuit_name = _circuit_arc_name(d, circuit, i) - (bus_limits.min <= ctl_limits.min <= ctl_limits.max <= bus_limits.max) || error("Bus limits for $bus_name disagree with control limits for circuit $circuit_name.") + (bus_limits.min <= ctl_limits.min <= ctl_limits.max <= bus_limits.max) || error( + "Bus limits for $bus_name disagree with control limits for circuit $circuit_name.", + ) lims = _voltage_limits(ctl_limits, network_model) vm = _voltage_magnitude(container, bus_name, network_model) @@ -2618,18 +2636,18 @@ function add_constraints!( flow = if use_slacks JuMP.@expression( - jump_model, - p[g.name, t] - slack_ub[g.name, t] + slack_lb[g.name, t] - ) + jump_model, + p[g.name, t] - slack_ub[g.name, t] + slack_lb[g.name, t] + ) else p[g.name, t] end cons[g.name, t] = if _tap_controlled(device_model, g) JuMP.@constraint( - jump_model, - flow * tap_var[g.name, t] == g.b_dc * g.adm.tap * angle - ) + jump_model, + flow * tap_var[g.name, t] == g.b_dc * g.adm.tap * angle + ) else JuMP.@constraint(jump_model, flow == g.b_dc * angle) end @@ -3043,9 +3061,9 @@ function add_constraints!( cons[g.name, t] = if _tap_controlled(device_model, g) JuMP.@constraint( - jump_model, - pft[g.name, t] * tap_var[g.name, t] == g.b_dc * g.adm.tap * angle - ) + jump_model, + pft[g.name, t] * tap_var[g.name, t] == g.b_dc * g.adm.tap * angle + ) else JuMP.@constraint(jump_model, pft[g.name, t] == g.b_dc * angle) end