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: 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 76a7a8d..859dc48 100644 --- a/src/ac_transmission_models/AC_branches.jl +++ b/src/ac_transmission_models/AC_branches.jl @@ -69,6 +69,45 @@ function get_default_time_series_names( return Dict{Type{<:TimeSeriesParameter}, String}() end +const ENABLE_CONTROLS_KEY = "enable_controls" + +_control_attribute( + ::Union{Type{PSY.TwoWindingTransformer}, Type{PSY.ThreeWindingTransformer}}, +) = (ENABLE_CONTROLS_KEY => false,) +_control_attribute(_) = () + +_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)) + +_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.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) +_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 @@ -84,6 +123,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 +134,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 +312,62 @@ 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)) + +_add_tap_control_variables!( + ::OptimizationContainer, + ::DeviceModel, + ::IS.FlattenIteratorWrapper, + ::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}, + devices::IS.FlattenIteratorWrapper{U}, + network_model::NetworkModel, +) where { + U <: _TRANSFORMERS, + 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)) + _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) + tap_var = add_variable_container!(container, TapRatioVariable, U, names, time_steps) + for (i, name) in enumerate(names), t in time_steps + 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, + upper_bound = bounds.max + ) + end + return +end + function _add_meta_flow_slack!( container::OptimizationContainer, ::Type{T}, @@ -1113,52 +1204,25 @@ 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 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 @@ -1172,7 +1236,46 @@ _dc_phase_shift(branch::PNM.AbstractReductionAggregate, nr::PNM.NetworkReduction _dc_phase_shift(branch::PSY.ACTransmission, ::PNM.NetworkReductionData) = PNM.get_series_phase_shift(branch) -function _branch_geometry( +_get_circuit(b::_TRANSFORMERS) = 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) = + 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) +_quantity_limits(c::PSY.TransformerCircuit) = PSY.get_controlled_quantity_limits(c) + +_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 +end + +function BranchGeometry( nr::PNM.NetworkReductionData, number_to_name::Dict{Int, String}, name::String, @@ -1181,7 +1284,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], @@ -1192,14 +1295,18 @@ function _branch_geometry( shift_dc = _dc_phase_shift(branch, nr), r_dc = PNM.arc_dc_resistance(nr, arc_tuple), direct = !_is_aggregate(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) @@ -1220,7 +1327,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, @@ -1382,14 +1489,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}, @@ -1411,30 +1510,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 +1651,104 @@ 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, + jump_model = get_jump_model(container) + vm = get_variable(container, VoltageMagnitude, PSY.ACBus) + va = get_variable(container, VoltageAngle, PSY.ACBus) + vmf, vmt = vm[from_bus, t], vm[to_bus, t] + vaf, vat = va[from_bus, t], va[to_bus, t] + return ( + 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)), ) - 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, + jump_model = get_jump_model(container) + vr = get_variable(container, VoltageReal, PSY.ACBus) + vi = get_variable(container, VoltageImaginary, PSY.ACBus) + 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 = 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), ) - 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} + 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_to = phi[from_bus, t], phi[to_bus, t] + return ( + 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(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 = 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, + ) +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) @@ -1694,28 +1766,40 @@ function add_constraints!( 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), + vp = _voltage_products(container, network_model, T, name, from_bus, to_bus, t) + 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( + jump_model, + pft[name, t] == + y.g11 * vp.v2_fr + y.g12 * vp.vv_cos + y.b12 * vp.vv_sin + + _slack_term(slacks.p_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_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_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), ) end @@ -1723,6 +1807,144 @@ function add_constraints!( return 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{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) + +function _add_voltage_control_constraints!( + container::OptimizationContainer, + sys::PSY.System, + devices::IS.FlattenIteratorWrapper{T}, + device_model::DeviceModel{T}, + network_model::NetworkModel{<:NativeACNetworkModel}, +) where {T <: _TRANSFORMERS} + _control_enabled(device_model) || return + + cons = add_constraints_container!( + container, + VoltageMagnitudeConstraint, + T, + String[], + Int[], + Int[]; + sparse = true, + ) + + 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 + + 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, 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) + cons[circuit_name, 2, t] = JuMP.@constraint(jump_model, vm[t] <= lims.max) + end + end + end + return +end + +_add_voltage_control_constraints!( + ::OptimizationContainer, + ::PSY.System, + ::IS.FlattenIteratorWrapper{T}, + ::DeviceModel{T}, + ::NetworkModel, +) where {T} = nothing + +function _add_reactive_control_constraints!( + container::OptimizationContainer, + devices::IS.FlattenIteratorWrapper{T}, + device_model::DeviceModel{T}, + ::NetworkModel{<:NativeACNetworkModel}, +) where {T <: _TRANSFORMERS} + _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 @@ -1869,106 +2091,23 @@ 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] +################################## IVRNetworkModel branch constraints ################## - # 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 +_branch_arc(d::PSY.ACTransmission) = PSY.get_arc(d) +_branch_arc(d::PSY.TwoWindingTransformer) = PSY.get_arc(PSY.get_circuit(d)) - 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 +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 -################################## IVRNetworkModel branch constraints ################## +# 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). @@ -2007,22 +2146,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}, @@ -2097,6 +2220,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, @@ -2160,6 +2284,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 @@ -2169,20 +2299,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] @@ -2467,15 +2597,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, @@ -2496,25 +2617,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 @@ -2769,65 +2905,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 @@ -2951,13 +3028,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, @@ -2979,14 +3049,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( - jump_model, - pft[g.name, t] == - g.b_dc * (va[g.from_name, t] - va[g.to_name, t] - g.shift_dc), - ) + 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] * 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 end return end diff --git a/src/ac_transmission_models/branch_constructor.jl b/src/ac_transmission_models/branch_constructor.jl index 270db4e..a2b625d 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_tap_control_variables!(container, device_model, devices, network_model) return end @@ -319,6 +320,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, ACPNetworkModel) add_constraint_dual!(container, sys, device_model) @@ -413,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 @@ -445,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) @@ -507,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, @@ -521,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 @@ -555,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) @@ -695,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 @@ -732,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) @@ -1123,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 @@ -1195,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/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/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/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/network_models/instantiate_network_model.jl b/src/network_models/instantiate_network_model.jl index cec315a..e519547 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,77 @@ 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 - -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 and their regulated buses. +function _controllable_transformers_irreducbile_buses(branch_models::BranchModelContainer) + irreducible_buses = Set{Int64}() + for model in values(branch_models) + _control_enabled(model) || continue + for transformer in get_device_cache(model) + for circuit in PSY.get_circuits(transformer) + _control_enabled(circuit) || continue + _push_component_buses!(irreducible_buses, circuit) + push!(irreducible_buses, PSY.get_regulated_bus_number(circuit)) + 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 +359,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 +380,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 +520,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/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, 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/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_native_transformer_tap.jl b/test/test_native_transformer_tap.jl deleted file mode 100644 index c7988ca..0000000 --- a/test/test_native_transformer_tap.jl +++ /dev/null @@ -1,114 +0,0 @@ -######################################################################################### -# 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/ -# 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. -######################################################################################### - -@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 "_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 - - # 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 - - # 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 -end 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 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..93b6986 --- /dev/null +++ b/test/test_transformer_controls.jl @@ -0,0 +1,755 @@ +######################################################################################### +# 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, + # 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") + 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, + formulation = StaticBranch, + kwargs..., +) + template = + get_thermal_dispatch_template_network(NetworkModel(network_formulation; kwargs...)) + set_device_model!( + template, + DeviceModel( + PSY.TwoWindingTransformer, + formulation; + attributes = _control_attributes(enable), + ), + ) + return template +end + +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 +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. + # 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) + @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.02, max = free_vm - 0.01) + + 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 + +############################ 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 "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) + + for network_formulation in ( + ACPNetworkModel, + ACRNetworkModel, + LPACCNetworkModel, + IVRNetworkModel, + ) + 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 + + 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 + +@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(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]] + 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) + _, _, _, 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 "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 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 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 ############################## + +@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. + # 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.TwoWindingTransformer(; + name = "parallel_to_Trans1", + 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) + 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 + +######################################################################################### +# 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_transformer_fixed_tap.jl b/test/test_transformer_fixed_tap.jl new file mode 100644 index 0000000..a363c80 --- /dev/null +++ b/test/test_transformer_fixed_tap.jl @@ -0,0 +1,99 @@ +######################################################################################### +# 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`). +# +# 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 + 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 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