From 085eb7703b384c8f309f11faa08f66b950d9f7e8 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Tue, 12 May 2026 14:52:09 -0400 Subject: [PATCH 1/7] Add HVDCTwoTerminalVSC two-terminal HVDC formulation Add two new formulations for PSY.TwoTerminalVSCLine: HVDCTwoTerminalVSC (NLP) and HVDCTwoTerminalVSCBin2 (SOS2 MILP). Both model independent active/reactive power control per terminal under a PQ capability circle, quadratic+linear converter losses, and explicit cable resistance via v_f - v_t = (1/g)*I. Reactive variables/constraints are added only on full AC networks. Refactor shared converter-loss helpers (loss expression builder, abs-value decomposition, linear-loss detection) out of mt_hvdc_models into src/common_models/quadratic_converter_loss.jl so the multi-terminal and two-terminal VSC paths reuse the same primitives. Rename the ConverterPositiveCurrent / NegativeCurrent / CurrentDirection variable types to generic PositiveCurrent / NegativeCurrent / CurrentDirection. Point Project.toml at the InfrastructureOptimizationModels ac/hvdc-vsc branch (which carries the per-device-bounds API used by these formulations) and restore runtests.jl to run the full test list. --- Project.toml | 3 +- src/PowerOperationsModels.jl | 13 +- .../branch_constructor.jl | 177 +++++++++ src/common_models/quadratic_converter_loss.jl | 124 ++++++ src/core/constraints.jl | 33 ++ src/core/formulations.jl | 23 +- src/core/variables.jl | 49 ++- src/mt_hvdc_models/HVDCsystems.jl | 88 +---- src/mt_hvdc_models/hvdcsystems_constructor.jl | 38 +- src/network_models/pm_translator.jl | 26 +- .../TwoTerminalDC_branches.jl | 359 ++++++++++++++++++ test/Project.toml | 2 +- test/runtests.jl | 6 + test/test_device_hvdc.jl | 146 ++++++- 14 files changed, 977 insertions(+), 110 deletions(-) create mode 100644 src/common_models/quadratic_converter_loss.jl diff --git a/Project.toml b/Project.toml index 0a828b7..0c2127c 100644 --- a/Project.toml +++ b/Project.toml @@ -26,7 +26,7 @@ PowerFlows = "94fada2c-0ca5-4b90-a1fb-4bc5b59ccfc7" [sources] InfrastructureSystems = {rev = "IS4", url = "https://github.com/Sienna-Platform/InfrastructureSystems.jl"} PowerSystems = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerSystems.jl"} -InfrastructureOptimizationModels = {rev = "main", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"} +InfrastructureOptimizationModels = {rev = "ac/hvdc-vsc", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"} [extensions] PowerFlowsExt = "PowerFlows" @@ -34,7 +34,6 @@ PowerFlowsExt = "PowerFlows" [compat] Dates = "1" DocStringExtensions = "~0.8, ~0.9" -InfrastructureOptimizationModels = "0.1.0" InfrastructureSystems = "3" InteractiveUtils = "1.11.0" JuMP = "^1.28" diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl index 738cc21..cef1246 100644 --- a/src/PowerOperationsModels.jl +++ b/src/PowerOperationsModels.jl @@ -213,6 +213,7 @@ include("common_models/add_to_expression.jl") include("common_models/add_parameters.jl") include("common_models/make_system_expressions.jl") include("common_models/reserve_range_constraints.jl") +include("common_models/quadratic_converter_loss.jl") # Market bid cost plumbing (PSY orchestration moved out of IOM). Must be included # before device-specific files that reference MBC_TYPES / IEC_TYPES. @@ -506,11 +507,15 @@ export PostContingencyActivePowerReserveDeploymentVariable export DCVoltage export DCLineCurrent export ConverterCurrent -export ConverterPositiveCurrent -export ConverterNegativeCurrent -export ConverterCurrentDirection +export PositiveCurrent +export NegativeCurrent +export CurrentDirection export HVDCFlowDirectionVariable export HVDCLosses +export HVDCFromDCVoltage +export HVDCToDCVoltage +export HVDCReactivePowerFromVariable +export HVDCReactivePowerToVariable # Load Variables export ShiftUpActivePowerVariable @@ -758,6 +763,8 @@ export HVDCTwoTerminalLossless export HVDCTwoTerminalDispatch export HVDCTwoTerminalPiecewiseLoss export HVDCTwoTerminalLCC +export HVDCTwoTerminalVSC +export HVDCTwoTerminalVSCBin2 # Converter Formulations export LosslessConverter diff --git a/src/ac_transmission_models/branch_constructor.jl b/src/ac_transmission_models/branch_constructor.jl index 568b072..3da7d40 100644 --- a/src/ac_transmission_models/branch_constructor.jl +++ b/src/ac_transmission_models/branch_constructor.jl @@ -1670,3 +1670,180 @@ function construct_device!( add_feedforward_constraints!(container, device_model, devices) return end + +############################################################################ +####################### Two-Terminal VSC HVDC Construct #################### +############################################################################ + +# Quadratic / bilinear approximation traits — same scheme used by the MT +# converter formulations. +_quad_config(::Type{HVDCTwoTerminalVSC}) = IOM.NoQuadApproxConfig() +_quad_config(::Type{HVDCTwoTerminalVSCBin2}) = + IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH) +_bilinear_config(::Type{HVDCTwoTerminalVSC}) = IOM.NoBilinearApproxConfig() +_bilinear_config(::Type{HVDCTwoTerminalVSCBin2}) = + IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)) + +# IOM's quadratic/bilinear approximation helpers take per-device bounds as +# `Vector{IOM.MinMax}` (one (min, max) pair per device, same order as `devices`). +function _vsc_v_from_bounds(devices) + n = length(devices) + bounds = Vector{IOM.MinMax}(undef, n) + for (k, d) in enumerate(devices) + lims = PSY.get_voltage_limits_from(d) + bounds[k] = IOM.MinMax((min = lims.min, max = lims.max)) + end + return bounds +end + +function _vsc_v_to_bounds(devices) + n = length(devices) + bounds = Vector{IOM.MinMax}(undef, n) + for (k, d) in enumerate(devices) + lims = PSY.get_voltage_limits_to(d) + bounds[k] = IOM.MinMax((min = lims.min, max = lims.max)) + end + return bounds +end + +function _vsc_i_bounds(devices) + n = length(devices) + bounds = Vector{IOM.MinMax}(undef, n) + for (k, d) in enumerate(devices) + i_max = _vsc_shared_i_max(d) + bounds[k] = IOM.MinMax((min = -i_max, max = i_max)) + end + return bounds +end + +function construct_device!( + container::OptimizationContainer, + sys::PSY.System, + ::ArgumentConstructStage, + device_model::DeviceModel{PSY.TwoTerminalVSCLine, F}, + network_model::NetworkModel{<:AbstractPowerModel}, +) where {F <: AbstractTwoTerminalVSCFormulation} + devices = get_available_components(device_model, sys) + + add_variables!(container, FlowActivePowerFromToVariable, devices, F) + add_variables!(container, FlowActivePowerToFromVariable, devices, F) + add_variables!(container, DCLineCurrentFlowVariable, devices, F) + add_variables!(container, HVDCFromDCVoltage, devices, F) + add_variables!(container, HVDCToDCVoltage, devices, F) + + _maybe_add_reactive_power_variables!(container, devices, device_model, network_model) + + if _use_linear_loss(F, device_model) + ll_devices = _devices_with_linear_loss(devices) + if isempty(ll_devices) + @warn "use_linear_loss is enabled but every TwoTerminalVSCLine has zero proportional loss terms; no linear-loss variables/constraints will be added." + else + _add_abs_value_decomposition_variables!(container, ll_devices, device_model) + end + end + + add_to_expression!( + container, ActivePowerBalance, FlowActivePowerFromToVariable, + devices, device_model, network_model, + ) + add_to_expression!( + container, ActivePowerBalance, FlowActivePowerToFromVariable, + devices, device_model, network_model, + ) + + add_feedforward_arguments!(container, device_model, devices) + return +end + +function construct_device!( + container::OptimizationContainer, + sys::PSY.System, + ::ModelConstructStage, + device_model::DeviceModel{PSY.TwoTerminalVSCLine, F}, + network_model::NetworkModel{<:AbstractPowerModel}, +) where {F <: AbstractTwoTerminalVSCFormulation} + devices = get_available_components(device_model, sys) + time_steps = get_time_steps(container) + line_names = [PSY.get_name(d) for d in devices] + + v_f_var = get_variable(container, HVDCFromDCVoltage, PSY.TwoTerminalVSCLine) + v_t_var = get_variable(container, HVDCToDCVoltage, PSY.TwoTerminalVSCLine) + i_var = get_variable(container, DCLineCurrentFlowVariable, PSY.TwoTerminalVSCLine) + + v_f_bounds = _vsc_v_from_bounds(devices) + v_t_bounds = _vsc_v_to_bounds(devices) + i_bounds = _vsc_i_bounds(devices) + + quad_cfg, bilin_cfg = _quad_config(F), _bilinear_config(F) + + v_f_sq_expr = IOM._add_quadratic_approx!( + quad_cfg, container, PSY.TwoTerminalVSCLine, + line_names, time_steps, v_f_var, v_f_bounds, "v_f_sq", + ) + v_t_sq_expr = IOM._add_quadratic_approx!( + quad_cfg, container, PSY.TwoTerminalVSCLine, + line_names, time_steps, v_t_var, v_t_bounds, "v_t_sq", + ) + i_sq_expr = IOM._add_quadratic_approx!( + quad_cfg, container, PSY.TwoTerminalVSCLine, + line_names, time_steps, i_var, i_bounds, "i_sq", + ) + + IOM._add_bilinear_approx!( + bilin_cfg, container, PSY.TwoTerminalVSCLine, + line_names, time_steps, + v_f_sq_expr, i_sq_expr, v_f_var, i_var, + v_f_bounds, i_bounds, "vi_ft", + ) + IOM._add_bilinear_approx!( + bilin_cfg, container, PSY.TwoTerminalVSCLine, + line_names, time_steps, + v_t_sq_expr, i_sq_expr, v_t_var, i_var, + v_t_bounds, i_bounds, "vi_tf", + ) + + add_constraints!( + container, HVDCCableOhmsLawConstraint, devices, device_model, network_model, + ) + add_constraints!( + container, HVDCVSCConverterPowerConstraint, devices, device_model, network_model, + ) + _maybe_add_reactive_power_constraints!(container, devices, device_model, network_model) + + if _use_linear_loss(F, device_model) + ll_devices = _devices_with_linear_loss(devices) + if !isempty(ll_devices) + _add_abs_value_decomposition_constraints!( + container, ll_devices, device_model, network_model, + DCLineCurrentFlowVariable, _vsc_shared_i_max, + ) + end + end + + add_constraint_dual!(container, sys, device_model) + add_feedforward_constraints!(container, device_model, devices) + return +end + +# AreaBalancePowerModel warning (consistent with other two-terminal formulations). +function construct_device!( + ::OptimizationContainer, + ::PSY.System, + ::ArgumentConstructStage, + ::DeviceModel{PSY.TwoTerminalVSCLine, <:AbstractTwoTerminalVSCFormulation}, + ::NetworkModel{AreaBalancePowerModel}, +) + @warn "AreaBalancePowerModel doesn't model individual line flows for PSY.TwoTerminalVSCLine. Arguments not built" + return +end + +function construct_device!( + ::OptimizationContainer, + ::PSY.System, + ::ModelConstructStage, + ::DeviceModel{PSY.TwoTerminalVSCLine, <:AbstractTwoTerminalVSCFormulation}, + ::NetworkModel{AreaBalancePowerModel}, +) + @warn "AreaBalancePowerModel doesn't model individual line flows for PSY.TwoTerminalVSCLine. Model not built" + return +end diff --git a/src/common_models/quadratic_converter_loss.jl b/src/common_models/quadratic_converter_loss.jl new file mode 100644 index 0000000..b2ad2eb --- /dev/null +++ b/src/common_models/quadratic_converter_loss.jl @@ -0,0 +1,124 @@ +# Shared helpers for quadratic / two-term converter losses +# loss(I) = a * I^2 + b * |I| + c +# Used by multi-terminal InterconnectingConverter formulations +# (Bin2QuadraticLossConverter, QuadraticLossConverter) and two-terminal +# HVDCTwoTerminalVSC formulations. + +######################################### +######## Loss-curve introspection ####### +######################################### + +_get_quadratic_term(loss_fn::PSY.QuadraticCurve) = PSY.get_quadratic_term(loss_fn) +_get_quadratic_term(loss_fn) = 0.0 + +# Whether the formulation wants the b*|I| linear-loss term (which requires +# decomposing the current into positive/negative parts with a direction binary). +_use_linear_loss(::Type{Bin2QuadraticLossConverter}, _) = true +_use_linear_loss(::Type{QuadraticLossConverter}, model) = + get_attribute(model, "use_linear_loss") +_use_linear_loss(::Type{HVDCTwoTerminalVSCBin2}, _) = true +_use_linear_loss(::Type{HVDCTwoTerminalVSC}, model) = + get_attribute(model, "use_linear_loss") + +# Per-device test: does this device have a nonzero linear loss term anywhere? +# Dispatched on device type because different PSY devices store loss curves on +# different fields. +_has_linear_loss(d::PSY.InterconnectingConverter) = + !iszero(PSY.get_proportional_term(PSY.get_loss_function(d))) +_has_linear_loss(d::PSY.TwoTerminalVSCLine) = + !iszero(PSY.get_proportional_term(PSY.get_converter_loss_from(d))) || + !iszero(PSY.get_proportional_term(PSY.get_converter_loss_to(d))) + +function _devices_with_linear_loss(devices) + return [d for d in devices if _has_linear_loss(d)] +end + +######################################### +######## Loss expression builder ######## +######################################### + +# Returns the JuMP expression a*i_sq + b*(i_pos + i_neg) + c +# for a single (device, time). The b*(i_pos+i_neg) term is included only when +# the formulation has opted into the linear-loss path AND b is nonzero for +# this specific device. +function _quadratic_converter_loss_expr( + a, b, c, i_sq_t, i_pos_t, i_neg_t; use_linear_loss::Bool, +) + loss = a * i_sq_t + c + if use_linear_loss && !iszero(b) + loss += b * (i_pos_t + i_neg_t) + end + return loss +end + +######################################### +####### Abs-value decomposition ######### +######################################### + +# Adds the three variables (PositiveCurrent, NegativeCurrent, CurrentDirection) +# that decompose a signed current variable into i = i^+ - i^- with a binary +# direction indicator. Called from the ArgumentConstructStage. +function _add_abs_value_decomposition_variables!( + container::OptimizationContainer, + devices, + ::DeviceModel{D, F}, +) where {D <: PSY.Device, F} + add_variables!(container, PositiveCurrent, devices, F) + add_variables!(container, NegativeCurrent, devices, F) + add_variables!(container, CurrentDirection, devices, F) + return +end + +# Adds the three constraints implementing i = i^+ - i^- with the big-M +# direction binary bounds. The CurrentAbsoluteValueConstraint container is +# created internally with three meta-tagged sub-containers ("", "pos_ub", +# "neg_ub"). Caller passes the parent current variable type and a function +# `d -> i_max` so device-specific bound lookups stay device-specific. +function _add_abs_value_decomposition_constraints!( + container::OptimizationContainer, + devices, + ::DeviceModel{D, F}, + ::NetworkModel{<:AbstractPowerModel}, + parent_var_type::Type{<:VariableType}, + i_max_getter::Function, +) where {D <: PSY.Device, F} + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + jump_model = get_jump_model(container) + i_var = get_variable(container, parent_var_type, D) + i_pos_var = get_variable(container, PositiveCurrent, D) + i_neg_var = get_variable(container, NegativeCurrent, D) + i_dir_var = get_variable(container, CurrentDirection, D) + + abs_val_const = add_constraints_container!( + container, CurrentAbsoluteValueConstraint, D, names, time_steps, + ) + pos_ub_const = add_constraints_container!( + container, CurrentAbsoluteValueConstraint, D, names, time_steps; + meta = "pos_ub", + ) + neg_ub_const = add_constraints_container!( + container, CurrentAbsoluteValueConstraint, D, names, time_steps; + meta = "neg_ub", + ) + + for d in devices + name = PSY.get_name(d) + i_max = i_max_getter(d) + for t in time_steps + abs_val_const[name, t] = JuMP.@constraint( + jump_model, + i_var[name, t] == i_pos_var[name, t] - i_neg_var[name, t], + ) + pos_ub_const[name, t] = JuMP.@constraint( + jump_model, + i_pos_var[name, t] <= i_max * i_dir_var[name, t], + ) + neg_ub_const[name, t] = JuMP.@constraint( + jump_model, + i_neg_var[name, t] <= i_max * (1 - i_dir_var[name, t]), + ) + end + end + return +end diff --git a/src/core/constraints.jl b/src/core/constraints.jl index 634e5c2..107a794 100644 --- a/src/core/constraints.jl +++ b/src/core/constraints.jl @@ -419,6 +419,39 @@ v_d^i = v_d^r - R_d I_d """ struct HVDCTransmissionDCLineConstraint <: ConstraintType end +""" +Per-terminal converter power-balance constraint for two-terminal VSC HVDC: + +```math +\\begin{aligned} +p_{ft} &= v_f \\cdot I - (a_f I^2 + b_f |I| + c_f) \\\\ +p_{tf} &= -v_t \\cdot I - (a_t I^2 + b_t |I| + c_t) +\\end{aligned} +``` +""" +struct HVDCVSCConverterPowerConstraint <: ConstraintType end + +""" +Cable Ohm's law for a two-terminal HVDC link with explicit DC resistance: + +```math +v_f - v_t = (1/g) \\cdot I +``` +""" +struct HVDCCableOhmsLawConstraint <: ConstraintType end + +""" +PQ capability constraint at each terminal of a two-terminal VSC HVDC, added only +on AC networks. For `HVDCTwoTerminalVSC` this is the exact circle: + +```math +p_k^2 + q_k^2 \\le (S_k^{\\max})^2 \\quad k \\in \\{f, t\\} +``` + +For `HVDCTwoTerminalVSCBin2` it is replaced by an inscribed polygon to stay MILP. +""" +struct HVDCVSCReactiveCapabilityConstraint <: ConstraintType end + abstract type PowerVariableLimitsConstraint <: ConstraintType end """ Struct to create the constraint to limit active power input expressions. diff --git a/src/core/formulations.jl b/src/core/formulations.jl index c37aec3..533f84e 100644 --- a/src/core/formulations.jl +++ b/src/core/formulations.jl @@ -184,8 +184,27 @@ Branch type to represent non-linear LCC (line commutated converter) model on two """ struct HVDCTwoTerminalLCC <: AbstractTwoTerminalDCLineFormulation end -# Not Implemented -# struct VoltageSourceDC <: AbstractTwoTerminalDCLineFormulation end +""" +Abstract supertype for two-terminal voltage-source converter (VSC) HVDC formulations. +Models per-terminal converters with quadratic / two-term losses +(``a I^2 + b |I| + c``), a shared signed cable current, an explicit DC-side +cable resistance (``v_f - v_t = (1/g) \\cdot I``), and (on AC networks) independent +reactive-power control bounded by per-terminal PQ capability. +""" +abstract type AbstractTwoTerminalVSCFormulation <: AbstractTwoTerminalDCLineFormulation end + +""" +Two-terminal VSC formulation that keeps the bilinear ``v \\cdot I`` and quadratic +``I^2`` terms exact. Requires an NLP-capable solver (e.g. Ipopt). +""" +struct HVDCTwoTerminalVSC <: AbstractTwoTerminalVSCFormulation end + +""" +Two-terminal VSC formulation that approximates the bilinear ``v \\cdot I`` and +quadratic ``I^2`` terms with SOS2-based piecewise-linear envelopes (the same +Bin2 scheme used by `Bin2QuadraticLossConverter`). Stays MILP. +""" +struct HVDCTwoTerminalVSCBin2 <: AbstractTwoTerminalVSCFormulation end ############################### AC/DC Converter Formulations ##################################### abstract type AbstractConverterFormulation <: AbstractDeviceFormulation end diff --git a/src/core/variables.jl b/src/core/variables.jl index 8a41d0d..44ab5bb 100644 --- a/src/core/variables.jl +++ b/src/core/variables.jl @@ -157,22 +157,25 @@ Docs abbreviation: ``i_c^{dc}`` struct ConverterCurrent <: VariableType end """ -Struct to dispatch the creation of DC Converter Positive Term Current Variables for DC formulations -Docs abbreviation: ``i_c^{+,dc}`` +Positive part of an absolute-value current decomposition (``i = i^+ - i^-``). +Used by both two-terminal HVDC links (cable current) and interconnecting converters. +Docs abbreviation: ``i^{+,dc}`` """ -struct ConverterPositiveCurrent <: VariableType end +struct PositiveCurrent <: VariableType end """ -Struct to dispatch the creation of DC Converter Negative Term Current Variables for DC formulations -Docs abbreviation: ``i_c^{-,dc}`` +Negative part of an absolute-value current decomposition (``i = i^+ - i^-``). +Used by both two-terminal HVDC links (cable current) and interconnecting converters. +Docs abbreviation: ``i^{-,dc}`` """ -struct ConverterNegativeCurrent <: VariableType end +struct NegativeCurrent <: VariableType end """ -Struct to dispatch the creation of DC Converter Binary for Absolute Value Current Variables for DC formulations -Docs abbreviation: `\\nu_c`` +Binary indicator selecting which sign the current takes in the absolute-value +decomposition (``i^+`` active when 1, ``i^-`` active when 0). +Docs abbreviation: ``\\nu`` """ -struct ConverterCurrentDirection <: VariableType end +struct CurrentDirection <: VariableType end ######################################################### ######################################################### @@ -384,6 +387,34 @@ Docs abbreviation: ``z`` """ struct HVDCPiecewiseBinaryLossVariable <: SparseVariableType end +""" +DC-side voltage at the from-terminal of a two-terminal HVDC link. +Used by `HVDCTwoTerminalVSC` formulations. +Docs abbreviation: ``v_f^{dc}`` +""" +struct HVDCFromDCVoltage <: VariableType end + +""" +DC-side voltage at the to-terminal of a two-terminal HVDC link. +Used by `HVDCTwoTerminalVSC` formulations. +Docs abbreviation: ``v_t^{dc}`` +""" +struct HVDCToDCVoltage <: VariableType end + +""" +Reactive power injected at the from-terminal AC bus by a two-terminal HVDC link. +Added only when the formulation runs on an AC network model. +Docs abbreviation: ``q_f`` +""" +struct HVDCReactivePowerFromVariable <: VariableType end + +""" +Reactive power injected at the to-terminal AC bus by a two-terminal HVDC link. +Added only when the formulation runs on an AC network model. +Docs abbreviation: ``q_t`` +""" +struct HVDCReactivePowerToVariable <: VariableType end + """ Struct to dispatch the creation of Interface Flow Slack Up variables diff --git a/src/mt_hvdc_models/HVDCsystems.jl b/src/mt_hvdc_models/HVDCsystems.jl index b3fcde8..3d197d8 100644 --- a/src/mt_hvdc_models/HVDCsystems.jl +++ b/src/mt_hvdc_models/HVDCsystems.jl @@ -120,22 +120,22 @@ end ## Binaries ### get_variable_binary(::Type{ConverterCurrent}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = false -get_variable_binary(::Type{ConverterPositiveCurrent}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = false -get_variable_binary(::Type{ConverterNegativeCurrent}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = false -get_variable_binary(::Type{ConverterCurrentDirection}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = true +get_variable_binary(::Type{PositiveCurrent}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = false +get_variable_binary(::Type{NegativeCurrent}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = false +get_variable_binary(::Type{CurrentDirection}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = true ### Warm Start ### get_variable_warm_start_value(::Type{ConverterCurrent}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = PSY.get_dc_current(d) ### Lower Bounds ### get_variable_lower_bound(::Type{ConverterCurrent}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = -PSY.get_max_dc_current(d) -get_variable_lower_bound(::Type{ConverterPositiveCurrent}, d::PSY.InterconnectingConverter,::Type{<:AbstractConverterFormulation}) = 0.0 -get_variable_lower_bound(::Type{ConverterNegativeCurrent}, d::PSY.InterconnectingConverter,::Type{<:AbstractConverterFormulation}) = 0.0 +get_variable_lower_bound(::Type{PositiveCurrent}, d::PSY.InterconnectingConverter,::Type{<:AbstractConverterFormulation}) = 0.0 +get_variable_lower_bound(::Type{NegativeCurrent}, d::PSY.InterconnectingConverter,::Type{<:AbstractConverterFormulation}) = 0.0 ### Upper Bounds ### get_variable_upper_bound(::Type{ConverterCurrent}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = PSY.get_max_dc_current(d) -get_variable_upper_bound(::Type{ConverterPositiveCurrent}, d::PSY.InterconnectingConverter,::Type{<:AbstractConverterFormulation}) = PSY.get_max_dc_current(d) -get_variable_upper_bound(::Type{ConverterNegativeCurrent}, d::PSY.InterconnectingConverter,::Type{<:AbstractConverterFormulation}) = PSY.get_max_dc_current(d) +get_variable_upper_bound(::Type{PositiveCurrent}, d::PSY.InterconnectingConverter,::Type{<:AbstractConverterFormulation}) = PSY.get_max_dc_current(d) +get_variable_upper_bound(::Type{NegativeCurrent}, d::PSY.InterconnectingConverter,::Type{<:AbstractConverterFormulation}) = PSY.get_max_dc_current(d) function get_default_attributes( @@ -540,20 +540,6 @@ end ############## Converters ################## -_get_quadratic_term(loss_fn::PSY.QuadraticCurve) = PSY.get_quadratic_term(loss_fn) -_get_quadratic_term(loss_fn) = 0.0 - -_use_linear_loss(::Type{Bin2QuadraticLossConverter}, _) = true -_use_linear_loss(::Type{QuadraticLossConverter}, model) = - get_attribute(model, "use_linear_loss") - -function _devices_with_linear_loss(devices) - return [ - d for d in devices if - !iszero(PSY.get_proportional_term(PSY.get_loss_function(d))) - ] -end - function add_constraints!( container::OptimizationContainer, ::Type{ConverterLossConstraint}, @@ -572,8 +558,8 @@ function add_constraints!( use_linear_loss = _use_linear_loss(V, model) && !isempty(_devices_with_linear_loss(devices)) if use_linear_loss - i_pos_var = get_variable(container, ConverterPositiveCurrent, U) - i_neg_var = get_variable(container, ConverterNegativeCurrent, U) + i_pos_var = get_variable(container, PositiveCurrent, U) + i_neg_var = get_variable(container, NegativeCurrent, U) end ipc_names = [PSY.get_name(d) for d in devices] @@ -589,13 +575,15 @@ function add_constraints!( b = PSY.get_proportional_term(loss_function) c = PSY.get_constant_term(loss_function) for t in time_steps - loss = a * i_sq_expr[name, t] + c - if use_linear_loss && !iszero(b) - loss += b * (i_pos_var[name, t] + i_neg_var[name, t]) - end + i_pos_t = use_linear_loss ? i_pos_var[name, t] : nothing + i_neg_t = use_linear_loss ? i_neg_var[name, t] : nothing + loss = _quadratic_converter_loss_expr( + a, b, c, i_sq_expr[name, t], i_pos_t, i_neg_t; + use_linear_loss = use_linear_loss, + ) loss_const[name, t] = JuMP.@constraint( jump_model, - P_ac_var[name, t] == vi_expr[name, t] - loss + P_ac_var[name, t] == vi_expr[name, t] - loss, ) end end @@ -606,52 +594,18 @@ function add_constraints!( container::OptimizationContainer, ::Type{T}, devices::W, - ::DeviceModel{U, V}, - ::NetworkModel{<:AbstractPowerModel}, + model::DeviceModel{U, V}, + network_model::NetworkModel{<:AbstractPowerModel}, ) where { T <: CurrentAbsoluteValueConstraint, U <: PSY.InterconnectingConverter, V <: AbstractQuadraticLossConverter, W <: Union{Vector{U}, IS.FlattenIteratorWrapper{U}}, } - time_steps = get_time_steps(container) - names = [PSY.get_name(d) for d in devices] - jump_model = get_jump_model(container) - i_var = get_variable(container, ConverterCurrent, U) - i_pos_var = get_variable(container, ConverterPositiveCurrent, U) - i_neg_var = get_variable(container, ConverterNegativeCurrent, U) - i_dir_var = get_variable(container, ConverterCurrentDirection, U) - - abs_val_const = add_constraints_container!( - container, CurrentAbsoluteValueConstraint, U, names, time_steps, - ) - pos_ub_const = add_constraints_container!( - container, CurrentAbsoluteValueConstraint, U, names, time_steps; - meta = "pos_ub", - ) - neg_ub_const = add_constraints_container!( - container, CurrentAbsoluteValueConstraint, U, names, time_steps; - meta = "neg_ub", + _add_abs_value_decomposition_constraints!( + container, devices, model, network_model, + ConverterCurrent, PSY.get_max_dc_current, ) - - for d in devices - name = PSY.get_name(d) - i_max = PSY.get_max_dc_current(d) - for t in time_steps - abs_val_const[name, t] = JuMP.@constraint( - jump_model, - i_var[name, t] == i_pos_var[name, t] - i_neg_var[name, t] - ) - pos_ub_const[name, t] = JuMP.@constraint( - jump_model, - i_pos_var[name, t] <= i_max * i_dir_var[name, t] - ) - neg_ub_const[name, t] = JuMP.@constraint( - jump_model, - i_neg_var[name, t] <= i_max * (1 - i_dir_var[name, t]) - ) - end - end return end diff --git a/src/mt_hvdc_models/hvdcsystems_constructor.jl b/src/mt_hvdc_models/hvdcsystems_constructor.jl index c1ec321..249d856 100644 --- a/src/mt_hvdc_models/hvdcsystems_constructor.jl +++ b/src/mt_hvdc_models/hvdcsystems_constructor.jl @@ -59,9 +59,9 @@ function construct_device!( if isempty(ll_devices) @warn "use_linear_loss is enabled but every InterconnectingConverter has a zero proportional loss term; no linear-loss variables/constraints will be added." else - add_variables!(container, ConverterPositiveCurrent, ll_devices, T) - add_variables!(container, ConverterNegativeCurrent, ll_devices, T) - add_variables!(container, ConverterCurrentDirection, ll_devices, T) + add_variables!(container, PositiveCurrent, ll_devices, T) + add_variables!(container, NegativeCurrent, ll_devices, T) + add_variables!(container, CurrentDirection, ll_devices, T) end end @@ -77,19 +77,19 @@ function construct_device!( return end - function _voltage_expr_per_converter( - container::OptimizationContainer, - devices, - ipc_names::Vector{String}, - time_steps, - ) - v_var = get_variable(container, DCVoltage, PSY.DCBus) - bus_names = [PSY.get_name(PSY.get_dc_bus(d)) for d in devices] - return JuMP.Containers.DenseAxisArray( - [v_var[b, t] for b in bus_names, t in time_steps], - ipc_names, time_steps, - ) - end +function _voltage_expr_per_converter( + container::OptimizationContainer, + devices, + ipc_names::Vector{String}, + time_steps, +) + v_var = get_variable(container, DCVoltage, PSY.DCBus) + bus_names = [PSY.get_name(PSY.get_dc_bus(d)) for d in devices] + return JuMP.Containers.DenseAxisArray( + [v_var[b, t] for b in bus_names, t in time_steps], + ipc_names, time_steps, + ) +end function _converter_vi_bounds(devices) n = length(devices) @@ -104,9 +104,11 @@ function _converter_vi_bounds(devices) return v_bounds, i_bounds end -_quad_config(::Type{Bin2QuadraticLossConverter}) = IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH) +_quad_config(::Type{Bin2QuadraticLossConverter}) = + IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH) _quad_config(::Type{QuadraticLossConverter}) = IOM.NoQuadApproxConfig() -_bilinear_config(::Type{Bin2QuadraticLossConverter}) = IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)) +_bilinear_config(::Type{Bin2QuadraticLossConverter}) = + IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)) _bilinear_config(::Type{QuadraticLossConverter}) = IOM.NoBilinearApproxConfig() function construct_device!( diff --git a/src/network_models/pm_translator.jl b/src/network_models/pm_translator.jl index 09f999e..46e8827 100644 --- a/src/network_models/pm_translator.jl +++ b/src/network_models/pm_translator.jl @@ -663,6 +663,27 @@ function get_branch_to_pm( return Dict{String, Any}() end +function get_branch_to_pm( + ix::Int, + branch::PSY.TwoTerminalVSCLine, + ::Type{<:AbstractTwoTerminalVSCFormulation}, + ::Type{<:AbstractPowerModel}, +) + return Dict{String, Any}() +end + +# Trait: should this DeviceModel be skipped when translating two-terminal HVDC +# branches into PowerModels? LCC and VSC formulations are handled by their own +# constraints in POM rather than via PM's built-in DC line model, so they're +# skipped here. Default is `false`; override per formulation via dispatch. +_skip_pm_two_terminal_translation(::DeviceModel) = false +_skip_pm_two_terminal_translation( + ::DeviceModel{PSY.TwoTerminalLCCLine, HVDCTwoTerminalLCC}, +) = true +_skip_pm_two_terminal_translation( + ::DeviceModel{PSY.TwoTerminalVSCLine, <:AbstractTwoTerminalVSCFormulation}, +) = true + function get_branches_to_pm( sys::PSY.System, network_model::NetworkModel{S}, @@ -720,10 +741,7 @@ function get_branches_to_pm( for (d, device_model) in branch_template comp_type = get_component_type(device_model) !(comp_type <: T) && continue - if comp_type <: PSY.TwoTerminalLCCLine && - get_formulation(device_model) <: HVDCTwoTerminalLCC - continue - end + _skip_pm_two_terminal_translation(device_model) && continue start_idx += length(PM_branches) for (i, branch) in enumerate(get_available_components(device_model, sys)) ix = i + start_idx diff --git a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl index 5422b6f..4fa388e 100644 --- a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl +++ b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl @@ -1433,3 +1433,362 @@ function add_constraints!( end return end + +############################################################################## +####################### Two-Terminal VSC Formulation ######################### +############################################################################## + +#! format: off + +# Variable trait methods for the shared cable current and DC voltages +get_variable_binary(::Type{DCLineCurrentFlowVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false +get_variable_binary(::Type{HVDCFromDCVoltage}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false +get_variable_binary(::Type{HVDCToDCVoltage}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false +get_variable_binary(::Type{HVDCReactivePowerFromVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false +get_variable_binary(::Type{HVDCReactivePowerToVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false +get_variable_binary(::Type{PositiveCurrent}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false +get_variable_binary(::Type{NegativeCurrent}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false +get_variable_binary(::Type{CurrentDirection}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = true +get_variable_binary(::Type{FlowActivePowerFromToVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false +get_variable_binary(::Type{FlowActivePowerToFromVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false + +# Warm starts +get_variable_warm_start_value(::Type{DCLineCurrentFlowVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_dc_current(d) +get_variable_warm_start_value(::Type{HVDCReactivePowerFromVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_reactive_power_from(d) +get_variable_warm_start_value(::Type{HVDCReactivePowerToVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_reactive_power_to(d) +get_variable_warm_start_value(::Type{FlowActivePowerFromToVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_active_power_flow(d) +get_variable_warm_start_value(::Type{FlowActivePowerToFromVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = -PSY.get_active_power_flow(d) + +# Active power flow bounds (per-terminal) +get_variable_lower_bound(::Type{FlowActivePowerFromToVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_active_power_limits_from(d).min +get_variable_upper_bound(::Type{FlowActivePowerFromToVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_active_power_limits_from(d).max +get_variable_lower_bound(::Type{FlowActivePowerToFromVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_active_power_limits_to(d).min +get_variable_upper_bound(::Type{FlowActivePowerToFromVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_active_power_limits_to(d).max + +# Reactive power bounds (per-terminal) +get_variable_lower_bound(::Type{HVDCReactivePowerFromVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_reactive_power_limits_from(d).min +get_variable_upper_bound(::Type{HVDCReactivePowerFromVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_reactive_power_limits_from(d).max +get_variable_lower_bound(::Type{HVDCReactivePowerToVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_reactive_power_limits_to(d).min +get_variable_upper_bound(::Type{HVDCReactivePowerToVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_reactive_power_limits_to(d).max + +# DC voltage bounds (per-terminal) +get_variable_lower_bound(::Type{HVDCFromDCVoltage}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_voltage_limits_from(d).min +get_variable_upper_bound(::Type{HVDCFromDCVoltage}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_voltage_limits_from(d).max +get_variable_lower_bound(::Type{HVDCToDCVoltage}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_voltage_limits_to(d).min +get_variable_upper_bound(::Type{HVDCToDCVoltage}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_voltage_limits_to(d).max + +# Shared cable current bounds — must respect BOTH terminals' I_max ratings +_vsc_shared_i_max(d::PSY.TwoTerminalVSCLine) = + min(PSY.get_max_dc_current_from(d), PSY.get_max_dc_current_to(d)) +get_variable_lower_bound(::Type{DCLineCurrentFlowVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = -_vsc_shared_i_max(d) +get_variable_upper_bound(::Type{DCLineCurrentFlowVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = _vsc_shared_i_max(d) + +# Positive/negative parts: each in [0, i_max] +get_variable_lower_bound(::Type{PositiveCurrent}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = 0.0 +get_variable_upper_bound(::Type{PositiveCurrent}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = _vsc_shared_i_max(d) +get_variable_lower_bound(::Type{NegativeCurrent}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = 0.0 +get_variable_upper_bound(::Type{NegativeCurrent}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = _vsc_shared_i_max(d) + +#! format: on + +####################### VSC reactive-power optional path ##################### +# The reactive power and PQ capability machinery is added only when the +# network actually models reactive power. On active-only networks these +# helpers are no-ops. + +function _maybe_add_reactive_power_variables!( + container::OptimizationContainer, + devices, + model::DeviceModel{PSY.TwoTerminalVSCLine, F}, + ::NetworkModel{<:AbstractPowerModel}, +) where {F <: AbstractTwoTerminalVSCFormulation} + add_variables!(container, HVDCReactivePowerFromVariable, devices, F) + add_variables!(container, HVDCReactivePowerToVariable, devices, F) + return +end + +_maybe_add_reactive_power_variables!( + ::OptimizationContainer, + _devices, + ::DeviceModel{PSY.TwoTerminalVSCLine, <:AbstractTwoTerminalVSCFormulation}, + ::NetworkModel{<:AbstractActivePowerModel}, +) = nothing + +function _maybe_add_reactive_power_constraints!( + container::OptimizationContainer, + devices, + model::DeviceModel{PSY.TwoTerminalVSCLine, F}, + network_model::NetworkModel{<:AbstractPowerModel}, +) where {F <: AbstractTwoTerminalVSCFormulation} + add_constraints!( + container, HVDCVSCReactiveCapabilityConstraint, + devices, model, network_model, + ) + return +end + +_maybe_add_reactive_power_constraints!( + ::OptimizationContainer, + _devices, + ::DeviceModel{PSY.TwoTerminalVSCLine, <:AbstractTwoTerminalVSCFormulation}, + ::NetworkModel{<:AbstractActivePowerModel}, +) = nothing + +####################### VSC core constraints ################################ + +# Cable Ohm's law: v_f - v_t = (1/g) * I +function add_constraints!( + container::OptimizationContainer, + ::Type{HVDCCableOhmsLawConstraint}, + devices::Union{ + Vector{PSY.TwoTerminalVSCLine}, + IS.FlattenIteratorWrapper{PSY.TwoTerminalVSCLine}, + }, + ::DeviceModel{PSY.TwoTerminalVSCLine, F}, + ::NetworkModel{<:AbstractPowerModel}, +) where {F <: AbstractTwoTerminalVSCFormulation} + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + jump_model = get_jump_model(container) + v_f = get_variable(container, HVDCFromDCVoltage, PSY.TwoTerminalVSCLine) + v_t = get_variable(container, HVDCToDCVoltage, PSY.TwoTerminalVSCLine) + i_var = get_variable(container, DCLineCurrentFlowVariable, PSY.TwoTerminalVSCLine) + + cons = add_constraints_container!( + container, HVDCCableOhmsLawConstraint, PSY.TwoTerminalVSCLine, + names, time_steps, + ) + + for d in devices + name = PSY.get_name(d) + g = PSY.get_g(d) + # If g == 0 we treat the cable as lossless (v_f == v_t). + # Otherwise R = 1/g. + for t in time_steps + cons[name, t] = if iszero(g) + JuMP.@constraint(jump_model, v_f[name, t] == v_t[name, t]) + else + JuMP.@constraint( + jump_model, + v_f[name, t] - v_t[name, t] == (1.0 / g) * i_var[name, t], + ) + end + end + end + return +end + +# Per-terminal converter power balance: +# p_ft == v_f * I + (a_f * I^2 + b_f * |I| + c_f) +# p_tf == -v_t * I + (a_t * I^2 + b_t * |I| + c_t) +# Sign convention: FlowActivePowerFromToVariable / ToFromVariable are positive +# when the corresponding AC bus is sourcing power into the converter (matches +# the existing add_to_expression! method's -1.0 multiplier). +function add_constraints!( + container::OptimizationContainer, + ::Type{HVDCVSCConverterPowerConstraint}, + devices::Union{ + Vector{PSY.TwoTerminalVSCLine}, + IS.FlattenIteratorWrapper{PSY.TwoTerminalVSCLine}, + }, + model::DeviceModel{PSY.TwoTerminalVSCLine, F}, + ::NetworkModel{<:AbstractPowerModel}, +) where {F <: AbstractTwoTerminalVSCFormulation} + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + jump_model = get_jump_model(container) + + p_ft = get_variable(container, FlowActivePowerFromToVariable, PSY.TwoTerminalVSCLine) + p_tf = get_variable(container, FlowActivePowerToFromVariable, PSY.TwoTerminalVSCLine) + vi_expr = get_expression( + container, + IOM.BilinearProductExpression, + PSY.TwoTerminalVSCLine, + "vi_ft", + ) + vi_expr_to = get_expression( + container, + IOM.BilinearProductExpression, + PSY.TwoTerminalVSCLine, + "vi_tf", + ) + i_sq_expr = + get_expression(container, IOM.QuadraticExpression, PSY.TwoTerminalVSCLine, "i_sq") + + use_linear_loss = + _use_linear_loss(F, model) && !isempty(_devices_with_linear_loss(devices)) + if use_linear_loss + i_pos_var = get_variable(container, PositiveCurrent, PSY.TwoTerminalVSCLine) + i_neg_var = get_variable(container, NegativeCurrent, PSY.TwoTerminalVSCLine) + end + + cons_ft = add_constraints_container!( + container, HVDCVSCConverterPowerConstraint, PSY.TwoTerminalVSCLine, + names, time_steps; meta = "ft", + ) + cons_tf = add_constraints_container!( + container, HVDCVSCConverterPowerConstraint, PSY.TwoTerminalVSCLine, + names, time_steps; meta = "tf", + ) + + for d in devices + name = PSY.get_name(d) + loss_from = PSY.get_converter_loss_from(d) + loss_to = PSY.get_converter_loss_to(d) + a_f = _get_quadratic_term(loss_from) + b_f = PSY.get_proportional_term(loss_from) + c_f = PSY.get_constant_term(loss_from) + a_t = _get_quadratic_term(loss_to) + b_t = PSY.get_proportional_term(loss_to) + c_t = PSY.get_constant_term(loss_to) + for t in time_steps + i_pos_t = use_linear_loss ? i_pos_var[name, t] : nothing + i_neg_t = use_linear_loss ? i_neg_var[name, t] : nothing + loss_ft = _quadratic_converter_loss_expr( + a_f, b_f, c_f, i_sq_expr[name, t], i_pos_t, i_neg_t; + use_linear_loss = use_linear_loss, + ) + loss_tf = _quadratic_converter_loss_expr( + a_t, b_t, c_t, i_sq_expr[name, t], i_pos_t, i_neg_t; + use_linear_loss = use_linear_loss, + ) + cons_ft[name, t] = JuMP.@constraint( + jump_model, + p_ft[name, t] == vi_expr[name, t] + loss_ft, + ) + cons_tf[name, t] = JuMP.@constraint( + jump_model, + p_tf[name, t] == -vi_expr_to[name, t] + loss_tf, + ) + end + end + return +end + +# PQ capability: p_k^2 + q_k^2 <= S_k^2 (NLP) or octagonal polygon (Bin2). +function add_constraints!( + container::OptimizationContainer, + ::Type{HVDCVSCReactiveCapabilityConstraint}, + devices::Union{ + Vector{PSY.TwoTerminalVSCLine}, + IS.FlattenIteratorWrapper{PSY.TwoTerminalVSCLine}, + }, + ::DeviceModel{PSY.TwoTerminalVSCLine, HVDCTwoTerminalVSC}, + ::NetworkModel{<:AbstractPowerModel}, +) + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + jump_model = get_jump_model(container) + + p_ft = get_variable(container, FlowActivePowerFromToVariable, PSY.TwoTerminalVSCLine) + p_tf = get_variable(container, FlowActivePowerToFromVariable, PSY.TwoTerminalVSCLine) + q_f = get_variable(container, HVDCReactivePowerFromVariable, PSY.TwoTerminalVSCLine) + q_t = get_variable(container, HVDCReactivePowerToVariable, PSY.TwoTerminalVSCLine) + + cons_f = add_constraints_container!( + container, HVDCVSCReactiveCapabilityConstraint, PSY.TwoTerminalVSCLine, + names, time_steps; meta = "from", + ) + cons_t = add_constraints_container!( + container, HVDCVSCReactiveCapabilityConstraint, PSY.TwoTerminalVSCLine, + names, time_steps; meta = "to", + ) + + for d in devices + name = PSY.get_name(d) + s_f = PSY.get_rating_from(d) + s_t = PSY.get_rating_to(d) + for t in time_steps + cons_f[name, t] = JuMP.@constraint( + jump_model, + p_ft[name, t]^2 + q_f[name, t]^2 <= s_f^2, + ) + cons_t[name, t] = JuMP.@constraint( + jump_model, + p_tf[name, t]^2 + q_t[name, t]^2 <= s_t^2, + ) + end + end + return +end + +# Bin2 variant: inscribed octagon in the rating circle (s/sqrt(2) half-diagonal). +# Eight linear constraints per terminal: ±p ± q ≤ s and ±p ≤ s, ±q ≤ s. +function add_constraints!( + container::OptimizationContainer, + ::Type{HVDCVSCReactiveCapabilityConstraint}, + devices::Union{ + Vector{PSY.TwoTerminalVSCLine}, + IS.FlattenIteratorWrapper{PSY.TwoTerminalVSCLine}, + }, + ::DeviceModel{PSY.TwoTerminalVSCLine, HVDCTwoTerminalVSCBin2}, + ::NetworkModel{<:AbstractPowerModel}, +) + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + jump_model = get_jump_model(container) + + p_ft = get_variable(container, FlowActivePowerFromToVariable, PSY.TwoTerminalVSCLine) + p_tf = get_variable(container, FlowActivePowerToFromVariable, PSY.TwoTerminalVSCLine) + q_f = get_variable(container, HVDCReactivePowerFromVariable, PSY.TwoTerminalVSCLine) + q_t = get_variable(container, HVDCReactivePowerToVariable, PSY.TwoTerminalVSCLine) + + cons = Dict{String, Any}() + for tag in ("from_pp", "from_pn", "from_np", "from_nn", + "to_pp", "to_pn", "to_np", "to_nn") + cons[tag] = add_constraints_container!( + container, HVDCVSCReactiveCapabilityConstraint, PSY.TwoTerminalVSCLine, + names, time_steps; meta = tag, + ) + end + + inv_sqrt2 = 1.0 / sqrt(2.0) + for d in devices + name = PSY.get_name(d) + s_f = PSY.get_rating_from(d) * inv_sqrt2 + s_t = PSY.get_rating_to(d) * inv_sqrt2 + for t in time_steps + # Octagon at the from terminal: project (p, q) onto 4 diagonal lines. + cons["from_pp"][name, t] = + JuMP.@constraint(jump_model, p_ft[name, t] + q_f[name, t] <= 2.0 * s_f) + cons["from_pn"][name, t] = + JuMP.@constraint(jump_model, p_ft[name, t] - q_f[name, t] <= 2.0 * s_f) + cons["from_np"][name, t] = + JuMP.@constraint(jump_model, -p_ft[name, t] + q_f[name, t] <= 2.0 * s_f) + cons["from_nn"][name, t] = + JuMP.@constraint(jump_model, -p_ft[name, t] - q_f[name, t] <= 2.0 * s_f) + cons["to_pp"][name, t] = + JuMP.@constraint(jump_model, p_tf[name, t] + q_t[name, t] <= 2.0 * s_t) + cons["to_pn"][name, t] = + JuMP.@constraint(jump_model, p_tf[name, t] - q_t[name, t] <= 2.0 * s_t) + cons["to_np"][name, t] = + JuMP.@constraint(jump_model, -p_tf[name, t] + q_t[name, t] <= 2.0 * s_t) + cons["to_nn"][name, t] = + JuMP.@constraint(jump_model, -p_tf[name, t] - q_t[name, t] <= 2.0 * s_t) + end + end + return +end + +####################### VSC defaults ######################################### + +function get_default_time_series_names( + ::Type{PSY.TwoTerminalVSCLine}, + ::Type{<:AbstractTwoTerminalVSCFormulation}, +) + return Dict{Type{<:TimeSeriesParameter}, String}() +end + +function get_default_attributes( + ::Type{PSY.TwoTerminalVSCLine}, + ::Type{HVDCTwoTerminalVSC}, +) + return Dict{String, Any}("use_linear_loss" => false) +end + +function get_default_attributes( + ::Type{PSY.TwoTerminalVSCLine}, + ::Type{HVDCTwoTerminalVSCBin2}, +) + return Dict{String, Any}() +end diff --git a/test/Project.toml b/test/Project.toml index a7425a5..43494b2 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -35,7 +35,7 @@ UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" InfrastructureSystems = {rev = "IS4", url = "https://github.com/Sienna-Platform/InfrastructureSystems.jl"} PowerSystemCaseBuilder = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerSystemCaseBuilder.jl"} PowerSystems = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerSystems.jl"} - +InfrastructureOptimizationModels = {rev = "ac/hvdc-vsc", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"} [compat] HiGHS = "1" diff --git a/test/runtests.jl b/test/runtests.jl index 7d8fc42..d290015 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -22,8 +22,14 @@ const DISABLED_TEST_FILES = [ # Can generate with ls -1 test | grep "test_.*.jl # "test_device_synchronous_condenser_constructors.jl", # "test_device_thermal_generation_constructors.jl", # "test_formulation_combinations.jl", +# "test_import_export_cost.jl", # "test_initialization_problem.jl", +# "test_is_time_variant_proportional.jl", +# "test_market_bid_cost.jl", +# "test_mbc_parameter_population.jl", # "test_model_decision.jl", +# "test_multi_interval.jl", +# "test_network_constructors.jl", # "test_network_constructors_with_dlr.jl", # "test_problem_template.jl", # "test_storage_device_models.jl", diff --git a/test/test_device_hvdc.jl b/test/test_device_hvdc.jl index 62652e8..08f6f63 100644 --- a/test/test_device_hvdc.jl +++ b/test/test_device_hvdc.jl @@ -208,9 +208,147 @@ end store_variable_names = true, optimizer = HiGHS_optimizer, ) - @test_logs( - (:warn, r"every InterconnectingConverter has a zero proportional loss"), - match_mode = :any, - build!(model; output_dir = mktempdir(; cleanup = true)), + # `build!` wraps its body in `Logging.with_logger(file_logger)`, which masks + # any `TestLogger` set up by `@test_logs`. The warning we want lands in the + # per-build log file instead, so read it back from there. + output_dir = mktempdir(; cleanup = true) + @test build!(model; output_dir = output_dir) == IOM.ModelBuildStatus.BUILT + log_path = joinpath(output_dir, IOM.PROBLEM_LOG_FILENAME) + @test occursin(r"linear[_ ]loss"i, read(log_path, String)) +end + +############################################################################## +################ Two-Terminal VSC HVDC tests ################################# +############################################################################## + +# Build a small AC test system and replace an AC line with a TwoTerminalVSCLine +# so we have a concrete VSC device to exercise the formulation against. +function _generate_test_vsc_sys(; + g = 50.0, + rating_from = 2.0, + rating_to = 2.0, + loss_a = 0.01, + loss_b = 0.0, + loss_c = 0.0, +) + sys = build_system(PSITestSystems, "c_sys5_uc"; force_build = true) + line = get_component(Line, sys, "1") + remove_component!(sys, line) + + vsc = TwoTerminalVSCLine(; + name = get_name(line), + available = true, + arc = get_arc(line), + active_power_flow = 0.0, + rating = max(rating_from, rating_to), + active_power_limits_from = (min = -rating_from, max = rating_from), + active_power_limits_to = (min = -rating_to, max = rating_to), + g = g, + dc_current = 0.0, + reactive_power_from = 0.0, + dc_voltage_control_from = true, + ac_voltage_control_from = true, + dc_setpoint_from = 0.0, + ac_setpoint_from = 1.0, + converter_loss_from = QuadraticCurve(loss_a, loss_b, loss_c), + max_dc_current_from = 5.0, + rating_from = rating_from, + reactive_power_limits_from = (min = -rating_from, max = rating_from), + power_factor_weighting_fraction_from = 1.0, + voltage_limits_from = (min = 0.95, max = 1.05), + reactive_power_to = 0.0, + dc_voltage_control_to = true, + ac_voltage_control_to = true, + dc_setpoint_to = 0.0, + ac_setpoint_to = 1.0, + converter_loss_to = QuadraticCurve(loss_a, loss_b, loss_c), + max_dc_current_to = 5.0, + rating_to = rating_to, + reactive_power_limits_to = (min = -rating_to, max = rating_to), + power_factor_weighting_fraction_to = 1.0, + voltage_limits_to = (min = 0.95, max = 1.05), + ) + add_component!(sys, vsc) + return sys +end + +function _build_vsc_model(formulation, network, optimizer; sys = _generate_test_vsc_sys()) + template = OperationsProblemTemplate(NetworkModel(network)) + set_device_model!(template, ThermalStandard, ThermalDispatchNoMin) + set_device_model!(template, RenewableDispatch, RenewableFullDispatch) + set_device_model!(template, PowerLoad, StaticPowerLoad) + set_device_model!(template, DeviceModel(Line, StaticBranch)) + set_device_model!(template, DeviceModel(TwoTerminalVSCLine, formulation)) + return DecisionModel( + template, sys; store_variable_names = true, optimizer = optimizer, ) end + +@testset "HVDC Two-Terminal VSC (Bin2) on DCP" begin + model = _build_vsc_model(HVDCTwoTerminalVSCBin2, DCPPowerModel, HiGHS_optimizer) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED +end + +@testset "HVDC Two-Terminal VSC (NLP) on DCP" begin + model = _build_vsc_model(HVDCTwoTerminalVSC, DCPPowerModel, ipopt_optimizer) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED +end + +@testset "HVDC Two-Terminal VSC on AC (NLP)" begin + model = _build_vsc_model(HVDCTwoTerminalVSC, ACPPowerModel, ipopt_optimizer) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED +end + +@testset "HVDC VSC Bin2 vs NLP objective agreement" begin + function _solve(formulation, optimizer) + sys = _generate_test_vsc_sys() + model = _build_vsc_model(formulation, DCPPowerModel, optimizer; sys = sys) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + return IOM.get_optimization_container(model).optimizer_stats.objective_value + end + bin2_obj = _solve(HVDCTwoTerminalVSCBin2, HiGHS_optimizer) + nlp_obj = _solve(HVDCTwoTerminalVSC, ipopt_optimizer) + # Bin2 PWL-approximates the same physics — objectives should agree closely. + @test isapprox(bin2_obj, nlp_obj; rtol = 0.05) +end + +@testset "HVDC VSC: higher cable resistance increases cost" begin + # Smaller g => larger R = 1/g => more losses => optimum should not improve. + function _solve_with_g(g_value) + sys = _generate_test_vsc_sys(; g = g_value) + model = _build_vsc_model( + HVDCTwoTerminalVSCBin2, DCPPowerModel, HiGHS_optimizer; sys = sys, + ) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + return IOM.get_optimization_container(model).optimizer_stats.objective_value + end + low_R_obj = _solve_with_g(100.0) # large g, small R + high_R_obj = _solve_with_g(20.0) # smaller g, larger R + @test high_R_obj >= low_R_obj - 1e-6 +end + +@testset "HVDC VSC: tighter PQ rating raises cost on AC" begin + function _solve_with_rating(s) + sys = _generate_test_vsc_sys(; rating_from = s, rating_to = s) + model = _build_vsc_model( + HVDCTwoTerminalVSC, ACPPowerModel, ipopt_optimizer; sys = sys, + ) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + return IOM.get_optimization_container(model).optimizer_stats.objective_value + end + looser = _solve_with_rating(2.0) + tighter = _solve_with_rating(1.0) + @test tighter >= looser - 1e-6 +end From 9d249b78786ae97722d55dd915b2d5169c16adb1 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Tue, 12 May 2026 17:35:56 -0400 Subject: [PATCH 2/7] Address PR #119 review on HVDCTwoTerminalVSC MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Rename Bin2 → MIP for VSC + IPC formulations and update tests - Make HVDCReactivePowerVariable parametric on From/To with const aliases - Wire HVDCReactivePower{From,To}Variable into ReactivePowerBalance on AC - Fix g==0 cable Ohm's law to force I==0 (open circuit, not v_f == v_t) - Add axis-aligned half-spaces to the MIP PQ-capability octagon - Default use_linear_loss=true for MIP and false for NLP via get_default_attributes; warn when NLP runs with use_linear_loss=true - Collapse abs-value decomposition into a single ModelConstructStage helper - Skip zero-coefficient terms in the quadratic converter loss expression - Fix docstring sign convention on HVDCVSCConverterPowerConstraint - Drop brittle log-file-regex linear-loss warning test Co-Authored-By: Claude Opus 4.7 --- src/PowerOperationsModels.jl | 4 +- .../branch_constructor.jl | 62 +++++---------- src/common_models/add_to_expression.jl | 72 +++++++++++++++++ src/common_models/quadratic_converter_loss.jl | 72 ++++++----------- src/core/constraints.jl | 6 +- src/core/formulations.jl | 10 +-- src/core/variables.jl | 18 ++--- src/mt_hvdc_models/HVDCsystems.jl | 30 ++------ src/mt_hvdc_models/hvdcsystems_constructor.jl | 35 +++------ .../TwoTerminalDC_branches.jl | 77 ++++++++++++------- test/test_device_hvdc.jl | 64 +++++---------- 11 files changed, 221 insertions(+), 229 deletions(-) diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl index cef1246..54380f2 100644 --- a/src/PowerOperationsModels.jl +++ b/src/PowerOperationsModels.jl @@ -764,13 +764,13 @@ export HVDCTwoTerminalDispatch export HVDCTwoTerminalPiecewiseLoss export HVDCTwoTerminalLCC export HVDCTwoTerminalVSC -export HVDCTwoTerminalVSCBin2 +export HVDCTwoTerminalVSCMIP # Converter Formulations export LosslessConverter export LinearLossConverter export AbstractQuadraticLossConverter -export Bin2QuadraticLossConverter +export MIPQuadraticLossConverter export QuadraticLossConverter # DC Line Formulations diff --git a/src/ac_transmission_models/branch_constructor.jl b/src/ac_transmission_models/branch_constructor.jl index 3da7d40..466e44c 100644 --- a/src/ac_transmission_models/branch_constructor.jl +++ b/src/ac_transmission_models/branch_constructor.jl @@ -1678,42 +1678,28 @@ end # Quadratic / bilinear approximation traits — same scheme used by the MT # converter formulations. _quad_config(::Type{HVDCTwoTerminalVSC}) = IOM.NoQuadApproxConfig() -_quad_config(::Type{HVDCTwoTerminalVSCBin2}) = +_quad_config(::Type{HVDCTwoTerminalVSCMIP}) = IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH) _bilinear_config(::Type{HVDCTwoTerminalVSC}) = IOM.NoBilinearApproxConfig() -_bilinear_config(::Type{HVDCTwoTerminalVSCBin2}) = +_bilinear_config(::Type{HVDCTwoTerminalVSCMIP}) = IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)) -# IOM's quadratic/bilinear approximation helpers take per-device bounds as -# `Vector{IOM.MinMax}` (one (min, max) pair per device, same order as `devices`). function _vsc_v_from_bounds(devices) - n = length(devices) - bounds = Vector{IOM.MinMax}(undef, n) - for (k, d) in enumerate(devices) - lims = PSY.get_voltage_limits_from(d) - bounds[k] = IOM.MinMax((min = lims.min, max = lims.max)) - end - return bounds + return [ + (min = PSY.get_voltage_limits_from(d).min, max = PSY.get_voltage_limits_from(d).max) + for d in devices + ] end function _vsc_v_to_bounds(devices) - n = length(devices) - bounds = Vector{IOM.MinMax}(undef, n) - for (k, d) in enumerate(devices) - lims = PSY.get_voltage_limits_to(d) - bounds[k] = IOM.MinMax((min = lims.min, max = lims.max)) - end - return bounds + return [ + (min = PSY.get_voltage_limits_to(d).min, max = PSY.get_voltage_limits_to(d).max) + for d in devices + ] end function _vsc_i_bounds(devices) - n = length(devices) - bounds = Vector{IOM.MinMax}(undef, n) - for (k, d) in enumerate(devices) - i_max = _vsc_shared_i_max(d) - bounds[k] = IOM.MinMax((min = -i_max, max = i_max)) - end - return bounds + return [(min = -_vsc_shared_i_max(d), max = _vsc_shared_i_max(d)) for d in devices] end function construct_device!( @@ -1733,15 +1719,6 @@ function construct_device!( _maybe_add_reactive_power_variables!(container, devices, device_model, network_model) - if _use_linear_loss(F, device_model) - ll_devices = _devices_with_linear_loss(devices) - if isempty(ll_devices) - @warn "use_linear_loss is enabled but every TwoTerminalVSCLine has zero proportional loss terms; no linear-loss variables/constraints will be added." - else - _add_abs_value_decomposition_variables!(container, ll_devices, device_model) - end - end - add_to_expression!( container, ActivePowerBalance, FlowActivePowerFromToVariable, devices, device_model, network_model, @@ -1810,14 +1787,15 @@ function construct_device!( ) _maybe_add_reactive_power_constraints!(container, devices, device_model, network_model) - if _use_linear_loss(F, device_model) - ll_devices = _devices_with_linear_loss(devices) - if !isempty(ll_devices) - _add_abs_value_decomposition_constraints!( - container, ll_devices, device_model, network_model, - DCLineCurrentFlowVariable, _vsc_shared_i_max, - ) - end + use_ll = get_attribute(device_model, "use_linear_loss") + if F === HVDCTwoTerminalVSC && use_ll + @warn "use_linear_loss = true on HVDCTwoTerminalVSC introduces a binary direction variable; the model is no longer a smooth NLP and requires a MINLP-capable solver." + end + if use_ll + _add_abs_value_decomposition!( + container, devices, device_model, network_model, + DCLineCurrentFlowVariable, _vsc_shared_i_max, + ) end add_constraint_dual!(container, sys, device_model) diff --git a/src/common_models/add_to_expression.jl b/src/common_models/add_to_expression.jl index e8e9c7e..6e00161 100644 --- a/src/common_models/add_to_expression.jl +++ b/src/common_models/add_to_expression.jl @@ -718,6 +718,78 @@ function add_to_expression!( return end +""" +Two-terminal VSC implementation to add ReactivePowerBalance for HVDCReactivePowerVariable{From} +""" +function add_to_expression!( + container::OptimizationContainer, + ::Type{T}, + ::Type{U}, + devices::IS.FlattenIteratorWrapper{V}, + ::DeviceModel{V, W}, + network_model::NetworkModel{X}, +) where { + T <: ReactivePowerBalance, + U <: HVDCReactivePowerVariable{From}, + V <: PSY.TwoTerminalVSCLine, + W <: AbstractTwoTerminalVSCFormulation, + X <: ACPPowerModel, +} + var = get_variable(container, U, V) + nodal_expr = get_expression(container, T, PSY.ACBus) + network_reduction = get_network_reduction(network_model) + time_steps = get_time_steps(container) + for d in devices + name = PSY.get_name(d) + bus_no_from = + PNM.get_mapped_bus_number(network_reduction, PSY.get_arc(d).from) + for t in time_steps + add_proportional_to_jump_expression!( + nodal_expr[bus_no_from, t], + var[name, t], + -1.0, + ) + end + end + return +end + +""" +Two-terminal VSC implementation to add ReactivePowerBalance for HVDCReactivePowerVariable{To} +""" +function add_to_expression!( + container::OptimizationContainer, + ::Type{T}, + ::Type{U}, + devices::IS.FlattenIteratorWrapper{V}, + ::DeviceModel{V, W}, + network_model::NetworkModel{X}, +) where { + T <: ReactivePowerBalance, + U <: HVDCReactivePowerVariable{To}, + V <: PSY.TwoTerminalVSCLine, + W <: AbstractTwoTerminalVSCFormulation, + X <: ACPPowerModel, +} + var = get_variable(container, U, V) + nodal_expr = get_expression(container, T, PSY.ACBus) + network_reduction = get_network_reduction(network_model) + time_steps = get_time_steps(container) + for d in devices + name = PSY.get_name(d) + bus_no_to = + PNM.get_mapped_bus_number(network_reduction, PSY.get_arc(d).to) + for t in time_steps + add_proportional_to_jump_expression!( + nodal_expr[bus_no_to, t], + var[name, t], + -1.0, + ) + end + end + return +end + """ PWL implementation to add FromTo branch variables to SystemBalanceExpressions """ diff --git a/src/common_models/quadratic_converter_loss.jl b/src/common_models/quadratic_converter_loss.jl index b2ad2eb..b9db0b8 100644 --- a/src/common_models/quadratic_converter_loss.jl +++ b/src/common_models/quadratic_converter_loss.jl @@ -1,7 +1,7 @@ # Shared helpers for quadratic / two-term converter losses # loss(I) = a * I^2 + b * |I| + c # Used by multi-terminal InterconnectingConverter formulations -# (Bin2QuadraticLossConverter, QuadraticLossConverter) and two-terminal +# (MIPQuadraticLossConverter, QuadraticLossConverter) and two-terminal # HVDCTwoTerminalVSC formulations. ######################################### @@ -11,79 +11,53 @@ _get_quadratic_term(loss_fn::PSY.QuadraticCurve) = PSY.get_quadratic_term(loss_fn) _get_quadratic_term(loss_fn) = 0.0 -# Whether the formulation wants the b*|I| linear-loss term (which requires -# decomposing the current into positive/negative parts with a direction binary). -_use_linear_loss(::Type{Bin2QuadraticLossConverter}, _) = true -_use_linear_loss(::Type{QuadraticLossConverter}, model) = - get_attribute(model, "use_linear_loss") -_use_linear_loss(::Type{HVDCTwoTerminalVSCBin2}, _) = true -_use_linear_loss(::Type{HVDCTwoTerminalVSC}, model) = - get_attribute(model, "use_linear_loss") - -# Per-device test: does this device have a nonzero linear loss term anywhere? -# Dispatched on device type because different PSY devices store loss curves on -# different fields. _has_linear_loss(d::PSY.InterconnectingConverter) = !iszero(PSY.get_proportional_term(PSY.get_loss_function(d))) _has_linear_loss(d::PSY.TwoTerminalVSCLine) = !iszero(PSY.get_proportional_term(PSY.get_converter_loss_from(d))) || !iszero(PSY.get_proportional_term(PSY.get_converter_loss_to(d))) -function _devices_with_linear_loss(devices) - return [d for d in devices if _has_linear_loss(d)] -end +_devices_with_linear_loss(devices) = filter(_has_linear_loss, devices) ######################################### ######## Loss expression builder ######## ######################################### -# Returns the JuMP expression a*i_sq + b*(i_pos + i_neg) + c -# for a single (device, time). The b*(i_pos+i_neg) term is included only when -# the formulation has opted into the linear-loss path AND b is nonzero for -# this specific device. function _quadratic_converter_loss_expr( - a, b, c, i_sq_t, i_pos_t, i_neg_t; use_linear_loss::Bool, + a::Float64, b::Float64, c::Float64, + i_sq_t, i_pos_t, i_neg_t; + use_linear_loss::Bool, ) - loss = a * i_sq_t + c - if use_linear_loss && !iszero(b) - loss += b * (i_pos_t + i_neg_t) - end - return loss + quad = iszero(a) ? 0 : a * i_sq_t + lin = (use_linear_loss && !iszero(b)) ? b * (i_pos_t + i_neg_t) : 0 + const_term = iszero(c) ? 0 : c + return quad + lin + const_term end ######################################### ####### Abs-value decomposition ######### ######################################### -# Adds the three variables (PositiveCurrent, NegativeCurrent, CurrentDirection) -# that decompose a signed current variable into i = i^+ - i^- with a binary -# direction indicator. Called from the ArgumentConstructStage. -function _add_abs_value_decomposition_variables!( - container::OptimizationContainer, - devices, - ::DeviceModel{D, F}, -) where {D <: PSY.Device, F} - add_variables!(container, PositiveCurrent, devices, F) - add_variables!(container, NegativeCurrent, devices, F) - add_variables!(container, CurrentDirection, devices, F) - return -end - -# Adds the three constraints implementing i = i^+ - i^- with the big-M -# direction binary bounds. The CurrentAbsoluteValueConstraint container is -# created internally with three meta-tagged sub-containers ("", "pos_ub", -# "neg_ub"). Caller passes the parent current variable type and a function -# `d -> i_max` so device-specific bound lookups stay device-specific. -function _add_abs_value_decomposition_constraints!( +function _add_abs_value_decomposition!( container::OptimizationContainer, devices, - ::DeviceModel{D, F}, + model::DeviceModel{D, F}, ::NetworkModel{<:AbstractPowerModel}, parent_var_type::Type{<:VariableType}, i_max_getter::Function, ) where {D <: PSY.Device, F} + ll_devices = _devices_with_linear_loss(devices) + if isempty(ll_devices) + @warn "use_linear_loss is enabled but no $(D) has a nonzero proportional loss term; no linear-loss variables/constraints will be added." + return + end + + add_variables!(container, PositiveCurrent, ll_devices, F) + add_variables!(container, NegativeCurrent, ll_devices, F) + add_variables!(container, CurrentDirection, ll_devices, F) + time_steps = get_time_steps(container) - names = [PSY.get_name(d) for d in devices] + names = [PSY.get_name(d) for d in ll_devices] jump_model = get_jump_model(container) i_var = get_variable(container, parent_var_type, D) i_pos_var = get_variable(container, PositiveCurrent, D) @@ -102,7 +76,7 @@ function _add_abs_value_decomposition_constraints!( meta = "neg_ub", ) - for d in devices + for d in ll_devices name = PSY.get_name(d) i_max = i_max_getter(d) for t in time_steps diff --git a/src/core/constraints.jl b/src/core/constraints.jl index 107a794..b0d54d5 100644 --- a/src/core/constraints.jl +++ b/src/core/constraints.jl @@ -424,8 +424,8 @@ Per-terminal converter power-balance constraint for two-terminal VSC HVDC: ```math \\begin{aligned} -p_{ft} &= v_f \\cdot I - (a_f I^2 + b_f |I| + c_f) \\\\ -p_{tf} &= -v_t \\cdot I - (a_t I^2 + b_t |I| + c_t) +p_{ft} &= v_f \\cdot I + (a_f I^2 + b_f |I| + c_f) \\\\ +p_{tf} &= -v_t \\cdot I + (a_t I^2 + b_t |I| + c_t) \\end{aligned} ``` """ @@ -448,7 +448,7 @@ on AC networks. For `HVDCTwoTerminalVSC` this is the exact circle: p_k^2 + q_k^2 \\le (S_k^{\\max})^2 \\quad k \\in \\{f, t\\} ``` -For `HVDCTwoTerminalVSCBin2` it is replaced by an inscribed polygon to stay MILP. +For `HVDCTwoTerminalVSCMIP` it is replaced by an inscribed polygon to stay MILP. """ struct HVDCVSCReactiveCapabilityConstraint <: ConstraintType end diff --git a/src/core/formulations.jl b/src/core/formulations.jl index 533f84e..e604e4e 100644 --- a/src/core/formulations.jl +++ b/src/core/formulations.jl @@ -202,9 +202,9 @@ struct HVDCTwoTerminalVSC <: AbstractTwoTerminalVSCFormulation end """ Two-terminal VSC formulation that approximates the bilinear ``v \\cdot I`` and quadratic ``I^2`` terms with SOS2-based piecewise-linear envelopes (the same -Bin2 scheme used by `Bin2QuadraticLossConverter`). Stays MILP. +scheme used by `MIPQuadraticLossConverter`). Stays MILP. """ -struct HVDCTwoTerminalVSCBin2 <: AbstractTwoTerminalVSCFormulation end +struct HVDCTwoTerminalVSCMIP <: AbstractTwoTerminalVSCFormulation end ############################### AC/DC Converter Formulations ##################################### abstract type AbstractConverterFormulation <: AbstractDeviceFormulation end @@ -225,10 +225,10 @@ Abstract supertype for InterconnectingConverter formulations with quadratic loss abstract type AbstractQuadraticLossConverter <: AbstractConverterFormulation end """ -Quadratic Loss InterconnectingConverter using the Bin2 separable bilinear approximation -(`v·i = ½((v+i)² − v² − i²)`) with a SOS2-based PWL approximation for x². +Quadratic Loss InterconnectingConverter using the separable bilinear approximation +(`v·i = ½((v+i)² − v² − i²)`) with a SOS2-based PWL approximation for x². Stays MILP. """ -struct Bin2QuadraticLossConverter <: AbstractQuadraticLossConverter end +struct MIPQuadraticLossConverter <: AbstractQuadraticLossConverter end """ Quadratic Loss InterconnectingConverter using exact bilinear (v·i) and quadratic (i²) diff --git a/src/core/variables.jl b/src/core/variables.jl index 44ab5bb..74448af 100644 --- a/src/core/variables.jl +++ b/src/core/variables.jl @@ -401,19 +401,19 @@ Docs abbreviation: ``v_t^{dc}`` """ struct HVDCToDCVoltage <: VariableType end -""" -Reactive power injected at the from-terminal AC bus by a two-terminal HVDC link. -Added only when the formulation runs on an AC network model. -Docs abbreviation: ``q_f`` -""" -struct HVDCReactivePowerFromVariable <: VariableType end +abstract type From end +abstract type To end """ -Reactive power injected at the to-terminal AC bus by a two-terminal HVDC link. +Reactive power injected at the from- (`HVDCReactivePowerVariable{From}`) or to- +(`HVDCReactivePowerVariable{To}`) terminal AC bus by a two-terminal HVDC link. Added only when the formulation runs on an AC network model. -Docs abbreviation: ``q_t`` +Docs abbreviation: ``q_f`` / ``q_t`` """ -struct HVDCReactivePowerToVariable <: VariableType end +struct HVDCReactivePowerVariable{Side} <: VariableType end + +const HVDCReactivePowerFromVariable = HVDCReactivePowerVariable{From} +const HVDCReactivePowerToVariable = HVDCReactivePowerVariable{To} """ Struct to dispatch the creation of Interface Flow Slack Up variables diff --git a/src/mt_hvdc_models/HVDCsystems.jl b/src/mt_hvdc_models/HVDCsystems.jl index 3d197d8..d2bc1ed 100644 --- a/src/mt_hvdc_models/HVDCsystems.jl +++ b/src/mt_hvdc_models/HVDCsystems.jl @@ -140,18 +140,16 @@ get_variable_upper_bound(::Type{NegativeCurrent}, d::PSY.InterconnectingConverte function get_default_attributes( ::Type{PSY.InterconnectingConverter}, - ::Type{Bin2QuadraticLossConverter}, + ::Type{MIPQuadraticLossConverter}, ) - return Dict{String, Any}() + return Dict{String, Any}("use_linear_loss" => true) end function get_default_attributes( ::Type{PSY.InterconnectingConverter}, ::Type{QuadraticLossConverter}, ) - return Dict{String, Any}( - "use_linear_loss" => false, - ) + return Dict{String, Any}("use_linear_loss" => false) end #! format: on @@ -556,7 +554,8 @@ function add_constraints!( vi_expr = get_expression(container, IOM.BilinearProductExpression, U, "vi") i_sq_expr = get_expression(container, IOM.QuadraticExpression, U, "i_sq") use_linear_loss = - _use_linear_loss(V, model) && !isempty(_devices_with_linear_loss(devices)) + get_attribute(model, "use_linear_loss") && + !isempty(_devices_with_linear_loss(devices)) if use_linear_loss i_pos_var = get_variable(container, PositiveCurrent, U) i_neg_var = get_variable(container, NegativeCurrent, U) @@ -590,25 +589,6 @@ function add_constraints!( return end -function add_constraints!( - container::OptimizationContainer, - ::Type{T}, - devices::W, - model::DeviceModel{U, V}, - network_model::NetworkModel{<:AbstractPowerModel}, -) where { - T <: CurrentAbsoluteValueConstraint, - U <: PSY.InterconnectingConverter, - V <: AbstractQuadraticLossConverter, - W <: Union{Vector{U}, IS.FlattenIteratorWrapper{U}}, -} - _add_abs_value_decomposition_constraints!( - container, devices, model, network_model, - ConverterCurrent, PSY.get_max_dc_current, - ) - return -end - ############################################ ########### Objective Function ############# ############################################ diff --git a/src/mt_hvdc_models/hvdcsystems_constructor.jl b/src/mt_hvdc_models/hvdcsystems_constructor.jl index 249d856..5c77594 100644 --- a/src/mt_hvdc_models/hvdcsystems_constructor.jl +++ b/src/mt_hvdc_models/hvdcsystems_constructor.jl @@ -54,16 +54,6 @@ function construct_device!( devices = get_available_components(model, sys) add_variables!(container, ActivePowerVariable, devices, T) add_variables!(container, ConverterCurrent, devices, T) - if _use_linear_loss(T, model) - ll_devices = _devices_with_linear_loss(devices) - if isempty(ll_devices) - @warn "use_linear_loss is enabled but every InterconnectingConverter has a zero proportional loss term; no linear-loss variables/constraints will be added." - else - add_variables!(container, PositiveCurrent, ll_devices, T) - add_variables!(container, NegativeCurrent, ll_devices, T) - add_variables!(container, CurrentDirection, ll_devices, T) - end - end add_to_expression!( container, ActivePowerBalance, ActivePowerVariable, @@ -104,10 +94,10 @@ function _converter_vi_bounds(devices) return v_bounds, i_bounds end -_quad_config(::Type{Bin2QuadraticLossConverter}) = +_quad_config(::Type{MIPQuadraticLossConverter}) = IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH) _quad_config(::Type{QuadraticLossConverter}) = IOM.NoQuadApproxConfig() -_bilinear_config(::Type{Bin2QuadraticLossConverter}) = +_bilinear_config(::Type{MIPQuadraticLossConverter}) = IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)) _bilinear_config(::Type{QuadraticLossConverter}) = IOM.NoBilinearApproxConfig() @@ -151,17 +141,16 @@ function construct_device!( ) add_constraints!(container, ConverterLossConstraint, devices, model, network_model) - if _use_linear_loss(T, model) - ll_devices = _devices_with_linear_loss(devices) - if !isempty(ll_devices) - add_constraints!( - container, - CurrentAbsoluteValueConstraint, - ll_devices, - model, - network_model, - ) - end + + use_ll = get_attribute(model, "use_linear_loss") + if T === QuadraticLossConverter && use_ll + @warn "use_linear_loss = true on QuadraticLossConverter introduces a binary direction variable; the model is no longer a smooth NLP and requires a MINLP-capable solver." + end + if use_ll + _add_abs_value_decomposition!( + container, devices, model, network_model, + ConverterCurrent, PSY.get_max_dc_current, + ) end add_feedforward_constraints!(container, model, devices) diff --git a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl index 4fa388e..22aae82 100644 --- a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl +++ b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl @@ -1444,8 +1444,7 @@ end get_variable_binary(::Type{DCLineCurrentFlowVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false get_variable_binary(::Type{HVDCFromDCVoltage}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false get_variable_binary(::Type{HVDCToDCVoltage}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false -get_variable_binary(::Type{HVDCReactivePowerFromVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false -get_variable_binary(::Type{HVDCReactivePowerToVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false +get_variable_binary(::Type{<:HVDCReactivePowerVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false get_variable_binary(::Type{PositiveCurrent}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false get_variable_binary(::Type{NegativeCurrent}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false get_variable_binary(::Type{CurrentDirection}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = true @@ -1500,10 +1499,18 @@ function _maybe_add_reactive_power_variables!( container::OptimizationContainer, devices, model::DeviceModel{PSY.TwoTerminalVSCLine, F}, - ::NetworkModel{<:AbstractPowerModel}, + network_model::NetworkModel{<:AbstractPowerModel}, ) where {F <: AbstractTwoTerminalVSCFormulation} add_variables!(container, HVDCReactivePowerFromVariable, devices, F) add_variables!(container, HVDCReactivePowerToVariable, devices, F) + add_to_expression!( + container, ReactivePowerBalance, HVDCReactivePowerFromVariable, + devices, model, network_model, + ) + add_to_expression!( + container, ReactivePowerBalance, HVDCReactivePowerToVariable, + devices, model, network_model, + ) return end @@ -1562,11 +1569,9 @@ function add_constraints!( for d in devices name = PSY.get_name(d) g = PSY.get_g(d) - # If g == 0 we treat the cable as lossless (v_f == v_t). - # Otherwise R = 1/g. for t in time_steps cons[name, t] = if iszero(g) - JuMP.@constraint(jump_model, v_f[name, t] == v_t[name, t]) + JuMP.@constraint(jump_model, i_var[name, t] == 0) else JuMP.@constraint( jump_model, @@ -1616,7 +1621,8 @@ function add_constraints!( get_expression(container, IOM.QuadraticExpression, PSY.TwoTerminalVSCLine, "i_sq") use_linear_loss = - _use_linear_loss(F, model) && !isempty(_devices_with_linear_loss(devices)) + get_attribute(model, "use_linear_loss") && + !isempty(_devices_with_linear_loss(devices)) if use_linear_loss i_pos_var = get_variable(container, PositiveCurrent, PSY.TwoTerminalVSCLine) i_neg_var = get_variable(container, NegativeCurrent, PSY.TwoTerminalVSCLine) @@ -1665,7 +1671,7 @@ function add_constraints!( return end -# PQ capability: p_k^2 + q_k^2 <= S_k^2 (NLP) or octagonal polygon (Bin2). +# PQ capability: p_k^2 + q_k^2 <= S_k^2 (NLP) or octagonal polygon (MIP). function add_constraints!( container::OptimizationContainer, ::Type{HVDCVSCReactiveCapabilityConstraint}, @@ -1712,8 +1718,6 @@ function add_constraints!( return end -# Bin2 variant: inscribed octagon in the rating circle (s/sqrt(2) half-diagonal). -# Eight linear constraints per terminal: ±p ± q ≤ s and ±p ≤ s, ±q ≤ s. function add_constraints!( container::OptimizationContainer, ::Type{HVDCVSCReactiveCapabilityConstraint}, @@ -1721,7 +1725,7 @@ function add_constraints!( Vector{PSY.TwoTerminalVSCLine}, IS.FlattenIteratorWrapper{PSY.TwoTerminalVSCLine}, }, - ::DeviceModel{PSY.TwoTerminalVSCLine, HVDCTwoTerminalVSCBin2}, + ::DeviceModel{PSY.TwoTerminalVSCLine, HVDCTwoTerminalVSCMIP}, ::NetworkModel{<:AbstractPowerModel}, ) time_steps = get_time_steps(container) @@ -1733,9 +1737,12 @@ function add_constraints!( q_f = get_variable(container, HVDCReactivePowerFromVariable, PSY.TwoTerminalVSCLine) q_t = get_variable(container, HVDCReactivePowerToVariable, PSY.TwoTerminalVSCLine) - cons = Dict{String, Any}() - for tag in ("from_pp", "from_pn", "from_np", "from_nn", + diag_tags = ("from_pp", "from_pn", "from_np", "from_nn", "to_pp", "to_pn", "to_np", "to_nn") + axis_tags = ("from_p_ub", "from_p_lb", "from_q_ub", "from_q_lb", + "to_p_ub", "to_p_lb", "to_q_ub", "to_q_lb") + cons = Dict{String, Any}() + for tag in (diag_tags..., axis_tags...) cons[tag] = add_constraints_container!( container, HVDCVSCReactiveCapabilityConstraint, PSY.TwoTerminalVSCLine, names, time_steps; meta = tag, @@ -1745,26 +1752,44 @@ function add_constraints!( inv_sqrt2 = 1.0 / sqrt(2.0) for d in devices name = PSY.get_name(d) - s_f = PSY.get_rating_from(d) * inv_sqrt2 - s_t = PSY.get_rating_to(d) * inv_sqrt2 + rating_f = PSY.get_rating_from(d) + rating_t = PSY.get_rating_to(d) + diag_f = rating_f * inv_sqrt2 * 2.0 + diag_t = rating_t * inv_sqrt2 * 2.0 for t in time_steps - # Octagon at the from terminal: project (p, q) onto 4 diagonal lines. cons["from_pp"][name, t] = - JuMP.@constraint(jump_model, p_ft[name, t] + q_f[name, t] <= 2.0 * s_f) + JuMP.@constraint(jump_model, p_ft[name, t] + q_f[name, t] <= diag_f) cons["from_pn"][name, t] = - JuMP.@constraint(jump_model, p_ft[name, t] - q_f[name, t] <= 2.0 * s_f) + JuMP.@constraint(jump_model, p_ft[name, t] - q_f[name, t] <= diag_f) cons["from_np"][name, t] = - JuMP.@constraint(jump_model, -p_ft[name, t] + q_f[name, t] <= 2.0 * s_f) + JuMP.@constraint(jump_model, -p_ft[name, t] + q_f[name, t] <= diag_f) cons["from_nn"][name, t] = - JuMP.@constraint(jump_model, -p_ft[name, t] - q_f[name, t] <= 2.0 * s_f) + JuMP.@constraint(jump_model, -p_ft[name, t] - q_f[name, t] <= diag_f) cons["to_pp"][name, t] = - JuMP.@constraint(jump_model, p_tf[name, t] + q_t[name, t] <= 2.0 * s_t) + JuMP.@constraint(jump_model, p_tf[name, t] + q_t[name, t] <= diag_t) cons["to_pn"][name, t] = - JuMP.@constraint(jump_model, p_tf[name, t] - q_t[name, t] <= 2.0 * s_t) + JuMP.@constraint(jump_model, p_tf[name, t] - q_t[name, t] <= diag_t) cons["to_np"][name, t] = - JuMP.@constraint(jump_model, -p_tf[name, t] + q_t[name, t] <= 2.0 * s_t) + JuMP.@constraint(jump_model, -p_tf[name, t] + q_t[name, t] <= diag_t) cons["to_nn"][name, t] = - JuMP.@constraint(jump_model, -p_tf[name, t] - q_t[name, t] <= 2.0 * s_t) + JuMP.@constraint(jump_model, -p_tf[name, t] - q_t[name, t] <= diag_t) + + cons["from_p_ub"][name, t] = + JuMP.@constraint(jump_model, p_ft[name, t] <= rating_f) + cons["from_p_lb"][name, t] = + JuMP.@constraint(jump_model, -p_ft[name, t] <= rating_f) + cons["from_q_ub"][name, t] = + JuMP.@constraint(jump_model, q_f[name, t] <= rating_f) + cons["from_q_lb"][name, t] = + JuMP.@constraint(jump_model, -q_f[name, t] <= rating_f) + cons["to_p_ub"][name, t] = + JuMP.@constraint(jump_model, p_tf[name, t] <= rating_t) + cons["to_p_lb"][name, t] = + JuMP.@constraint(jump_model, -p_tf[name, t] <= rating_t) + cons["to_q_ub"][name, t] = + JuMP.@constraint(jump_model, q_t[name, t] <= rating_t) + cons["to_q_lb"][name, t] = + JuMP.@constraint(jump_model, -q_t[name, t] <= rating_t) end end return @@ -1788,7 +1813,7 @@ end function get_default_attributes( ::Type{PSY.TwoTerminalVSCLine}, - ::Type{HVDCTwoTerminalVSCBin2}, + ::Type{HVDCTwoTerminalVSCMIP}, ) - return Dict{String, Any}() + return Dict{String, Any}("use_linear_loss" => true) end diff --git a/test/test_device_hvdc.jl b/test/test_device_hvdc.jl index 08f6f63..27c1f7b 100644 --- a/test/test_device_hvdc.jl +++ b/test/test_device_hvdc.jl @@ -97,7 +97,7 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED end -@testset "HVDC System with Losses Network (Bin2QuadraticLossConverter)" begin +@testset "HVDC System with Losses Network (MIPQuadraticLossConverter)" begin sys = _generate_test_hvdc_sys() template = OperationsProblemTemplate() set_device_model!(template, ThermalStandard, ThermalDispatchNoMin) @@ -106,7 +106,7 @@ end set_device_model!(template, TModelHVDCLine, DCLossyLine) ipc_model = DeviceModel( InterconnectingConverter, - Bin2QuadraticLossConverter, + MIPQuadraticLossConverter, ) set_device_model!(template, ipc_model) set_hvdc_network_model!(template, VoltageDispatchHVDCNetworkModel) @@ -147,7 +147,7 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED end -@testset "HVDC Bin2 vs NLP QuadraticLossConverter agreement" begin +@testset "HVDC MIP vs NLP QuadraticLossConverter agreement" begin function _build_and_solve(formulation, optimizer) sys = _generate_test_hvdc_sys() template = OperationsProblemTemplate() @@ -176,45 +176,13 @@ end return model end - bin2_model = _build_and_solve(Bin2QuadraticLossConverter, HiGHS_optimizer) + mip_model = _build_and_solve(MIPQuadraticLossConverter, HiGHS_optimizer) nlp_model = _build_and_solve(QuadraticLossConverter, ipopt_optimizer) - bin2_obj = IOM.get_objective_value(OptimizationProblemOutputs(bin2_model)) + mip_obj = IOM.get_objective_value(OptimizationProblemOutputs(mip_model)) nlp_obj = IOM.get_objective_value(OptimizationProblemOutputs(nlp_model)) - # Bin2 is a relaxation/PWL approximation of the exact NLP. The two objectives - # should agree to within a few percent on this small system. - @test isapprox(bin2_obj, nlp_obj; rtol = 0.05) -end - -@testset "HVDC linear-loss warning when all converters have b=0" begin - sys = _generate_test_hvdc_sys() - for ipc in get_components(InterconnectingConverter, sys) - set_loss_function!(ipc, QuadraticCurve(0.01, 0.0, 0.0)) - end - template = OperationsProblemTemplate() - set_device_model!(template, ThermalStandard, ThermalDispatchNoMin) - set_device_model!(template, PowerLoad, StaticPowerLoad) - set_device_model!(template, DeviceModel(Line, StaticBranch)) - set_device_model!(template, TModelHVDCLine, DCLossyLine) - set_device_model!( - template, - DeviceModel(InterconnectingConverter, Bin2QuadraticLossConverter), - ) - set_hvdc_network_model!(template, VoltageDispatchHVDCNetworkModel) - model = DecisionModel( - template, - sys; - store_variable_names = true, - optimizer = HiGHS_optimizer, - ) - # `build!` wraps its body in `Logging.with_logger(file_logger)`, which masks - # any `TestLogger` set up by `@test_logs`. The warning we want lands in the - # per-build log file instead, so read it back from there. - output_dir = mktempdir(; cleanup = true) - @test build!(model; output_dir = output_dir) == IOM.ModelBuildStatus.BUILT - log_path = joinpath(output_dir, IOM.PROBLEM_LOG_FILENAME) - @test occursin(r"linear[_ ]loss"i, read(log_path, String)) + @test isapprox(mip_obj, nlp_obj; rtol = 0.05) end ############################################################################## @@ -284,8 +252,8 @@ function _build_vsc_model(formulation, network, optimizer; sys = _generate_test_ ) end -@testset "HVDC Two-Terminal VSC (Bin2) on DCP" begin - model = _build_vsc_model(HVDCTwoTerminalVSCBin2, DCPPowerModel, HiGHS_optimizer) +@testset "HVDC Two-Terminal VSC (MIP) on DCP" begin + model = _build_vsc_model(HVDCTwoTerminalVSCMIP, DCPPowerModel, HiGHS_optimizer) @test build!(model; output_dir = mktempdir(; cleanup = true)) == IOM.ModelBuildStatus.BUILT @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED @@ -305,7 +273,14 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED end -@testset "HVDC VSC Bin2 vs NLP objective agreement" begin +@testset "HVDC Two-Terminal VSC (MIP) on AC" begin + model = _build_vsc_model(HVDCTwoTerminalVSCMIP, ACPPowerModel, HiGHS_optimizer) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED +end + +@testset "HVDC VSC MIP vs NLP objective agreement" begin function _solve(formulation, optimizer) sys = _generate_test_vsc_sys() model = _build_vsc_model(formulation, DCPPowerModel, optimizer; sys = sys) @@ -314,10 +289,9 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED return IOM.get_optimization_container(model).optimizer_stats.objective_value end - bin2_obj = _solve(HVDCTwoTerminalVSCBin2, HiGHS_optimizer) + mip_obj = _solve(HVDCTwoTerminalVSCMIP, HiGHS_optimizer) nlp_obj = _solve(HVDCTwoTerminalVSC, ipopt_optimizer) - # Bin2 PWL-approximates the same physics — objectives should agree closely. - @test isapprox(bin2_obj, nlp_obj; rtol = 0.05) + @test isapprox(mip_obj, nlp_obj; rtol = 0.05) end @testset "HVDC VSC: higher cable resistance increases cost" begin @@ -325,7 +299,7 @@ end function _solve_with_g(g_value) sys = _generate_test_vsc_sys(; g = g_value) model = _build_vsc_model( - HVDCTwoTerminalVSCBin2, DCPPowerModel, HiGHS_optimizer; sys = sys, + HVDCTwoTerminalVSCMIP, DCPPowerModel, HiGHS_optimizer; sys = sys, ) @test build!(model; output_dir = mktempdir(; cleanup = true)) == IOM.ModelBuildStatus.BUILT From b5649b57944478b38c1c3b63d3c13d9019ba6199 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Tue, 12 May 2026 21:46:08 -0400 Subject: [PATCH 3/7] Address PR #119 Review 2 on HVDCTwoTerminalVSC - Add HVDCTwoTerminalVSCLP formulation that owns the octagonal linear outer-approximation of the per-terminal PQ disk; document the octagon geometry and its tightness tradeoffs vs the PWL path. - Refactor HVDCTwoTerminalVSC (NLP) and HVDCTwoTerminalVSCMIP (PWL) to share a single PQ-capability add_constraints! method that writes `p_sq + q_sq <= s^2` over IOM._add_quadratic_approx! expressions. - Collapse the HVDCReactivePowerVariable{From}/{To} add_to_expression methods into one parametric method via a new `_vsc_arc_bus` helper; collapse `_vsc_v_from_bounds`/`_vsc_v_to_bounds` similarly. - Dispatch `_quadratic_converter_loss_expr` on i_sq_t (QuadExpr vs AffExpr) and accumulate via JuMP.add_to_expression!; degrade to AffExpr when a == 0 even on the QuadExpr branch. - Drop the explicit-type-check MINLP warnings on both VSC and MT HVDC paths; leave brief comments in their place. - Document the use_linear_loss defaults on all four get_default_attributes methods. - Reorder construct_device! to add the abs-value decomposition variables before the constraints that consume them (PositiveCurrent / NegativeCurrent were previously requested before they were added). - Replace `filter(_has_linear_loss, devices)` with a comprehension so it works on FlattenIteratorWrapper as well as Vector. - Omit the VSCMIP/VSCLP on ACPPowerModel build tests: HiGHS cannot handle the ACP network's trigonometric branch constraints. Co-Authored-By: Claude Opus 4.7 --- src/PowerOperationsModels.jl | 1 + .../branch_constructor.jl | 55 +++-- src/common_models/add_to_expression.jl | 51 +--- src/common_models/quadratic_converter_loss.jl | 41 +++- src/core/formulations.jl | 12 + src/mt_hvdc_models/HVDCsystems.jl | 6 + src/mt_hvdc_models/hvdcsystems_constructor.jl | 14 +- .../TwoTerminalDC_branches.jl | 232 ++++++++++++++---- test/test_device_hvdc.jl | 13 +- 9 files changed, 294 insertions(+), 131 deletions(-) diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl index 54380f2..af4387f 100644 --- a/src/PowerOperationsModels.jl +++ b/src/PowerOperationsModels.jl @@ -765,6 +765,7 @@ export HVDCTwoTerminalPiecewiseLoss export HVDCTwoTerminalLCC export HVDCTwoTerminalVSC export HVDCTwoTerminalVSCMIP +export HVDCTwoTerminalVSCLP # Converter Formulations export LosslessConverter diff --git a/src/ac_transmission_models/branch_constructor.jl b/src/ac_transmission_models/branch_constructor.jl index 466e44c..8535765 100644 --- a/src/ac_transmission_models/branch_constructor.jl +++ b/src/ac_transmission_models/branch_constructor.jl @@ -1680,22 +1680,24 @@ end _quad_config(::Type{HVDCTwoTerminalVSC}) = IOM.NoQuadApproxConfig() _quad_config(::Type{HVDCTwoTerminalVSCMIP}) = IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH) +_quad_config(::Type{HVDCTwoTerminalVSCLP}) = + IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH) _bilinear_config(::Type{HVDCTwoTerminalVSC}) = IOM.NoBilinearApproxConfig() _bilinear_config(::Type{HVDCTwoTerminalVSCMIP}) = IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)) +_bilinear_config(::Type{HVDCTwoTerminalVSCLP}) = + IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)) -function _vsc_v_from_bounds(devices) - return [ - (min = PSY.get_voltage_limits_from(d).min, max = PSY.get_voltage_limits_from(d).max) - for d in devices - ] -end +# True when the PQ-capability constraint is `p_sq + q_sq <= s^2` using IOM's +# quadratic-approx machinery (NLP exact + MIP PWL). False for the LP path, +# which uses an octagonal linear outer-approximation instead and therefore +# does not need p_sq / q_sq expressions. +_uses_pq_sq_approx(::Type{HVDCTwoTerminalVSC}) = true +_uses_pq_sq_approx(::Type{HVDCTwoTerminalVSCMIP}) = true +_uses_pq_sq_approx(::Type{HVDCTwoTerminalVSCLP}) = false -function _vsc_v_to_bounds(devices) - return [ - (min = PSY.get_voltage_limits_to(d).min, max = PSY.get_voltage_limits_to(d).max) - for d in devices - ] +function _vsc_v_bounds(devices, ::Type{S}) where {S <: Union{From, To}} + return [(min = _vsc_v_limits(d, S).min, max = _vsc_v_limits(d, S).max) for d in devices] end function _vsc_i_bounds(devices) @@ -1747,8 +1749,8 @@ function construct_device!( v_t_var = get_variable(container, HVDCToDCVoltage, PSY.TwoTerminalVSCLine) i_var = get_variable(container, DCLineCurrentFlowVariable, PSY.TwoTerminalVSCLine) - v_f_bounds = _vsc_v_from_bounds(devices) - v_t_bounds = _vsc_v_to_bounds(devices) + v_f_bounds = _vsc_v_bounds(devices, From) + v_t_bounds = _vsc_v_bounds(devices, To) i_bounds = _vsc_i_bounds(devices) quad_cfg, bilin_cfg = _quad_config(F), _bilinear_config(F) @@ -1779,18 +1781,19 @@ function construct_device!( v_t_bounds, i_bounds, "vi_tf", ) - add_constraints!( - container, HVDCCableOhmsLawConstraint, devices, device_model, network_model, + _maybe_add_pq_sq_approx!( + container, devices, line_names, time_steps, quad_cfg, F, network_model, ) - add_constraints!( - container, HVDCVSCConverterPowerConstraint, devices, device_model, network_model, - ) - _maybe_add_reactive_power_constraints!(container, devices, device_model, network_model) + # The abs-value decomposition adds PositiveCurrent / NegativeCurrent + # variables that the HVDCVSCConverterPowerConstraint reads, so it must + # run before that constraint is added. + # + # Note: use_linear_loss = true on HVDCTwoTerminalVSC adds a binary + # direction variable (CurrentDirection), making the model MINLP rather + # than smooth NLP. Caller is responsible for supplying a MINLP-capable + # solver in that case. use_ll = get_attribute(device_model, "use_linear_loss") - if F === HVDCTwoTerminalVSC && use_ll - @warn "use_linear_loss = true on HVDCTwoTerminalVSC introduces a binary direction variable; the model is no longer a smooth NLP and requires a MINLP-capable solver." - end if use_ll _add_abs_value_decomposition!( container, devices, device_model, network_model, @@ -1798,6 +1801,14 @@ function construct_device!( ) end + add_constraints!( + container, HVDCCableOhmsLawConstraint, devices, device_model, network_model, + ) + add_constraints!( + container, HVDCVSCConverterPowerConstraint, devices, device_model, network_model, + ) + _maybe_add_reactive_power_constraints!(container, devices, device_model, network_model) + add_constraint_dual!(container, sys, device_model) add_feedforward_constraints!(container, device_model, devices) return diff --git a/src/common_models/add_to_expression.jl b/src/common_models/add_to_expression.jl index 6e00161..c40e5d6 100644 --- a/src/common_models/add_to_expression.jl +++ b/src/common_models/add_to_expression.jl @@ -719,69 +719,34 @@ function add_to_expression!( end """ -Two-terminal VSC implementation to add ReactivePowerBalance for HVDCReactivePowerVariable{From} +Two-terminal VSC implementation to add ReactivePowerBalance for either side +(From/To) of HVDCReactivePowerVariable. Dispatches on the variable's `Side` +parameter via the `_vsc_arc_bus` helper. """ function add_to_expression!( container::OptimizationContainer, ::Type{T}, - ::Type{U}, - devices::IS.FlattenIteratorWrapper{V}, - ::DeviceModel{V, W}, - network_model::NetworkModel{X}, -) where { - T <: ReactivePowerBalance, - U <: HVDCReactivePowerVariable{From}, - V <: PSY.TwoTerminalVSCLine, - W <: AbstractTwoTerminalVSCFormulation, - X <: ACPPowerModel, -} - var = get_variable(container, U, V) - nodal_expr = get_expression(container, T, PSY.ACBus) - network_reduction = get_network_reduction(network_model) - time_steps = get_time_steps(container) - for d in devices - name = PSY.get_name(d) - bus_no_from = - PNM.get_mapped_bus_number(network_reduction, PSY.get_arc(d).from) - for t in time_steps - add_proportional_to_jump_expression!( - nodal_expr[bus_no_from, t], - var[name, t], - -1.0, - ) - end - end - return -end - -""" -Two-terminal VSC implementation to add ReactivePowerBalance for HVDCReactivePowerVariable{To} -""" -function add_to_expression!( - container::OptimizationContainer, - ::Type{T}, - ::Type{U}, + ::Type{HVDCReactivePowerVariable{S}}, devices::IS.FlattenIteratorWrapper{V}, ::DeviceModel{V, W}, network_model::NetworkModel{X}, ) where { T <: ReactivePowerBalance, - U <: HVDCReactivePowerVariable{To}, + S <: Union{From, To}, V <: PSY.TwoTerminalVSCLine, W <: AbstractTwoTerminalVSCFormulation, X <: ACPPowerModel, } - var = get_variable(container, U, V) + var = get_variable(container, HVDCReactivePowerVariable{S}, V) nodal_expr = get_expression(container, T, PSY.ACBus) network_reduction = get_network_reduction(network_model) time_steps = get_time_steps(container) for d in devices name = PSY.get_name(d) - bus_no_to = - PNM.get_mapped_bus_number(network_reduction, PSY.get_arc(d).to) + bus_no = PNM.get_mapped_bus_number(network_reduction, _vsc_arc_bus(d, S)) for t in time_steps add_proportional_to_jump_expression!( - nodal_expr[bus_no_to, t], + nodal_expr[bus_no, t], var[name, t], -1.0, ) diff --git a/src/common_models/quadratic_converter_loss.jl b/src/common_models/quadratic_converter_loss.jl index b9db0b8..ed8bac6 100644 --- a/src/common_models/quadratic_converter_loss.jl +++ b/src/common_models/quadratic_converter_loss.jl @@ -17,21 +17,50 @@ _has_linear_loss(d::PSY.TwoTerminalVSCLine) = !iszero(PSY.get_proportional_term(PSY.get_converter_loss_from(d))) || !iszero(PSY.get_proportional_term(PSY.get_converter_loss_to(d))) -_devices_with_linear_loss(devices) = filter(_has_linear_loss, devices) +# `filter` has no method for FlattenIteratorWrapper, so use a comprehension +# that works on any iterable. +_devices_with_linear_loss(devices) = [d for d in devices if _has_linear_loss(d)] ######################################### ######## Loss expression builder ######## ######################################### +# Dispatched on the JuMP type of `i_sq_t`: +# - JuMP.QuadExpr: exact i^2 product (NoQuadApproxConfig / NLP path). When +# `a == 0` the result has no quadratic term, so we degrade to AffExpr. +# - JuMP.AffExpr: PWL approximation of i^2 (SOS2QuadConfig / MIP path). +# In both cases we accumulate with `add_to_expression!` to avoid the +# intermediate JuMP expressions the previous `+` chain produced. function _quadratic_converter_loss_expr( a::Float64, b::Float64, c::Float64, - i_sq_t, i_pos_t, i_neg_t; + i_sq_t::JuMP.QuadExpr, i_pos_t, i_neg_t; use_linear_loss::Bool, ) - quad = iszero(a) ? 0 : a * i_sq_t - lin = (use_linear_loss && !iszero(b)) ? b * (i_pos_t + i_neg_t) : 0 - const_term = iszero(c) ? 0 : c - return quad + lin + const_term + if iszero(a) + expr = JuMP.AffExpr(c) + else + expr = JuMP.QuadExpr(JuMP.AffExpr(c)) + JuMP.add_to_expression!(expr, a, i_sq_t) + end + if use_linear_loss && !iszero(b) + JuMP.add_to_expression!(expr, b, i_pos_t) + JuMP.add_to_expression!(expr, b, i_neg_t) + end + return expr +end + +function _quadratic_converter_loss_expr( + a::Float64, b::Float64, c::Float64, + i_sq_t::JuMP.AffExpr, i_pos_t, i_neg_t; + use_linear_loss::Bool, +) + expr = JuMP.AffExpr(c) + iszero(a) || JuMP.add_to_expression!(expr, a, i_sq_t) + if use_linear_loss && !iszero(b) + JuMP.add_to_expression!(expr, b, i_pos_t) + JuMP.add_to_expression!(expr, b, i_neg_t) + end + return expr end ######################################### diff --git a/src/core/formulations.jl b/src/core/formulations.jl index e604e4e..ada69d4 100644 --- a/src/core/formulations.jl +++ b/src/core/formulations.jl @@ -206,6 +206,18 @@ scheme used by `MIPQuadraticLossConverter`). Stays MILP. """ struct HVDCTwoTerminalVSCMIP <: AbstractTwoTerminalVSCFormulation end +""" +Two-terminal VSC formulation that shares the SOS2 PWL loss model of +`HVDCTwoTerminalVSCMIP` but enforces the per-terminal PQ capability via an +octagonal linear outer-approximation of the disk ``p^2 + q^2 \\le \\text{rating}^2`` +(8 linear constraints per terminal per timestep, no extra binaries beyond +those required by the loss model). The "LP" name refers to the PQ-capability +constraint shape — it is strictly linear and is guaranteed to be a relaxation +of the disk, in contrast to the SOS2 PWL approximation `HVDCTwoTerminalVSCMIP` +uses for PQ which is tighter but not strictly outer- or inner-bounding. +""" +struct HVDCTwoTerminalVSCLP <: AbstractTwoTerminalVSCFormulation end + ############################### AC/DC Converter Formulations ##################################### abstract type AbstractConverterFormulation <: AbstractDeviceFormulation end diff --git a/src/mt_hvdc_models/HVDCsystems.jl b/src/mt_hvdc_models/HVDCsystems.jl index d2bc1ed..23982a8 100644 --- a/src/mt_hvdc_models/HVDCsystems.jl +++ b/src/mt_hvdc_models/HVDCsystems.jl @@ -138,6 +138,9 @@ get_variable_upper_bound(::Type{PositiveCurrent}, d::PSY.InterconnectingConverte get_variable_upper_bound(::Type{NegativeCurrent}, d::PSY.InterconnectingConverter,::Type{<:AbstractConverterFormulation}) = PSY.get_max_dc_current(d) +# Default `use_linear_loss = true`: the SOS2 PWL approximations already make the +# model MIP, so the extra binary direction variable from the linear-loss +# decomposition is free. function get_default_attributes( ::Type{PSY.InterconnectingConverter}, ::Type{MIPQuadraticLossConverter}, @@ -145,6 +148,9 @@ function get_default_attributes( return Dict{String, Any}("use_linear_loss" => true) end +# Default `use_linear_loss = false`: the formulation is a smooth NLP solvable by +# Ipopt; enabling the linear-loss term introduces a binary direction variable +# that pushes the model to MINLP, which pure NLP solvers cannot handle. function get_default_attributes( ::Type{PSY.InterconnectingConverter}, ::Type{QuadraticLossConverter}, diff --git a/src/mt_hvdc_models/hvdcsystems_constructor.jl b/src/mt_hvdc_models/hvdcsystems_constructor.jl index 5c77594..5872af8 100644 --- a/src/mt_hvdc_models/hvdcsystems_constructor.jl +++ b/src/mt_hvdc_models/hvdcsystems_constructor.jl @@ -140,12 +140,14 @@ function construct_device!( "vi", ) - add_constraints!(container, ConverterLossConstraint, devices, model, network_model) - + # The abs-value decomposition adds PositiveCurrent / NegativeCurrent + # variables that the ConverterLossConstraint reads, so it must run first. + # + # Note: use_linear_loss = true on QuadraticLossConverter adds a binary + # direction variable (CurrentDirection), making the model MINLP rather + # than smooth NLP. Caller is responsible for supplying a MINLP-capable + # solver in that case. use_ll = get_attribute(model, "use_linear_loss") - if T === QuadraticLossConverter && use_ll - @warn "use_linear_loss = true on QuadraticLossConverter introduces a binary direction variable; the model is no longer a smooth NLP and requires a MINLP-capable solver." - end if use_ll _add_abs_value_decomposition!( container, devices, model, network_model, @@ -153,6 +155,8 @@ function construct_device!( ) end + add_constraints!(container, ConverterLossConstraint, devices, model, network_model) + add_feedforward_constraints!(container, model, devices) add_to_objective_function!( container, devices, model, get_network_formulation(network_model), diff --git a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl index 22aae82..06566a0 100644 --- a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl +++ b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl @@ -1476,6 +1476,14 @@ get_variable_upper_bound(::Type{HVDCFromDCVoltage}, d::PSY.TwoTerminalVSCLine, : get_variable_lower_bound(::Type{HVDCToDCVoltage}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_voltage_limits_to(d).min get_variable_upper_bound(::Type{HVDCToDCVoltage}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_voltage_limits_to(d).max +# Side-parametric accessors so per-terminal methods (From/To) can be written once. +_vsc_arc_bus(d::PSY.TwoTerminalVSCLine, ::Type{From}) = PSY.get_arc(d).from +_vsc_arc_bus(d::PSY.TwoTerminalVSCLine, ::Type{To}) = PSY.get_arc(d).to +_vsc_v_limits(d::PSY.TwoTerminalVSCLine, ::Type{From}) = PSY.get_voltage_limits_from(d) +_vsc_v_limits(d::PSY.TwoTerminalVSCLine, ::Type{To}) = PSY.get_voltage_limits_to(d) +_vsc_rating(d::PSY.TwoTerminalVSCLine, ::Type{From}) = PSY.get_rating_from(d) +_vsc_rating(d::PSY.TwoTerminalVSCLine, ::Type{To}) = PSY.get_rating_to(d) + # Shared cable current bounds — must respect BOTH terminals' I_max ratings _vsc_shared_i_max(d::PSY.TwoTerminalVSCLine) = min(PSY.get_max_dc_current_from(d), PSY.get_max_dc_current_to(d)) @@ -1541,6 +1549,72 @@ _maybe_add_reactive_power_constraints!( ::NetworkModel{<:AbstractActivePowerModel}, ) = nothing +# Register the four `p_*_sq` / `q_*_sq` IOM quadratic-approx expressions used +# by the unified PQ-capability constraint (`p_sq + q_sq ≤ s²`). Only fired +# when the network actually models reactive power AND the formulation opts +# into the IOM-PWL PQ path via `_uses_pq_sq_approx`. The LP formulation +# (octagon) bypasses this entirely. +function _maybe_add_pq_sq_approx!( + container::OptimizationContainer, + devices, + line_names, + time_steps, + quad_cfg, + ::Type{F}, + ::NetworkModel{<:AbstractPowerModel}, +) where {F <: AbstractTwoTerminalVSCFormulation} + _uses_pq_sq_approx(F) || return + p_ft_var = get_variable(container, FlowActivePowerFromToVariable, PSY.TwoTerminalVSCLine) + p_tf_var = get_variable(container, FlowActivePowerToFromVariable, PSY.TwoTerminalVSCLine) + q_f_var = get_variable(container, HVDCReactivePowerFromVariable, PSY.TwoTerminalVSCLine) + q_t_var = get_variable(container, HVDCReactivePowerToVariable, PSY.TwoTerminalVSCLine) + + p_ft_bounds = [ + (min = PSY.get_active_power_limits_from(d).min, + max = PSY.get_active_power_limits_from(d).max) for d in devices + ] + p_tf_bounds = [ + (min = PSY.get_active_power_limits_to(d).min, + max = PSY.get_active_power_limits_to(d).max) for d in devices + ] + q_f_bounds = [ + (min = PSY.get_reactive_power_limits_from(d).min, + max = PSY.get_reactive_power_limits_from(d).max) for d in devices + ] + q_t_bounds = [ + (min = PSY.get_reactive_power_limits_to(d).min, + max = PSY.get_reactive_power_limits_to(d).max) for d in devices + ] + + IOM._add_quadratic_approx!( + quad_cfg, container, PSY.TwoTerminalVSCLine, + line_names, time_steps, p_ft_var, p_ft_bounds, "p_ft_sq", + ) + IOM._add_quadratic_approx!( + quad_cfg, container, PSY.TwoTerminalVSCLine, + line_names, time_steps, p_tf_var, p_tf_bounds, "p_tf_sq", + ) + IOM._add_quadratic_approx!( + quad_cfg, container, PSY.TwoTerminalVSCLine, + line_names, time_steps, q_f_var, q_f_bounds, "q_f_sq", + ) + IOM._add_quadratic_approx!( + quad_cfg, container, PSY.TwoTerminalVSCLine, + line_names, time_steps, q_t_var, q_t_bounds, "q_t_sq", + ) + return +end + +_maybe_add_pq_sq_approx!( + ::OptimizationContainer, + _devices, + _line_names, + _time_steps, + _quad_cfg, + ::Type{<:AbstractTwoTerminalVSCFormulation}, + ::NetworkModel{<:AbstractActivePowerModel}, +) = nothing + ####################### VSC core constraints ################################ # Cable Ohm's law: v_f - v_t = (1/g) * I @@ -1671,7 +1745,22 @@ function add_constraints!( return end -# PQ capability: p_k^2 + q_k^2 <= S_k^2 (NLP) or octagonal polygon (MIP). +# PQ capability for the exact-NLP (HVDCTwoTerminalVSC) and SOS2-PWL +# (HVDCTwoTerminalVSCMIP) paths. +# +# Constraint: p_sq + q_sq ≤ s² per terminal, where p_sq / q_sq are the +# expressions produced by `IOM._add_quadratic_approx!` (registered as +# `IOM.QuadraticExpression` with metas `p_*_sq` / `q_*_sq` by +# `_maybe_add_pq_sq_approx!`). +# +# • HVDCTwoTerminalVSC uses `NoQuadApproxConfig`, so `p_sq[name,t]` IS the +# exact `p_ft[name,t]^2` (JuMP.QuadExpr). The constraint reduces to +# `p² + q² ≤ s²` and the model stays a smooth NLP. +# • HVDCTwoTerminalVSCMIP uses `SolverSOS2QuadConfig(K)`, so `p_sq[name,t]` +# is an SOS2-based PWL approximation of `p²` (JuMP.AffExpr). Tightness +# grows with `K`; the approximation is not guaranteed to be a strict +# outer- or inner-bound of the disk — its position depends on PWL +# direction. function add_constraints!( container::OptimizationContainer, ::Type{HVDCVSCReactiveCapabilityConstraint}, @@ -1679,17 +1768,25 @@ function add_constraints!( Vector{PSY.TwoTerminalVSCLine}, IS.FlattenIteratorWrapper{PSY.TwoTerminalVSCLine}, }, - ::DeviceModel{PSY.TwoTerminalVSCLine, HVDCTwoTerminalVSC}, + ::DeviceModel{PSY.TwoTerminalVSCLine, F}, ::NetworkModel{<:AbstractPowerModel}, -) +) where {F <: Union{HVDCTwoTerminalVSC, HVDCTwoTerminalVSCMIP}} time_steps = get_time_steps(container) names = [PSY.get_name(d) for d in devices] jump_model = get_jump_model(container) - p_ft = get_variable(container, FlowActivePowerFromToVariable, PSY.TwoTerminalVSCLine) - p_tf = get_variable(container, FlowActivePowerToFromVariable, PSY.TwoTerminalVSCLine) - q_f = get_variable(container, HVDCReactivePowerFromVariable, PSY.TwoTerminalVSCLine) - q_t = get_variable(container, HVDCReactivePowerToVariable, PSY.TwoTerminalVSCLine) + p_ft_sq = get_expression( + container, IOM.QuadraticExpression, PSY.TwoTerminalVSCLine, "p_ft_sq", + ) + p_tf_sq = get_expression( + container, IOM.QuadraticExpression, PSY.TwoTerminalVSCLine, "p_tf_sq", + ) + q_f_sq = get_expression( + container, IOM.QuadraticExpression, PSY.TwoTerminalVSCLine, "q_f_sq", + ) + q_t_sq = get_expression( + container, IOM.QuadraticExpression, PSY.TwoTerminalVSCLine, "q_t_sq", + ) cons_f = add_constraints_container!( container, HVDCVSCReactiveCapabilityConstraint, PSY.TwoTerminalVSCLine, @@ -1702,22 +1799,46 @@ function add_constraints!( for d in devices name = PSY.get_name(d) - s_f = PSY.get_rating_from(d) - s_t = PSY.get_rating_to(d) + s_f2 = PSY.get_rating_from(d)^2 + s_t2 = PSY.get_rating_to(d)^2 for t in time_steps cons_f[name, t] = JuMP.@constraint( - jump_model, - p_ft[name, t]^2 + q_f[name, t]^2 <= s_f^2, + jump_model, p_ft_sq[name, t] + q_f_sq[name, t] <= s_f2, ) cons_t[name, t] = JuMP.@constraint( - jump_model, - p_tf[name, t]^2 + q_t[name, t]^2 <= s_t^2, + jump_model, p_tf_sq[name, t] + q_t_sq[name, t] <= s_t2, ) end end return end +# PQ capability for the linear (HVDCTwoTerminalVSCLP) path. +# +# The disk `p² + q² ≤ rating²` is approximated by a regular octagon that +# *contains* the disk: the intersection of two squares circumscribing the +# disk, where each square is tangent to the disk at four points. +# +# • Axis-aligned square: |p| ≤ rating, |q| ≤ rating +# (tangent at (±rating, 0) and (0, ±rating)) +# • 45°-rotated square: |p| + |q| ≤ rating·√2 +# (tangent at the four diagonal points) +# +# Their intersection is a regular octagon. Eight linear constraints per +# terminal per timestep (four diagonal + four axis-aligned). +# +# Tightness: the octagon is loose by at most ≈8.2% in area (octagon-to-disk +# area ratio = 8·tan(π/8)/π ≈ 1.082). It is *guaranteed* to be an outer +# approximation: any solution feasible to the real disk is feasible to the +# octagon, so the LP relaxation never blocks a physically realizable +# operating point. +# +# Tradeoff vs the unified PWL path used by `HVDCTwoTerminalVSC` / +# `HVDCTwoTerminalVSCMIP` (above): the SOS2 PWL on `p²` and `q²` separately +# can be made arbitrarily tight by increasing the breakpoint count, but the +# resulting feasible region is not guaranteed to be a strict outer- or +# inner-approximation of the disk. The octagon is the simplest always- +# conservative linear envelope and introduces no extra variables. function add_constraints!( container::OptimizationContainer, ::Type{HVDCVSCReactiveCapabilityConstraint}, @@ -1725,7 +1846,7 @@ function add_constraints!( Vector{PSY.TwoTerminalVSCLine}, IS.FlattenIteratorWrapper{PSY.TwoTerminalVSCLine}, }, - ::DeviceModel{PSY.TwoTerminalVSCLine, HVDCTwoTerminalVSCMIP}, + ::DeviceModel{PSY.TwoTerminalVSCLine, HVDCTwoTerminalVSCLP}, ::NetworkModel{<:AbstractPowerModel}, ) time_steps = get_time_steps(container) @@ -1749,47 +1870,37 @@ function add_constraints!( ) end - inv_sqrt2 = 1.0 / sqrt(2.0) + # Run the same octagon block twice (from/to) using a per-side spec, so + # the eight inner constraints aren't textually duplicated. + side_specs = ( + (prefix = "from", p_var = p_ft, q_var = q_f, side = From), + (prefix = "to", p_var = p_tf, q_var = q_t, side = To), + ) for d in devices name = PSY.get_name(d) - rating_f = PSY.get_rating_from(d) - rating_t = PSY.get_rating_to(d) - diag_f = rating_f * inv_sqrt2 * 2.0 - diag_t = rating_t * inv_sqrt2 * 2.0 - for t in time_steps - cons["from_pp"][name, t] = - JuMP.@constraint(jump_model, p_ft[name, t] + q_f[name, t] <= diag_f) - cons["from_pn"][name, t] = - JuMP.@constraint(jump_model, p_ft[name, t] - q_f[name, t] <= diag_f) - cons["from_np"][name, t] = - JuMP.@constraint(jump_model, -p_ft[name, t] + q_f[name, t] <= diag_f) - cons["from_nn"][name, t] = - JuMP.@constraint(jump_model, -p_ft[name, t] - q_f[name, t] <= diag_f) - cons["to_pp"][name, t] = - JuMP.@constraint(jump_model, p_tf[name, t] + q_t[name, t] <= diag_t) - cons["to_pn"][name, t] = - JuMP.@constraint(jump_model, p_tf[name, t] - q_t[name, t] <= diag_t) - cons["to_np"][name, t] = - JuMP.@constraint(jump_model, -p_tf[name, t] + q_t[name, t] <= diag_t) - cons["to_nn"][name, t] = - JuMP.@constraint(jump_model, -p_tf[name, t] - q_t[name, t] <= diag_t) - - cons["from_p_ub"][name, t] = - JuMP.@constraint(jump_model, p_ft[name, t] <= rating_f) - cons["from_p_lb"][name, t] = - JuMP.@constraint(jump_model, -p_ft[name, t] <= rating_f) - cons["from_q_ub"][name, t] = - JuMP.@constraint(jump_model, q_f[name, t] <= rating_f) - cons["from_q_lb"][name, t] = - JuMP.@constraint(jump_model, -q_f[name, t] <= rating_f) - cons["to_p_ub"][name, t] = - JuMP.@constraint(jump_model, p_tf[name, t] <= rating_t) - cons["to_p_lb"][name, t] = - JuMP.@constraint(jump_model, -p_tf[name, t] <= rating_t) - cons["to_q_ub"][name, t] = - JuMP.@constraint(jump_model, q_t[name, t] <= rating_t) - cons["to_q_lb"][name, t] = - JuMP.@constraint(jump_model, -q_t[name, t] <= rating_t) + for spec in side_specs + rating = _vsc_rating(d, spec.side) + diag = rating * sqrt(2.0) + prefix = spec.prefix + p_var, q_var = spec.p_var, spec.q_var + for t in time_steps + cons[prefix * "_pp"][name, t] = + JuMP.@constraint(jump_model, p_var[name, t] + q_var[name, t] <= diag) + cons[prefix * "_pn"][name, t] = + JuMP.@constraint(jump_model, p_var[name, t] - q_var[name, t] <= diag) + cons[prefix * "_np"][name, t] = + JuMP.@constraint(jump_model, -p_var[name, t] + q_var[name, t] <= diag) + cons[prefix * "_nn"][name, t] = + JuMP.@constraint(jump_model, -p_var[name, t] - q_var[name, t] <= diag) + cons[prefix * "_p_ub"][name, t] = + JuMP.@constraint(jump_model, p_var[name, t] <= rating) + cons[prefix * "_p_lb"][name, t] = + JuMP.@constraint(jump_model, -p_var[name, t] <= rating) + cons[prefix * "_q_ub"][name, t] = + JuMP.@constraint(jump_model, q_var[name, t] <= rating) + cons[prefix * "_q_lb"][name, t] = + JuMP.@constraint(jump_model, -q_var[name, t] <= rating) + end end end return @@ -1804,6 +1915,9 @@ function get_default_time_series_names( return Dict{Type{<:TimeSeriesParameter}, String}() end +# Default `use_linear_loss = false`: HVDCTwoTerminalVSC is a smooth NLP solvable +# by Ipopt; enabling the linear-loss term introduces a binary direction variable +# that pushes the model to MINLP, which pure NLP solvers cannot handle. function get_default_attributes( ::Type{PSY.TwoTerminalVSCLine}, ::Type{HVDCTwoTerminalVSC}, @@ -1811,9 +1925,21 @@ function get_default_attributes( return Dict{String, Any}("use_linear_loss" => false) end +# Default `use_linear_loss = true`: HVDCTwoTerminalVSCMIP already uses SOS2 PWL +# (MIP), so the extra binary direction variable from the linear-loss +# decomposition is free. function get_default_attributes( ::Type{PSY.TwoTerminalVSCLine}, ::Type{HVDCTwoTerminalVSCMIP}, ) return Dict{String, Any}("use_linear_loss" => true) end + +# Default `use_linear_loss = true`: same MIP rationale as VSCMIP — the SOS2 PWL +# loss model already makes the formulation MIP-compatible. +function get_default_attributes( + ::Type{PSY.TwoTerminalVSCLine}, + ::Type{HVDCTwoTerminalVSCLP}, +) + return Dict{String, Any}("use_linear_loss" => true) +end diff --git a/test/test_device_hvdc.jl b/test/test_device_hvdc.jl index 27c1f7b..a4a8c5e 100644 --- a/test/test_device_hvdc.jl +++ b/test/test_device_hvdc.jl @@ -273,13 +273,22 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED end -@testset "HVDC Two-Terminal VSC (MIP) on AC" begin - model = _build_vsc_model(HVDCTwoTerminalVSCMIP, ACPPowerModel, HiGHS_optimizer) +# Omitted: HVDCTwoTerminalVSCMIP on ACPPowerModel cannot be solved by HiGHS +# because ACPPowerModel introduces trigonometric (cos/sin) constraints in the +# branch ohms law that require an NLP-capable solver. The MIP loss model +# (SOS2 PWL) is independent of the network's nonlinearity. To exercise this +# combination we'd need a MINLP solver (e.g. Gurobi). + +@testset "HVDC Two-Terminal VSC (LP) on DCP" begin + model = _build_vsc_model(HVDCTwoTerminalVSCLP, DCPPowerModel, HiGHS_optimizer) @test build!(model; output_dir = mktempdir(; cleanup = true)) == IOM.ModelBuildStatus.BUILT @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED end +# Omitted: HVDCTwoTerminalVSCLP on ACPPowerModel — same reason as VSCMIP/AC +# above (HiGHS cannot solve the ACP network's trigonometric branch constraints). + @testset "HVDC VSC MIP vs NLP objective agreement" begin function _solve(formulation, optimizer) sys = _generate_test_vsc_sys() From c6ca90a997e03d6de639e39380b22a0479b2905a Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Wed, 13 May 2026 14:12:15 -0400 Subject: [PATCH 4/7] Address PR #119 Review 3 on HVDCTwoTerminalVSC - Drop parametric From/To types in favor of two concrete reactive-power variable structs (HVDCReactivePowerFromVariable / ToVariable); unroll the side-parametric _vsc_* accessors into two named methods each. - Rename HVDCTwoTerminalVSCMIP -> HVDCTwoTerminalVSCMILP (forward-looking for future MINLP support) and HVDCVSCReactiveCapabilityConstraint -> HVDCVSCApparentPowerLimitConstraint with a docstring covering NLP/MILP/LP paths accurately. - Split _add_abs_value_decomposition! into separate variables (ArgumentStage) and constraints (ModelStage) helpers; move the variable-creation call to ArgumentConstructStage in both the VSC and MT converter constructors. - Move _maybe_add_reactive_power_variables! / _constraints! to a new common_models/network_conditional.jl as generic helpers parameterized by variable types and constraint type, addressing issue #120. Co-Authored-By: Claude Opus 4.7 --- src/PowerOperationsModels.jl | 3 +- .../branch_constructor.jl | 49 ++++---- src/common_models/add_to_expression.jl | 47 +++++-- src/common_models/network_conditional.jl | 52 ++++++++ src/common_models/quadratic_converter_loss.jl | 24 +++- src/core/constraints.jl | 24 ++-- src/core/formulations.jl | 6 +- src/core/variables.jl | 18 +-- src/mt_hvdc_models/hvdcsystems_constructor.jl | 25 ++-- .../TwoTerminalDC_branches.jl | 115 ++++++------------ test/test_device_hvdc.jl | 8 +- 11 files changed, 223 insertions(+), 148 deletions(-) create mode 100644 src/common_models/network_conditional.jl diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl index af4387f..ade659e 100644 --- a/src/PowerOperationsModels.jl +++ b/src/PowerOperationsModels.jl @@ -214,6 +214,7 @@ include("common_models/add_parameters.jl") include("common_models/make_system_expressions.jl") include("common_models/reserve_range_constraints.jl") include("common_models/quadratic_converter_loss.jl") +include("common_models/network_conditional.jl") # Market bid cost plumbing (PSY orchestration moved out of IOM). Must be included # before device-specific files that reference MBC_TYPES / IEC_TYPES. @@ -764,7 +765,7 @@ export HVDCTwoTerminalDispatch export HVDCTwoTerminalPiecewiseLoss export HVDCTwoTerminalLCC export HVDCTwoTerminalVSC -export HVDCTwoTerminalVSCMIP +export HVDCTwoTerminalVSCMILP export HVDCTwoTerminalVSCLP # Converter Formulations diff --git a/src/ac_transmission_models/branch_constructor.jl b/src/ac_transmission_models/branch_constructor.jl index 8535765..133a487 100644 --- a/src/ac_transmission_models/branch_constructor.jl +++ b/src/ac_transmission_models/branch_constructor.jl @@ -1678,12 +1678,12 @@ end # Quadratic / bilinear approximation traits — same scheme used by the MT # converter formulations. _quad_config(::Type{HVDCTwoTerminalVSC}) = IOM.NoQuadApproxConfig() -_quad_config(::Type{HVDCTwoTerminalVSCMIP}) = +_quad_config(::Type{HVDCTwoTerminalVSCMILP}) = IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH) _quad_config(::Type{HVDCTwoTerminalVSCLP}) = IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH) _bilinear_config(::Type{HVDCTwoTerminalVSC}) = IOM.NoBilinearApproxConfig() -_bilinear_config(::Type{HVDCTwoTerminalVSCMIP}) = +_bilinear_config(::Type{HVDCTwoTerminalVSCMILP}) = IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)) _bilinear_config(::Type{HVDCTwoTerminalVSCLP}) = IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)) @@ -1693,12 +1693,13 @@ _bilinear_config(::Type{HVDCTwoTerminalVSCLP}) = # which uses an octagonal linear outer-approximation instead and therefore # does not need p_sq / q_sq expressions. _uses_pq_sq_approx(::Type{HVDCTwoTerminalVSC}) = true -_uses_pq_sq_approx(::Type{HVDCTwoTerminalVSCMIP}) = true +_uses_pq_sq_approx(::Type{HVDCTwoTerminalVSCMILP}) = true _uses_pq_sq_approx(::Type{HVDCTwoTerminalVSCLP}) = false -function _vsc_v_bounds(devices, ::Type{S}) where {S <: Union{From, To}} - return [(min = _vsc_v_limits(d, S).min, max = _vsc_v_limits(d, S).max) for d in devices] -end +_vsc_v_bounds_from(devices) = + [(min = _vsc_v_limits_from(d).min, max = _vsc_v_limits_from(d).max) for d in devices] +_vsc_v_bounds_to(devices) = + [(min = _vsc_v_limits_to(d).min, max = _vsc_v_limits_to(d).max) for d in devices] function _vsc_i_bounds(devices) return [(min = -_vsc_shared_i_max(d), max = _vsc_shared_i_max(d)) for d in devices] @@ -1719,7 +1720,19 @@ function construct_device!( add_variables!(container, HVDCFromDCVoltage, devices, F) add_variables!(container, HVDCToDCVoltage, devices, F) - _maybe_add_reactive_power_variables!(container, devices, device_model, network_model) + _maybe_add_reactive_power_variables!( + container, devices, device_model, network_model, + (HVDCReactivePowerFromVariable, HVDCReactivePowerToVariable), + ) + + # use_linear_loss = true adds a binary direction variable (CurrentDirection), + # making the model MINLP rather than smooth NLP. Caller is responsible for + # supplying a MINLP-capable solver in that case. + if get_attribute(device_model, "use_linear_loss") + _add_abs_value_decomposition_variables!( + container, devices, device_model, network_model, + ) + end add_to_expression!( container, ActivePowerBalance, FlowActivePowerFromToVariable, @@ -1749,8 +1762,8 @@ function construct_device!( v_t_var = get_variable(container, HVDCToDCVoltage, PSY.TwoTerminalVSCLine) i_var = get_variable(container, DCLineCurrentFlowVariable, PSY.TwoTerminalVSCLine) - v_f_bounds = _vsc_v_bounds(devices, From) - v_t_bounds = _vsc_v_bounds(devices, To) + v_f_bounds = _vsc_v_bounds_from(devices) + v_t_bounds = _vsc_v_bounds_to(devices) i_bounds = _vsc_i_bounds(devices) quad_cfg, bilin_cfg = _quad_config(F), _bilinear_config(F) @@ -1785,17 +1798,8 @@ function construct_device!( container, devices, line_names, time_steps, quad_cfg, F, network_model, ) - # The abs-value decomposition adds PositiveCurrent / NegativeCurrent - # variables that the HVDCVSCConverterPowerConstraint reads, so it must - # run before that constraint is added. - # - # Note: use_linear_loss = true on HVDCTwoTerminalVSC adds a binary - # direction variable (CurrentDirection), making the model MINLP rather - # than smooth NLP. Caller is responsible for supplying a MINLP-capable - # solver in that case. - use_ll = get_attribute(device_model, "use_linear_loss") - if use_ll - _add_abs_value_decomposition!( + if get_attribute(device_model, "use_linear_loss") + _add_abs_value_decomposition_constraints!( container, devices, device_model, network_model, DCLineCurrentFlowVariable, _vsc_shared_i_max, ) @@ -1807,7 +1811,10 @@ function construct_device!( add_constraints!( container, HVDCVSCConverterPowerConstraint, devices, device_model, network_model, ) - _maybe_add_reactive_power_constraints!(container, devices, device_model, network_model) + _maybe_add_reactive_power_constraints!( + container, devices, device_model, network_model, + HVDCVSCApparentPowerLimitConstraint, + ) add_constraint_dual!(container, sys, device_model) add_feedforward_constraints!(container, device_model, devices) diff --git a/src/common_models/add_to_expression.jl b/src/common_models/add_to_expression.jl index c40e5d6..0d1c9e8 100644 --- a/src/common_models/add_to_expression.jl +++ b/src/common_models/add_to_expression.jl @@ -719,31 +719,64 @@ function add_to_expression!( end """ -Two-terminal VSC implementation to add ReactivePowerBalance for either side -(From/To) of HVDCReactivePowerVariable. Dispatches on the variable's `Side` -parameter via the `_vsc_arc_bus` helper. +Two-terminal VSC: add `HVDCReactivePowerFromVariable` to `ReactivePowerBalance` +at the from-terminal AC bus. """ function add_to_expression!( container::OptimizationContainer, ::Type{T}, - ::Type{HVDCReactivePowerVariable{S}}, + ::Type{HVDCReactivePowerFromVariable}, devices::IS.FlattenIteratorWrapper{V}, ::DeviceModel{V, W}, network_model::NetworkModel{X}, ) where { T <: ReactivePowerBalance, - S <: Union{From, To}, V <: PSY.TwoTerminalVSCLine, W <: AbstractTwoTerminalVSCFormulation, X <: ACPPowerModel, } - var = get_variable(container, HVDCReactivePowerVariable{S}, V) + var = get_variable(container, HVDCReactivePowerFromVariable, V) nodal_expr = get_expression(container, T, PSY.ACBus) network_reduction = get_network_reduction(network_model) time_steps = get_time_steps(container) for d in devices name = PSY.get_name(d) - bus_no = PNM.get_mapped_bus_number(network_reduction, _vsc_arc_bus(d, S)) + bus_no = PNM.get_mapped_bus_number(network_reduction, PSY.get_arc(d).from) + for t in time_steps + add_proportional_to_jump_expression!( + nodal_expr[bus_no, t], + var[name, t], + -1.0, + ) + end + end + return +end + +""" +Two-terminal VSC: add `HVDCReactivePowerToVariable` to `ReactivePowerBalance` +at the to-terminal AC bus. +""" +function add_to_expression!( + container::OptimizationContainer, + ::Type{T}, + ::Type{HVDCReactivePowerToVariable}, + devices::IS.FlattenIteratorWrapper{V}, + ::DeviceModel{V, W}, + network_model::NetworkModel{X}, +) where { + T <: ReactivePowerBalance, + V <: PSY.TwoTerminalVSCLine, + W <: AbstractTwoTerminalVSCFormulation, + X <: ACPPowerModel, +} + var = get_variable(container, HVDCReactivePowerToVariable, V) + nodal_expr = get_expression(container, T, PSY.ACBus) + network_reduction = get_network_reduction(network_model) + time_steps = get_time_steps(container) + for d in devices + name = PSY.get_name(d) + bus_no = PNM.get_mapped_bus_number(network_reduction, PSY.get_arc(d).to) for t in time_steps add_proportional_to_jump_expression!( nodal_expr[bus_no, t], diff --git a/src/common_models/network_conditional.jl b/src/common_models/network_conditional.jl new file mode 100644 index 0000000..cc5dbb9 --- /dev/null +++ b/src/common_models/network_conditional.jl @@ -0,0 +1,52 @@ +# Helpers for variables/constraints that should only appear when the network +# model actually represents the relevant physical quantity (e.g. reactive +# power on AC networks). Generic over device type, variable types, and +# constraint type so any formulation can adopt the pattern. + +""" +Add reactive-power variables for a device and register them in the system's +`ReactivePowerBalance` expression. Only fires on AC networks +(`NetworkModel{<:AbstractPowerModel}`); no-op on active-power-only networks. + +`var_types` is a tuple/iterable of `VariableType` subtypes; each is added via +`add_variables!` and then linked into `ReactivePowerBalance` via +`add_to_expression!`. The caller's device-specific `add_to_expression!` +methods are responsible for the actual bus mapping and sign convention. +""" +function _maybe_add_reactive_power_variables!( + container::OptimizationContainer, + devices, + model::DeviceModel{D, F}, + network_model::NetworkModel{N}, + var_types, +) where {D <: PSY.Device, F, N <: AbstractPowerModel} + if N <: AbstractActivePowerModel + return + end + for V in var_types + add_variables!(container, V, devices, F) + add_to_expression!( + container, ReactivePowerBalance, V, + devices, model, network_model, + ) + end + return +end + +""" +Add a reactive-power-related constraint for a device only on AC networks +(`NetworkModel{<:AbstractPowerModel}`). No-op on active-power-only networks. +""" +function _maybe_add_reactive_power_constraints!( + container::OptimizationContainer, + devices, + model::DeviceModel{D, F}, + network_model::NetworkModel{N}, + constraint_type::Type{<:ConstraintType}, +) where {D <: PSY.Device, F, N <: AbstractPowerModel} + if N <: AbstractActivePowerModel + return + end + add_constraints!(container, constraint_type, devices, model, network_model) + return +end diff --git a/src/common_models/quadratic_converter_loss.jl b/src/common_models/quadratic_converter_loss.jl index ed8bac6..fe9fa5b 100644 --- a/src/common_models/quadratic_converter_loss.jl +++ b/src/common_models/quadratic_converter_loss.jl @@ -67,23 +67,37 @@ end ####### Abs-value decomposition ######### ######################################### -function _add_abs_value_decomposition!( +# Split into two stages so the variables are added in ArgumentConstructStage +# and the constraints in ModelConstructStage, matching the IOM construction +# convention. Both no-op when no device has a nonzero proportional loss term; +# the variables-stage helper additionally emits a configuration warning. +function _add_abs_value_decomposition_variables!( container::OptimizationContainer, devices, - model::DeviceModel{D, F}, + ::DeviceModel{D, F}, ::NetworkModel{<:AbstractPowerModel}, - parent_var_type::Type{<:VariableType}, - i_max_getter::Function, ) where {D <: PSY.Device, F} ll_devices = _devices_with_linear_loss(devices) if isempty(ll_devices) @warn "use_linear_loss is enabled but no $(D) has a nonzero proportional loss term; no linear-loss variables/constraints will be added." return end - add_variables!(container, PositiveCurrent, ll_devices, F) add_variables!(container, NegativeCurrent, ll_devices, F) add_variables!(container, CurrentDirection, ll_devices, F) + return +end + +function _add_abs_value_decomposition_constraints!( + container::OptimizationContainer, + devices, + ::DeviceModel{D, F}, + ::NetworkModel{<:AbstractPowerModel}, + parent_var_type::Type{<:VariableType}, + i_max_getter::Function, +) where {D <: PSY.Device, F} + ll_devices = _devices_with_linear_loss(devices) + isempty(ll_devices) && return time_steps = get_time_steps(container) names = [PSY.get_name(d) for d in ll_devices] diff --git a/src/core/constraints.jl b/src/core/constraints.jl index b0d54d5..6389173 100644 --- a/src/core/constraints.jl +++ b/src/core/constraints.jl @@ -441,16 +441,20 @@ v_f - v_t = (1/g) \\cdot I struct HVDCCableOhmsLawConstraint <: ConstraintType end """ -PQ capability constraint at each terminal of a two-terminal VSC HVDC, added only -on AC networks. For `HVDCTwoTerminalVSC` this is the exact circle: - -```math -p_k^2 + q_k^2 \\le (S_k^{\\max})^2 \\quad k \\in \\{f, t\\} -``` - -For `HVDCTwoTerminalVSCMIP` it is replaced by an inscribed polygon to stay MILP. -""" -struct HVDCVSCReactiveCapabilityConstraint <: ConstraintType end +Apparent-power limit at each terminal of a two-terminal VSC HVDC, added only on +AC networks. Enforces ``|S_k| \\le S_k^{\\max}`` for ``k \\in \\{f, t\\}`` via +one of three formulation-specific shapes: + +- `HVDCTwoTerminalVSC` (NLP): exact disk ``p_k^2 + q_k^2 \\le (S_k^{\\max})^2``. +- `HVDCTwoTerminalVSCMILP` (MILP): same `p_sq + q_sq ≤ s²` constraint, but with + `p_sq` and `q_sq` replaced by SOS2-based piecewise-linear approximations of + ``p^2`` and ``q^2`` (tightness grows with breakpoint count; not guaranteed to + be a strict outer- or inner-approximation). +- `HVDCTwoTerminalVSCLP` (LP): octagonal outer-approximation of the disk + (8 linear constraints per terminal per timestep, guaranteed outer-bound; + loose by at most ≈8.2% in area). +""" +struct HVDCVSCApparentPowerLimitConstraint <: ConstraintType end abstract type PowerVariableLimitsConstraint <: ConstraintType end """ diff --git a/src/core/formulations.jl b/src/core/formulations.jl index ada69d4..e58e136 100644 --- a/src/core/formulations.jl +++ b/src/core/formulations.jl @@ -204,16 +204,16 @@ Two-terminal VSC formulation that approximates the bilinear ``v \\cdot I`` and quadratic ``I^2`` terms with SOS2-based piecewise-linear envelopes (the same scheme used by `MIPQuadraticLossConverter`). Stays MILP. """ -struct HVDCTwoTerminalVSCMIP <: AbstractTwoTerminalVSCFormulation end +struct HVDCTwoTerminalVSCMILP <: AbstractTwoTerminalVSCFormulation end """ Two-terminal VSC formulation that shares the SOS2 PWL loss model of -`HVDCTwoTerminalVSCMIP` but enforces the per-terminal PQ capability via an +`HVDCTwoTerminalVSCMILP` but enforces the per-terminal PQ capability via an octagonal linear outer-approximation of the disk ``p^2 + q^2 \\le \\text{rating}^2`` (8 linear constraints per terminal per timestep, no extra binaries beyond those required by the loss model). The "LP" name refers to the PQ-capability constraint shape — it is strictly linear and is guaranteed to be a relaxation -of the disk, in contrast to the SOS2 PWL approximation `HVDCTwoTerminalVSCMIP` +of the disk, in contrast to the SOS2 PWL approximation `HVDCTwoTerminalVSCMILP` uses for PQ which is tighter but not strictly outer- or inner-bounding. """ struct HVDCTwoTerminalVSCLP <: AbstractTwoTerminalVSCFormulation end diff --git a/src/core/variables.jl b/src/core/variables.jl index 74448af..44ab5bb 100644 --- a/src/core/variables.jl +++ b/src/core/variables.jl @@ -401,19 +401,19 @@ Docs abbreviation: ``v_t^{dc}`` """ struct HVDCToDCVoltage <: VariableType end -abstract type From end -abstract type To end - """ -Reactive power injected at the from- (`HVDCReactivePowerVariable{From}`) or to- -(`HVDCReactivePowerVariable{To}`) terminal AC bus by a two-terminal HVDC link. +Reactive power injected at the from-terminal AC bus by a two-terminal HVDC link. Added only when the formulation runs on an AC network model. -Docs abbreviation: ``q_f`` / ``q_t`` +Docs abbreviation: ``q_f`` """ -struct HVDCReactivePowerVariable{Side} <: VariableType end +struct HVDCReactivePowerFromVariable <: VariableType end -const HVDCReactivePowerFromVariable = HVDCReactivePowerVariable{From} -const HVDCReactivePowerToVariable = HVDCReactivePowerVariable{To} +""" +Reactive power injected at the to-terminal AC bus by a two-terminal HVDC link. +Added only when the formulation runs on an AC network model. +Docs abbreviation: ``q_t`` +""" +struct HVDCReactivePowerToVariable <: VariableType end """ Struct to dispatch the creation of Interface Flow Slack Up variables diff --git a/src/mt_hvdc_models/hvdcsystems_constructor.jl b/src/mt_hvdc_models/hvdcsystems_constructor.jl index 5872af8..4b5779b 100644 --- a/src/mt_hvdc_models/hvdcsystems_constructor.jl +++ b/src/mt_hvdc_models/hvdcsystems_constructor.jl @@ -55,6 +55,15 @@ function construct_device!( add_variables!(container, ActivePowerVariable, devices, T) add_variables!(container, ConverterCurrent, devices, T) + # use_linear_loss = true on QuadraticLossConverter adds a binary direction + # variable (CurrentDirection), making the model MINLP rather than smooth + # NLP. Caller is responsible for supplying a MINLP-capable solver. + if get_attribute(model, "use_linear_loss") + _add_abs_value_decomposition_variables!( + container, devices, model, network_model, + ) + end + add_to_expression!( container, ActivePowerBalance, ActivePowerVariable, devices, model, network_model, @@ -140,16 +149,12 @@ function construct_device!( "vi", ) - # The abs-value decomposition adds PositiveCurrent / NegativeCurrent - # variables that the ConverterLossConstraint reads, so it must run first. - # - # Note: use_linear_loss = true on QuadraticLossConverter adds a binary - # direction variable (CurrentDirection), making the model MINLP rather - # than smooth NLP. Caller is responsible for supplying a MINLP-capable - # solver in that case. - use_ll = get_attribute(model, "use_linear_loss") - if use_ll - _add_abs_value_decomposition!( + # PositiveCurrent / NegativeCurrent / CurrentDirection variables were added + # in the ArgumentConstructStage; only the decomposition constraints need + # the JuMP model now. ConverterLossConstraint reads these variables, so + # the decomposition constraints must be added first. + if get_attribute(model, "use_linear_loss") + _add_abs_value_decomposition_constraints!( container, devices, model, network_model, ConverterCurrent, PSY.get_max_dc_current, ) diff --git a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl index 06566a0..41548fe 100644 --- a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl +++ b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl @@ -1444,7 +1444,8 @@ end get_variable_binary(::Type{DCLineCurrentFlowVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false get_variable_binary(::Type{HVDCFromDCVoltage}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false get_variable_binary(::Type{HVDCToDCVoltage}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false -get_variable_binary(::Type{<:HVDCReactivePowerVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false +get_variable_binary(::Type{HVDCReactivePowerFromVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false +get_variable_binary(::Type{HVDCReactivePowerToVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false get_variable_binary(::Type{PositiveCurrent}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false get_variable_binary(::Type{NegativeCurrent}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false get_variable_binary(::Type{CurrentDirection}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = true @@ -1476,13 +1477,12 @@ get_variable_upper_bound(::Type{HVDCFromDCVoltage}, d::PSY.TwoTerminalVSCLine, : get_variable_lower_bound(::Type{HVDCToDCVoltage}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_voltage_limits_to(d).min get_variable_upper_bound(::Type{HVDCToDCVoltage}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_voltage_limits_to(d).max -# Side-parametric accessors so per-terminal methods (From/To) can be written once. -_vsc_arc_bus(d::PSY.TwoTerminalVSCLine, ::Type{From}) = PSY.get_arc(d).from -_vsc_arc_bus(d::PSY.TwoTerminalVSCLine, ::Type{To}) = PSY.get_arc(d).to -_vsc_v_limits(d::PSY.TwoTerminalVSCLine, ::Type{From}) = PSY.get_voltage_limits_from(d) -_vsc_v_limits(d::PSY.TwoTerminalVSCLine, ::Type{To}) = PSY.get_voltage_limits_to(d) -_vsc_rating(d::PSY.TwoTerminalVSCLine, ::Type{From}) = PSY.get_rating_from(d) -_vsc_rating(d::PSY.TwoTerminalVSCLine, ::Type{To}) = PSY.get_rating_to(d) +# Per-terminal accessors. Two concrete methods rather than one side-parametric +# method to keep the public type surface narrow (no From/To tag types). +_vsc_v_limits_from(d::PSY.TwoTerminalVSCLine) = PSY.get_voltage_limits_from(d) +_vsc_v_limits_to(d::PSY.TwoTerminalVSCLine) = PSY.get_voltage_limits_to(d) +_vsc_rating_from(d::PSY.TwoTerminalVSCLine) = PSY.get_rating_from(d) +_vsc_rating_to(d::PSY.TwoTerminalVSCLine) = PSY.get_rating_to(d) # Shared cable current bounds — must respect BOTH terminals' I_max ratings _vsc_shared_i_max(d::PSY.TwoTerminalVSCLine) = @@ -1498,56 +1498,13 @@ get_variable_upper_bound(::Type{NegativeCurrent}, d::PSY.TwoTerminalVSCLine, ::T #! format: on -####################### VSC reactive-power optional path ##################### -# The reactive power and PQ capability machinery is added only when the -# network actually models reactive power. On active-only networks these -# helpers are no-ops. - -function _maybe_add_reactive_power_variables!( - container::OptimizationContainer, - devices, - model::DeviceModel{PSY.TwoTerminalVSCLine, F}, - network_model::NetworkModel{<:AbstractPowerModel}, -) where {F <: AbstractTwoTerminalVSCFormulation} - add_variables!(container, HVDCReactivePowerFromVariable, devices, F) - add_variables!(container, HVDCReactivePowerToVariable, devices, F) - add_to_expression!( - container, ReactivePowerBalance, HVDCReactivePowerFromVariable, - devices, model, network_model, - ) - add_to_expression!( - container, ReactivePowerBalance, HVDCReactivePowerToVariable, - devices, model, network_model, - ) - return -end - -_maybe_add_reactive_power_variables!( - ::OptimizationContainer, - _devices, - ::DeviceModel{PSY.TwoTerminalVSCLine, <:AbstractTwoTerminalVSCFormulation}, - ::NetworkModel{<:AbstractActivePowerModel}, -) = nothing - -function _maybe_add_reactive_power_constraints!( - container::OptimizationContainer, - devices, - model::DeviceModel{PSY.TwoTerminalVSCLine, F}, - network_model::NetworkModel{<:AbstractPowerModel}, -) where {F <: AbstractTwoTerminalVSCFormulation} - add_constraints!( - container, HVDCVSCReactiveCapabilityConstraint, - devices, model, network_model, - ) - return -end - -_maybe_add_reactive_power_constraints!( - ::OptimizationContainer, - _devices, - ::DeviceModel{PSY.TwoTerminalVSCLine, <:AbstractTwoTerminalVSCFormulation}, - ::NetworkModel{<:AbstractActivePowerModel}, -) = nothing +####################### VSC PQ-approx registration ########################### +# The generic `_maybe_add_reactive_power_variables!` / +# `_maybe_add_reactive_power_constraints!` helpers live in +# `common_models/network_conditional.jl` and are called directly from the VSC +# constructor. This section only defines the VSC-specific +# `_maybe_add_pq_sq_approx!` since it registers four named expressions that +# are specific to the VSC variables (`p_ft_sq`, `p_tf_sq`, `q_f_sq`, `q_t_sq`). # Register the four `p_*_sq` / `q_*_sq` IOM quadratic-approx expressions used # by the unified PQ-capability constraint (`p_sq + q_sq ≤ s²`). Only fired @@ -1564,26 +1521,28 @@ function _maybe_add_pq_sq_approx!( ::NetworkModel{<:AbstractPowerModel}, ) where {F <: AbstractTwoTerminalVSCFormulation} _uses_pq_sq_approx(F) || return - p_ft_var = get_variable(container, FlowActivePowerFromToVariable, PSY.TwoTerminalVSCLine) - p_tf_var = get_variable(container, FlowActivePowerToFromVariable, PSY.TwoTerminalVSCLine) + p_ft_var = + get_variable(container, FlowActivePowerFromToVariable, PSY.TwoTerminalVSCLine) + p_tf_var = + get_variable(container, FlowActivePowerToFromVariable, PSY.TwoTerminalVSCLine) q_f_var = get_variable(container, HVDCReactivePowerFromVariable, PSY.TwoTerminalVSCLine) q_t_var = get_variable(container, HVDCReactivePowerToVariable, PSY.TwoTerminalVSCLine) p_ft_bounds = [ (min = PSY.get_active_power_limits_from(d).min, - max = PSY.get_active_power_limits_from(d).max) for d in devices + max = PSY.get_active_power_limits_from(d).max) for d in devices ] p_tf_bounds = [ (min = PSY.get_active_power_limits_to(d).min, - max = PSY.get_active_power_limits_to(d).max) for d in devices + max = PSY.get_active_power_limits_to(d).max) for d in devices ] q_f_bounds = [ (min = PSY.get_reactive_power_limits_from(d).min, - max = PSY.get_reactive_power_limits_from(d).max) for d in devices + max = PSY.get_reactive_power_limits_from(d).max) for d in devices ] q_t_bounds = [ (min = PSY.get_reactive_power_limits_to(d).min, - max = PSY.get_reactive_power_limits_to(d).max) for d in devices + max = PSY.get_reactive_power_limits_to(d).max) for d in devices ] IOM._add_quadratic_approx!( @@ -1746,7 +1705,7 @@ function add_constraints!( end # PQ capability for the exact-NLP (HVDCTwoTerminalVSC) and SOS2-PWL -# (HVDCTwoTerminalVSCMIP) paths. +# (HVDCTwoTerminalVSCMILP) paths. # # Constraint: p_sq + q_sq ≤ s² per terminal, where p_sq / q_sq are the # expressions produced by `IOM._add_quadratic_approx!` (registered as @@ -1756,21 +1715,21 @@ end # • HVDCTwoTerminalVSC uses `NoQuadApproxConfig`, so `p_sq[name,t]` IS the # exact `p_ft[name,t]^2` (JuMP.QuadExpr). The constraint reduces to # `p² + q² ≤ s²` and the model stays a smooth NLP. -# • HVDCTwoTerminalVSCMIP uses `SolverSOS2QuadConfig(K)`, so `p_sq[name,t]` +# • HVDCTwoTerminalVSCMILP uses `SolverSOS2QuadConfig(K)`, so `p_sq[name,t]` # is an SOS2-based PWL approximation of `p²` (JuMP.AffExpr). Tightness # grows with `K`; the approximation is not guaranteed to be a strict # outer- or inner-bound of the disk — its position depends on PWL # direction. function add_constraints!( container::OptimizationContainer, - ::Type{HVDCVSCReactiveCapabilityConstraint}, + ::Type{HVDCVSCApparentPowerLimitConstraint}, devices::Union{ Vector{PSY.TwoTerminalVSCLine}, IS.FlattenIteratorWrapper{PSY.TwoTerminalVSCLine}, }, ::DeviceModel{PSY.TwoTerminalVSCLine, F}, ::NetworkModel{<:AbstractPowerModel}, -) where {F <: Union{HVDCTwoTerminalVSC, HVDCTwoTerminalVSCMIP}} +) where {F <: Union{HVDCTwoTerminalVSC, HVDCTwoTerminalVSCMILP}} time_steps = get_time_steps(container) names = [PSY.get_name(d) for d in devices] jump_model = get_jump_model(container) @@ -1789,11 +1748,11 @@ function add_constraints!( ) cons_f = add_constraints_container!( - container, HVDCVSCReactiveCapabilityConstraint, PSY.TwoTerminalVSCLine, + container, HVDCVSCApparentPowerLimitConstraint, PSY.TwoTerminalVSCLine, names, time_steps; meta = "from", ) cons_t = add_constraints_container!( - container, HVDCVSCReactiveCapabilityConstraint, PSY.TwoTerminalVSCLine, + container, HVDCVSCApparentPowerLimitConstraint, PSY.TwoTerminalVSCLine, names, time_steps; meta = "to", ) @@ -1834,14 +1793,14 @@ end # operating point. # # Tradeoff vs the unified PWL path used by `HVDCTwoTerminalVSC` / -# `HVDCTwoTerminalVSCMIP` (above): the SOS2 PWL on `p²` and `q²` separately +# `HVDCTwoTerminalVSCMILP` (above): the SOS2 PWL on `p²` and `q²` separately # can be made arbitrarily tight by increasing the breakpoint count, but the # resulting feasible region is not guaranteed to be a strict outer- or # inner-approximation of the disk. The octagon is the simplest always- # conservative linear envelope and introduces no extra variables. function add_constraints!( container::OptimizationContainer, - ::Type{HVDCVSCReactiveCapabilityConstraint}, + ::Type{HVDCVSCApparentPowerLimitConstraint}, devices::Union{ Vector{PSY.TwoTerminalVSCLine}, IS.FlattenIteratorWrapper{PSY.TwoTerminalVSCLine}, @@ -1865,7 +1824,7 @@ function add_constraints!( cons = Dict{String, Any}() for tag in (diag_tags..., axis_tags...) cons[tag] = add_constraints_container!( - container, HVDCVSCReactiveCapabilityConstraint, PSY.TwoTerminalVSCLine, + container, HVDCVSCApparentPowerLimitConstraint, PSY.TwoTerminalVSCLine, names, time_steps; meta = tag, ) end @@ -1873,13 +1832,13 @@ function add_constraints!( # Run the same octagon block twice (from/to) using a per-side spec, so # the eight inner constraints aren't textually duplicated. side_specs = ( - (prefix = "from", p_var = p_ft, q_var = q_f, side = From), - (prefix = "to", p_var = p_tf, q_var = q_t, side = To), + (prefix = "from", p_var = p_ft, q_var = q_f, rating_getter = _vsc_rating_from), + (prefix = "to", p_var = p_tf, q_var = q_t, rating_getter = _vsc_rating_to), ) for d in devices name = PSY.get_name(d) for spec in side_specs - rating = _vsc_rating(d, spec.side) + rating = spec.rating_getter(d) diag = rating * sqrt(2.0) prefix = spec.prefix p_var, q_var = spec.p_var, spec.q_var @@ -1925,12 +1884,12 @@ function get_default_attributes( return Dict{String, Any}("use_linear_loss" => false) end -# Default `use_linear_loss = true`: HVDCTwoTerminalVSCMIP already uses SOS2 PWL +# Default `use_linear_loss = true`: HVDCTwoTerminalVSCMILP already uses SOS2 PWL # (MIP), so the extra binary direction variable from the linear-loss # decomposition is free. function get_default_attributes( ::Type{PSY.TwoTerminalVSCLine}, - ::Type{HVDCTwoTerminalVSCMIP}, + ::Type{HVDCTwoTerminalVSCMILP}, ) return Dict{String, Any}("use_linear_loss" => true) end diff --git a/test/test_device_hvdc.jl b/test/test_device_hvdc.jl index a4a8c5e..2417240 100644 --- a/test/test_device_hvdc.jl +++ b/test/test_device_hvdc.jl @@ -253,7 +253,7 @@ function _build_vsc_model(formulation, network, optimizer; sys = _generate_test_ end @testset "HVDC Two-Terminal VSC (MIP) on DCP" begin - model = _build_vsc_model(HVDCTwoTerminalVSCMIP, DCPPowerModel, HiGHS_optimizer) + model = _build_vsc_model(HVDCTwoTerminalVSCMILP, DCPPowerModel, HiGHS_optimizer) @test build!(model; output_dir = mktempdir(; cleanup = true)) == IOM.ModelBuildStatus.BUILT @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED @@ -273,7 +273,7 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED end -# Omitted: HVDCTwoTerminalVSCMIP on ACPPowerModel cannot be solved by HiGHS +# Omitted: HVDCTwoTerminalVSCMILP on ACPPowerModel cannot be solved by HiGHS # because ACPPowerModel introduces trigonometric (cos/sin) constraints in the # branch ohms law that require an NLP-capable solver. The MIP loss model # (SOS2 PWL) is independent of the network's nonlinearity. To exercise this @@ -298,7 +298,7 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED return IOM.get_optimization_container(model).optimizer_stats.objective_value end - mip_obj = _solve(HVDCTwoTerminalVSCMIP, HiGHS_optimizer) + mip_obj = _solve(HVDCTwoTerminalVSCMILP, HiGHS_optimizer) nlp_obj = _solve(HVDCTwoTerminalVSC, ipopt_optimizer) @test isapprox(mip_obj, nlp_obj; rtol = 0.05) end @@ -308,7 +308,7 @@ end function _solve_with_g(g_value) sys = _generate_test_vsc_sys(; g = g_value) model = _build_vsc_model( - HVDCTwoTerminalVSCMIP, DCPPowerModel, HiGHS_optimizer; sys = sys, + HVDCTwoTerminalVSCMILP, DCPPowerModel, HiGHS_optimizer; sys = sys, ) @test build!(model; output_dir = mktempdir(; cleanup = true)) == IOM.ModelBuildStatus.BUILT From c4b50e6a73aed9eb2648e62b391c1efa79e7b3b6 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Thu, 14 May 2026 12:50:57 -0400 Subject: [PATCH 5/7] Address PR #119 Review 3 follow-up on HVDCTwoTerminalVSC MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Drop the HVDCTwoTerminalVSCMILP formulation: the PWL surrogate for the loss model is already MILP-grade, so a separate MILP-on-disk path didn't buy anything we can solve. The LP path now carries a `use_octagon` attribute (default `true`) that toggles the four diagonal half-planes on top of the always-on axis-aligned box, so users can pick either the octagon (≤8.2% loose linear outer-approximation of the disk) or just the box without adding new formulations. Refactors driven by the review feedback: - Rename MIPQuadraticLossConverter -> MILPQuadraticLossConverter (forward-looking for potential MINLP variants). - Delete _maybe_add_pq_sq_approx! and the per-formulation _uses_pq_sq_approx trait; the p_sq/q_sq IOM expressions are now registered by _register_pq_sq_expressions!, resolved by dispatch on (formulation, network) so only HVDCTwoTerminalVSC on an AC network registers them. - Drop the VSC PSY wrappers (_vsc_v_limits_*, _vsc_rating_*, _vsc_v_bounds_*, _vsc_i_bounds) and call PSY.get_voltage_limits_* / get_rating_* directly via broadcast. - Consolidate _vsc_shared_i_max into _linear_loss_i_max with method overloads for TwoTerminalVSCLine and InterconnectingConverter. - Collapse the From/To HVDCReactivePower* add_to_expression! pair into one method dispatching on the variable type via _vsc_q_terminal_bus, and flip the coefficient from -1.0 to +1.0 (VSC injects/consumes Q freely, unlike LCC which is strictly a consumer). - Replace the runtime `if N <: AbstractActivePowerModel` guards in network_conditional.jl with explicit no-op methods dispatched on NetworkModel{<:AbstractActivePowerModel}. - Combine the two _quadratic_converter_loss_expr methods via a JuMP-type-dispatched _loss_seed helper. - Drop the i_max_getter::Function argument from _add_abs_value_decomposition_constraints!; the getter is now picked by dispatch on the device type. - Switch the four VSC add_constraints! signatures to the parameterized Union{Vector{U}, IS.FlattenIteratorWrapper{U}} where {U <: ...} form used elsewhere in the file. Tests: - Delete the HVDCTwoTerminalVSCMILP testsets and rename the MILPQuadraticLossConverter testset accordingly. - Switch the MIP-vs-NLP agreement test to LP-vs-NLP isapprox at 5% tolerance (on DCP the PQ disk is inactive, so the LP and NLP differ only by the i² loss surrogate; the SOS2 PWL can be slightly above the exact i², so a strict ordering doesn't hold). - Add a new testset that asserts use_octagon=true on the LP path gives an objective ≥ the box-only case, pinning the new attribute's semantics. Co-Authored-By: Claude Opus 4.7 --- src/PowerOperationsModels.jl | 3 +- .../branch_constructor.jl | 36 +- src/common_models/add_to_expression.jl | 57 +-- src/common_models/network_conditional.jl | 53 +-- src/common_models/quadratic_converter_loss.jl | 46 +-- src/core/constraints.jl | 17 +- src/core/formulations.jl | 26 +- src/mt_hvdc_models/HVDCsystems.jl | 2 +- src/mt_hvdc_models/hvdcsystems_constructor.jl | 7 +- .../TwoTerminalDC_branches.jl | 350 +++++++----------- test/test_device_hvdc.jl | 75 ++-- 11 files changed, 283 insertions(+), 389 deletions(-) diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl index ade659e..9f9195e 100644 --- a/src/PowerOperationsModels.jl +++ b/src/PowerOperationsModels.jl @@ -765,14 +765,13 @@ export HVDCTwoTerminalDispatch export HVDCTwoTerminalPiecewiseLoss export HVDCTwoTerminalLCC export HVDCTwoTerminalVSC -export HVDCTwoTerminalVSCMILP export HVDCTwoTerminalVSCLP # Converter Formulations export LosslessConverter export LinearLossConverter export AbstractQuadraticLossConverter -export MIPQuadraticLossConverter +export MILPQuadraticLossConverter export QuadraticLossConverter # DC Line Formulations diff --git a/src/ac_transmission_models/branch_constructor.jl b/src/ac_transmission_models/branch_constructor.jl index 133a487..7ec35b1 100644 --- a/src/ac_transmission_models/branch_constructor.jl +++ b/src/ac_transmission_models/branch_constructor.jl @@ -1678,33 +1678,12 @@ end # Quadratic / bilinear approximation traits — same scheme used by the MT # converter formulations. _quad_config(::Type{HVDCTwoTerminalVSC}) = IOM.NoQuadApproxConfig() -_quad_config(::Type{HVDCTwoTerminalVSCMILP}) = - IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH) _quad_config(::Type{HVDCTwoTerminalVSCLP}) = IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH) _bilinear_config(::Type{HVDCTwoTerminalVSC}) = IOM.NoBilinearApproxConfig() -_bilinear_config(::Type{HVDCTwoTerminalVSCMILP}) = - IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)) _bilinear_config(::Type{HVDCTwoTerminalVSCLP}) = IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)) -# True when the PQ-capability constraint is `p_sq + q_sq <= s^2` using IOM's -# quadratic-approx machinery (NLP exact + MIP PWL). False for the LP path, -# which uses an octagonal linear outer-approximation instead and therefore -# does not need p_sq / q_sq expressions. -_uses_pq_sq_approx(::Type{HVDCTwoTerminalVSC}) = true -_uses_pq_sq_approx(::Type{HVDCTwoTerminalVSCMILP}) = true -_uses_pq_sq_approx(::Type{HVDCTwoTerminalVSCLP}) = false - -_vsc_v_bounds_from(devices) = - [(min = _vsc_v_limits_from(d).min, max = _vsc_v_limits_from(d).max) for d in devices] -_vsc_v_bounds_to(devices) = - [(min = _vsc_v_limits_to(d).min, max = _vsc_v_limits_to(d).max) for d in devices] - -function _vsc_i_bounds(devices) - return [(min = -_vsc_shared_i_max(d), max = _vsc_shared_i_max(d)) for d in devices] -end - function construct_device!( container::OptimizationContainer, sys::PSY.System, @@ -1762,9 +1741,11 @@ function construct_device!( v_t_var = get_variable(container, HVDCToDCVoltage, PSY.TwoTerminalVSCLine) i_var = get_variable(container, DCLineCurrentFlowVariable, PSY.TwoTerminalVSCLine) - v_f_bounds = _vsc_v_bounds_from(devices) - v_t_bounds = _vsc_v_bounds_to(devices) - i_bounds = _vsc_i_bounds(devices) + v_f_bounds = PSY.get_voltage_limits_from.(devices) + v_t_bounds = PSY.get_voltage_limits_to.(devices) + i_bounds = [ + (min = -_linear_loss_i_max(d), max = _linear_loss_i_max(d)) for d in devices + ] quad_cfg, bilin_cfg = _quad_config(F), _bilinear_config(F) @@ -1794,14 +1775,15 @@ function construct_device!( v_t_bounds, i_bounds, "vi_tf", ) - _maybe_add_pq_sq_approx!( - container, devices, line_names, time_steps, quad_cfg, F, network_model, + _register_pq_sq_expressions!( + container, devices, line_names, time_steps, quad_cfg, device_model, + network_model, ) if get_attribute(device_model, "use_linear_loss") _add_abs_value_decomposition_constraints!( container, devices, device_model, network_model, - DCLineCurrentFlowVariable, _vsc_shared_i_max, + DCLineCurrentFlowVariable, ) end diff --git a/src/common_models/add_to_expression.jl b/src/common_models/add_to_expression.jl index 0d1c9e8..576f056 100644 --- a/src/common_models/add_to_expression.jl +++ b/src/common_models/add_to_expression.jl @@ -718,70 +718,39 @@ function add_to_expression!( return end -""" -Two-terminal VSC: add `HVDCReactivePowerFromVariable` to `ReactivePowerBalance` -at the from-terminal AC bus. -""" -function add_to_expression!( - container::OptimizationContainer, - ::Type{T}, - ::Type{HVDCReactivePowerFromVariable}, - devices::IS.FlattenIteratorWrapper{V}, - ::DeviceModel{V, W}, - network_model::NetworkModel{X}, -) where { - T <: ReactivePowerBalance, - V <: PSY.TwoTerminalVSCLine, - W <: AbstractTwoTerminalVSCFormulation, - X <: ACPPowerModel, -} - var = get_variable(container, HVDCReactivePowerFromVariable, V) - nodal_expr = get_expression(container, T, PSY.ACBus) - network_reduction = get_network_reduction(network_model) - time_steps = get_time_steps(container) - for d in devices - name = PSY.get_name(d) - bus_no = PNM.get_mapped_bus_number(network_reduction, PSY.get_arc(d).from) - for t in time_steps - add_proportional_to_jump_expression!( - nodal_expr[bus_no, t], - var[name, t], - -1.0, - ) - end - end - return -end +# A VSC terminal can inject or consume Q freely, so the variable enters +# `ReactivePowerBalance` as a signed injection (+1.0) rather than a load (−1.0). +# Side selection picks the from- or to-terminal bus via dispatch on the +# variable type, so the body is written once. +_vsc_q_terminal_bus(d, ::Type{HVDCReactivePowerFromVariable}) = PSY.get_arc(d).from +_vsc_q_terminal_bus(d, ::Type{HVDCReactivePowerToVariable}) = PSY.get_arc(d).to -""" -Two-terminal VSC: add `HVDCReactivePowerToVariable` to `ReactivePowerBalance` -at the to-terminal AC bus. -""" function add_to_expression!( container::OptimizationContainer, ::Type{T}, - ::Type{HVDCReactivePowerToVariable}, + ::Type{U}, devices::IS.FlattenIteratorWrapper{V}, ::DeviceModel{V, W}, network_model::NetworkModel{X}, ) where { T <: ReactivePowerBalance, + U <: Union{HVDCReactivePowerFromVariable, HVDCReactivePowerToVariable}, V <: PSY.TwoTerminalVSCLine, W <: AbstractTwoTerminalVSCFormulation, X <: ACPPowerModel, } - var = get_variable(container, HVDCReactivePowerToVariable, V) + var = get_variable(container, U, V) nodal_expr = get_expression(container, T, PSY.ACBus) network_reduction = get_network_reduction(network_model) time_steps = get_time_steps(container) for d in devices name = PSY.get_name(d) - bus_no = PNM.get_mapped_bus_number(network_reduction, PSY.get_arc(d).to) + bus_no = PNM.get_mapped_bus_number( + network_reduction, _vsc_q_terminal_bus(d, U), + ) for t in time_steps add_proportional_to_jump_expression!( - nodal_expr[bus_no, t], - var[name, t], - -1.0, + nodal_expr[bus_no, t], var[name, t], 1.0, ) end end diff --git a/src/common_models/network_conditional.jl b/src/common_models/network_conditional.jl index cc5dbb9..e11a5c6 100644 --- a/src/common_models/network_conditional.jl +++ b/src/common_models/network_conditional.jl @@ -1,52 +1,59 @@ # Helpers for variables/constraints that should only appear when the network # model actually represents the relevant physical quantity (e.g. reactive -# power on AC networks). Generic over device type, variable types, and -# constraint type so any formulation can adopt the pattern. +# power on AC networks). Each helper has a no-op method that dispatches on +# `NetworkModel{<:AbstractActivePowerModel}`; Julia's method resolution picks +# the more specific no-op over the AC body for DC formulations. """ Add reactive-power variables for a device and register them in the system's -`ReactivePowerBalance` expression. Only fires on AC networks -(`NetworkModel{<:AbstractPowerModel}`); no-op on active-power-only networks. - -`var_types` is a tuple/iterable of `VariableType` subtypes; each is added via -`add_variables!` and then linked into `ReactivePowerBalance` via -`add_to_expression!`. The caller's device-specific `add_to_expression!` -methods are responsible for the actual bus mapping and sign convention. +`ReactivePowerBalance` expression. `var_types` is a tuple/iterable of +`VariableType` subtypes; each is added via `add_variables!` and then linked +into `ReactivePowerBalance` via `add_to_expression!`. The caller's +device-specific `add_to_expression!` methods are responsible for the actual +bus mapping and sign convention. """ function _maybe_add_reactive_power_variables!( container::OptimizationContainer, devices, model::DeviceModel{D, F}, - network_model::NetworkModel{N}, + network_model::NetworkModel{<:AbstractPowerModel}, var_types, -) where {D <: PSY.Device, F, N <: AbstractPowerModel} - if N <: AbstractActivePowerModel - return - end +) where {D <: PSY.Device, F} for V in var_types add_variables!(container, V, devices, F) add_to_expression!( - container, ReactivePowerBalance, V, - devices, model, network_model, + container, ReactivePowerBalance, V, devices, model, network_model, ) end return end +_maybe_add_reactive_power_variables!( + ::OptimizationContainer, + _devices, + ::DeviceModel{D, F}, + ::NetworkModel{<:AbstractActivePowerModel}, + _var_types, +) where {D <: PSY.Device, F} = nothing + """ -Add a reactive-power-related constraint for a device only on AC networks -(`NetworkModel{<:AbstractPowerModel}`). No-op on active-power-only networks. +Add a reactive-power-related constraint for a device on AC networks. """ function _maybe_add_reactive_power_constraints!( container::OptimizationContainer, devices, model::DeviceModel{D, F}, - network_model::NetworkModel{N}, + network_model::NetworkModel{<:AbstractPowerModel}, constraint_type::Type{<:ConstraintType}, -) where {D <: PSY.Device, F, N <: AbstractPowerModel} - if N <: AbstractActivePowerModel - return - end +) where {D <: PSY.Device, F} add_constraints!(container, constraint_type, devices, model, network_model) return end + +_maybe_add_reactive_power_constraints!( + ::OptimizationContainer, + _devices, + ::DeviceModel{D, F}, + ::NetworkModel{<:AbstractActivePowerModel}, + ::Type{<:ConstraintType}, +) where {D <: PSY.Device, F} = nothing diff --git a/src/common_models/quadratic_converter_loss.jl b/src/common_models/quadratic_converter_loss.jl index fe9fa5b..9d301b1 100644 --- a/src/common_models/quadratic_converter_loss.jl +++ b/src/common_models/quadratic_converter_loss.jl @@ -1,7 +1,7 @@ # Shared helpers for quadratic / two-term converter losses # loss(I) = a * I^2 + b * |I| + c # Used by multi-terminal InterconnectingConverter formulations -# (MIPQuadraticLossConverter, QuadraticLossConverter) and two-terminal +# (MILPQuadraticLossConverter, QuadraticLossConverter) and two-terminal # HVDCTwoTerminalVSC formulations. ######################################### @@ -25,36 +25,19 @@ _devices_with_linear_loss(devices) = [d for d in devices if _has_linear_loss(d)] ######## Loss expression builder ######## ######################################### -# Dispatched on the JuMP type of `i_sq_t`: -# - JuMP.QuadExpr: exact i^2 product (NoQuadApproxConfig / NLP path). When -# `a == 0` the result has no quadratic term, so we degrade to AffExpr. -# - JuMP.AffExpr: PWL approximation of i^2 (SOS2QuadConfig / MIP path). -# In both cases we accumulate with `add_to_expression!` to avoid the -# intermediate JuMP expressions the previous `+` chain produced. -function _quadratic_converter_loss_expr( - a::Float64, b::Float64, c::Float64, - i_sq_t::JuMP.QuadExpr, i_pos_t, i_neg_t; - use_linear_loss::Bool, -) - if iszero(a) - expr = JuMP.AffExpr(c) - else - expr = JuMP.QuadExpr(JuMP.AffExpr(c)) - JuMP.add_to_expression!(expr, a, i_sq_t) - end - if use_linear_loss && !iszero(b) - JuMP.add_to_expression!(expr, b, i_pos_t) - JuMP.add_to_expression!(expr, b, i_neg_t) - end - return expr -end +# `_loss_seed` picks the right JuMP expression flavor for the I^2 term up +# front: QuadExpr when `i_sq_t` is the exact bilinear i*i (NLP path, +# NoQuadApproxConfig), AffExpr when it's the SOS2 PWL surrogate (MIP path). +# With `a == 0` the I^2 term drops, so we degrade to AffExpr in either case. +_loss_seed(c::Float64, ::JuMP.QuadExpr) = JuMP.QuadExpr(JuMP.AffExpr(c)) +_loss_seed(c::Float64, ::JuMP.AffExpr) = JuMP.AffExpr(c) function _quadratic_converter_loss_expr( a::Float64, b::Float64, c::Float64, - i_sq_t::JuMP.AffExpr, i_pos_t, i_neg_t; + i_sq_t, i_pos_t, i_neg_t; use_linear_loss::Bool, ) - expr = JuMP.AffExpr(c) + expr = iszero(a) ? JuMP.AffExpr(c) : _loss_seed(c, i_sq_t) iszero(a) || JuMP.add_to_expression!(expr, a, i_sq_t) if use_linear_loss && !iszero(b) JuMP.add_to_expression!(expr, b, i_pos_t) @@ -88,13 +71,20 @@ function _add_abs_value_decomposition_variables!( return end +# I_max comes from different PSY accessors depending on the device, so we +# dispatch on the device type rather than asking the caller to thread a getter +# through. Two terminals on a VSC line share the same DC current variable, so +# we take the binding (min) rating; MT converters expose a single getter. +_linear_loss_i_max(d::PSY.TwoTerminalVSCLine) = + min(PSY.get_max_dc_current_from(d), PSY.get_max_dc_current_to(d)) +_linear_loss_i_max(d::PSY.InterconnectingConverter) = PSY.get_max_dc_current(d) + function _add_abs_value_decomposition_constraints!( container::OptimizationContainer, devices, ::DeviceModel{D, F}, ::NetworkModel{<:AbstractPowerModel}, parent_var_type::Type{<:VariableType}, - i_max_getter::Function, ) where {D <: PSY.Device, F} ll_devices = _devices_with_linear_loss(devices) isempty(ll_devices) && return @@ -121,7 +111,7 @@ function _add_abs_value_decomposition_constraints!( for d in ll_devices name = PSY.get_name(d) - i_max = i_max_getter(d) + i_max = _linear_loss_i_max(d) for t in time_steps abs_val_const[name, t] = JuMP.@constraint( jump_model, diff --git a/src/core/constraints.jl b/src/core/constraints.jl index 6389173..a6a405b 100644 --- a/src/core/constraints.jl +++ b/src/core/constraints.jl @@ -442,17 +442,16 @@ struct HVDCCableOhmsLawConstraint <: ConstraintType end """ Apparent-power limit at each terminal of a two-terminal VSC HVDC, added only on -AC networks. Enforces ``|S_k| \\le S_k^{\\max}`` for ``k \\in \\{f, t\\}`` via -one of three formulation-specific shapes: +AC networks. Enforces ``|S_k| \\le S_k^{\\max}`` for ``k \\in \\{f, t\\}`` via one +of two formulation-specific shapes: - `HVDCTwoTerminalVSC` (NLP): exact disk ``p_k^2 + q_k^2 \\le (S_k^{\\max})^2``. -- `HVDCTwoTerminalVSCMILP` (MILP): same `p_sq + q_sq ≤ s²` constraint, but with - `p_sq` and `q_sq` replaced by SOS2-based piecewise-linear approximations of - ``p^2`` and ``q^2`` (tightness grows with breakpoint count; not guaranteed to - be a strict outer- or inner-approximation). -- `HVDCTwoTerminalVSCLP` (LP): octagonal outer-approximation of the disk - (8 linear constraints per terminal per timestep, guaranteed outer-bound; - loose by at most ≈8.2% in area). +- `HVDCTwoTerminalVSCLP` (LP): linear outer-approximation. The axis-aligned box + ``|p_k|, |q_k| \\le S_k^{\\max}`` is added unconditionally; the four diagonal + half-planes ``|p_k| \\pm q_k \\le S_k^{\\max}\\sqrt{2}`` are added when the + device-model attribute `use_octagon` (default `true`) is on, in which case + the intersection is a regular octagon circumscribing the disk + (loose by at most ≈8.2% in area). """ struct HVDCVSCApparentPowerLimitConstraint <: ConstraintType end diff --git a/src/core/formulations.jl b/src/core/formulations.jl index e58e136..37b20e1 100644 --- a/src/core/formulations.jl +++ b/src/core/formulations.jl @@ -200,21 +200,15 @@ Two-terminal VSC formulation that keeps the bilinear ``v \\cdot I`` and quadrati struct HVDCTwoTerminalVSC <: AbstractTwoTerminalVSCFormulation end """ -Two-terminal VSC formulation that approximates the bilinear ``v \\cdot I`` and -quadratic ``I^2`` terms with SOS2-based piecewise-linear envelopes (the same -scheme used by `MIPQuadraticLossConverter`). Stays MILP. -""" -struct HVDCTwoTerminalVSCMILP <: AbstractTwoTerminalVSCFormulation end - -""" -Two-terminal VSC formulation that shares the SOS2 PWL loss model of -`HVDCTwoTerminalVSCMILP` but enforces the per-terminal PQ capability via an -octagonal linear outer-approximation of the disk ``p^2 + q^2 \\le \\text{rating}^2`` -(8 linear constraints per terminal per timestep, no extra binaries beyond -those required by the loss model). The "LP" name refers to the PQ-capability -constraint shape — it is strictly linear and is guaranteed to be a relaxation -of the disk, in contrast to the SOS2 PWL approximation `HVDCTwoTerminalVSCMILP` -uses for PQ which is tighter but not strictly outer- or inner-bounding. +Two-terminal VSC formulation that uses SOS2 piecewise-linear surrogates for the +bilinear ``v \\cdot I`` and quadratic ``I^2`` terms (so the loss model itself is +mixed-integer linear) and enforces the per-terminal PQ capability via a linear +outer-approximation of the disk ``p^2 + q^2 \\le \\text{rating}^2``: axis-aligned +box constraints ``|p|, |q| \\le \\text{rating}`` always, plus four diagonal +constraints ``|p| \\pm q \\le \\text{rating}\\sqrt{2}`` when the device-model +attribute `use_octagon` (default `true`) is on. With the diagonals in place the +feasible region is a regular octagon circumscribing the disk; turning them off +leaves only the box. """ struct HVDCTwoTerminalVSCLP <: AbstractTwoTerminalVSCFormulation end @@ -240,7 +234,7 @@ abstract type AbstractQuadraticLossConverter <: AbstractConverterFormulation end Quadratic Loss InterconnectingConverter using the separable bilinear approximation (`v·i = ½((v+i)² − v² − i²)`) with a SOS2-based PWL approximation for x². Stays MILP. """ -struct MIPQuadraticLossConverter <: AbstractQuadraticLossConverter end +struct MILPQuadraticLossConverter <: AbstractQuadraticLossConverter end """ Quadratic Loss InterconnectingConverter using exact bilinear (v·i) and quadratic (i²) diff --git a/src/mt_hvdc_models/HVDCsystems.jl b/src/mt_hvdc_models/HVDCsystems.jl index 23982a8..37dd418 100644 --- a/src/mt_hvdc_models/HVDCsystems.jl +++ b/src/mt_hvdc_models/HVDCsystems.jl @@ -143,7 +143,7 @@ get_variable_upper_bound(::Type{NegativeCurrent}, d::PSY.InterconnectingConverte # decomposition is free. function get_default_attributes( ::Type{PSY.InterconnectingConverter}, - ::Type{MIPQuadraticLossConverter}, + ::Type{MILPQuadraticLossConverter}, ) return Dict{String, Any}("use_linear_loss" => true) end diff --git a/src/mt_hvdc_models/hvdcsystems_constructor.jl b/src/mt_hvdc_models/hvdcsystems_constructor.jl index 4b5779b..dd15cfa 100644 --- a/src/mt_hvdc_models/hvdcsystems_constructor.jl +++ b/src/mt_hvdc_models/hvdcsystems_constructor.jl @@ -103,10 +103,10 @@ function _converter_vi_bounds(devices) return v_bounds, i_bounds end -_quad_config(::Type{MIPQuadraticLossConverter}) = +_quad_config(::Type{MILPQuadraticLossConverter}) = IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH) _quad_config(::Type{QuadraticLossConverter}) = IOM.NoQuadApproxConfig() -_bilinear_config(::Type{MIPQuadraticLossConverter}) = +_bilinear_config(::Type{MILPQuadraticLossConverter}) = IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)) _bilinear_config(::Type{QuadraticLossConverter}) = IOM.NoBilinearApproxConfig() @@ -155,8 +155,7 @@ function construct_device!( # the decomposition constraints must be added first. if get_attribute(model, "use_linear_loss") _add_abs_value_decomposition_constraints!( - container, devices, model, network_model, - ConverterCurrent, PSY.get_max_dc_current, + container, devices, model, network_model, ConverterCurrent, ) end diff --git a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl index 41548fe..2a7eaf2 100644 --- a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl +++ b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl @@ -1477,100 +1477,73 @@ get_variable_upper_bound(::Type{HVDCFromDCVoltage}, d::PSY.TwoTerminalVSCLine, : get_variable_lower_bound(::Type{HVDCToDCVoltage}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_voltage_limits_to(d).min get_variable_upper_bound(::Type{HVDCToDCVoltage}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_voltage_limits_to(d).max -# Per-terminal accessors. Two concrete methods rather than one side-parametric -# method to keep the public type surface narrow (no From/To tag types). -_vsc_v_limits_from(d::PSY.TwoTerminalVSCLine) = PSY.get_voltage_limits_from(d) -_vsc_v_limits_to(d::PSY.TwoTerminalVSCLine) = PSY.get_voltage_limits_to(d) -_vsc_rating_from(d::PSY.TwoTerminalVSCLine) = PSY.get_rating_from(d) -_vsc_rating_to(d::PSY.TwoTerminalVSCLine) = PSY.get_rating_to(d) - -# Shared cable current bounds — must respect BOTH terminals' I_max ratings -_vsc_shared_i_max(d::PSY.TwoTerminalVSCLine) = - min(PSY.get_max_dc_current_from(d), PSY.get_max_dc_current_to(d)) -get_variable_lower_bound(::Type{DCLineCurrentFlowVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = -_vsc_shared_i_max(d) -get_variable_upper_bound(::Type{DCLineCurrentFlowVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = _vsc_shared_i_max(d) +# Shared cable current bounds — must respect BOTH terminals' I_max ratings. +# `_linear_loss_i_max(::PSY.TwoTerminalVSCLine)` lives in +# common_models/quadratic_converter_loss.jl and returns the same `min(...)`. +get_variable_lower_bound(::Type{DCLineCurrentFlowVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = -_linear_loss_i_max(d) +get_variable_upper_bound(::Type{DCLineCurrentFlowVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = _linear_loss_i_max(d) # Positive/negative parts: each in [0, i_max] get_variable_lower_bound(::Type{PositiveCurrent}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = 0.0 -get_variable_upper_bound(::Type{PositiveCurrent}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = _vsc_shared_i_max(d) +get_variable_upper_bound(::Type{PositiveCurrent}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = _linear_loss_i_max(d) get_variable_lower_bound(::Type{NegativeCurrent}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = 0.0 -get_variable_upper_bound(::Type{NegativeCurrent}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = _vsc_shared_i_max(d) +get_variable_upper_bound(::Type{NegativeCurrent}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = _linear_loss_i_max(d) #! format: on ####################### VSC PQ-approx registration ########################### -# The generic `_maybe_add_reactive_power_variables!` / -# `_maybe_add_reactive_power_constraints!` helpers live in -# `common_models/network_conditional.jl` and are called directly from the VSC -# constructor. This section only defines the VSC-specific -# `_maybe_add_pq_sq_approx!` since it registers four named expressions that -# are specific to the VSC variables (`p_ft_sq`, `p_tf_sq`, `q_f_sq`, `q_t_sq`). - -# Register the four `p_*_sq` / `q_*_sq` IOM quadratic-approx expressions used -# by the unified PQ-capability constraint (`p_sq + q_sq ≤ s²`). Only fired -# when the network actually models reactive power AND the formulation opts -# into the IOM-PWL PQ path via `_uses_pq_sq_approx`. The LP formulation -# (octagon) bypasses this entirely. -function _maybe_add_pq_sq_approx!( + +# Register the four IOM `QuadraticExpression` handles (`p_ft_sq`, `p_tf_sq`, +# `q_f_sq`, `q_t_sq`) that the disk constraint reads. Only the NLP formulation +# on an AC network actually emits the disk constraint, so only that combo +# registers anything; the LP path and active-power-only networks no-op. +function _register_pq_sq_expressions!( container::OptimizationContainer, devices, line_names, time_steps, quad_cfg, - ::Type{F}, + ::DeviceModel{PSY.TwoTerminalVSCLine, HVDCTwoTerminalVSC}, ::NetworkModel{<:AbstractPowerModel}, -) where {F <: AbstractTwoTerminalVSCFormulation} - _uses_pq_sq_approx(F) || return - p_ft_var = - get_variable(container, FlowActivePowerFromToVariable, PSY.TwoTerminalVSCLine) - p_tf_var = - get_variable(container, FlowActivePowerToFromVariable, PSY.TwoTerminalVSCLine) - q_f_var = get_variable(container, HVDCReactivePowerFromVariable, PSY.TwoTerminalVSCLine) - q_t_var = get_variable(container, HVDCReactivePowerToVariable, PSY.TwoTerminalVSCLine) - - p_ft_bounds = [ - (min = PSY.get_active_power_limits_from(d).min, - max = PSY.get_active_power_limits_from(d).max) for d in devices - ] - p_tf_bounds = [ - (min = PSY.get_active_power_limits_to(d).min, - max = PSY.get_active_power_limits_to(d).max) for d in devices - ] - q_f_bounds = [ - (min = PSY.get_reactive_power_limits_from(d).min, - max = PSY.get_reactive_power_limits_from(d).max) for d in devices - ] - q_t_bounds = [ - (min = PSY.get_reactive_power_limits_to(d).min, - max = PSY.get_reactive_power_limits_to(d).max) for d in devices - ] - +) + p_ft = get_variable(container, FlowActivePowerFromToVariable, PSY.TwoTerminalVSCLine) + p_tf = get_variable(container, FlowActivePowerToFromVariable, PSY.TwoTerminalVSCLine) + q_f = get_variable(container, HVDCReactivePowerFromVariable, PSY.TwoTerminalVSCLine) + q_t = get_variable(container, HVDCReactivePowerToVariable, PSY.TwoTerminalVSCLine) + p_ft_bounds = PSY.get_active_power_limits_from.(devices) + p_tf_bounds = PSY.get_active_power_limits_to.(devices) + q_f_bounds = PSY.get_reactive_power_limits_from.(devices) + q_t_bounds = PSY.get_reactive_power_limits_to.(devices) IOM._add_quadratic_approx!( quad_cfg, container, PSY.TwoTerminalVSCLine, - line_names, time_steps, p_ft_var, p_ft_bounds, "p_ft_sq", + line_names, time_steps, p_ft, p_ft_bounds, "p_ft_sq", ) IOM._add_quadratic_approx!( quad_cfg, container, PSY.TwoTerminalVSCLine, - line_names, time_steps, p_tf_var, p_tf_bounds, "p_tf_sq", + line_names, time_steps, p_tf, p_tf_bounds, "p_tf_sq", ) IOM._add_quadratic_approx!( quad_cfg, container, PSY.TwoTerminalVSCLine, - line_names, time_steps, q_f_var, q_f_bounds, "q_f_sq", + line_names, time_steps, q_f, q_f_bounds, "q_f_sq", ) IOM._add_quadratic_approx!( quad_cfg, container, PSY.TwoTerminalVSCLine, - line_names, time_steps, q_t_var, q_t_bounds, "q_t_sq", + line_names, time_steps, q_t, q_t_bounds, "q_t_sq", ) return end -_maybe_add_pq_sq_approx!( - ::OptimizationContainer, - _devices, - _line_names, - _time_steps, - _quad_cfg, - ::Type{<:AbstractTwoTerminalVSCFormulation}, +# LP path: no disk constraint, so no p_sq/q_sq are needed. +_register_pq_sq_expressions!( + ::OptimizationContainer, _devices, _names, _times, _cfg, + ::DeviceModel{PSY.TwoTerminalVSCLine, HVDCTwoTerminalVSCLP}, + ::NetworkModel, +) = nothing + +# Active-power-only networks don't carry reactive variables at all. +_register_pq_sq_expressions!( + ::OptimizationContainer, _devices, _names, _times, _cfg, + ::DeviceModel{PSY.TwoTerminalVSCLine, HVDCTwoTerminalVSC}, ::NetworkModel{<:AbstractActivePowerModel}, ) = nothing @@ -1580,23 +1553,19 @@ _maybe_add_pq_sq_approx!( function add_constraints!( container::OptimizationContainer, ::Type{HVDCCableOhmsLawConstraint}, - devices::Union{ - Vector{PSY.TwoTerminalVSCLine}, - IS.FlattenIteratorWrapper{PSY.TwoTerminalVSCLine}, - }, - ::DeviceModel{PSY.TwoTerminalVSCLine, F}, + devices::Union{Vector{U}, IS.FlattenIteratorWrapper{U}}, + ::DeviceModel{U, F}, ::NetworkModel{<:AbstractPowerModel}, -) where {F <: AbstractTwoTerminalVSCFormulation} +) where {U <: PSY.TwoTerminalVSCLine, F <: AbstractTwoTerminalVSCFormulation} time_steps = get_time_steps(container) names = [PSY.get_name(d) for d in devices] jump_model = get_jump_model(container) - v_f = get_variable(container, HVDCFromDCVoltage, PSY.TwoTerminalVSCLine) - v_t = get_variable(container, HVDCToDCVoltage, PSY.TwoTerminalVSCLine) - i_var = get_variable(container, DCLineCurrentFlowVariable, PSY.TwoTerminalVSCLine) + v_f = get_variable(container, HVDCFromDCVoltage, U) + v_t = get_variable(container, HVDCToDCVoltage, U) + i_var = get_variable(container, DCLineCurrentFlowVariable, U) cons = add_constraints_container!( - container, HVDCCableOhmsLawConstraint, PSY.TwoTerminalVSCLine, - names, time_steps, + container, HVDCCableOhmsLawConstraint, U, names, time_steps, ) for d in devices @@ -1619,55 +1588,40 @@ end # Per-terminal converter power balance: # p_ft == v_f * I + (a_f * I^2 + b_f * |I| + c_f) # p_tf == -v_t * I + (a_t * I^2 + b_t * |I| + c_t) -# Sign convention: FlowActivePowerFromToVariable / ToFromVariable are positive -# when the corresponding AC bus is sourcing power into the converter (matches -# the existing add_to_expression! method's -1.0 multiplier). +# Active power enters `ActivePowerBalance` with a -1 multiplier (the AC bus +# sees the converter as a load drawing p_ft / p_tf), so positive values of +# `FlowActivePower*Variable` correspond to power flowing AC → DC at the +# respective terminal. function add_constraints!( container::OptimizationContainer, ::Type{HVDCVSCConverterPowerConstraint}, - devices::Union{ - Vector{PSY.TwoTerminalVSCLine}, - IS.FlattenIteratorWrapper{PSY.TwoTerminalVSCLine}, - }, - model::DeviceModel{PSY.TwoTerminalVSCLine, F}, + devices::Union{Vector{U}, IS.FlattenIteratorWrapper{U}}, + model::DeviceModel{U, F}, ::NetworkModel{<:AbstractPowerModel}, -) where {F <: AbstractTwoTerminalVSCFormulation} +) where {U <: PSY.TwoTerminalVSCLine, F <: AbstractTwoTerminalVSCFormulation} time_steps = get_time_steps(container) names = [PSY.get_name(d) for d in devices] jump_model = get_jump_model(container) - p_ft = get_variable(container, FlowActivePowerFromToVariable, PSY.TwoTerminalVSCLine) - p_tf = get_variable(container, FlowActivePowerToFromVariable, PSY.TwoTerminalVSCLine) - vi_expr = get_expression( - container, - IOM.BilinearProductExpression, - PSY.TwoTerminalVSCLine, - "vi_ft", - ) - vi_expr_to = get_expression( - container, - IOM.BilinearProductExpression, - PSY.TwoTerminalVSCLine, - "vi_tf", - ) - i_sq_expr = - get_expression(container, IOM.QuadraticExpression, PSY.TwoTerminalVSCLine, "i_sq") + p_ft = get_variable(container, FlowActivePowerFromToVariable, U) + p_tf = get_variable(container, FlowActivePowerToFromVariable, U) + vi_expr_ft = get_expression(container, IOM.BilinearProductExpression, U, "vi_ft") + vi_expr_tf = get_expression(container, IOM.BilinearProductExpression, U, "vi_tf") + i_sq_expr = get_expression(container, IOM.QuadraticExpression, U, "i_sq") use_linear_loss = get_attribute(model, "use_linear_loss") && !isempty(_devices_with_linear_loss(devices)) if use_linear_loss - i_pos_var = get_variable(container, PositiveCurrent, PSY.TwoTerminalVSCLine) - i_neg_var = get_variable(container, NegativeCurrent, PSY.TwoTerminalVSCLine) + i_pos_var = get_variable(container, PositiveCurrent, U) + i_neg_var = get_variable(container, NegativeCurrent, U) end cons_ft = add_constraints_container!( - container, HVDCVSCConverterPowerConstraint, PSY.TwoTerminalVSCLine, - names, time_steps; meta = "ft", + container, HVDCVSCConverterPowerConstraint, U, names, time_steps; meta = "ft", ) cons_tf = add_constraints_container!( - container, HVDCVSCConverterPowerConstraint, PSY.TwoTerminalVSCLine, - names, time_steps; meta = "tf", + container, HVDCVSCConverterPowerConstraint, U, names, time_steps; meta = "tf", ) for d in devices @@ -1693,67 +1647,45 @@ function add_constraints!( ) cons_ft[name, t] = JuMP.@constraint( jump_model, - p_ft[name, t] == vi_expr[name, t] + loss_ft, + p_ft[name, t] == vi_expr_ft[name, t] + loss_ft, ) cons_tf[name, t] = JuMP.@constraint( jump_model, - p_tf[name, t] == -vi_expr_to[name, t] + loss_tf, + p_tf[name, t] == -vi_expr_tf[name, t] + loss_tf, ) end end return end -# PQ capability for the exact-NLP (HVDCTwoTerminalVSC) and SOS2-PWL -# (HVDCTwoTerminalVSCMILP) paths. -# -# Constraint: p_sq + q_sq ≤ s² per terminal, where p_sq / q_sq are the -# expressions produced by `IOM._add_quadratic_approx!` (registered as -# `IOM.QuadraticExpression` with metas `p_*_sq` / `q_*_sq` by -# `_maybe_add_pq_sq_approx!`). -# -# • HVDCTwoTerminalVSC uses `NoQuadApproxConfig`, so `p_sq[name,t]` IS the -# exact `p_ft[name,t]^2` (JuMP.QuadExpr). The constraint reduces to -# `p² + q² ≤ s²` and the model stays a smooth NLP. -# • HVDCTwoTerminalVSCMILP uses `SolverSOS2QuadConfig(K)`, so `p_sq[name,t]` -# is an SOS2-based PWL approximation of `p²` (JuMP.AffExpr). Tightness -# grows with `K`; the approximation is not guaranteed to be a strict -# outer- or inner-bound of the disk — its position depends on PWL -# direction. +# PQ capability — exact disk for the NLP formulation. `p_*_sq` / `q_*_sq` are +# the `IOM.QuadraticExpression` handles registered by +# `_register_pq_sq_expressions!`. Under `NoQuadApproxConfig` (what +# `HVDCTwoTerminalVSC` uses) they are exact QuadExprs, so the constraint is +# the smooth `p² + q² ≤ s²` and the model stays an NLP. function add_constraints!( container::OptimizationContainer, ::Type{HVDCVSCApparentPowerLimitConstraint}, - devices::Union{ - Vector{PSY.TwoTerminalVSCLine}, - IS.FlattenIteratorWrapper{PSY.TwoTerminalVSCLine}, - }, - ::DeviceModel{PSY.TwoTerminalVSCLine, F}, + devices::Union{Vector{U}, IS.FlattenIteratorWrapper{U}}, + ::DeviceModel{U, HVDCTwoTerminalVSC}, ::NetworkModel{<:AbstractPowerModel}, -) where {F <: Union{HVDCTwoTerminalVSC, HVDCTwoTerminalVSCMILP}} +) where {U <: PSY.TwoTerminalVSCLine} time_steps = get_time_steps(container) names = [PSY.get_name(d) for d in devices] jump_model = get_jump_model(container) - p_ft_sq = get_expression( - container, IOM.QuadraticExpression, PSY.TwoTerminalVSCLine, "p_ft_sq", - ) - p_tf_sq = get_expression( - container, IOM.QuadraticExpression, PSY.TwoTerminalVSCLine, "p_tf_sq", - ) - q_f_sq = get_expression( - container, IOM.QuadraticExpression, PSY.TwoTerminalVSCLine, "q_f_sq", - ) - q_t_sq = get_expression( - container, IOM.QuadraticExpression, PSY.TwoTerminalVSCLine, "q_t_sq", - ) + p_ft_sq = get_expression(container, IOM.QuadraticExpression, U, "p_ft_sq") + p_tf_sq = get_expression(container, IOM.QuadraticExpression, U, "p_tf_sq") + q_f_sq = get_expression(container, IOM.QuadraticExpression, U, "q_f_sq") + q_t_sq = get_expression(container, IOM.QuadraticExpression, U, "q_t_sq") cons_f = add_constraints_container!( - container, HVDCVSCApparentPowerLimitConstraint, PSY.TwoTerminalVSCLine, - names, time_steps; meta = "from", + container, HVDCVSCApparentPowerLimitConstraint, U, names, time_steps; + meta = "from", ) cons_t = add_constraints_container!( - container, HVDCVSCApparentPowerLimitConstraint, PSY.TwoTerminalVSCLine, - names, time_steps; meta = "to", + container, HVDCVSCApparentPowerLimitConstraint, U, names, time_steps; + meta = "to", ) for d in devices @@ -1772,68 +1704,56 @@ function add_constraints!( return end -# PQ capability for the linear (HVDCTwoTerminalVSCLP) path. -# -# The disk `p² + q² ≤ rating²` is approximated by a regular octagon that -# *contains* the disk: the intersection of two squares circumscribing the -# disk, where each square is tangent to the disk at four points. +# PQ capability — linear outer-approximation for the LP formulation. # -# • Axis-aligned square: |p| ≤ rating, |q| ≤ rating -# (tangent at (±rating, 0) and (0, ±rating)) -# • 45°-rotated square: |p| + |q| ≤ rating·√2 -# (tangent at the four diagonal points) +# We always add the axis-aligned box |p|, |q| ≤ rating. When the +# device-model attribute `use_octagon` (default `true`) is on, we also add +# the four 45°-rotated diagonals |p| ± q ≤ rating·√2 ; their intersection +# with the box is a regular octagon circumscribing the disk +# p² + q² ≤ rating². # -# Their intersection is a regular octagon. Eight linear constraints per -# terminal per timestep (four diagonal + four axis-aligned). -# -# Tightness: the octagon is loose by at most ≈8.2% in area (octagon-to-disk -# area ratio = 8·tan(π/8)/π ≈ 1.082). It is *guaranteed* to be an outer -# approximation: any solution feasible to the real disk is feasible to the -# octagon, so the LP relaxation never blocks a physically realizable -# operating point. -# -# Tradeoff vs the unified PWL path used by `HVDCTwoTerminalVSC` / -# `HVDCTwoTerminalVSCMILP` (above): the SOS2 PWL on `p²` and `q²` separately -# can be made arbitrarily tight by increasing the breakpoint count, but the -# resulting feasible region is not guaranteed to be a strict outer- or -# inner-approximation of the disk. The octagon is the simplest always- -# conservative linear envelope and introduces no extra variables. +# Outer-approximation proof: for any (p, q) on the disk, p² ≤ p² + q² ≤ r² +# gives |p|, |q| ≤ r, and Cauchy–Schwarz gives (|p|+|q|)² ≤ 2(p²+q²) ≤ 2r² +# so |p|+|q| ≤ r√2. Both half-plane families contain the disk, and so does +# their intersection. The octagon is loose by at most ≈8.2% in area +# (octagon-to-disk area ratio 8·tan(π/8)/π ≈ 1.082). function add_constraints!( container::OptimizationContainer, ::Type{HVDCVSCApparentPowerLimitConstraint}, - devices::Union{ - Vector{PSY.TwoTerminalVSCLine}, - IS.FlattenIteratorWrapper{PSY.TwoTerminalVSCLine}, - }, - ::DeviceModel{PSY.TwoTerminalVSCLine, HVDCTwoTerminalVSCLP}, + devices::Union{Vector{U}, IS.FlattenIteratorWrapper{U}}, + model::DeviceModel{U, HVDCTwoTerminalVSCLP}, ::NetworkModel{<:AbstractPowerModel}, -) +) where {U <: PSY.TwoTerminalVSCLine} time_steps = get_time_steps(container) names = [PSY.get_name(d) for d in devices] jump_model = get_jump_model(container) - p_ft = get_variable(container, FlowActivePowerFromToVariable, PSY.TwoTerminalVSCLine) - p_tf = get_variable(container, FlowActivePowerToFromVariable, PSY.TwoTerminalVSCLine) - q_f = get_variable(container, HVDCReactivePowerFromVariable, PSY.TwoTerminalVSCLine) - q_t = get_variable(container, HVDCReactivePowerToVariable, PSY.TwoTerminalVSCLine) + p_ft = get_variable(container, FlowActivePowerFromToVariable, U) + p_tf = get_variable(container, FlowActivePowerToFromVariable, U) + q_f = get_variable(container, HVDCReactivePowerFromVariable, U) + q_t = get_variable(container, HVDCReactivePowerToVariable, U) - diag_tags = ("from_pp", "from_pn", "from_np", "from_nn", + use_octagon = get_attribute(model, "use_octagon") + side_tags = if use_octagon + ("from_p_ub", "from_p_lb", "from_q_ub", "from_q_lb", + "to_p_ub", "to_p_lb", "to_q_ub", "to_q_lb", + "from_pp", "from_pn", "from_np", "from_nn", "to_pp", "to_pn", "to_np", "to_nn") - axis_tags = ("from_p_ub", "from_p_lb", "from_q_ub", "from_q_lb", + else + ("from_p_ub", "from_p_lb", "from_q_ub", "from_q_lb", "to_p_ub", "to_p_lb", "to_q_ub", "to_q_lb") + end cons = Dict{String, Any}() - for tag in (diag_tags..., axis_tags...) + for tag in side_tags cons[tag] = add_constraints_container!( - container, HVDCVSCApparentPowerLimitConstraint, PSY.TwoTerminalVSCLine, + container, HVDCVSCApparentPowerLimitConstraint, U, names, time_steps; meta = tag, ) end - # Run the same octagon block twice (from/to) using a per-side spec, so - # the eight inner constraints aren't textually duplicated. side_specs = ( - (prefix = "from", p_var = p_ft, q_var = q_f, rating_getter = _vsc_rating_from), - (prefix = "to", p_var = p_tf, q_var = q_t, rating_getter = _vsc_rating_to), + (prefix = "from", p_var = p_ft, q_var = q_f, rating_getter = PSY.get_rating_from), + (prefix = "to", p_var = p_tf, q_var = q_t, rating_getter = PSY.get_rating_to), ) for d in devices name = PSY.get_name(d) @@ -1843,14 +1763,6 @@ function add_constraints!( prefix = spec.prefix p_var, q_var = spec.p_var, spec.q_var for t in time_steps - cons[prefix * "_pp"][name, t] = - JuMP.@constraint(jump_model, p_var[name, t] + q_var[name, t] <= diag) - cons[prefix * "_pn"][name, t] = - JuMP.@constraint(jump_model, p_var[name, t] - q_var[name, t] <= diag) - cons[prefix * "_np"][name, t] = - JuMP.@constraint(jump_model, -p_var[name, t] + q_var[name, t] <= diag) - cons[prefix * "_nn"][name, t] = - JuMP.@constraint(jump_model, -p_var[name, t] - q_var[name, t] <= diag) cons[prefix * "_p_ub"][name, t] = JuMP.@constraint(jump_model, p_var[name, t] <= rating) cons[prefix * "_p_lb"][name, t] = @@ -1859,6 +1771,28 @@ function add_constraints!( JuMP.@constraint(jump_model, q_var[name, t] <= rating) cons[prefix * "_q_lb"][name, t] = JuMP.@constraint(jump_model, -q_var[name, t] <= rating) + if use_octagon + cons[prefix * "_pp"][name, t] = + JuMP.@constraint( + jump_model, + p_var[name, t] + q_var[name, t] <= diag + ) + cons[prefix * "_pn"][name, t] = + JuMP.@constraint( + jump_model, + p_var[name, t] - q_var[name, t] <= diag + ) + cons[prefix * "_np"][name, t] = + JuMP.@constraint( + jump_model, + -p_var[name, t] + q_var[name, t] <= diag + ) + cons[prefix * "_nn"][name, t] = + JuMP.@constraint( + jump_model, + -p_var[name, t] - q_var[name, t] <= diag + ) + end end end end @@ -1884,21 +1818,17 @@ function get_default_attributes( return Dict{String, Any}("use_linear_loss" => false) end -# Default `use_linear_loss = true`: HVDCTwoTerminalVSCMILP already uses SOS2 PWL -# (MIP), so the extra binary direction variable from the linear-loss -# decomposition is free. -function get_default_attributes( - ::Type{PSY.TwoTerminalVSCLine}, - ::Type{HVDCTwoTerminalVSCMILP}, -) - return Dict{String, Any}("use_linear_loss" => true) -end - -# Default `use_linear_loss = true`: same MIP rationale as VSCMIP — the SOS2 PWL -# loss model already makes the formulation MIP-compatible. +# `use_linear_loss = true`: the SOS2 PWL loss model already pushes the +# formulation into MILP territory, so the extra binary direction variable from +# the linear-loss decomposition is free. +# `use_octagon = true`: adds the four diagonals |p| ± q ≤ rating·√2 on top of +# the axis-aligned box |p|, |q| ≤ rating. The intersection is a regular octagon +# circumscribing the disk p² + q² ≤ rating² and is a guaranteed outer +# approximation (loose by at most ≈8.2% in area). Setting it to false leaves +# only the box, which is cheaper but a looser linear envelope of the disk. function get_default_attributes( ::Type{PSY.TwoTerminalVSCLine}, ::Type{HVDCTwoTerminalVSCLP}, ) - return Dict{String, Any}("use_linear_loss" => true) + return Dict{String, Any}("use_linear_loss" => true, "use_octagon" => true) end diff --git a/test/test_device_hvdc.jl b/test/test_device_hvdc.jl index 2417240..f36f60a 100644 --- a/test/test_device_hvdc.jl +++ b/test/test_device_hvdc.jl @@ -97,7 +97,7 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED end -@testset "HVDC System with Losses Network (MIPQuadraticLossConverter)" begin +@testset "HVDC System with Losses Network (MILPQuadraticLossConverter)" begin sys = _generate_test_hvdc_sys() template = OperationsProblemTemplate() set_device_model!(template, ThermalStandard, ThermalDispatchNoMin) @@ -106,7 +106,7 @@ end set_device_model!(template, TModelHVDCLine, DCLossyLine) ipc_model = DeviceModel( InterconnectingConverter, - MIPQuadraticLossConverter, + MILPQuadraticLossConverter, ) set_device_model!(template, ipc_model) set_hvdc_network_model!(template, VoltageDispatchHVDCNetworkModel) @@ -147,7 +147,7 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED end -@testset "HVDC MIP vs NLP QuadraticLossConverter agreement" begin +@testset "HVDC MILP vs NLP QuadraticLossConverter agreement" begin function _build_and_solve(formulation, optimizer) sys = _generate_test_hvdc_sys() template = OperationsProblemTemplate() @@ -176,13 +176,13 @@ end return model end - mip_model = _build_and_solve(MIPQuadraticLossConverter, HiGHS_optimizer) + milp_model = _build_and_solve(MILPQuadraticLossConverter, HiGHS_optimizer) nlp_model = _build_and_solve(QuadraticLossConverter, ipopt_optimizer) - mip_obj = IOM.get_objective_value(OptimizationProblemOutputs(mip_model)) + milp_obj = IOM.get_objective_value(OptimizationProblemOutputs(milp_model)) nlp_obj = IOM.get_objective_value(OptimizationProblemOutputs(nlp_model)) - @test isapprox(mip_obj, nlp_obj; rtol = 0.05) + @test isapprox(milp_obj, nlp_obj; rtol = 0.05) end ############################################################################## @@ -252,13 +252,6 @@ function _build_vsc_model(formulation, network, optimizer; sys = _generate_test_ ) end -@testset "HVDC Two-Terminal VSC (MIP) on DCP" begin - model = _build_vsc_model(HVDCTwoTerminalVSCMILP, DCPPowerModel, HiGHS_optimizer) - @test build!(model; output_dir = mktempdir(; cleanup = true)) == - IOM.ModelBuildStatus.BUILT - @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED -end - @testset "HVDC Two-Terminal VSC (NLP) on DCP" begin model = _build_vsc_model(HVDCTwoTerminalVSC, DCPPowerModel, ipopt_optimizer) @test build!(model; output_dir = mktempdir(; cleanup = true)) == @@ -273,12 +266,6 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED end -# Omitted: HVDCTwoTerminalVSCMILP on ACPPowerModel cannot be solved by HiGHS -# because ACPPowerModel introduces trigonometric (cos/sin) constraints in the -# branch ohms law that require an NLP-capable solver. The MIP loss model -# (SOS2 PWL) is independent of the network's nonlinearity. To exercise this -# combination we'd need a MINLP solver (e.g. Gurobi). - @testset "HVDC Two-Terminal VSC (LP) on DCP" begin model = _build_vsc_model(HVDCTwoTerminalVSCLP, DCPPowerModel, HiGHS_optimizer) @test build!(model; output_dir = mktempdir(; cleanup = true)) == @@ -286,10 +273,15 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED end -# Omitted: HVDCTwoTerminalVSCLP on ACPPowerModel — same reason as VSCMIP/AC -# above (HiGHS cannot solve the ACP network's trigonometric branch constraints). +# Omitted: HVDCTwoTerminalVSCLP on ACPPowerModel — HiGHS cannot solve the +# ACP network's trigonometric (cos/sin) branch ohms-law constraints; we'd +# need an MINLP-capable solver (e.g. Gurobi) for that combination. -@testset "HVDC VSC MIP vs NLP objective agreement" begin +@testset "HVDC VSC LP vs NLP objective agreement" begin + # On a DC network the PQ disk constraint is inactive (no reactive + # variables exist), so the LP and NLP differ only by the i² loss model + # (SOS2 PWL vs exact). For a smooth convex loss curve the two should agree + # within a few percent. function _solve(formulation, optimizer) sys = _generate_test_vsc_sys() model = _build_vsc_model(formulation, DCPPowerModel, optimizer; sys = sys) @@ -298,9 +290,42 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED return IOM.get_optimization_container(model).optimizer_stats.objective_value end - mip_obj = _solve(HVDCTwoTerminalVSCMILP, HiGHS_optimizer) + lp_obj = _solve(HVDCTwoTerminalVSCLP, HiGHS_optimizer) nlp_obj = _solve(HVDCTwoTerminalVSC, ipopt_optimizer) - @test isapprox(mip_obj, nlp_obj; rtol = 0.05) + @test isapprox(lp_obj, nlp_obj; rtol = 0.05) +end + +@testset "HVDC VSC LP: octagon is at least as tight as box-only" begin + # With `use_octagon = true` the LP adds four diagonals on top of the box, + # shrinking the feasible region (or leaving it unchanged if the diagonals + # are non-binding). The min-cost objective should therefore be ≥ the + # box-only case. + function _solve_lp(use_octagon) + sys = _generate_test_vsc_sys() + template = OperationsProblemTemplate(NetworkModel(DCPPowerModel)) + set_device_model!(template, ThermalStandard, ThermalDispatchNoMin) + set_device_model!(template, RenewableDispatch, RenewableFullDispatch) + set_device_model!(template, PowerLoad, StaticPowerLoad) + set_device_model!(template, DeviceModel(Line, StaticBranch)) + set_device_model!( + template, + DeviceModel( + TwoTerminalVSCLine, + HVDCTwoTerminalVSCLP; + attributes = Dict("use_linear_loss" => true, "use_octagon" => use_octagon), + ), + ) + model = DecisionModel( + template, sys; store_variable_names = true, optimizer = HiGHS_optimizer, + ) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + return IOM.get_optimization_container(model).optimizer_stats.objective_value + end + box_only_obj = _solve_lp(false) + octagon_obj = _solve_lp(true) + @test octagon_obj >= box_only_obj - 1e-6 end @testset "HVDC VSC: higher cable resistance increases cost" begin @@ -308,7 +333,7 @@ end function _solve_with_g(g_value) sys = _generate_test_vsc_sys(; g = g_value) model = _build_vsc_model( - HVDCTwoTerminalVSCMILP, DCPPowerModel, HiGHS_optimizer; sys = sys, + HVDCTwoTerminalVSCLP, DCPPowerModel, HiGHS_optimizer; sys = sys, ) @test build!(model; output_dir = mktempdir(; cleanup = true)) == IOM.ModelBuildStatus.BUILT From 45eb412352ce1ae9341c301413fc8d72632c24ce Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Thu, 14 May 2026 14:06:45 -0400 Subject: [PATCH 6/7] formatting --- src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl index 2a7eaf2..f26df4a 100644 --- a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl +++ b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl @@ -1736,12 +1736,12 @@ function add_constraints!( use_octagon = get_attribute(model, "use_octagon") side_tags = if use_octagon ("from_p_ub", "from_p_lb", "from_q_ub", "from_q_lb", - "to_p_ub", "to_p_lb", "to_q_ub", "to_q_lb", - "from_pp", "from_pn", "from_np", "from_nn", - "to_pp", "to_pn", "to_np", "to_nn") + "to_p_ub", "to_p_lb", "to_q_ub", "to_q_lb", + "from_pp", "from_pn", "from_np", "from_nn", + "to_pp", "to_pn", "to_np", "to_nn") else ("from_p_ub", "from_p_lb", "from_q_ub", "from_q_lb", - "to_p_ub", "to_p_lb", "to_q_ub", "to_q_lb") + "to_p_ub", "to_p_lb", "to_q_ub", "to_q_lb") end cons = Dict{String, Any}() for tag in side_tags From d176fb011c62da2d9c667bf039732bdb24fe1018 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Thu, 14 May 2026 14:48:32 -0400 Subject: [PATCH 7/7] Address PR #119 Review 4 on HVDCTwoTerminalVSC --- src/PowerOperationsModels.jl | 2 +- .../branch_constructor.jl | 15 +++++++++------ src/common_models/quadratic_converter_loss.jl | 2 +- src/core/constraints.jl | 2 +- src/core/formulations.jl | 2 +- src/core/variables.jl | 4 ++-- src/mt_hvdc_models/hvdcsystems_constructor.jl | 9 ++++++--- .../TwoTerminalDC_branches.jl | 19 ++++++++++--------- test/test_device_hvdc.jl | 8 ++++---- 9 files changed, 35 insertions(+), 28 deletions(-) diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl index 9f9195e..6da666b 100644 --- a/src/PowerOperationsModels.jl +++ b/src/PowerOperationsModels.jl @@ -764,7 +764,7 @@ export HVDCTwoTerminalLossless export HVDCTwoTerminalDispatch export HVDCTwoTerminalPiecewiseLoss export HVDCTwoTerminalLCC -export HVDCTwoTerminalVSC +export HVDCTwoTerminalVSCNLP export HVDCTwoTerminalVSCLP # Converter Formulations diff --git a/src/ac_transmission_models/branch_constructor.jl b/src/ac_transmission_models/branch_constructor.jl index 7ec35b1..fcc96db 100644 --- a/src/ac_transmission_models/branch_constructor.jl +++ b/src/ac_transmission_models/branch_constructor.jl @@ -1677,10 +1677,10 @@ end # Quadratic / bilinear approximation traits — same scheme used by the MT # converter formulations. -_quad_config(::Type{HVDCTwoTerminalVSC}) = IOM.NoQuadApproxConfig() +_quad_config(::Type{HVDCTwoTerminalVSCNLP}) = IOM.NoQuadApproxConfig() _quad_config(::Type{HVDCTwoTerminalVSCLP}) = IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH) -_bilinear_config(::Type{HVDCTwoTerminalVSC}) = IOM.NoBilinearApproxConfig() +_bilinear_config(::Type{HVDCTwoTerminalVSCNLP}) = IOM.NoBilinearApproxConfig() _bilinear_config(::Type{HVDCTwoTerminalVSCLP}) = IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)) @@ -1704,9 +1704,12 @@ function construct_device!( (HVDCReactivePowerFromVariable, HVDCReactivePowerToVariable), ) - # use_linear_loss = true adds a binary direction variable (CurrentDirection), - # making the model MINLP rather than smooth NLP. Caller is responsible for - # supplying a MINLP-capable solver in that case. + # use_linear_loss = true adds a binary direction variable (CurrentDirection). + # Both HVDCTwoTerminalVSCNLP and HVDCTwoTerminalVSCLP dispatch through here: + # on the NLP formulation this pushes the model from a smooth NLP to MINLP + # (caller must supply a MINLP-capable solver); on the LP formulation the + # SOS2 PWL approximations already make the model MILP, so the extra binary + # is free. if get_attribute(device_model, "use_linear_loss") _add_abs_value_decomposition_variables!( container, devices, device_model, network_model, @@ -1776,7 +1779,7 @@ function construct_device!( ) _register_pq_sq_expressions!( - container, devices, line_names, time_steps, quad_cfg, device_model, + container, devices, line_names, time_steps, device_model, network_model, ) diff --git a/src/common_models/quadratic_converter_loss.jl b/src/common_models/quadratic_converter_loss.jl index 9d301b1..9ad09d7 100644 --- a/src/common_models/quadratic_converter_loss.jl +++ b/src/common_models/quadratic_converter_loss.jl @@ -2,7 +2,7 @@ # loss(I) = a * I^2 + b * |I| + c # Used by multi-terminal InterconnectingConverter formulations # (MILPQuadraticLossConverter, QuadraticLossConverter) and two-terminal -# HVDCTwoTerminalVSC formulations. +# HVDCTwoTerminalVSCNLP formulations. ######################################### ######## Loss-curve introspection ####### diff --git a/src/core/constraints.jl b/src/core/constraints.jl index a6a405b..1eb459d 100644 --- a/src/core/constraints.jl +++ b/src/core/constraints.jl @@ -445,7 +445,7 @@ Apparent-power limit at each terminal of a two-terminal VSC HVDC, added only on AC networks. Enforces ``|S_k| \\le S_k^{\\max}`` for ``k \\in \\{f, t\\}`` via one of two formulation-specific shapes: -- `HVDCTwoTerminalVSC` (NLP): exact disk ``p_k^2 + q_k^2 \\le (S_k^{\\max})^2``. +- `HVDCTwoTerminalVSCNLP` (NLP): exact disk ``p_k^2 + q_k^2 \\le (S_k^{\\max})^2``. - `HVDCTwoTerminalVSCLP` (LP): linear outer-approximation. The axis-aligned box ``|p_k|, |q_k| \\le S_k^{\\max}`` is added unconditionally; the four diagonal half-planes ``|p_k| \\pm q_k \\le S_k^{\\max}\\sqrt{2}`` are added when the diff --git a/src/core/formulations.jl b/src/core/formulations.jl index 37b20e1..32201ae 100644 --- a/src/core/formulations.jl +++ b/src/core/formulations.jl @@ -197,7 +197,7 @@ abstract type AbstractTwoTerminalVSCFormulation <: AbstractTwoTerminalDCLineForm Two-terminal VSC formulation that keeps the bilinear ``v \\cdot I`` and quadratic ``I^2`` terms exact. Requires an NLP-capable solver (e.g. Ipopt). """ -struct HVDCTwoTerminalVSC <: AbstractTwoTerminalVSCFormulation end +struct HVDCTwoTerminalVSCNLP <: AbstractTwoTerminalVSCFormulation end """ Two-terminal VSC formulation that uses SOS2 piecewise-linear surrogates for the diff --git a/src/core/variables.jl b/src/core/variables.jl index 44ab5bb..bb3ce26 100644 --- a/src/core/variables.jl +++ b/src/core/variables.jl @@ -389,14 +389,14 @@ struct HVDCPiecewiseBinaryLossVariable <: SparseVariableType end """ DC-side voltage at the from-terminal of a two-terminal HVDC link. -Used by `HVDCTwoTerminalVSC` formulations. +Used by `HVDCTwoTerminalVSCNLP` formulations. Docs abbreviation: ``v_f^{dc}`` """ struct HVDCFromDCVoltage <: VariableType end """ DC-side voltage at the to-terminal of a two-terminal HVDC link. -Used by `HVDCTwoTerminalVSC` formulations. +Used by `HVDCTwoTerminalVSCNLP` formulations. Docs abbreviation: ``v_t^{dc}`` """ struct HVDCToDCVoltage <: VariableType end diff --git a/src/mt_hvdc_models/hvdcsystems_constructor.jl b/src/mt_hvdc_models/hvdcsystems_constructor.jl index dd15cfa..e8d00a3 100644 --- a/src/mt_hvdc_models/hvdcsystems_constructor.jl +++ b/src/mt_hvdc_models/hvdcsystems_constructor.jl @@ -55,9 +55,12 @@ function construct_device!( add_variables!(container, ActivePowerVariable, devices, T) add_variables!(container, ConverterCurrent, devices, T) - # use_linear_loss = true on QuadraticLossConverter adds a binary direction - # variable (CurrentDirection), making the model MINLP rather than smooth - # NLP. Caller is responsible for supplying a MINLP-capable solver. + # use_linear_loss = true adds a binary direction variable (CurrentDirection). + # Both MILPQuadraticLossConverter and QuadraticLossConverter dispatch through + # here: on the NLP QuadraticLossConverter this pushes the model from a smooth + # NLP to MINLP (caller must supply a MINLP-capable solver); on + # MILPQuadraticLossConverter the SOS2 PWL approximations already make the + # model MILP, so the extra binary is free. if get_attribute(model, "use_linear_loss") _add_abs_value_decomposition_variables!( container, devices, model, network_model, diff --git a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl index f26df4a..d320236 100644 --- a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl +++ b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl @@ -1502,10 +1502,11 @@ function _register_pq_sq_expressions!( devices, line_names, time_steps, - quad_cfg, - ::DeviceModel{PSY.TwoTerminalVSCLine, HVDCTwoTerminalVSC}, + ::DeviceModel{PSY.TwoTerminalVSCLine, HVDCTwoTerminalVSCNLP}, ::NetworkModel{<:AbstractPowerModel}, ) + # This dispatch is the NLP path, so the quad config is fixed. + quad_cfg = IOM.NoQuadApproxConfig() p_ft = get_variable(container, FlowActivePowerFromToVariable, PSY.TwoTerminalVSCLine) p_tf = get_variable(container, FlowActivePowerToFromVariable, PSY.TwoTerminalVSCLine) q_f = get_variable(container, HVDCReactivePowerFromVariable, PSY.TwoTerminalVSCLine) @@ -1535,15 +1536,15 @@ end # LP path: no disk constraint, so no p_sq/q_sq are needed. _register_pq_sq_expressions!( - ::OptimizationContainer, _devices, _names, _times, _cfg, + ::OptimizationContainer, _devices, _names, _times, ::DeviceModel{PSY.TwoTerminalVSCLine, HVDCTwoTerminalVSCLP}, ::NetworkModel, ) = nothing # Active-power-only networks don't carry reactive variables at all. _register_pq_sq_expressions!( - ::OptimizationContainer, _devices, _names, _times, _cfg, - ::DeviceModel{PSY.TwoTerminalVSCLine, HVDCTwoTerminalVSC}, + ::OptimizationContainer, _devices, _names, _times, + ::DeviceModel{PSY.TwoTerminalVSCLine, HVDCTwoTerminalVSCNLP}, ::NetworkModel{<:AbstractActivePowerModel}, ) = nothing @@ -1661,13 +1662,13 @@ end # PQ capability — exact disk for the NLP formulation. `p_*_sq` / `q_*_sq` are # the `IOM.QuadraticExpression` handles registered by # `_register_pq_sq_expressions!`. Under `NoQuadApproxConfig` (what -# `HVDCTwoTerminalVSC` uses) they are exact QuadExprs, so the constraint is +# `HVDCTwoTerminalVSCNLP` uses) they are exact QuadExprs, so the constraint is # the smooth `p² + q² ≤ s²` and the model stays an NLP. function add_constraints!( container::OptimizationContainer, ::Type{HVDCVSCApparentPowerLimitConstraint}, devices::Union{Vector{U}, IS.FlattenIteratorWrapper{U}}, - ::DeviceModel{U, HVDCTwoTerminalVSC}, + ::DeviceModel{U, HVDCTwoTerminalVSCNLP}, ::NetworkModel{<:AbstractPowerModel}, ) where {U <: PSY.TwoTerminalVSCLine} time_steps = get_time_steps(container) @@ -1808,12 +1809,12 @@ function get_default_time_series_names( return Dict{Type{<:TimeSeriesParameter}, String}() end -# Default `use_linear_loss = false`: HVDCTwoTerminalVSC is a smooth NLP solvable +# Default `use_linear_loss = false`: HVDCTwoTerminalVSCNLP is a smooth NLP solvable # by Ipopt; enabling the linear-loss term introduces a binary direction variable # that pushes the model to MINLP, which pure NLP solvers cannot handle. function get_default_attributes( ::Type{PSY.TwoTerminalVSCLine}, - ::Type{HVDCTwoTerminalVSC}, + ::Type{HVDCTwoTerminalVSCNLP}, ) return Dict{String, Any}("use_linear_loss" => false) end diff --git a/test/test_device_hvdc.jl b/test/test_device_hvdc.jl index f36f60a..2632e26 100644 --- a/test/test_device_hvdc.jl +++ b/test/test_device_hvdc.jl @@ -253,14 +253,14 @@ function _build_vsc_model(formulation, network, optimizer; sys = _generate_test_ end @testset "HVDC Two-Terminal VSC (NLP) on DCP" begin - model = _build_vsc_model(HVDCTwoTerminalVSC, DCPPowerModel, ipopt_optimizer) + model = _build_vsc_model(HVDCTwoTerminalVSCNLP, DCPPowerModel, ipopt_optimizer) @test build!(model; output_dir = mktempdir(; cleanup = true)) == IOM.ModelBuildStatus.BUILT @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED end @testset "HVDC Two-Terminal VSC on AC (NLP)" begin - model = _build_vsc_model(HVDCTwoTerminalVSC, ACPPowerModel, ipopt_optimizer) + model = _build_vsc_model(HVDCTwoTerminalVSCNLP, ACPPowerModel, ipopt_optimizer) @test build!(model; output_dir = mktempdir(; cleanup = true)) == IOM.ModelBuildStatus.BUILT @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED @@ -291,7 +291,7 @@ end return IOM.get_optimization_container(model).optimizer_stats.objective_value end lp_obj = _solve(HVDCTwoTerminalVSCLP, HiGHS_optimizer) - nlp_obj = _solve(HVDCTwoTerminalVSC, ipopt_optimizer) + nlp_obj = _solve(HVDCTwoTerminalVSCNLP, ipopt_optimizer) @test isapprox(lp_obj, nlp_obj; rtol = 0.05) end @@ -349,7 +349,7 @@ end function _solve_with_rating(s) sys = _generate_test_vsc_sys(; rating_from = s, rating_to = s) model = _build_vsc_model( - HVDCTwoTerminalVSC, ACPPowerModel, ipopt_optimizer; sys = sys, + HVDCTwoTerminalVSCNLP, ACPPowerModel, ipopt_optimizer; sys = sys, ) @test build!(model; output_dir = mktempdir(; cleanup = true)) == IOM.ModelBuildStatus.BUILT