From 35145a814a6c577c234a448c69a01891e5fb5442 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Wed, 29 Apr 2026 15:50:17 -0600 Subject: [PATCH 1/9] claude hvdc --- Project.toml | 12 +- src/PowerOperationsModels.jl | 16 +- src/core/constraints.jl | 73 --- src/core/formulations.jl | 22 +- src/core/variables.jl | 79 ---- src/mt_hvdc_models/HVDCsystems.jl | 437 ++++-------------- src/mt_hvdc_models/hvdcsystems_constructor.jl | 240 ++++------ test/Project.toml | 8 +- test/test_device_hvdc.jl | 35 +- 9 files changed, 245 insertions(+), 677 deletions(-) diff --git a/Project.toml b/Project.toml index e749a78..0a828b7 100644 --- a/Project.toml +++ b/Project.toml @@ -23,18 +23,18 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [weakdeps] 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"} + [extensions] PowerFlowsExt = "PowerFlows" -[sources] -InfrastructureSystems = {url = "https://github.com/NREL-Sienna/InfrastructureSystems.jl", rev = "IS4"} -PowerSystems = {url = "https://github.com/NREL-Sienna/PowerSystems.jl", rev = "psy6"} -InfrastructureOptimizationModels = {url = "https://github.com/NREL-Sienna/InfrastructureOptimizationModels.jl", rev = "lk/pom-test-fixes"} - [compat] Dates = "1" DocStringExtensions = "~0.8, ~0.9" -InfrastructureOptimizationModels = "0.1" +InfrastructureOptimizationModels = "0.1.0" InfrastructureSystems = "3" InteractiveUtils = "1.11.0" JuMP = "^1.28" diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl index 27411c1..738cc21 100644 --- a/src/PowerOperationsModels.jl +++ b/src/PowerOperationsModels.jl @@ -505,24 +505,12 @@ export PostContingencyActivePowerReserveDeploymentVariable # HVDC Variables export DCVoltage export DCLineCurrent -export ConverterPowerDirection export ConverterCurrent -export SquaredConverterCurrent -export InterpolationSquaredCurrentVariable -export InterpolationBinarySquaredCurrentVariable export ConverterPositiveCurrent export ConverterNegativeCurrent -export SquaredDCVoltage -export InterpolationSquaredVoltageVariable -export InterpolationBinarySquaredVoltageVariable -export AuxBilinearConverterVariable -export AuxBilinearSquaredConverterVariable -export InterpolationSquaredBilinearVariable -export InterpolationBinarySquaredBilinearVariable +export ConverterCurrentDirection export HVDCFlowDirectionVariable export HVDCLosses -export ConverterDCPower -export ConverterCurrentDirection # Load Variables export ShiftUpActivePowerVariable @@ -774,6 +762,8 @@ export HVDCTwoTerminalLCC # Converter Formulations export LosslessConverter export LinearLossConverter +export AbstractQuadraticLossConverter +export Bin2QuadraticLossConverter export QuadraticLossConverter # DC Line Formulations diff --git a/src/core/constraints.jl b/src/core/constraints.jl index 2dc7a17..634e5c2 100644 --- a/src/core/constraints.jl +++ b/src/core/constraints.jl @@ -577,19 +577,6 @@ struct DCLineCurrentConstraint <: ConstraintType end struct NodalBalanceCurrentConstraint <: ConstraintType end -""" -Struct to create the constraints that compute the converter DC power based on current and voltage. - -The specified constraints are formulated as: -```math -\\begin{align*} -& p_c = 0.5 * (γ^sq - v^sq - i^sq), \\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ -& γ_c = v_c + i_c, \\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ -\\end{align*} -``` -""" -struct ConverterPowerCalculationConstraint <: ConstraintType end - """ Struct to create the constraints that decide the balance of AC and DC power of the converter. @@ -603,66 +590,6 @@ The specified constraints are formulated as: """ struct ConverterLossConstraint <: ConstraintType end -""" -Struct to create the McCormick envelopes constraints that decide the bounds on the DC active power. - -The specified constraints are formulated as: -```math -\\begin{align*} -& p_c >= V^{min} i_c + v_c I^{min} - I^{min}V^{min}, \\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ -& p_c >= V^{max} i_c + v_c I^{max} - I^{max}V^{max}, \\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ -& p_c <= V^{max} i_c + v_c I^{min} - I^{min}V^{max}, \\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ -& p_c <= V^{min} i_c + v_c I^{max} - I^{max}V^{min}, \\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ -\\end{align*} -``` -""" -struct ConverterMcCormickEnvelopes <: ConstraintType end - -""" -Struct to create the Quadratic PWL interpolation constraints that decide square value of the voltage. -In this case x = voltage and y = squared_voltage. -The specified constraints are formulated as: -```math -\\begin{align*} -& x = x_0 + \\sum_{k=1}^K (x_{k} - x_{k-1}) \\delta_k, \\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ -& y = y_0 + \\sum_{k=1}^K (x_{k} - x_{k-1}) \\delta_k, \\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ -& z_k \\le \\delta_k, \\quad \\forall t \\in \\{1,\\dots, T\\}, \\forall k \\in \\{1,\\dots, K-1\\} \\\\ -& z_k \\ge \\delta_{k+1}, \\quad \\forall t \\in \\{1,\\dots, T\\}, \\forall k \\in \\{1,\\dots, K-1\\} \\\\ -\\end{align*} -``` -""" -struct InterpolationVoltageConstraints <: ConstraintType end - -""" -Struct to create the Quadratic PWL interpolation constraints that decide square value of the current. -In this case x = current and y = squared_current. -The specified constraints are formulated as: -```math -\\begin{align*} -& x = x_0 + \\sum_{k=1}^K (x_{k} - x_{k-1}) \\delta_k, \\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ -& y = y_0 + \\sum_{k=1}^K (x_{k} - x_{k-1}) \\delta_k, \\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ -& z_k \\le \\delta_k, \\quad \\forall t \\in \\{1,\\dots, T\\}, \\forall k \\in \\{1,\\dots, K-1\\} \\\\ -& z_k \\ge \\delta_{k+1}, \\quad \\forall t \\in \\{1,\\dots, T\\}, \\forall k \\in \\{1,\\dots, K-1\\} \\\\ -\\end{align*} -``` -""" -struct InterpolationCurrentConstraints <: ConstraintType end - -""" -Struct to create the Quadratic PWL interpolation constraints that decide square value of the bilinear variable γ. -In this case x = γ and y = squared_γ. -The specified constraints are formulated as: -```math -\\begin{align*} -& x = x_0 + \\sum_{k=1}^K (x_{k} - x_{k-1}) \\delta_k, \\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ -& y = y_0 + \\sum_{k=1}^K (x_{k} - x_{k-1}) \\delta_k, \\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ -& z_k \\le \\delta_k, \\quad \\forall t \\in \\{1,\\dots, T\\}, \\forall k \\in \\{1,\\dots, K-1\\} \\\\ -& z_k \\ge \\delta_{k+1}, \\quad \\forall t \\in \\{1,\\dots, T\\}, \\forall k \\in \\{1,\\dots, K-1\\} \\\\ -\\end{align*} -``` -""" -struct InterpolationBilinearConstraints <: ConstraintType end - """ Struct to create the constraints that set the absolute value for the current to use in losses through a lossy Interconnecting Power Converter. The specified constraint is formulated as: diff --git a/src/core/formulations.jl b/src/core/formulations.jl index 1159327..c37aec3 100644 --- a/src/core/formulations.jl +++ b/src/core/formulations.jl @@ -201,9 +201,27 @@ Linear Loss InterconnectingConverter Model struct LinearLossConverter <: AbstractConverterFormulation end """ -Quadratic Loss InterconnectingConverter Model +Abstract supertype for InterconnectingConverter formulations with quadratic losses. """ -struct QuadraticLossConverter <: AbstractConverterFormulation end +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². +""" +struct Bin2QuadraticLossConverter <: AbstractQuadraticLossConverter end + +""" +Quadratic Loss InterconnectingConverter using exact bilinear (v·i) and quadratic (i²) +products. Requires an NLP-capable solver (e.g., Ipopt). + +Attributes (set on `DeviceModel`): +- `use_linear_loss`: Include the linear loss term `b·(i⁺ + i⁻)` (default false). The + linear term requires an absolute-value reformulation that introduces a binary + variable, so this is only solvable by MINLP solvers (e.g., Gurobi); pure NLP + solvers like Ipopt cannot solve the model with this option enabled. +""" +struct QuadraticLossConverter <: AbstractQuadraticLossConverter end ############################## HVDC Lines Formulations ################################## abstract type AbstractDCLineFormulation <: AbstractBranchFormulation end diff --git a/src/core/variables.jl b/src/core/variables.jl index 660ef5e..8a41d0d 100644 --- a/src/core/variables.jl +++ b/src/core/variables.jl @@ -150,30 +150,12 @@ Docs abbreviation: ``i_l^{dc}`` """ struct DCLineCurrent <: VariableType end -""" -Struct to dispatch the creation of Squared Voltage Variables for DC formulations -Docs abbreviation: ``v^{sq,dc}`` -""" -struct SquaredDCVoltage <: VariableType end - """ Struct to dispatch the creation of DC Converter Current Variables for DC formulations Docs abbreviation: ``i_c^{dc}`` """ struct ConverterCurrent <: VariableType end -""" -Struct to dispatch the creation of DC Converter Power Variables for DC formulations -Docs abbreviation: ``p_c^{dc}`` -""" -struct ConverterDCPower <: VariableType end - -""" -Struct to dispatch the creation of Squared DC Converter Current Variables for DC formulations -Docs abbreviation: ``i_c^{sq,dc}`` -""" -struct SquaredConverterCurrent <: VariableType end - """ Struct to dispatch the creation of DC Converter Positive Term Current Variables for DC formulations Docs abbreviation: ``i_c^{+,dc}`` @@ -192,67 +174,6 @@ Docs abbreviation: `\\nu_c`` """ struct ConverterCurrentDirection <: VariableType end -""" -Struct to dispatch the creation of Binary Variable for Converter Power Direction -Docs abbreviation: ``\\kappa_c^{dc}`` -""" -struct ConverterPowerDirection <: VariableType end - -""" -Struct to dispatch the creation of Auxiliary Variable for Converter Bilinear term: v * i -Docs abbreviation: ``\\gamma_c^{dc}`` -""" -struct AuxBilinearConverterVariable <: VariableType end - -""" -Struct to dispatch the creation of Auxiliary Variable for Squared Converter Bilinear term: v * i - -Docs abbreviation: ``\\gamma_c^{sq,dc}`` -""" -struct AuxBilinearSquaredConverterVariable <: VariableType end - -""" -Struct to dispatch the creation of Continuous Interpolation Variable for Squared Converter Voltage - -Docs abbreviation: ``\\delta_c^{v}`` -""" -struct InterpolationSquaredVoltageVariable <: InterpolationVariableType end - -""" -Struct to dispatch the creation of Binary Interpolation Variable for Squared Converter Voltage - -Docs abbreviation: ``z_c^{v}`` -""" -struct InterpolationBinarySquaredVoltageVariable <: BinaryInterpolationVariableType end - -""" -Struct to dispatch the creation of Continuous Interpolation Variable for Squared Converter Current - -Docs abbreviation: ``\\delta_c^{i}`` -""" -struct InterpolationSquaredCurrentVariable <: InterpolationVariableType end - -""" -Struct to dispatch the creation of Binary Interpolation Variable for Squared Converter Current - -Docs abbreviation: ``z_c^{i}`` -""" -struct InterpolationBinarySquaredCurrentVariable <: BinaryInterpolationVariableType end - -""" -Struct to dispatch the creation of Continuous Interpolation Variable for Squared Converter AuxVar - -Docs abbreviation: ``\\delta_c^{\\gamma}`` -""" -struct InterpolationSquaredBilinearVariable <: InterpolationVariableType end - -""" -Struct to dispatch the creation of Binary Interpolation Variable for Squared Converter AuxVar - -Docs abbreviation: ``z_c^{\\gamma}`` -""" -struct InterpolationBinarySquaredBilinearVariable <: BinaryInterpolationVariableType end - ######################################################### ######################################################### diff --git a/src/mt_hvdc_models/HVDCsystems.jl b/src/mt_hvdc_models/HVDCsystems.jl index cfa40e1..f374d29 100644 --- a/src/mt_hvdc_models/HVDCsystems.jl +++ b/src/mt_hvdc_models/HVDCsystems.jl @@ -119,63 +119,40 @@ end ############################################ ## Binaries ### -get_variable_binary(::Type{ConverterDCPower}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = false -get_variable_binary(::Type{ConverterPowerDirection}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = true 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{SquaredConverterCurrent}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = false -get_variable_binary(::Type{SquaredDCVoltage}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = false -get_variable_binary(::Type{AuxBilinearConverterVariable}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = false -get_variable_binary(::Type{AuxBilinearSquaredConverterVariable}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = false -function get_variable_binary( - ::Type{W}, - ::Type{PSY.InterconnectingConverter}, - ::Type{<:AbstractConverterFormulation} -) where W <: InterpolationVariableType - return false -end -function get_variable_binary( - ::Type{W}, - ::Type{PSY.InterconnectingConverter}, - ::Type{<:AbstractConverterFormulation} -) where W <: BinaryInterpolationVariableType - return true -end - ### 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{ConverterDCPower}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = PSY.get_active_power_limits(d).min get_variable_lower_bound(::Type{ConverterCurrent}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = -PSY.get_max_dc_current(d) -get_variable_lower_bound(::Type{SquaredConverterCurrent}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = 0.0 -get_variable_lower_bound(::Type{SquaredDCVoltage}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = PSY.get_voltage_limits(d.dc_bus).min^2 -get_variable_lower_bound(::Type{<:InterpolationVariableType}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = 0.0 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 ### Upper Bounds ### -get_variable_upper_bound(::Type{ConverterDCPower}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = PSY.get_active_power_limits(d).max get_variable_upper_bound(::Type{ConverterCurrent}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = PSY.get_max_dc_current(d) -get_variable_upper_bound(::Type{SquaredConverterCurrent}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = PSY.get_max_dc_current(d)^2 -get_variable_upper_bound(::Type{SquaredDCVoltage}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = PSY.get_voltage_limits(d.dc_bus).max^2 -get_variable_upper_bound(::Type{<:InterpolationVariableType}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = 1.0 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) +function get_default_attributes( + ::Type{PSY.InterconnectingConverter}, + ::Type{Bin2QuadraticLossConverter}, +) + return Dict{String, Any}( + "use_linear_loss" => true + ) +end + function get_default_attributes( ::Type{PSY.InterconnectingConverter}, ::Type{QuadraticLossConverter}, ) return Dict{String, Any}( - "voltage_segments" => 3, - "current_segments" => 6, - "bilinear_segments" => 10, - "use_linear_loss" => true, + "use_linear_loss" => false, ) end @@ -422,7 +399,7 @@ function add_to_expression!( T <: ActivePowerBalance, U <: ActivePowerVariable, V <: PSY.InterconnectingConverter, - W <: QuadraticLossConverter, + W <: AbstractQuadraticLossConverter, } variable = get_variable(container, U, V) sys_expr = get_expression(container, T, PSY.System) @@ -452,7 +429,7 @@ function add_to_expression!( T <: DCCurrentBalance, U <: ConverterCurrent, V <: PSY.InterconnectingConverter, - W <: QuadraticLossConverter, + W <: AbstractQuadraticLossConverter, } variable = get_variable(container, U, V) expression_dc = get_expression(container, T, PSY.DCBus) @@ -564,174 +541,85 @@ function add_constraints!( end ############## Converters ################## -function add_constraints!( - container::OptimizationContainer, - ::Type{ConverterPowerCalculationConstraint}, - devices::IS.FlattenIteratorWrapper{U}, - model::DeviceModel{U, V}, - network_model::NetworkModel{X}, -) where { - U <: PSY.InterconnectingConverter, - V <: QuadraticLossConverter, - X <: AbstractActivePowerModel, -} - time_steps = get_time_steps(container) - varcurrent = get_variable(container, ConverterCurrent, U) - var_dcvoltage = get_variable(container, DCVoltage, PSY.DCBus) - var_sq_current = get_variable(container, SquaredConverterCurrent, U) - var_sq_voltage = get_variable(container, SquaredDCVoltage, U) - var_bilinear = get_variable(container, AuxBilinearConverterVariable, U) - var_sq_bilinear = get_variable(container, AuxBilinearSquaredConverterVariable, U) - var_dc_power = get_variable(container, ConverterDCPower, U) - ipc_names = axes(varcurrent, 1) - constraint = - add_constraints_container!(container, ConverterPowerCalculationConstraint, - U, - ipc_names, - time_steps, - ) - constraint_aux = - add_constraints_container!(container, ConverterPowerCalculationConstraint, - U, - ipc_names, - time_steps; - meta = "aux", - ) - for device in devices - name = PSY.get_name(device) - dc_bus_name = PSY.get_name(PSY.get_dc_bus(device)) +# Builds, per converter, an AffExpr container of the DC-bus voltage variable indexed +# by converter name. IOM bilinear/quadratic approximations consume an x_var indexed by +# component name; DCVoltage is stored per DC bus, so we wrap it for compatibility. +function _voltage_expr_per_converter( + container::OptimizationContainer, + devices, + ipc_names::Vector{String}, + time_steps, +) + v_var = get_variable(container, DCVoltage, PSY.DCBus) + v_expr = JuMP.Containers.DenseAxisArray{JuMP.AffExpr}(undef, ipc_names, time_steps) + for d in devices + name = PSY.get_name(d) + dc_bus_name = PSY.get_name(PSY.get_dc_bus(d)) for t in time_steps - # p_dc = v_dc * i_dc = 0.5 * (bilinear - v_dc^2 - i_dc^2) - constraint[name, t] = JuMP.@constraint( - get_jump_model(container), - var_dc_power[name, t] == - 0.5 * ( - var_sq_bilinear[name, t] - var_sq_voltage[name, t] - - var_sq_current[name, t] - ) - ) - constraint_aux[name, t] = JuMP.@constraint( - get_jump_model(container), - var_bilinear[name, t] == - var_dcvoltage[dc_bus_name, t] + varcurrent[name, t] - ) + ae = JuMP.AffExpr(0.0) + add_proportional_to_jump_expression!(ae, v_var[dc_bus_name, t], 1.0) + v_expr[name, t] = ae end end - return + return v_expr +end + +function _converter_vi_bounds(devices) + n = length(devices) + v_bounds = Vector{IOM.MinMax}(undef, n) + i_bounds = Vector{IOM.MinMax}(undef, n) + for (k, d) in enumerate(devices) + v_min, v_max = PSY.get_voltage_limits(PSY.get_dc_bus(d)) + i_max = PSY.get_max_dc_current(d) + v_bounds[k] = IOM.MinMax((min = v_min, max = v_max)) + i_bounds[k] = IOM.MinMax((min = -i_max, max = i_max)) + end + return v_bounds, i_bounds end function add_constraints!( container::OptimizationContainer, - ::Type{ConverterMcCormickEnvelopes}, + ::Type{ConverterLossConstraint}, devices::IS.FlattenIteratorWrapper{U}, model::DeviceModel{U, V}, - network_model::NetworkModel{X}, + ::NetworkModel{X}, ) where { U <: PSY.InterconnectingConverter, - V <: QuadraticLossConverter, + V <: AbstractQuadraticLossConverter, X <: AbstractActivePowerModel, } - time_steps = get_time_steps(container) - varcurrent = get_variable(container, ConverterCurrent, U) - var_dcvoltage = get_variable(container, DCVoltage, PSY.DCBus) - var_dc_power = get_variable(container, ConverterDCPower, U) - ipc_names = axes(varcurrent, 1) - constraint1_under = - add_constraints_container!(container, ConverterMcCormickEnvelopes, - U, - ipc_names, - time_steps; - meta = "under_1", - ) - constraint2_under = - add_constraints_container!(container, ConverterMcCormickEnvelopes, - U, - ipc_names, - time_steps; - meta = "under_2", - ) - constraint1_over = - add_constraints_container!(container, ConverterMcCormickEnvelopes, - U, - ipc_names, - time_steps; - meta = "over_1", - ) - constraint2_over = - add_constraints_container!(container, ConverterMcCormickEnvelopes, - U, - ipc_names, - time_steps; - meta = "over_2", - ) - - for device in devices - name = PSY.get_name(device) - dc_bus = PSY.get_dc_bus(device) - dc_bus_name = PSY.get_name(dc_bus) - V_min, V_max = PSY.get_voltage_limits(dc_bus) - I_max = PSY.get_max_dc_current(device) - I_min = -I_max - for t in time_steps - constraint1_under[name, t] = JuMP.@constraint( - get_jump_model(container), - var_dc_power[name, t] >= - V_min * varcurrent[name, t] + var_dcvoltage[dc_bus_name, t] * I_min - - I_min * V_min - ) - constraint2_under[name, t] = JuMP.@constraint( - get_jump_model(container), - var_dc_power[name, t] >= - V_max * varcurrent[name, t] + var_dcvoltage[dc_bus_name, t] * I_max - - I_max * V_max - ) - constraint1_over[name, t] = JuMP.@constraint( - get_jump_model(container), - var_dc_power[name, t] <= - V_max * varcurrent[name, t] + var_dcvoltage[dc_bus_name, t] * I_min - - I_min * V_max - ) - constraint2_over[name, t] = JuMP.@constraint( - get_jump_model(container), - var_dc_power[name, t] <= - V_min * varcurrent[name, t] + var_dcvoltage[dc_bus_name, t] * I_max - - I_max * V_min - ) - end - end + _add_converter_loss_constraint!( + container, devices, model; + use_linear_loss = get_attribute(model, "use_linear_loss"), + ) return end -function add_constraints!( +function _add_converter_loss_constraint!( container::OptimizationContainer, - ::Type{ConverterLossConstraint}, devices::IS.FlattenIteratorWrapper{U}, - model::DeviceModel{U, V}, - network_model::NetworkModel{X}, + model::DeviceModel{U, V}; + use_linear_loss::Bool, ) where { U <: PSY.InterconnectingConverter, - V <: QuadraticLossConverter, - X <: AbstractActivePowerModel, + V <: AbstractQuadraticLossConverter, } time_steps = get_time_steps(container) - var_sq_current = get_variable(container, SquaredConverterCurrent, U) - var_ac_power = get_variable(container, ActivePowerVariable, U) - var_dc_power = get_variable(container, ConverterDCPower, U) - ipc_names = axes(var_sq_current, 1) - constraint = - add_constraints_container!(container, ConverterLossConstraint, - U, - ipc_names, - time_steps, - ) - - use_linear_loss = get_attribute(model, "use_linear_loss") + P_ac_var = get_variable(container, ActivePowerVariable, U) + vi_expr = get_expression(container, IOM.BilinearProductExpression, U, "vi") + i_sq_expr = get_expression(container, IOM.QuadraticExpression, U, "i_sq") if use_linear_loss - pos_current = get_variable(container, ConverterPositiveCurrent, U) - neg_current = get_variable(container, ConverterNegativeCurrent, U) + i_pos_var = get_variable(container, ConverterPositiveCurrent, U) + i_neg_var = get_variable(container, ConverterNegativeCurrent, U) end + ipc_names = [PSY.get_name(d) for d in devices] + loss_const = add_constraints_container!( + container, ConverterLossConstraint, U, ipc_names, time_steps, + ) + + jump_model = get_jump_model(container) for device in devices name = PSY.get_name(device) loss_function = PSY.get_loss_function(device) @@ -745,16 +633,13 @@ function add_constraints!( c = PSY.get_constant_term(loss_function) end for t in time_steps + loss = a * i_sq_expr[name, t] + c if use_linear_loss - loss = - a * var_sq_current[name, t] + - b * (pos_current[name, t] + neg_current[name, t]) + c - else - loss = a * var_sq_current[name, t] + c + loss += b * (i_pos_var[name, t] + i_neg_var[name, t]) end - constraint[name, t] = JuMP.@constraint( - get_jump_model(container), - var_ac_power[name, t] == var_dc_power[name, t] - loss + loss_const[name, t] = JuMP.@constraint( + jump_model, + P_ac_var[name, t] == vi_expr[name, t] - loss ) end end @@ -770,175 +655,49 @@ function add_constraints!( ) where { T <: CurrentAbsoluteValueConstraint, U <: PSY.InterconnectingConverter, - V <: QuadraticLossConverter, + V <: AbstractQuadraticLossConverter, } time_steps = get_time_steps(container) names = [PSY.get_name(d) for d in devices] - JuMPmodel = get_jump_model(container) - # current vars # - current_var = get_variable(container, ConverterCurrent, U) # From direction - current_var_pos = get_variable(container, ConverterPositiveCurrent, U) # From direction - current_var_neg = get_variable(container, ConverterNegativeCurrent, U) # From direction - current_dir = get_variable(container, ConverterCurrentDirection, U) - - constraint = - add_constraints_container!(container, CurrentAbsoluteValueConstraint, - U, - names, - time_steps, - ) - constraint_pos_ub = - add_constraints_container!(container, CurrentAbsoluteValueConstraint, - U, - names, - time_steps; - meta = "pos_ub", - ) - constraint_neg_ub = - add_constraints_container!(container, CurrentAbsoluteValueConstraint, - U, - names, - time_steps; - meta = "neg_ub", - ) + 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", + ) for d in devices name = PSY.get_name(d) - I_max = PSY.get_max_dc_current(d) + i_max = PSY.get_max_dc_current(d) for t in time_steps - constraint[name, t] = JuMP.@constraint( - JuMPmodel, - current_var[name, t] == current_var_pos[name, t] - current_var_neg[name, t] + abs_val_const[name, t] = JuMP.@constraint( + jump_model, + i_var[name, t] == i_pos_var[name, t] - i_neg_var[name, t] ) - constraint_pos_ub[name, t] = JuMP.@constraint( - JuMPmodel, - current_var_pos[name, t] <= I_max * current_dir[name, t] + pos_ub_const[name, t] = JuMP.@constraint( + jump_model, + i_pos_var[name, t] <= i_max * i_dir_var[name, t] ) - constraint_neg_ub[name, t] = JuMP.@constraint( - JuMPmodel, - current_var_neg[name, t] <= I_max * (1 - current_dir[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 -function add_constraints!( - container::OptimizationContainer, - ::Type{T}, - devices::IS.FlattenIteratorWrapper{U}, - model::DeviceModel{U, V}, - ::NetworkModel{<:AbstractPowerModel}, -) where { - T <: InterpolationVoltageConstraints, - U <: PSY.InterconnectingConverter, - V <: QuadraticLossConverter, -} - dic_var_bkpts = Dict{String, Vector{Float64}}() - dic_function_bkpts = Dict{String, Vector{Float64}}() - num_segments = get_attribute(model, "voltage_segments") - for d in devices - name = PSY.get_name(d) - vmin, vmax = PSY.get_voltage_limits(d.dc_bus) - var_bkpts, function_bkpts = - _get_breakpoints_for_pwl_function(vmin, vmax, x -> x^2; num_segments) - dic_var_bkpts[name] = var_bkpts - dic_function_bkpts[name] = function_bkpts - end - - _add_generic_incremental_interpolation_constraint!( - container, - DCVoltage, - SquaredDCVoltage, - InterpolationSquaredVoltageVariable, - InterpolationBinarySquaredVoltageVariable, - InterpolationVoltageConstraints, - devices, - dic_var_bkpts, - dic_function_bkpts, - ) - return -end - -function add_constraints!( - container::OptimizationContainer, - ::Type{T}, - devices::IS.FlattenIteratorWrapper{U}, - model::DeviceModel{U, V}, - ::NetworkModel{<:AbstractPowerModel}, -) where { - T <: InterpolationCurrentConstraints, - U <: PSY.InterconnectingConverter, - V <: QuadraticLossConverter, -} - dic_var_bkpts = Dict{String, Vector{Float64}}() - dic_function_bkpts = Dict{String, Vector{Float64}}() - num_segments = get_attribute(model, "current_segments") - for d in devices - name = PSY.get_name(d) - Imax = PSY.get_max_dc_current(d) - Imin = -Imax - var_bkpts, function_bkpts = - _get_breakpoints_for_pwl_function(Imin, Imax, x -> x^2; num_segments) - dic_var_bkpts[name] = var_bkpts - dic_function_bkpts[name] = function_bkpts - end - - _add_generic_incremental_interpolation_constraint!( - container, - ConverterCurrent, - SquaredConverterCurrent, - InterpolationSquaredCurrentVariable, - InterpolationBinarySquaredCurrentVariable, - InterpolationCurrentConstraints, - devices, - dic_var_bkpts, - dic_function_bkpts, - ) - return -end - -function add_constraints!( - container::OptimizationContainer, - ::Type{T}, - devices::IS.FlattenIteratorWrapper{U}, - model::DeviceModel{U, V}, - ::NetworkModel{<:AbstractPowerModel}, -) where { - T <: InterpolationBilinearConstraints, - U <: PSY.InterconnectingConverter, - V <: QuadraticLossConverter, -} - dic_var_bkpts = Dict{String, Vector{Float64}}() - dic_function_bkpts = Dict{String, Vector{Float64}}() - num_segments = get_attribute(model, "bilinear_segments") - for d in devices - name = PSY.get_name(d) - vmin, vmax = PSY.get_voltage_limits(d.dc_bus) - Imax = PSY.get_max_dc_current(d) - Imin = -Imax - γ_min = vmin * Imin - γ_max = vmax * Imax - var_bkpts, function_bkpts = - _get_breakpoints_for_pwl_function(γ_min, γ_max, x -> x^2; num_segments) - dic_var_bkpts[name] = var_bkpts - dic_function_bkpts[name] = function_bkpts - end - - _add_generic_incremental_interpolation_constraint!( - container, - AuxBilinearConverterVariable, - AuxBilinearSquaredConverterVariable, - InterpolationSquaredBilinearVariable, - InterpolationBinarySquaredBilinearVariable, - InterpolationBilinearConstraints, - devices, - dic_var_bkpts, - dic_function_bkpts, - ) - return -end - ############################################ ########### Objective Function ############# ############################################ diff --git a/src/mt_hvdc_models/hvdcsystems_constructor.jl b/src/mt_hvdc_models/hvdcsystems_constructor.jl index 138eada..0713519 100644 --- a/src/mt_hvdc_models/hvdcsystems_constructor.jl +++ b/src/mt_hvdc_models/hvdcsystems_constructor.jl @@ -48,105 +48,25 @@ function construct_device!( container::OptimizationContainer, sys::PSY.System, ::ArgumentConstructStage, - model::DeviceModel{PSY.InterconnectingConverter, QuadraticLossConverter}, + model::DeviceModel{PSY.InterconnectingConverter, T}, network_model::NetworkModel{<:AbstractActivePowerModel}, -) - devices = get_available_components( - model, - sys, - ) - ##################### - ##### Variables ##### - ##################### - - # Add Power Variable - add_variables!(container, ActivePowerVariable, devices, QuadraticLossConverter) # p_c^{ac} - add_variables!(container, ConverterDCPower, devices, QuadraticLossConverter) # p_c - # Add Current Variables: i, i+, i- - add_variables!(container, ConverterCurrent, devices, QuadraticLossConverter) # i - add_variables!(container, SquaredConverterCurrent, devices, QuadraticLossConverter) # i^sq - use_linear_loss = get_attribute(model, "use_linear_loss") - if use_linear_loss - add_variables!( - container, - ConverterPositiveCurrent, - devices, - QuadraticLossConverter, - ) # i^+ - add_variables!( - container, - ConverterNegativeCurrent, - devices, - QuadraticLossConverter, - ) # i^- - add_variables!( - container, - ConverterCurrentDirection, - devices, - QuadraticLossConverter, - ) # ν +) where {T <: AbstractQuadraticLossConverter} + devices = get_available_components(model, sys) + add_variables!(container, ActivePowerVariable, devices, T) + add_variables!(container, ConverterCurrent, devices, T) + if get_attribute(model, "use_linear_loss") + add_variables!(container, ConverterPositiveCurrent, devices, T) + add_variables!(container, ConverterNegativeCurrent, devices, T) + add_variables!(container, ConverterCurrentDirection, devices, T) end - # Add Voltage Variables: v^sq - add_variables!(container, SquaredDCVoltage, devices, QuadraticLossConverter) - # Add Bilinear Variables: γ, γ^{sq} - add_variables!( - container, - AuxBilinearConverterVariable, - devices, - QuadraticLossConverter, - ) # γ - add_variables!( - container, - AuxBilinearSquaredConverterVariable, - devices, - QuadraticLossConverter, - ) # γ^{sq} - - #### Add Interpolation Variables #### - - v_segments = get_attribute(model, "voltage_segments") - i_segments = get_attribute(model, "current_segments") - γ_segments = get_attribute(model, "bilinear_segments") - - vars_vector = [ - # Voltage v # - (InterpolationSquaredVoltageVariable, v_segments), # δ^v - (InterpolationBinarySquaredVoltageVariable, v_segments), # z^v - # Current i # - (InterpolationSquaredCurrentVariable, i_segments), # δ^i - (InterpolationBinarySquaredCurrentVariable, i_segments), # z^i - # Bilinear γ # - (InterpolationSquaredBilinearVariable, γ_segments), # δ^γ - (InterpolationBinarySquaredBilinearVariable, γ_segments), # z^γ - ] - - for (T, len_segments) in vars_vector - add_sparse_pwl_interpolation_variables!(container, T, - devices, - model, - len_segments, - ) - end - - ##################### - #### Expressions #### - ##################### add_to_expression!( - container, - ActivePowerBalance, - ActivePowerVariable, - devices, - model, - network_model, + container, ActivePowerBalance, ActivePowerVariable, + devices, model, network_model, ) add_to_expression!( - container, - DCCurrentBalance, - ConverterCurrent, - devices, - model, - network_model, + container, DCCurrentBalance, ConverterCurrent, + devices, model, network_model, ) add_feedforward_arguments!(container, model, devices) return @@ -156,75 +76,89 @@ function construct_device!( container::OptimizationContainer, sys::PSY.System, ::ModelConstructStage, - model::DeviceModel{PSY.InterconnectingConverter, QuadraticLossConverter}, + model::DeviceModel{PSY.InterconnectingConverter, Bin2QuadraticLossConverter}, network_model::NetworkModel{<:AbstractActivePowerModel}, ) - devices = get_available_components( - model, - sys, - ) + devices = get_available_components(model, sys) + time_steps = get_time_steps(container) + ipc_names = [PSY.get_name(d) for d in devices] + v_bounds, i_bounds = _converter_vi_bounds(devices) + v_expr = _voltage_expr_per_converter(container, devices, ipc_names, time_steps) + i_var = get_variable(container, ConverterCurrent, PSY.InterconnectingConverter) + + v_sq_expr = IOM._add_quadratic_approx!( + IOM.ManualSOS2QuadConfig(IOM.DEFAULT_INTERPOLATION_LENGTH), + container, PSY.InterconnectingConverter, + ipc_names, time_steps, + v_expr, v_bounds, + "v_sq", + ) + i_sq_expr = IOM._add_quadratic_approx!( + IOM.ManualSOS2QuadConfig(IOM.DEFAULT_INTERPOLATION_LENGTH), + container, PSY.InterconnectingConverter, + ipc_names, time_steps, + i_var, i_bounds, + "i_sq", + ) + IOM._add_bilinear_approx!( + IOM.Bin2Config(IOM.ManualSOS2QuadConfig(IOM.DEFAULT_INTERPOLATION_LENGTH)), + container, PSY.InterconnectingConverter, + ipc_names, time_steps, + v_sq_expr, i_sq_expr, + v_expr, i_var, + v_bounds, i_bounds, + "vi", + ) + + add_constraints!(container, ConverterLossConstraint, devices, model, network_model) + add_constraints!(container, CurrentAbsoluteValueConstraint, devices, model, network_model) - add_constraints!( - container, - ConverterPowerCalculationConstraint, - devices, - model, - network_model, - ) - add_constraints!( - container, - ConverterMcCormickEnvelopes, - devices, - model, - network_model, - ) - 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), ) - use_linear_loss = get_attribute(model, "use_linear_loss") - if use_linear_loss - add_constraints!( - container, - CurrentAbsoluteValueConstraint, - devices, - model, - network_model, - ) + return +end + +function construct_device!( + container::OptimizationContainer, + sys::PSY.System, + ::ModelConstructStage, + model::DeviceModel{PSY.InterconnectingConverter, QuadraticLossConverter}, + network_model::NetworkModel{<:AbstractActivePowerModel}, +) + devices = get_available_components(model, sys) + time_steps = get_time_steps(container) + ipc_names = [PSY.get_name(d) for d in devices] + v_bounds, i_bounds = _converter_vi_bounds(devices) + v_expr = _voltage_expr_per_converter(container, devices, ipc_names, time_steps) + i_var = get_variable(container, ConverterCurrent, PSY.InterconnectingConverter) + + IOM._add_quadratic_approx!( + IOM.NoQuadApproxConfig(), + container, PSY.InterconnectingConverter, + ipc_names, time_steps, + i_var, i_bounds, + "i_sq", + ) + IOM._add_bilinear_approx!( + IOM.NoBilinearApproxConfig(), + container, PSY.InterconnectingConverter, + ipc_names, time_steps, + v_expr, i_var, + v_bounds, i_bounds, + "vi", + ) + + add_constraints!(container, ConverterLossConstraint, devices, model, network_model) + if get_attribute(model, "use_linear_loss") + add_constraints!(container, CurrentAbsoluteValueConstraint, devices, model, network_model) end - add_constraints!( - container, - InterpolationVoltageConstraints, - devices, - model, - network_model, - ) - add_constraints!( - container, - InterpolationCurrentConstraints, - devices, - model, - network_model, - ) - add_constraints!( - container, - InterpolationBilinearConstraints, - devices, - model, - network_model, - ) add_feedforward_constraints!(container, model, devices) add_to_objective_function!( - container, - devices, - model, - get_network_formulation(network_model), + container, devices, model, get_network_formulation(network_model), ) - #add_constraint_dual!(container, sys, model) return end diff --git a/test/Project.toml b/test/Project.toml index 9445668..a7425a5 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -32,10 +32,10 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" [sources] -InfrastructureSystems = {url = "https://github.com/NREL-Sienna/InfrastructureSystems.jl", rev = "IS4"} -PowerSystems = {url = "https://github.com/NREL-Sienna/PowerSystems.jl", rev = "psy6"} -InfrastructureOptimizationModels = {url = "https://github.com/NREL-Sienna/InfrastructureOptimizationModels.jl", rev = "lk/pom-test-fixes"} -PowerSystemCaseBuilder = {url = "https://github.com/NREL-Sienna/PowerSystemCaseBuilder.jl", rev = "psy6"} +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"} + [compat] HiGHS = "1" diff --git a/test/test_device_hvdc.jl b/test/test_device_hvdc.jl index 2b158b3..5452cdf 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" begin +@testset "HVDC System with Losses Network (Bin2QuadraticLossConverter)" begin sys = _generate_test_hvdc_sys() template = OperationsProblemTemplate() set_device_model!(template, ThermalStandard, ThermalDispatchNoMin) @@ -106,13 +106,7 @@ end set_device_model!(template, TModelHVDCLine, DCLossyLine) ipc_model = DeviceModel( InterconnectingConverter, - QuadraticLossConverter; - attributes = Dict( - "voltage_segments" => 3, - "current_segments" => 3, - "bilinear_segments" => 3, - "use_linear_loss" => true, - ), + Bin2QuadraticLossConverter, ) set_device_model!(template, ipc_model) set_hvdc_network_model!(template, VoltageDispatchHVDCNetworkModel) @@ -127,3 +121,28 @@ end IOM.ModelBuildStatus.BUILT @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED end + +@testset "HVDC System with Losses Network (QuadraticLossConverter NLP)" begin + sys = _generate_test_hvdc_sys() + 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) + ipc_model = DeviceModel( + InterconnectingConverter, + QuadraticLossConverter, + ) + set_device_model!(template, ipc_model) + set_hvdc_network_model!(template, VoltageDispatchHVDCNetworkModel) + model = + DecisionModel( + template, + sys; + store_variable_names = true, + optimizer = ipopt_optimizer, + ) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED +end From b39c63d4470c3070cb4cdaa5a018b09c9054b916 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Wed, 29 Apr 2026 16:43:17 -0600 Subject: [PATCH 2/9] formatting --- src/mt_hvdc_models/hvdcsystems_constructor.jl | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/mt_hvdc_models/hvdcsystems_constructor.jl b/src/mt_hvdc_models/hvdcsystems_constructor.jl index 0713519..7d6e712 100644 --- a/src/mt_hvdc_models/hvdcsystems_constructor.jl +++ b/src/mt_hvdc_models/hvdcsystems_constructor.jl @@ -111,7 +111,13 @@ function construct_device!( ) add_constraints!(container, ConverterLossConstraint, devices, model, network_model) - add_constraints!(container, CurrentAbsoluteValueConstraint, devices, model, network_model) + add_constraints!( + container, + CurrentAbsoluteValueConstraint, + devices, + model, + network_model, + ) add_feedforward_constraints!(container, model, devices) add_to_objective_function!( @@ -152,7 +158,13 @@ function construct_device!( add_constraints!(container, ConverterLossConstraint, devices, model, network_model) if get_attribute(model, "use_linear_loss") - add_constraints!(container, CurrentAbsoluteValueConstraint, devices, model, network_model) + add_constraints!( + container, + CurrentAbsoluteValueConstraint, + devices, + model, + network_model, + ) end add_feedforward_constraints!(container, model, devices) From b7457919664706083d1fa92e4480fe7b95f88f4b Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Wed, 29 Apr 2026 17:40:01 -0600 Subject: [PATCH 3/9] fix linear loss adding;fix hdf output --- src/mt_hvdc_models/HVDCsystems.jl | 35 ++++-------------- src/mt_hvdc_models/hvdcsystems_constructor.jl | 22 ++++++----- src/operation/decision_model.jl | 4 +- src/operation/emulation_model.jl | 4 +- test/test_device_hvdc.jl | 37 +++++++++++++++++++ 5 files changed, 61 insertions(+), 41 deletions(-) diff --git a/src/mt_hvdc_models/HVDCsystems.jl b/src/mt_hvdc_models/HVDCsystems.jl index f374d29..02754e8 100644 --- a/src/mt_hvdc_models/HVDCsystems.jl +++ b/src/mt_hvdc_models/HVDCsystems.jl @@ -578,6 +578,9 @@ function _converter_vi_bounds(devices) return v_bounds, i_bounds end +_get_quadratic_term(loss_fn::PSY.QuadraticCurve) = PSY.get_quadratic_term(loss_fn) +_get_quadratic_term(loss_fn) = 0.0 + function add_constraints!( container::OptimizationContainer, ::Type{ConverterLossConstraint}, @@ -588,28 +591,12 @@ function add_constraints!( U <: PSY.InterconnectingConverter, V <: AbstractQuadraticLossConverter, X <: AbstractActivePowerModel, -} - _add_converter_loss_constraint!( - container, devices, model; - use_linear_loss = get_attribute(model, "use_linear_loss"), - ) - return -end - -function _add_converter_loss_constraint!( - container::OptimizationContainer, - devices::IS.FlattenIteratorWrapper{U}, - model::DeviceModel{U, V}; - use_linear_loss::Bool, -) where { - U <: PSY.InterconnectingConverter, - V <: AbstractQuadraticLossConverter, } time_steps = get_time_steps(container) P_ac_var = get_variable(container, ActivePowerVariable, U) vi_expr = get_expression(container, IOM.BilinearProductExpression, U, "vi") i_sq_expr = get_expression(container, IOM.QuadraticExpression, U, "i_sq") - if use_linear_loss + if get_attribute(model, "use_linear_loss") i_pos_var = get_variable(container, ConverterPositiveCurrent, U) i_neg_var = get_variable(container, ConverterNegativeCurrent, U) end @@ -623,18 +610,12 @@ function _add_converter_loss_constraint!( for device in devices name = PSY.get_name(device) loss_function = PSY.get_loss_function(device) - if isa(loss_function, PSY.QuadraticCurve) - a = PSY.get_quadratic_term(loss_function) - b = PSY.get_proportional_term(loss_function) - c = PSY.get_constant_term(loss_function) - else - a = 0.0 - b = PSY.get_proportional_term(loss_function) - c = PSY.get_constant_term(loss_function) - end + a = _get_quadratic_term(loss_function) + 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 + if get_attribute(model, "use_linear_loss") loss += b * (i_pos_var[name, t] + i_neg_var[name, t]) end loss_const[name, t] = JuMP.@constraint( diff --git a/src/mt_hvdc_models/hvdcsystems_constructor.jl b/src/mt_hvdc_models/hvdcsystems_constructor.jl index 7d6e712..4df2d00 100644 --- a/src/mt_hvdc_models/hvdcsystems_constructor.jl +++ b/src/mt_hvdc_models/hvdcsystems_constructor.jl @@ -87,21 +87,21 @@ function construct_device!( i_var = get_variable(container, ConverterCurrent, PSY.InterconnectingConverter) v_sq_expr = IOM._add_quadratic_approx!( - IOM.ManualSOS2QuadConfig(IOM.DEFAULT_INTERPOLATION_LENGTH), + IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH), container, PSY.InterconnectingConverter, ipc_names, time_steps, v_expr, v_bounds, "v_sq", ) i_sq_expr = IOM._add_quadratic_approx!( - IOM.ManualSOS2QuadConfig(IOM.DEFAULT_INTERPOLATION_LENGTH), + IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH), container, PSY.InterconnectingConverter, ipc_names, time_steps, i_var, i_bounds, "i_sq", ) IOM._add_bilinear_approx!( - IOM.Bin2Config(IOM.ManualSOS2QuadConfig(IOM.DEFAULT_INTERPOLATION_LENGTH)), + IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)), container, PSY.InterconnectingConverter, ipc_names, time_steps, v_sq_expr, i_sq_expr, @@ -111,13 +111,15 @@ function construct_device!( ) add_constraints!(container, ConverterLossConstraint, devices, model, network_model) - add_constraints!( - container, - CurrentAbsoluteValueConstraint, - devices, - model, - network_model, - ) + if get_attribute(model, "use_linear_loss") + add_constraints!( + container, + CurrentAbsoluteValueConstraint, + devices, + model, + network_model, + ) + end add_feedforward_constraints!(container, model, devices) add_to_objective_function!( diff --git a/src/operation/decision_model.jl b/src/operation/decision_model.jl index b6ce1b9..be9cadc 100644 --- a/src/operation/decision_model.jl +++ b/src/operation/decision_model.jl @@ -67,7 +67,7 @@ function build!( IOM.register_recorders!(model, file_mode) logger = IS.configure_logging(get_internal(model), IOM.PROBLEM_LOG_FILENAME, file_mode) if store_system_in_results - @warn "store_system_in_results for $(model) is set to true. This will do nothing unless a Simulation is being built." + @warn "store_system_in_results for $(model.name) is set to true. This will do nothing unless a Simulation is being built." end try Logging.with_logger(logger) do @@ -151,7 +151,7 @@ function solve!( kwargs..., ) if store_system_in_results - @warn "store_system_in_results for $(model) is set to true. This will do nothing unless a Simulation is being built." + @warn "store_system_in_results for $(model.name) is set to true. This will do nothing unless a Simulation is being built." end build_if_not_already_built!( model; diff --git a/src/operation/emulation_model.jl b/src/operation/emulation_model.jl index 9241729..6a4d6f3 100644 --- a/src/operation/emulation_model.jl +++ b/src/operation/emulation_model.jl @@ -69,7 +69,7 @@ function build!( file_mode, ) if store_system_in_results - @warn "store_system_in_results for $(model) is set to true. This will do nothing unless a Simulation is being built." + @warn "store_system_in_results for $(model.name) is set to true. This will do nothing unless a Simulation is being built." end try Logging.with_logger(logger) do @@ -197,7 +197,7 @@ function run!( kwargs..., ) if store_system_in_results - @warn "store_system_in_results for $(model) is set to true. This will do nothing unless a Simulation is being built." + @warn "store_system_in_results for $(model.name) is set to true. This will do nothing unless a Simulation is being built." end build_if_not_already_built!( model; diff --git a/test/test_device_hvdc.jl b/test/test_device_hvdc.jl index 5452cdf..2250bd3 100644 --- a/test/test_device_hvdc.jl +++ b/test/test_device_hvdc.jl @@ -146,3 +146,40 @@ end IOM.ModelBuildStatus.BUILT @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED end + +@testset "HVDC Bin2 vs NLP QuadraticLossConverter agreement" begin + function _build_and_solve(formulation, optimizer) + sys = _generate_test_hvdc_sys() + 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, formulation); + attributes = Dict("use_linear_loss" => false) + ) + set_hvdc_network_model!(template, VoltageDispatchHVDCNetworkModel) + model = DecisionModel( + template, + sys; + store_variable_names = true, + optimizer = optimizer, + ) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + return model + end + + bin2_model = _build_and_solve(Bin2QuadraticLossConverter, HiGHS_optimizer) + nlp_model = _build_and_solve(QuadraticLossConverter, ipopt_optimizer) + + bin2_obj = IOM.get_objective_value(OptimizationProblemOutputs(bin2_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 From 4c8c8fee51d11da86173c1657bf12e3871bed410 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Thu, 30 Apr 2026 09:33:44 -0600 Subject: [PATCH 4/9] formatting --- test/test_device_hvdc.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_device_hvdc.jl b/test/test_device_hvdc.jl index 2250bd3..2f55602 100644 --- a/test/test_device_hvdc.jl +++ b/test/test_device_hvdc.jl @@ -158,7 +158,7 @@ end set_device_model!( template, DeviceModel(InterconnectingConverter, formulation); - attributes = Dict("use_linear_loss" => false) + attributes = Dict("use_linear_loss" => false), ) set_hvdc_network_model!(template, VoltageDispatchHVDCNetworkModel) model = DecisionModel( From abbc7929788425b2fb7c71128f7eab85ed7f6fe9 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Thu, 30 Apr 2026 15:04:37 -0600 Subject: [PATCH 5/9] move _voltage helpers to constructor combine ModelConstruct stages for both formulations only add linear loss formulations for devices with non-zero linear loss --- src/mt_hvdc_models/HVDCsystems.jl | 60 +++------ src/mt_hvdc_models/hvdcsystems_constructor.jl | 123 +++++++++--------- test/test_device_hvdc.jl | 35 ++++- 3 files changed, 109 insertions(+), 109 deletions(-) diff --git a/src/mt_hvdc_models/HVDCsystems.jl b/src/mt_hvdc_models/HVDCsystems.jl index 02754e8..b3fcde8 100644 --- a/src/mt_hvdc_models/HVDCsystems.jl +++ b/src/mt_hvdc_models/HVDCsystems.jl @@ -142,9 +142,7 @@ function get_default_attributes( ::Type{PSY.InterconnectingConverter}, ::Type{Bin2QuadraticLossConverter}, ) - return Dict{String, Any}( - "use_linear_loss" => true - ) + return Dict{String, Any}() end function get_default_attributes( @@ -542,45 +540,20 @@ end ############## Converters ################## -# Builds, per converter, an AffExpr container of the DC-bus voltage variable indexed -# by converter name. IOM bilinear/quadratic approximations consume an x_var indexed by -# component name; DCVoltage is stored per DC bus, so we wrap it for compatibility. -function _voltage_expr_per_converter( - container::OptimizationContainer, - devices, - ipc_names::Vector{String}, - time_steps, -) - v_var = get_variable(container, DCVoltage, PSY.DCBus) - v_expr = JuMP.Containers.DenseAxisArray{JuMP.AffExpr}(undef, ipc_names, time_steps) - for d in devices - name = PSY.get_name(d) - dc_bus_name = PSY.get_name(PSY.get_dc_bus(d)) - for t in time_steps - ae = JuMP.AffExpr(0.0) - add_proportional_to_jump_expression!(ae, v_var[dc_bus_name, t], 1.0) - v_expr[name, t] = ae - end - end - return v_expr -end - -function _converter_vi_bounds(devices) - n = length(devices) - v_bounds = Vector{IOM.MinMax}(undef, n) - i_bounds = Vector{IOM.MinMax}(undef, n) - for (k, d) in enumerate(devices) - v_min, v_max = PSY.get_voltage_limits(PSY.get_dc_bus(d)) - i_max = PSY.get_max_dc_current(d) - v_bounds[k] = IOM.MinMax((min = v_min, max = v_max)) - i_bounds[k] = IOM.MinMax((min = -i_max, max = i_max)) - end - return v_bounds, i_bounds -end - _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}, @@ -596,7 +569,9 @@ function add_constraints!( P_ac_var = get_variable(container, ActivePowerVariable, U) vi_expr = get_expression(container, IOM.BilinearProductExpression, U, "vi") i_sq_expr = get_expression(container, IOM.QuadraticExpression, U, "i_sq") - if get_attribute(model, "use_linear_loss") + 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) end @@ -615,7 +590,7 @@ function add_constraints!( c = PSY.get_constant_term(loss_function) for t in time_steps loss = a * i_sq_expr[name, t] + c - if get_attribute(model, "use_linear_loss") + if use_linear_loss && !iszero(b) loss += b * (i_pos_var[name, t] + i_neg_var[name, t]) end loss_const[name, t] = JuMP.@constraint( @@ -630,13 +605,14 @@ end function add_constraints!( container::OptimizationContainer, ::Type{T}, - devices::IS.FlattenIteratorWrapper{U}, + devices::W, ::DeviceModel{U, V}, ::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] diff --git a/src/mt_hvdc_models/hvdcsystems_constructor.jl b/src/mt_hvdc_models/hvdcsystems_constructor.jl index 4df2d00..c1ec321 100644 --- a/src/mt_hvdc_models/hvdcsystems_constructor.jl +++ b/src/mt_hvdc_models/hvdcsystems_constructor.jl @@ -54,10 +54,15 @@ function construct_device!( devices = get_available_components(model, sys) add_variables!(container, ActivePowerVariable, devices, T) add_variables!(container, ConverterCurrent, devices, T) - if get_attribute(model, "use_linear_loss") - add_variables!(container, ConverterPositiveCurrent, devices, T) - add_variables!(container, ConverterNegativeCurrent, devices, T) - add_variables!(container, ConverterCurrentDirection, 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, ConverterPositiveCurrent, ll_devices, T) + add_variables!(container, ConverterNegativeCurrent, ll_devices, T) + add_variables!(container, ConverterCurrentDirection, ll_devices, T) + end end add_to_expression!( @@ -72,13 +77,45 @@ 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 _converter_vi_bounds(devices) + n = length(devices) + v_bounds = Vector{IOM.MinMax}(undef, n) + i_bounds = Vector{IOM.MinMax}(undef, n) + for (k, d) in enumerate(devices) + v_min, v_max = PSY.get_voltage_limits(PSY.get_dc_bus(d)) + i_max = PSY.get_max_dc_current(d) + v_bounds[k] = IOM.MinMax((min = v_min, max = v_max)) + i_bounds[k] = IOM.MinMax((min = -i_max, max = i_max)) + end + return v_bounds, i_bounds +end + +_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{QuadraticLossConverter}) = IOM.NoBilinearApproxConfig() + function construct_device!( container::OptimizationContainer, sys::PSY.System, ::ModelConstructStage, - model::DeviceModel{PSY.InterconnectingConverter, Bin2QuadraticLossConverter}, + model::DeviceModel{PSY.InterconnectingConverter, T}, network_model::NetworkModel{<:AbstractActivePowerModel}, -) +) where {T <: AbstractQuadraticLossConverter} devices = get_available_components(model, sys) time_steps = get_time_steps(container) ipc_names = [PSY.get_name(d) for d in devices] @@ -86,22 +123,23 @@ function construct_device!( v_expr = _voltage_expr_per_converter(container, devices, ipc_names, time_steps) i_var = get_variable(container, ConverterCurrent, PSY.InterconnectingConverter) + quad_cfg, bilin_cfg = _quad_config(T), _bilinear_config(T) v_sq_expr = IOM._add_quadratic_approx!( - IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH), + quad_cfg, container, PSY.InterconnectingConverter, ipc_names, time_steps, v_expr, v_bounds, "v_sq", ) i_sq_expr = IOM._add_quadratic_approx!( - IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH), + quad_cfg, container, PSY.InterconnectingConverter, ipc_names, time_steps, i_var, i_bounds, "i_sq", ) IOM._add_bilinear_approx!( - IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)), + bilin_cfg, container, PSY.InterconnectingConverter, ipc_names, time_steps, v_sq_expr, i_sq_expr, @@ -111,62 +149,17 @@ function construct_device!( ) add_constraints!(container, ConverterLossConstraint, devices, model, network_model) - if get_attribute(model, "use_linear_loss") - add_constraints!( - container, - CurrentAbsoluteValueConstraint, - devices, - model, - network_model, - ) - end - - add_feedforward_constraints!(container, model, devices) - add_to_objective_function!( - container, devices, model, get_network_formulation(network_model), - ) - return -end - -function construct_device!( - container::OptimizationContainer, - sys::PSY.System, - ::ModelConstructStage, - model::DeviceModel{PSY.InterconnectingConverter, QuadraticLossConverter}, - network_model::NetworkModel{<:AbstractActivePowerModel}, -) - devices = get_available_components(model, sys) - time_steps = get_time_steps(container) - ipc_names = [PSY.get_name(d) for d in devices] - v_bounds, i_bounds = _converter_vi_bounds(devices) - v_expr = _voltage_expr_per_converter(container, devices, ipc_names, time_steps) - i_var = get_variable(container, ConverterCurrent, PSY.InterconnectingConverter) - - IOM._add_quadratic_approx!( - IOM.NoQuadApproxConfig(), - container, PSY.InterconnectingConverter, - ipc_names, time_steps, - i_var, i_bounds, - "i_sq", - ) - IOM._add_bilinear_approx!( - IOM.NoBilinearApproxConfig(), - container, PSY.InterconnectingConverter, - ipc_names, time_steps, - v_expr, i_var, - v_bounds, i_bounds, - "vi", - ) - - add_constraints!(container, ConverterLossConstraint, devices, model, network_model) - if get_attribute(model, "use_linear_loss") - add_constraints!( - container, - CurrentAbsoluteValueConstraint, - 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 end add_feedforward_constraints!(container, model, devices) diff --git a/test/test_device_hvdc.jl b/test/test_device_hvdc.jl index 2f55602..62652e8 100644 --- a/test/test_device_hvdc.jl +++ b/test/test_device_hvdc.jl @@ -157,8 +157,11 @@ end set_device_model!(template, TModelHVDCLine, DCLossyLine) set_device_model!( template, - DeviceModel(InterconnectingConverter, formulation); - attributes = Dict("use_linear_loss" => false), + DeviceModel( + InterconnectingConverter, + formulation; + attributes = Dict("use_linear_loss" => false), + ), ) set_hvdc_network_model!(template, VoltageDispatchHVDCNetworkModel) model = DecisionModel( @@ -183,3 +186,31 @@ end # 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, + ) + @test_logs( + (:warn, r"every InterconnectingConverter has a zero proportional loss"), + match_mode = :any, + build!(model; output_dir = mktempdir(; cleanup = true)), + ) +end From baa0288041dd6e3f5968b537192c862fb1fb7538 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli <8549957+acostarelli@users.noreply.github.com> Date: Mon, 18 May 2026 00:52:39 -0400 Subject: [PATCH 6/9] Add HVDCTwoTerminalVSC two-terminal HVDC formulation (#119) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * 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. * Address PR #119 review on HVDCTwoTerminalVSC - 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 * 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 * 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 * Address PR #119 Review 3 follow-up on HVDCTwoTerminalVSC 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 * formatting * Address PR #119 Review 4 on HVDCTwoTerminalVSC --------- Co-authored-by: Anthony Costarelli Co-authored-by: Claude Opus 4.7 --- Project.toml | 3 +- src/PowerOperationsModels.jl | 16 +- .../branch_constructor.jl | 158 +++++++ src/common_models/add_to_expression.jl | 39 ++ src/common_models/network_conditional.jl | 59 +++ src/common_models/quadratic_converter_loss.jl | 131 ++++++ src/core/constraints.jl | 36 ++ src/core/formulations.jl | 35 +- src/core/variables.jl | 49 ++- src/mt_hvdc_models/HVDCsystems.jl | 114 ++--- src/mt_hvdc_models/hvdcsystems_constructor.jl | 74 ++-- src/network_models/pm_translator.jl | 26 +- .../TwoTerminalDC_branches.jl | 400 ++++++++++++++++++ test/Project.toml | 2 +- test/runtests.jl | 6 + test/test_device_hvdc.jl | 206 +++++++-- 16 files changed, 1176 insertions(+), 178 deletions(-) create mode 100644 src/common_models/network_conditional.jl 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..6da666b 100644 --- a/src/PowerOperationsModels.jl +++ b/src/PowerOperationsModels.jl @@ -213,6 +213,8 @@ 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") +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. @@ -506,11 +508,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,12 +764,14 @@ export HVDCTwoTerminalLossless export HVDCTwoTerminalDispatch export HVDCTwoTerminalPiecewiseLoss export HVDCTwoTerminalLCC +export HVDCTwoTerminalVSCNLP +export HVDCTwoTerminalVSCLP # Converter Formulations export LosslessConverter export LinearLossConverter export AbstractQuadraticLossConverter -export Bin2QuadraticLossConverter +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 568b072..fcc96db 100644 --- a/src/ac_transmission_models/branch_constructor.jl +++ b/src/ac_transmission_models/branch_constructor.jl @@ -1670,3 +1670,161 @@ 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{HVDCTwoTerminalVSCNLP}) = IOM.NoQuadApproxConfig() +_quad_config(::Type{HVDCTwoTerminalVSCLP}) = + IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH) +_bilinear_config(::Type{HVDCTwoTerminalVSCNLP}) = IOM.NoBilinearApproxConfig() +_bilinear_config(::Type{HVDCTwoTerminalVSCLP}) = + IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)) + +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, + (HVDCReactivePowerFromVariable, HVDCReactivePowerToVariable), + ) + + # 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, + ) + 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 = 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) + + 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", + ) + + _register_pq_sq_expressions!( + container, devices, line_names, time_steps, device_model, + network_model, + ) + + if get_attribute(device_model, "use_linear_loss") + _add_abs_value_decomposition_constraints!( + container, devices, device_model, network_model, + DCLineCurrentFlowVariable, + ) + 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, + HVDCVSCApparentPowerLimitConstraint, + ) + + 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/add_to_expression.jl b/src/common_models/add_to_expression.jl index e8e9c7e..576f056 100644 --- a/src/common_models/add_to_expression.jl +++ b/src/common_models/add_to_expression.jl @@ -718,6 +718,45 @@ function add_to_expression!( 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 + +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 <: Union{HVDCReactivePowerFromVariable, HVDCReactivePowerToVariable}, + 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 = 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, + ) + end + end + return +end + """ PWL implementation to add FromTo branch variables to SystemBalanceExpressions """ diff --git a/src/common_models/network_conditional.jl b/src/common_models/network_conditional.jl new file mode 100644 index 0000000..e11a5c6 --- /dev/null +++ b/src/common_models/network_conditional.jl @@ -0,0 +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). 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. `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{<:AbstractPowerModel}, + var_types, +) 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, + ) + 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 on AC networks. +""" +function _maybe_add_reactive_power_constraints!( + container::OptimizationContainer, + devices, + model::DeviceModel{D, F}, + network_model::NetworkModel{<:AbstractPowerModel}, + constraint_type::Type{<:ConstraintType}, +) 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 new file mode 100644 index 0000000..9ad09d7 --- /dev/null +++ b/src/common_models/quadratic_converter_loss.jl @@ -0,0 +1,131 @@ +# Shared helpers for quadratic / two-term converter losses +# loss(I) = a * I^2 + b * |I| + c +# Used by multi-terminal InterconnectingConverter formulations +# (MILPQuadraticLossConverter, QuadraticLossConverter) and two-terminal +# HVDCTwoTerminalVSCNLP formulations. + +######################################### +######## Loss-curve introspection ####### +######################################### + +_get_quadratic_term(loss_fn::PSY.QuadraticCurve) = PSY.get_quadratic_term(loss_fn) +_get_quadratic_term(loss_fn) = 0.0 + +_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))) + +# `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 ######## +######################################### + +# `_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, i_pos_t, i_neg_t; + use_linear_loss::Bool, +) + 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) + JuMP.add_to_expression!(expr, b, i_neg_t) + end + return expr +end + +######################################### +####### 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, + ::DeviceModel{D, F}, + ::NetworkModel{<:AbstractPowerModel}, +) 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 + +# 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}, +) 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] + 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 ll_devices + name = PSY.get_name(d) + i_max = _linear_loss_i_max(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..1eb459d 100644 --- a/src/core/constraints.jl +++ b/src/core/constraints.jl @@ -419,6 +419,42 @@ 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 + +""" +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: + +- `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 + 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 + 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..32201ae 100644 --- a/src/core/formulations.jl +++ b/src/core/formulations.jl @@ -184,8 +184,33 @@ 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 HVDCTwoTerminalVSCNLP <: AbstractTwoTerminalVSCFormulation end + +""" +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 ############################### AC/DC Converter Formulations ##################################### abstract type AbstractConverterFormulation <: AbstractDeviceFormulation end @@ -206,10 +231,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 MILPQuadraticLossConverter <: 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 8a41d0d..bb3ce26 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 `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 `HVDCTwoTerminalVSCNLP` 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..37dd418 100644 --- a/src/mt_hvdc_models/HVDCsystems.jl +++ b/src/mt_hvdc_models/HVDCsystems.jl @@ -120,38 +120,42 @@ 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) +# 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{Bin2QuadraticLossConverter}, + ::Type{MILPQuadraticLossConverter}, ) - return Dict{String, Any}() + 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}, ) - return Dict{String, Any}( - "use_linear_loss" => false, - ) + return Dict{String, Any}("use_linear_loss" => false) end #! format: on @@ -540,20 +544,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}, @@ -570,10 +560,11 @@ 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, 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,66 +580,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 - loss_const[name, t] = JuMP.@constraint( - jump_model, - P_ac_var[name, t] == vi_expr[name, t] - loss - ) - end - end - return -end - -function add_constraints!( - container::OptimizationContainer, - ::Type{T}, - devices::W, - ::DeviceModel{U, V}, - ::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", - ) - - 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] + 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, ) - neg_ub_const[name, t] = JuMP.@constraint( + loss_const[name, t] = JuMP.@constraint( jump_model, - i_neg_var[name, t] <= i_max * (1 - i_dir_var[name, t]) + P_ac_var[name, t] == vi_expr[name, t] - loss, ) end end diff --git a/src/mt_hvdc_models/hvdcsystems_constructor.jl b/src/mt_hvdc_models/hvdcsystems_constructor.jl index c1ec321..e8d00a3 100644 --- a/src/mt_hvdc_models/hvdcsystems_constructor.jl +++ b/src/mt_hvdc_models/hvdcsystems_constructor.jl @@ -54,15 +54,17 @@ 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, ConverterPositiveCurrent, ll_devices, T) - add_variables!(container, ConverterNegativeCurrent, ll_devices, T) - add_variables!(container, ConverterCurrentDirection, ll_devices, T) - end + + # 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, + ) end add_to_expression!( @@ -77,19 +79,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 +106,11 @@ function _converter_vi_bounds(devices) return v_bounds, i_bounds end -_quad_config(::Type{Bin2QuadraticLossConverter}) = IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH) +_quad_config(::Type{MILPQuadraticLossConverter}) = + 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{MILPQuadraticLossConverter}) = + IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)) _bilinear_config(::Type{QuadraticLossConverter}) = IOM.NoBilinearApproxConfig() function construct_device!( @@ -148,20 +152,18 @@ function construct_device!( "vi", ) - 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 + # 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, + ) 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/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..d320236 100644 --- a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl +++ b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl @@ -1433,3 +1433,403 @@ 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. +# `_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}) = _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}) = _linear_loss_i_max(d) + +#! format: on + +####################### VSC PQ-approx registration ########################### + +# 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, + ::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) + 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, p_ft_bounds, "p_ft_sq", + ) + IOM._add_quadratic_approx!( + quad_cfg, container, PSY.TwoTerminalVSCLine, + 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, q_f_bounds, "q_f_sq", + ) + IOM._add_quadratic_approx!( + quad_cfg, container, PSY.TwoTerminalVSCLine, + line_names, time_steps, q_t, q_t_bounds, "q_t_sq", + ) + return +end + +# LP path: no disk constraint, so no p_sq/q_sq are needed. +_register_pq_sq_expressions!( + ::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, + ::DeviceModel{PSY.TwoTerminalVSCLine, HVDCTwoTerminalVSCNLP}, + ::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{U}, IS.FlattenIteratorWrapper{U}}, + ::DeviceModel{U, F}, + ::NetworkModel{<:AbstractPowerModel}, +) 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, U) + v_t = get_variable(container, HVDCToDCVoltage, U) + i_var = get_variable(container, DCLineCurrentFlowVariable, U) + + cons = add_constraints_container!( + container, HVDCCableOhmsLawConstraint, U, names, time_steps, + ) + + for d in devices + name = PSY.get_name(d) + g = PSY.get_g(d) + for t in time_steps + cons[name, t] = if iszero(g) + JuMP.@constraint(jump_model, i_var[name, t] == 0) + 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) +# 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{U}, IS.FlattenIteratorWrapper{U}}, + model::DeviceModel{U, F}, + ::NetworkModel{<:AbstractPowerModel}, +) 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, 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, U) + i_neg_var = get_variable(container, NegativeCurrent, U) + end + + cons_ft = add_constraints_container!( + container, HVDCVSCConverterPowerConstraint, U, names, time_steps; meta = "ft", + ) + cons_tf = add_constraints_container!( + container, HVDCVSCConverterPowerConstraint, U, 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_ft[name, t] + loss_ft, + ) + cons_tf[name, t] = JuMP.@constraint( + jump_model, + p_tf[name, t] == -vi_expr_tf[name, t] + loss_tf, + ) + end + end + return +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 +# `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, HVDCTwoTerminalVSCNLP}, + ::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_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, U, names, time_steps; + meta = "from", + ) + cons_t = add_constraints_container!( + container, HVDCVSCApparentPowerLimitConstraint, U, names, time_steps; + meta = "to", + ) + + for d in devices + name = PSY.get_name(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_sq[name, t] + q_f_sq[name, t] <= s_f2, + ) + cons_t[name, t] = JuMP.@constraint( + jump_model, p_tf_sq[name, t] + q_t_sq[name, t] <= s_t2, + ) + end + end + return +end + +# PQ capability — linear outer-approximation for the LP formulation. +# +# 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². +# +# 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{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, U) + p_tf = get_variable(container, FlowActivePowerToFromVariable, U) + q_f = get_variable(container, HVDCReactivePowerFromVariable, U) + q_t = get_variable(container, HVDCReactivePowerToVariable, U) + + 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") + 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 side_tags + cons[tag] = add_constraints_container!( + container, HVDCVSCApparentPowerLimitConstraint, U, + names, time_steps; meta = tag, + ) + end + + side_specs = ( + (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) + for spec in side_specs + rating = spec.rating_getter(d) + 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 * "_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) + 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 + return +end + +####################### VSC defaults ######################################### + +function get_default_time_series_names( + ::Type{PSY.TwoTerminalVSCLine}, + ::Type{<:AbstractTwoTerminalVSCFormulation}, +) + return Dict{Type{<:TimeSeriesParameter}, String}() +end + +# 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{HVDCTwoTerminalVSCNLP}, +) + return Dict{String, Any}("use_linear_loss" => false) +end + +# `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, "use_octagon" => true) +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..2632e26 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 (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, - Bin2QuadraticLossConverter, + 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 Bin2 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,41 +176,187 @@ end return model end - bin2_model = _build_and_solve(Bin2QuadraticLossConverter, HiGHS_optimizer) + milp_model = _build_and_solve(MILPQuadraticLossConverter, HiGHS_optimizer) nlp_model = _build_and_solve(QuadraticLossConverter, ipopt_optimizer) - bin2_obj = IOM.get_objective_value(OptimizationProblemOutputs(bin2_model)) + milp_obj = IOM.get_objective_value(OptimizationProblemOutputs(milp_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) + @test isapprox(milp_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() +############################################################################## +################ 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, 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, - ) - @test_logs( - (:warn, r"every InterconnectingConverter has a zero proportional loss"), - match_mode = :any, - build!(model; output_dir = mktempdir(; cleanup = true)), + set_device_model!(template, DeviceModel(TwoTerminalVSCLine, formulation)) + return DecisionModel( + template, sys; store_variable_names = true, optimizer = optimizer, ) end + +@testset "HVDC Two-Terminal VSC (NLP) on DCP" begin + 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(HVDCTwoTerminalVSCNLP, ACPPowerModel, 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 (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 — 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 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) + @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 + lp_obj = _solve(HVDCTwoTerminalVSCLP, HiGHS_optimizer) + nlp_obj = _solve(HVDCTwoTerminalVSCNLP, ipopt_optimizer) + @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 + # 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( + HVDCTwoTerminalVSCLP, 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( + HVDCTwoTerminalVSCNLP, 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 f0a590af6897b3660ef2213b5f7ea7e3d8f0b319 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Mon, 18 May 2026 16:38:40 -0400 Subject: [PATCH 7/9] Simplify HVDC converter loss modeling; consolidate HVDC tests MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Replace the abs-value decomposition (PositiveCurrent + NegativeCurrent + CurrentDirection binary, gated by use_linear_loss attribute) with a single linear AbsoluteValueCurrent variable bounded by abs_i >= ±i. The LP relaxation pins abs_i = |i| at the optimum because the loss term b*abs_i is minimized via the generation-cost objective, so no binary or complementarity constraint is needed. With the binary gone, NLP solvers (Ipopt) can also handle the linear loss term, so use_linear_loss collapses into the default codepath and is deleted on all four formulations. _quadratic_converter_loss_expr becomes two type-dispatched methods on typeof(i_sq_t): AffExpr for the MILP/PWL path, QuadExpr for the NLP/exact bilinear path. Each method is type-stable and returns the natural expression type for its formulation. Rename MILPQuadraticLossConverter -> QuadraticLossConverterMILP and QuadraticLossConverter -> QuadraticLossConverterNLP for NLP/MILP suffix consistency. Consolidate test/test_device_hvdc.jl: drop two redundant standalone MT build+solve tests, three redundant standalone VSC build+solve tests, and the octagon test (which silently passed on DCPPowerModel because HVDCVSCApparentPowerLimitConstraint is a no-op on active-power-only networks). Add one new test that verifies AbsoluteValueCurrent ≈ |ConverterCurrent| at the MILP optimum — direct evidence that the binary-free abs-value formulation is tight. Also addresses PR #103 review comments: - Drop misleading "NLP-only" wording from HVDCFromDCVoltage / HVDCToDCVoltage docstrings. - New header docstring on quadratic_converter_loss.jl explains both NLP and LP paths. - Mixed-loss-device KeyError bug (Copilot's comment) disappears by construction: AbsoluteValueCurrent is now registered for every device. - Silent-drop bug (use_linear_loss=false dropping nonzero b) disappears too: there's no toggle and b*abs_i is always modeled. Co-Authored-By: Claude Opus 4.7 --- src/PowerOperationsModels.jl | 8 +- .../branch_constructor.jl | 24 +-- src/common_models/quadratic_converter_loss.jl | 118 ++++------- src/core/formulations.jl | 10 +- src/core/variables.jl | 28 +-- src/mt_hvdc_models/HVDCsystems.jl | 44 +--- src/mt_hvdc_models/hvdcsystems_constructor.jl | 36 +--- .../TwoTerminalDC_branches.jl | 52 ++--- test/test_device_hvdc.jl | 191 +++++++----------- 9 files changed, 162 insertions(+), 349 deletions(-) diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl index 6da666b..0460cd7 100644 --- a/src/PowerOperationsModels.jl +++ b/src/PowerOperationsModels.jl @@ -508,9 +508,7 @@ export PostContingencyActivePowerReserveDeploymentVariable export DCVoltage export DCLineCurrent export ConverterCurrent -export PositiveCurrent -export NegativeCurrent -export CurrentDirection +export AbsoluteValueCurrent export HVDCFlowDirectionVariable export HVDCLosses export HVDCFromDCVoltage @@ -771,8 +769,8 @@ export HVDCTwoTerminalVSCLP export LosslessConverter export LinearLossConverter export AbstractQuadraticLossConverter -export MILPQuadraticLossConverter -export QuadraticLossConverter +export QuadraticLossConverterMILP +export QuadraticLossConverterNLP # DC Line Formulations export DCLosslessLine diff --git a/src/ac_transmission_models/branch_constructor.jl b/src/ac_transmission_models/branch_constructor.jl index fcc96db..f1ac408 100644 --- a/src/ac_transmission_models/branch_constructor.jl +++ b/src/ac_transmission_models/branch_constructor.jl @@ -1704,17 +1704,7 @@ function construct_device!( (HVDCReactivePowerFromVariable, HVDCReactivePowerToVariable), ) - # 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, - ) - end + _add_abs_value_variables!(container, devices, device_model, network_model) add_to_expression!( container, ActivePowerBalance, FlowActivePowerFromToVariable, @@ -1747,7 +1737,7 @@ function construct_device!( 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 + (min = -_vsc_cable_i_max(d), max = _vsc_cable_i_max(d)) for d in devices ] quad_cfg, bilin_cfg = _quad_config(F), _bilinear_config(F) @@ -1783,12 +1773,10 @@ function construct_device!( network_model, ) - if get_attribute(device_model, "use_linear_loss") - _add_abs_value_decomposition_constraints!( - container, devices, device_model, network_model, - DCLineCurrentFlowVariable, - ) - end + _add_abs_value_constraints!( + container, devices, device_model, network_model, + DCLineCurrentFlowVariable, + ) add_constraints!( container, HVDCCableOhmsLawConstraint, devices, device_model, network_model, diff --git a/src/common_models/quadratic_converter_loss.jl b/src/common_models/quadratic_converter_loss.jl index 9ad09d7..7e09294 100644 --- a/src/common_models/quadratic_converter_loss.jl +++ b/src/common_models/quadratic_converter_loss.jl @@ -1,8 +1,14 @@ # Shared helpers for quadratic / two-term converter losses # loss(I) = a * I^2 + b * |I| + c # Used by multi-terminal InterconnectingConverter formulations -# (MILPQuadraticLossConverter, QuadraticLossConverter) and two-terminal -# HVDCTwoTerminalVSCNLP formulations. +# (QuadraticLossConverterMILP, QuadraticLossConverterNLP) and two-terminal +# TwoTerminalVSCLine formulations (HVDCTwoTerminalVSCLP, HVDCTwoTerminalVSCNLP). +# +# `|I|` is represented by an LP surrogate: a single non-negative variable +# `AbsoluteValueCurrent` bounded below by both `i` and `-i`. The optimum +# pins it to `|i|` because the loss term `b · abs_i` is being minimized +# via the generation-cost objective; no binary or complementarity constraint +# is required. ######################################### ######## Loss-curve introspection ####### @@ -11,119 +17,83 @@ _get_quadratic_term(loss_fn::PSY.QuadraticCurve) = PSY.get_quadratic_term(loss_fn) _get_quadratic_term(loss_fn) = 0.0 -_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))) - -# `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 ######## ######################################### -# `_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) +# Two methods dispatched on `typeof(i_sq_t)` so each call site is fully +# type-stable. `i_sq_t` is `AffExpr` on the MILP/PWL path (SOS2 surrogate +# under `SolverSOS2QuadConfig`) and `QuadExpr` on the NLP/bilinear path +# (under `NoQuadApproxConfig`). MILP returns AffExpr → HiGHS gets a linear +# constraint; NLP returns QuadExpr → Ipopt gets the genuine quadratic +# term in the right field. function _quadratic_converter_loss_expr( a::Float64, b::Float64, c::Float64, - i_sq_t, i_pos_t, i_neg_t; - use_linear_loss::Bool, + i_sq_t::JuMP.AffExpr, + abs_i_t, ) - expr = iszero(a) ? JuMP.AffExpr(c) : _loss_seed(c, i_sq_t) + 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 + iszero(b) || JuMP.add_to_expression!(expr, b, abs_i_t) + return expr +end + +function _quadratic_converter_loss_expr( + a::Float64, b::Float64, c::Float64, + i_sq_t::JuMP.QuadExpr, + abs_i_t, +) + expr = JuMP.QuadExpr(JuMP.AffExpr(c)) + iszero(a) || JuMP.add_to_expression!(expr, a, i_sq_t) + iszero(b) || JuMP.add_to_expression!(expr, b, abs_i_t) return expr end ######################################### -####### Abs-value decomposition ######### +######## Absolute-value surrogate ####### ######################################### -# 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!( +function _add_abs_value_variables!( container::OptimizationContainer, devices, ::DeviceModel{D, F}, ::NetworkModel{<:AbstractPowerModel}, ) 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) + add_variables!(container, AbsoluteValueCurrent, devices, F) 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!( +function _add_abs_value_constraints!( container::OptimizationContainer, devices, ::DeviceModel{D, F}, ::NetworkModel{<:AbstractPowerModel}, parent_var_type::Type{<:VariableType}, ) 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] + 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_i_var = get_variable(container, AbsoluteValueCurrent, D) - abs_val_const = add_constraints_container!( - container, CurrentAbsoluteValueConstraint, D, names, time_steps, - ) - pos_ub_const = add_constraints_container!( + lower_const = add_constraints_container!( container, CurrentAbsoluteValueConstraint, D, names, time_steps; - meta = "pos_ub", + meta = "ge_pos", ) - neg_ub_const = add_constraints_container!( + upper_const = add_constraints_container!( container, CurrentAbsoluteValueConstraint, D, names, time_steps; - meta = "neg_ub", + meta = "ge_neg", ) - for d in ll_devices + for d in devices name = PSY.get_name(d) - i_max = _linear_loss_i_max(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], + lower_const[name, t] = JuMP.@constraint( + jump_model, abs_i_var[name, t] >= i_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]), + upper_const[name, t] = JuMP.@constraint( + jump_model, abs_i_var[name, t] >= -i_var[name, t], ) end end diff --git a/src/core/formulations.jl b/src/core/formulations.jl index 32201ae..97bda90 100644 --- a/src/core/formulations.jl +++ b/src/core/formulations.jl @@ -234,19 +234,13 @@ 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 MILPQuadraticLossConverter <: AbstractQuadraticLossConverter end +struct QuadraticLossConverterMILP <: AbstractQuadraticLossConverter end """ Quadratic Loss InterconnectingConverter using exact bilinear (v·i) and quadratic (i²) products. Requires an NLP-capable solver (e.g., Ipopt). - -Attributes (set on `DeviceModel`): -- `use_linear_loss`: Include the linear loss term `b·(i⁺ + i⁻)` (default false). The - linear term requires an absolute-value reformulation that introduces a binary - variable, so this is only solvable by MINLP solvers (e.g., Gurobi); pure NLP - solvers like Ipopt cannot solve the model with this option enabled. """ -struct QuadraticLossConverter <: AbstractQuadraticLossConverter end +struct QuadraticLossConverterNLP <: AbstractQuadraticLossConverter end ############################## HVDC Lines Formulations ################################## abstract type AbstractDCLineFormulation <: AbstractBranchFormulation end diff --git a/src/core/variables.jl b/src/core/variables.jl index bb3ce26..e8b9c2d 100644 --- a/src/core/variables.jl +++ b/src/core/variables.jl @@ -157,25 +157,13 @@ Docs abbreviation: ``i_c^{dc}`` struct ConverterCurrent <: VariableType end """ -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}`` +Non-negative LP surrogate for the absolute value of a DC current variable +(`abs_i ≥ i`, `abs_i ≥ -i`, `abs_i ≥ 0`; tight at optimum because the loss +term `b · abs_i` is minimized via the generation-cost objective). Used by +both two-terminal HVDC links (cable current) and interconnecting converters. +Docs abbreviation: ``|i|^{dc}`` """ -struct PositiveCurrent <: VariableType end - -""" -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 NegativeCurrent <: VariableType end - -""" -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 CurrentDirection <: VariableType end +struct AbsoluteValueCurrent <: VariableType end ######################################################### ######################################################### @@ -389,14 +377,14 @@ struct HVDCPiecewiseBinaryLossVariable <: SparseVariableType end """ DC-side voltage at the from-terminal of a two-terminal HVDC link. -Used by `HVDCTwoTerminalVSCNLP` formulations. +Used by both `HVDCTwoTerminalVSCNLP` and `HVDCTwoTerminalVSCLP` 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 `HVDCTwoTerminalVSCNLP` formulations. +Used by both `HVDCTwoTerminalVSCNLP` and `HVDCTwoTerminalVSCLP` formulations. Docs abbreviation: ``v_t^{dc}`` """ struct HVDCToDCVoltage <: VariableType end diff --git a/src/mt_hvdc_models/HVDCsystems.jl b/src/mt_hvdc_models/HVDCsystems.jl index 37dd418..07267a1 100644 --- a/src/mt_hvdc_models/HVDCsystems.jl +++ b/src/mt_hvdc_models/HVDCsystems.jl @@ -120,43 +120,18 @@ end ## Binaries ### get_variable_binary(::Type{ConverterCurrent}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = false -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 +get_variable_binary(::Type{AbsoluteValueCurrent}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = false ### 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{PositiveCurrent}, d::PSY.InterconnectingConverter,::Type{<:AbstractConverterFormulation}) = 0.0 -get_variable_lower_bound(::Type{NegativeCurrent}, d::PSY.InterconnectingConverter,::Type{<:AbstractConverterFormulation}) = 0.0 +get_variable_lower_bound(::Type{AbsoluteValueCurrent}, 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{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) - - -# 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{MILPQuadraticLossConverter}, -) - 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}, -) - return Dict{String, Any}("use_linear_loss" => false) -end +get_variable_upper_bound(::Type{AbsoluteValueCurrent}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = PSY.get_max_dc_current(d) #! format: on @@ -559,13 +534,7 @@ function add_constraints!( P_ac_var = get_variable(container, ActivePowerVariable, U) vi_expr = get_expression(container, IOM.BilinearProductExpression, U, "vi") 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, U) - i_neg_var = get_variable(container, NegativeCurrent, U) - end + abs_i_var = get_variable(container, AbsoluteValueCurrent, U) ipc_names = [PSY.get_name(d) for d in devices] loss_const = add_constraints_container!( @@ -580,11 +549,8 @@ function add_constraints!( b = PSY.get_proportional_term(loss_function) c = PSY.get_constant_term(loss_function) 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 = _quadratic_converter_loss_expr( - a, b, c, i_sq_expr[name, t], i_pos_t, i_neg_t; - use_linear_loss = use_linear_loss, + a, b, c, i_sq_expr[name, t], abs_i_var[name, t], ) loss_const[name, t] = JuMP.@constraint( jump_model, diff --git a/src/mt_hvdc_models/hvdcsystems_constructor.jl b/src/mt_hvdc_models/hvdcsystems_constructor.jl index e8d00a3..6f792d1 100644 --- a/src/mt_hvdc_models/hvdcsystems_constructor.jl +++ b/src/mt_hvdc_models/hvdcsystems_constructor.jl @@ -54,18 +54,7 @@ function construct_device!( devices = get_available_components(model, sys) add_variables!(container, ActivePowerVariable, devices, T) add_variables!(container, ConverterCurrent, devices, T) - - # 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, - ) - end + _add_abs_value_variables!(container, devices, model, network_model) add_to_expression!( container, ActivePowerBalance, ActivePowerVariable, @@ -106,12 +95,12 @@ function _converter_vi_bounds(devices) return v_bounds, i_bounds end -_quad_config(::Type{MILPQuadraticLossConverter}) = +_quad_config(::Type{QuadraticLossConverterMILP}) = IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH) -_quad_config(::Type{QuadraticLossConverter}) = IOM.NoQuadApproxConfig() -_bilinear_config(::Type{MILPQuadraticLossConverter}) = +_quad_config(::Type{QuadraticLossConverterNLP}) = IOM.NoQuadApproxConfig() +_bilinear_config(::Type{QuadraticLossConverterMILP}) = IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)) -_bilinear_config(::Type{QuadraticLossConverter}) = IOM.NoBilinearApproxConfig() +_bilinear_config(::Type{QuadraticLossConverterNLP}) = IOM.NoBilinearApproxConfig() function construct_device!( container::OptimizationContainer, @@ -152,15 +141,12 @@ function construct_device!( "vi", ) - # 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, - ) - end + # AbsoluteValueCurrent variable was added in the ArgumentConstructStage; + # only the `abs_i ≥ ±i` constraints need the JuMP model now. + # ConverterLossConstraint reads `abs_i`, so its constraints must come first. + _add_abs_value_constraints!( + container, devices, model, network_model, ConverterCurrent, + ) add_constraints!(container, ConverterLossConstraint, devices, model, network_model) diff --git a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl index d320236..3d1b17d 100644 --- a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl +++ b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl @@ -1446,9 +1446,7 @@ get_variable_binary(::Type{HVDCFromDCVoltage}, ::Type{PSY.TwoTerminalVSCLine}, : 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{AbsoluteValueCurrent}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false get_variable_binary(::Type{FlowActivePowerFromToVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false get_variable_binary(::Type{FlowActivePowerToFromVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false @@ -1478,16 +1476,14 @@ get_variable_lower_bound(::Type{HVDCToDCVoltage}, d::PSY.TwoTerminalVSCLine, ::T 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. -# `_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) +_vsc_cable_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_cable_i_max(d) +get_variable_upper_bound(::Type{DCLineCurrentFlowVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = _vsc_cable_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}) = _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}) = _linear_loss_i_max(d) +# AbsoluteValueCurrent: 0 ≤ abs_i ≤ I_max (LP surrogate for |i|) +get_variable_lower_bound(::Type{AbsoluteValueCurrent}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = 0.0 +get_variable_upper_bound(::Type{AbsoluteValueCurrent}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = _vsc_cable_i_max(d) #! format: on @@ -1610,13 +1606,7 @@ function add_constraints!( 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, U) - i_neg_var = get_variable(container, NegativeCurrent, U) - end + abs_i_var = get_variable(container, AbsoluteValueCurrent, U) cons_ft = add_constraints_container!( container, HVDCVSCConverterPowerConstraint, U, names, time_steps; meta = "ft", @@ -1636,15 +1626,12 @@ function add_constraints!( 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 + abs_i_t = abs_i_var[name, t] 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, + a_f, b_f, c_f, i_sq_expr[name, t], abs_i_t, ) 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, + a_t, b_t, c_t, i_sq_expr[name, t], abs_i_t, ) cons_ft[name, t] = JuMP.@constraint( jump_model, @@ -1809,19 +1796,6 @@ function get_default_time_series_names( return Dict{Type{<:TimeSeriesParameter}, String}() end -# 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{HVDCTwoTerminalVSCNLP}, -) - return Dict{String, Any}("use_linear_loss" => false) -end - -# `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 @@ -1831,5 +1805,5 @@ function get_default_attributes( ::Type{PSY.TwoTerminalVSCLine}, ::Type{HVDCTwoTerminalVSCLP}, ) - return Dict{String, Any}("use_linear_loss" => true, "use_octagon" => true) + return Dict{String, Any}("use_octagon" => true) end diff --git a/test/test_device_hvdc.jl b/test/test_device_hvdc.jl index 2632e26..0b5e9c2 100644 --- a/test/test_device_hvdc.jl +++ b/test/test_device_hvdc.jl @@ -97,59 +97,10 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED end -@testset "HVDC System with Losses Network (MILPQuadraticLossConverter)" begin - sys = _generate_test_hvdc_sys() - 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) - ipc_model = DeviceModel( - InterconnectingConverter, - MILPQuadraticLossConverter, - ) - set_device_model!(template, ipc_model) - set_hvdc_network_model!(template, VoltageDispatchHVDCNetworkModel) - 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 -end - -@testset "HVDC System with Losses Network (QuadraticLossConverter NLP)" begin - sys = _generate_test_hvdc_sys() - 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) - ipc_model = DeviceModel( - InterconnectingConverter, - QuadraticLossConverter, - ) - set_device_model!(template, ipc_model) - set_hvdc_network_model!(template, VoltageDispatchHVDCNetworkModel) - model = - DecisionModel( - template, - sys; - store_variable_names = true, - optimizer = ipopt_optimizer, - ) - @test build!(model; output_dir = mktempdir(; cleanup = true)) == - IOM.ModelBuildStatus.BUILT - @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED -end - @testset "HVDC MILP vs NLP QuadraticLossConverter agreement" begin - function _build_and_solve(formulation, optimizer) - sys = _generate_test_hvdc_sys() + # Build both formulations on the same system; compare objective values + # (Rodrigo's "same order of magnitude" ask from PR #103). + function _build_and_solve(sys, formulation, optimizer) template = OperationsProblemTemplate() set_device_model!(template, ThermalStandard, ThermalDispatchNoMin) set_device_model!(template, PowerLoad, StaticPowerLoad) @@ -157,18 +108,12 @@ end set_device_model!(template, TModelHVDCLine, DCLossyLine) set_device_model!( template, - DeviceModel( - InterconnectingConverter, - formulation; - attributes = Dict("use_linear_loss" => false), - ), + DeviceModel(InterconnectingConverter, formulation), ) set_hvdc_network_model!(template, VoltageDispatchHVDCNetworkModel) model = DecisionModel( - template, - sys; - store_variable_names = true, - optimizer = optimizer, + template, sys; + store_variable_names = true, optimizer = optimizer, ) @test build!(model; output_dir = mktempdir(; cleanup = true)) == IOM.ModelBuildStatus.BUILT @@ -176,15 +121,61 @@ end return model end - milp_model = _build_and_solve(MILPQuadraticLossConverter, HiGHS_optimizer) - nlp_model = _build_and_solve(QuadraticLossConverter, ipopt_optimizer) + sys = _generate_test_hvdc_sys() + milp_model = _build_and_solve(sys, QuadraticLossConverterMILP, HiGHS_optimizer) + nlp_model = _build_and_solve(sys, QuadraticLossConverterNLP, ipopt_optimizer) + # Objective is the right level of strictness for "same order of magnitude" + # (Rodrigo's ask on the PR #103 review). Per-converter or system-total + # current/power comparisons fail unpredictably on this fixture because the + # MT-HVDC fleet carries essentially no power either way (the loss term + # drives it toward zero on both sides), so the MILP's SOS2 PWL surrogate + # vs the NLP's exact bilinear leave residuals at very different + # magnitudes — both still tiny in absolute terms, just not within a + # rtol-comparable factor of each other. milp_obj = IOM.get_objective_value(OptimizationProblemOutputs(milp_model)) nlp_obj = IOM.get_objective_value(OptimizationProblemOutputs(nlp_model)) - @test isapprox(milp_obj, nlp_obj; rtol = 0.05) end +@testset "HVDC AbsoluteValueCurrent matches |ConverterCurrent| at MILP optimum" begin + # Direct evidence that the binary-free LP abs-value formulation is tight: + # the loss objective drives abs_i down to exactly |i| at the optimum. + sys = _generate_test_hvdc_sys() + 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, QuadraticLossConverterMILP), + ) + set_hvdc_network_model!(template, VoltageDispatchHVDCNetworkModel) + 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 + + container = IOM.get_optimization_container(model) + i_vals = + JuMP.value.( + IOM.get_variable(container, ConverterCurrent, InterconnectingConverter).data, + ) + abs_i_vals = + JuMP.value.( + IOM.get_variable( + container, + AbsoluteValueCurrent, + InterconnectingConverter, + ).data, + ) + @test isapprox(abs_i_vals, abs.(i_vals); atol = 1e-6) +end + ############################################################################## ################ Two-Terminal VSC HVDC tests ################################# ############################################################################## @@ -252,30 +243,21 @@ function _build_vsc_model(formulation, network, optimizer; sys = _generate_test_ ) end -@testset "HVDC Two-Terminal VSC (NLP) on DCP" begin - 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(HVDCTwoTerminalVSCNLP, ACPPowerModel, 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 (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 — 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. +# Standalone build+solve smoke tests for each (formulation, network) combo +# are covered by the agreement / property tests further down: +# - HVDCTwoTerminalVSCNLP on DCP → "HVDC VSC LP vs NLP objective agreement" +# - HVDCTwoTerminalVSCLP on DCP → same agreement test + cable-resistance test +# - HVDCTwoTerminalVSCNLP on AC → "HVDC VSC: tighter PQ rating raises cost on AC" +# HVDCTwoTerminalVSCLP on ACPPowerModel is omitted: HiGHS can't solve the +# ACP network's trig (cos/sin) branch ohms-law constraints, and no MINLP +# solver with trigonometric support is wired into the test deps. +# +# TODO: Re-add an `octagon vs box-only` LP property test once an MINLP solver +# with trig support is available. The previous version of that test ran on +# `DCPPowerModel`, which never adds `HVDCVSCApparentPowerLimitConstraint` +# (the constraint is gated by `_maybe_add_reactive_power_constraints!`, a +# no-op on `AbstractActivePowerModel`), so it asserted the same model against +# itself. @testset "HVDC VSC LP vs NLP objective agreement" begin # On a DC network the PQ disk constraint is inactive (no reactive @@ -295,39 +277,6 @@ end @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 # Smaller g => larger R = 1/g => more losses => optimum should not improve. function _solve_with_g(g_value) From 79200d9edf018ab84ae9bed1565ebb9245be9e73 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Tue, 19 May 2026 10:46:48 -0400 Subject: [PATCH 8/9] Tidy quadratic converter loss helpers MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Rename AbsoluteValueCurrent → CurrentAbsoluteValueVariable to match its sister CurrentAbsoluteValueConstraint. Inline the one-line _add_abs_value_variables! wrapper at its two call sites. Add a docstring to _quadratic_converter_loss_expr explaining why the iszero guards drop zero coefficient terms, and annotate abs_i_t::JuMP.VariableRef so the methods don't infer Any. Co-Authored-By: Claude Opus 4.7 --- src/PowerOperationsModels.jl | 2 +- .../branch_constructor.jl | 2 +- src/common_models/quadratic_converter_loss.jl | 40 ++++++++----------- src/core/variables.jl | 2 +- src/mt_hvdc_models/HVDCsystems.jl | 8 ++-- src/mt_hvdc_models/hvdcsystems_constructor.jl | 4 +- .../TwoTerminalDC_branches.jl | 10 ++--- test/test_device_hvdc.jl | 4 +- 8 files changed, 33 insertions(+), 39 deletions(-) diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl index 0460cd7..d0bdfaf 100644 --- a/src/PowerOperationsModels.jl +++ b/src/PowerOperationsModels.jl @@ -508,7 +508,7 @@ export PostContingencyActivePowerReserveDeploymentVariable export DCVoltage export DCLineCurrent export ConverterCurrent -export AbsoluteValueCurrent +export CurrentAbsoluteValueVariable export HVDCFlowDirectionVariable export HVDCLosses export HVDCFromDCVoltage diff --git a/src/ac_transmission_models/branch_constructor.jl b/src/ac_transmission_models/branch_constructor.jl index f1ac408..0e93c03 100644 --- a/src/ac_transmission_models/branch_constructor.jl +++ b/src/ac_transmission_models/branch_constructor.jl @@ -1704,7 +1704,7 @@ function construct_device!( (HVDCReactivePowerFromVariable, HVDCReactivePowerToVariable), ) - _add_abs_value_variables!(container, devices, device_model, network_model) + add_variables!(container, CurrentAbsoluteValueVariable, devices, F) add_to_expression!( container, ActivePowerBalance, FlowActivePowerFromToVariable, diff --git a/src/common_models/quadratic_converter_loss.jl b/src/common_models/quadratic_converter_loss.jl index 7e09294..06dffab 100644 --- a/src/common_models/quadratic_converter_loss.jl +++ b/src/common_models/quadratic_converter_loss.jl @@ -5,10 +5,10 @@ # TwoTerminalVSCLine formulations (HVDCTwoTerminalVSCLP, HVDCTwoTerminalVSCNLP). # # `|I|` is represented by an LP surrogate: a single non-negative variable -# `AbsoluteValueCurrent` bounded below by both `i` and `-i`. The optimum -# pins it to `|i|` because the loss term `b · abs_i` is being minimized -# via the generation-cost objective; no binary or complementarity constraint -# is required. +# `CurrentAbsoluteValueVariable` bounded below by both `i` and `-i`. The +# optimum pins it to `|i|` because the loss term `b · abs_i` is being +# minimized via the generation-cost objective; no binary or complementarity +# constraint is required. ######################################### ######## Loss-curve introspection ####### @@ -21,17 +21,21 @@ _get_quadratic_term(loss_fn) = 0.0 ######## Loss expression builder ######## ######################################### -# Two methods dispatched on `typeof(i_sq_t)` so each call site is fully -# type-stable. `i_sq_t` is `AffExpr` on the MILP/PWL path (SOS2 surrogate -# under `SolverSOS2QuadConfig`) and `QuadExpr` on the NLP/bilinear path -# (under `NoQuadApproxConfig`). MILP returns AffExpr → HiGHS gets a linear -# constraint; NLP returns QuadExpr → Ipopt gets the genuine quadratic -# term in the right field. +""" + _quadratic_converter_loss_expr(a, b, c, i_sq_t, abs_i_t) +Build the per-timestep converter loss expression `a·I² + b·|I| + c`. + +The `iszero(a)` / `iszero(b)` guards skip the corresponding +`add_to_expression!` call when the coefficient is exactly zero — this +avoids dragging a `0·I²` or `0·|I|` term into the JuMP expression, which +would otherwise allocate and force the wrong expression type (e.g., a +`QuadExpr` with a zero quadratic part on the MILP path). +""" function _quadratic_converter_loss_expr( a::Float64, b::Float64, c::Float64, i_sq_t::JuMP.AffExpr, - abs_i_t, + abs_i_t::JuMP.VariableRef, ) expr = JuMP.AffExpr(c) iszero(a) || JuMP.add_to_expression!(expr, a, i_sq_t) @@ -42,7 +46,7 @@ end function _quadratic_converter_loss_expr( a::Float64, b::Float64, c::Float64, i_sq_t::JuMP.QuadExpr, - abs_i_t, + abs_i_t::JuMP.VariableRef, ) expr = JuMP.QuadExpr(JuMP.AffExpr(c)) iszero(a) || JuMP.add_to_expression!(expr, a, i_sq_t) @@ -54,16 +58,6 @@ end ######## Absolute-value surrogate ####### ######################################### -function _add_abs_value_variables!( - container::OptimizationContainer, - devices, - ::DeviceModel{D, F}, - ::NetworkModel{<:AbstractPowerModel}, -) where {D <: PSY.Device, F} - add_variables!(container, AbsoluteValueCurrent, devices, F) - return -end - function _add_abs_value_constraints!( container::OptimizationContainer, devices, @@ -75,7 +69,7 @@ function _add_abs_value_constraints!( names = [PSY.get_name(d) for d in devices] jump_model = get_jump_model(container) i_var = get_variable(container, parent_var_type, D) - abs_i_var = get_variable(container, AbsoluteValueCurrent, D) + abs_i_var = get_variable(container, CurrentAbsoluteValueVariable, D) lower_const = add_constraints_container!( container, CurrentAbsoluteValueConstraint, D, names, time_steps; diff --git a/src/core/variables.jl b/src/core/variables.jl index e8b9c2d..1667837 100644 --- a/src/core/variables.jl +++ b/src/core/variables.jl @@ -163,7 +163,7 @@ term `b · abs_i` is minimized via the generation-cost objective). Used by both two-terminal HVDC links (cable current) and interconnecting converters. Docs abbreviation: ``|i|^{dc}`` """ -struct AbsoluteValueCurrent <: VariableType end +struct CurrentAbsoluteValueVariable <: VariableType end ######################################################### ######################################################### diff --git a/src/mt_hvdc_models/HVDCsystems.jl b/src/mt_hvdc_models/HVDCsystems.jl index 07267a1..8cf9cc9 100644 --- a/src/mt_hvdc_models/HVDCsystems.jl +++ b/src/mt_hvdc_models/HVDCsystems.jl @@ -120,18 +120,18 @@ end ## Binaries ### get_variable_binary(::Type{ConverterCurrent}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = false -get_variable_binary(::Type{AbsoluteValueCurrent}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = false +get_variable_binary(::Type{CurrentAbsoluteValueVariable}, ::Type{PSY.InterconnectingConverter}, ::Type{<:AbstractConverterFormulation}) = false ### 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{AbsoluteValueCurrent}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = 0.0 +get_variable_lower_bound(::Type{CurrentAbsoluteValueVariable}, 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{AbsoluteValueCurrent}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = PSY.get_max_dc_current(d) +get_variable_upper_bound(::Type{CurrentAbsoluteValueVariable}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = PSY.get_max_dc_current(d) #! format: on @@ -534,7 +534,7 @@ function add_constraints!( P_ac_var = get_variable(container, ActivePowerVariable, U) vi_expr = get_expression(container, IOM.BilinearProductExpression, U, "vi") i_sq_expr = get_expression(container, IOM.QuadraticExpression, U, "i_sq") - abs_i_var = get_variable(container, AbsoluteValueCurrent, U) + abs_i_var = get_variable(container, CurrentAbsoluteValueVariable, U) ipc_names = [PSY.get_name(d) for d in devices] loss_const = add_constraints_container!( diff --git a/src/mt_hvdc_models/hvdcsystems_constructor.jl b/src/mt_hvdc_models/hvdcsystems_constructor.jl index 6f792d1..743be6a 100644 --- a/src/mt_hvdc_models/hvdcsystems_constructor.jl +++ b/src/mt_hvdc_models/hvdcsystems_constructor.jl @@ -54,7 +54,7 @@ function construct_device!( devices = get_available_components(model, sys) add_variables!(container, ActivePowerVariable, devices, T) add_variables!(container, ConverterCurrent, devices, T) - _add_abs_value_variables!(container, devices, model, network_model) + add_variables!(container, CurrentAbsoluteValueVariable, devices, T) add_to_expression!( container, ActivePowerBalance, ActivePowerVariable, @@ -141,7 +141,7 @@ function construct_device!( "vi", ) - # AbsoluteValueCurrent variable was added in the ArgumentConstructStage; + # CurrentAbsoluteValueVariable was added in the ArgumentConstructStage; # only the `abs_i ≥ ±i` constraints need the JuMP model now. # ConverterLossConstraint reads `abs_i`, so its constraints must come first. _add_abs_value_constraints!( diff --git a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl index 3d1b17d..b71ed14 100644 --- a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl +++ b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl @@ -1446,7 +1446,7 @@ get_variable_binary(::Type{HVDCFromDCVoltage}, ::Type{PSY.TwoTerminalVSCLine}, : 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{AbsoluteValueCurrent}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false +get_variable_binary(::Type{CurrentAbsoluteValueVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false get_variable_binary(::Type{FlowActivePowerFromToVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false get_variable_binary(::Type{FlowActivePowerToFromVariable}, ::Type{PSY.TwoTerminalVSCLine}, ::Type{<:AbstractTwoTerminalVSCFormulation}) = false @@ -1481,9 +1481,9 @@ _vsc_cable_i_max(d::PSY.TwoTerminalVSCLine) = get_variable_lower_bound(::Type{DCLineCurrentFlowVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = -_vsc_cable_i_max(d) get_variable_upper_bound(::Type{DCLineCurrentFlowVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = _vsc_cable_i_max(d) -# AbsoluteValueCurrent: 0 ≤ abs_i ≤ I_max (LP surrogate for |i|) -get_variable_lower_bound(::Type{AbsoluteValueCurrent}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = 0.0 -get_variable_upper_bound(::Type{AbsoluteValueCurrent}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = _vsc_cable_i_max(d) +# CurrentAbsoluteValueVariable: 0 ≤ abs_i ≤ I_max (LP surrogate for |i|) +get_variable_lower_bound(::Type{CurrentAbsoluteValueVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = 0.0 +get_variable_upper_bound(::Type{CurrentAbsoluteValueVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = _vsc_cable_i_max(d) #! format: on @@ -1606,7 +1606,7 @@ function add_constraints!( vi_expr_tf = get_expression(container, IOM.BilinearProductExpression, U, "vi_tf") i_sq_expr = get_expression(container, IOM.QuadraticExpression, U, "i_sq") - abs_i_var = get_variable(container, AbsoluteValueCurrent, U) + abs_i_var = get_variable(container, CurrentAbsoluteValueVariable, U) cons_ft = add_constraints_container!( container, HVDCVSCConverterPowerConstraint, U, names, time_steps; meta = "ft", diff --git a/test/test_device_hvdc.jl b/test/test_device_hvdc.jl index 0b5e9c2..30acb05 100644 --- a/test/test_device_hvdc.jl +++ b/test/test_device_hvdc.jl @@ -138,7 +138,7 @@ end @test isapprox(milp_obj, nlp_obj; rtol = 0.05) end -@testset "HVDC AbsoluteValueCurrent matches |ConverterCurrent| at MILP optimum" begin +@testset "HVDC CurrentAbsoluteValueVariable matches |ConverterCurrent| at MILP optimum" begin # Direct evidence that the binary-free LP abs-value formulation is tight: # the loss objective drives abs_i down to exactly |i| at the optimum. sys = _generate_test_hvdc_sys() @@ -169,7 +169,7 @@ end JuMP.value.( IOM.get_variable( container, - AbsoluteValueCurrent, + CurrentAbsoluteValueVariable, InterconnectingConverter, ).data, ) From 046c508858ea4232748adf1d2ee3ef69e42b2985 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Tue, 19 May 2026 10:53:59 -0400 Subject: [PATCH 9/9] use jump utils; fix docstring --- src/common_models/quadratic_converter_loss.jl | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/common_models/quadratic_converter_loss.jl b/src/common_models/quadratic_converter_loss.jl index 06dffab..06ffd6c 100644 --- a/src/common_models/quadratic_converter_loss.jl +++ b/src/common_models/quadratic_converter_loss.jl @@ -26,11 +26,9 @@ _get_quadratic_term(loss_fn) = 0.0 Build the per-timestep converter loss expression `a·I² + b·|I| + c`. -The `iszero(a)` / `iszero(b)` guards skip the corresponding -`add_to_expression!` call when the coefficient is exactly zero — this -avoids dragging a `0·I²` or `0·|I|` term into the JuMP expression, which -would otherwise allocate and force the wrong expression type (e.g., a -`QuadExpr` with a zero quadratic part on the MILP path). +In MILP formulations the i_sq_t is still an AffExpr. + +The `iszero` guards avoid adding 0s to the JuMP expression which might slightly hurt the solver. """ function _quadratic_converter_loss_expr( a::Float64, b::Float64, c::Float64, @@ -38,8 +36,8 @@ function _quadratic_converter_loss_expr( abs_i_t::JuMP.VariableRef, ) expr = JuMP.AffExpr(c) - iszero(a) || JuMP.add_to_expression!(expr, a, i_sq_t) - iszero(b) || JuMP.add_to_expression!(expr, b, abs_i_t) + iszero(a) || add_proportional_to_jump_expression!(expr, i_sq_t, a) + iszero(b) || add_proportional_to_jump_expression!(expr, abs_i_t, b) return expr end @@ -49,8 +47,8 @@ function _quadratic_converter_loss_expr( abs_i_t::JuMP.VariableRef, ) expr = JuMP.QuadExpr(JuMP.AffExpr(c)) - iszero(a) || JuMP.add_to_expression!(expr, a, i_sq_t) - iszero(b) || JuMP.add_to_expression!(expr, b, abs_i_t) + iszero(a) || add_proportional_to_jump_expression!(expr, i_sq_t, a) + iszero(b) || add_proportional_to_jump_expression!(expr, abs_i_t, b) return expr end