Add TwoTerminalVSC formulation

Project.toml

@@ -23,18 +23,17 @@ 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 = "ac/hvdc-vsc", 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"
 InfrastructureSystems = "3"
 InteractiveUtils = "1.11.0"
 JuMP = "^1.28"

src/PowerOperationsModels.jl

@@ -213,6 +213,7 @@ include("common_models/add_to_expression.jl")
 include("common_models/add_parameters.jl")
 include("common_models/make_system_expressions.jl")
 include("common_models/reserve_range_constraints.jl")
+include("common_models/quadratic_converter_loss.jl")
 
 # Market bid cost plumbing (PSY orchestration moved out of IOM). Must be included
 # before device-specific files that reference MBC_TYPES / IEC_TYPES.
@@ -505,24 +506,16 @@ 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 PositiveCurrent
+export NegativeCurrent
+export CurrentDirection
 export HVDCFlowDirectionVariable
 export HVDCLosses
-export ConverterDCPower
-export ConverterCurrentDirection
+export HVDCFromDCVoltage
+export HVDCToDCVoltage
+export HVDCReactivePowerFromVariable
+export HVDCReactivePowerToVariable
 
 # Load Variables
 export ShiftUpActivePowerVariable
@@ -770,10 +763,14 @@ export HVDCTwoTerminalLossless
 export HVDCTwoTerminalDispatch
 export HVDCTwoTerminalPiecewiseLoss
 export HVDCTwoTerminalLCC
+export HVDCTwoTerminalVSC
+export HVDCTwoTerminalVSCBin2
 
 # Converter Formulations
 export LosslessConverter
 export LinearLossConverter
+export AbstractQuadraticLossConverter
+export Bin2QuadraticLossConverter
 export QuadraticLossConverter
 
 # DC Line Formulations

src/ac_transmission_models/branch_constructor.jl

@@ -1670,3 +1670,180 @@ function construct_device!(
     add_feedforward_constraints!(container, device_model, devices)
     return
 end
+
+############################################################################
+####################### Two-Terminal VSC HVDC Construct ####################
+############################################################################
+
+# Quadratic / bilinear approximation traits — same scheme used by the MT
+# converter formulations.
+_quad_config(::Type{HVDCTwoTerminalVSC}) = IOM.NoQuadApproxConfig()
+_quad_config(::Type{HVDCTwoTerminalVSCBin2}) =
+    IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)
+_bilinear_config(::Type{HVDCTwoTerminalVSC}) = IOM.NoBilinearApproxConfig()
+_bilinear_config(::Type{HVDCTwoTerminalVSCBin2}) =
+    IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH))
+
+# IOM's quadratic/bilinear approximation helpers take per-device bounds as
+# `Vector{IOM.MinMax}` (one (min, max) pair per device, same order as `devices`).
+function _vsc_v_from_bounds(devices)
+    n = length(devices)
+    bounds = Vector{IOM.MinMax}(undef, n)
+    for (k, d) in enumerate(devices)
+        lims = PSY.get_voltage_limits_from(d)
+        bounds[k] = IOM.MinMax((min = lims.min, max = lims.max))
+    end
+    return bounds
+end
+
+function _vsc_v_to_bounds(devices)
+    n = length(devices)
+    bounds = Vector{IOM.MinMax}(undef, n)
+    for (k, d) in enumerate(devices)
+        lims = PSY.get_voltage_limits_to(d)
+        bounds[k] = IOM.MinMax((min = lims.min, max = lims.max))
+    end
+    return bounds
+end
+
+function _vsc_i_bounds(devices)
+    n = length(devices)
+    bounds = Vector{IOM.MinMax}(undef, n)
+    for (k, d) in enumerate(devices)
+        i_max = _vsc_shared_i_max(d)
+        bounds[k] = IOM.MinMax((min = -i_max, max = i_max))
+    end
+    return bounds
+end
+
+function construct_device!(
+    container::OptimizationContainer,
+    sys::PSY.System,
+    ::ArgumentConstructStage,
+    device_model::DeviceModel{PSY.TwoTerminalVSCLine, F},
+    network_model::NetworkModel{<:AbstractPowerModel},
+) where {F <: AbstractTwoTerminalVSCFormulation}
+    devices = get_available_components(device_model, sys)
+
+    add_variables!(container, FlowActivePowerFromToVariable, devices, F)
+    add_variables!(container, FlowActivePowerToFromVariable, devices, F)
+    add_variables!(container, DCLineCurrentFlowVariable, devices, F)
+    add_variables!(container, HVDCFromDCVoltage, devices, F)
+    add_variables!(container, HVDCToDCVoltage, devices, F)
+
+    _maybe_add_reactive_power_variables!(container, devices, device_model, network_model)
+
+    if _use_linear_loss(F, device_model)
+        ll_devices = _devices_with_linear_loss(devices)
+        if isempty(ll_devices)
+            @warn "use_linear_loss is enabled but every TwoTerminalVSCLine has zero proportional loss terms; no linear-loss variables/constraints will be added."
+        else
+            _add_abs_value_decomposition_variables!(container, ll_devices, device_model)
+        end
+    end
+
+    add_to_expression!(
+        container, ActivePowerBalance, FlowActivePowerFromToVariable,
+        devices, device_model, network_model,
+    )
+    add_to_expression!(
+        container, ActivePowerBalance, FlowActivePowerToFromVariable,
+        devices, device_model, network_model,
+    )
+
+    add_feedforward_arguments!(container, device_model, devices)
+    return
+end
+
+function construct_device!(
+    container::OptimizationContainer,
+    sys::PSY.System,
+    ::ModelConstructStage,
+    device_model::DeviceModel{PSY.TwoTerminalVSCLine, F},
+    network_model::NetworkModel{<:AbstractPowerModel},
+) where {F <: AbstractTwoTerminalVSCFormulation}
+    devices = get_available_components(device_model, sys)
+    time_steps = get_time_steps(container)
+    line_names = [PSY.get_name(d) for d in devices]
+
+    v_f_var = get_variable(container, HVDCFromDCVoltage, PSY.TwoTerminalVSCLine)
+    v_t_var = get_variable(container, HVDCToDCVoltage, PSY.TwoTerminalVSCLine)
+    i_var = get_variable(container, DCLineCurrentFlowVariable, PSY.TwoTerminalVSCLine)
+
+    v_f_bounds = _vsc_v_from_bounds(devices)
+    v_t_bounds = _vsc_v_to_bounds(devices)
+    i_bounds = _vsc_i_bounds(devices)
+
+    quad_cfg, bilin_cfg = _quad_config(F), _bilinear_config(F)
+
+    v_f_sq_expr = IOM._add_quadratic_approx!(
+        quad_cfg, container, PSY.TwoTerminalVSCLine,
+        line_names, time_steps, v_f_var, v_f_bounds, "v_f_sq",
+    )
+    v_t_sq_expr = IOM._add_quadratic_approx!(
+        quad_cfg, container, PSY.TwoTerminalVSCLine,
+        line_names, time_steps, v_t_var, v_t_bounds, "v_t_sq",
+    )
+    i_sq_expr = IOM._add_quadratic_approx!(
+        quad_cfg, container, PSY.TwoTerminalVSCLine,
+        line_names, time_steps, i_var, i_bounds, "i_sq",
+    )
+
+    IOM._add_bilinear_approx!(
+        bilin_cfg, container, PSY.TwoTerminalVSCLine,
+        line_names, time_steps,
+        v_f_sq_expr, i_sq_expr, v_f_var, i_var,
+        v_f_bounds, i_bounds, "vi_ft",
+    )
+    IOM._add_bilinear_approx!(
+        bilin_cfg, container, PSY.TwoTerminalVSCLine,
+        line_names, time_steps,
+        v_t_sq_expr, i_sq_expr, v_t_var, i_var,
+        v_t_bounds, i_bounds, "vi_tf",
+    )
+
+    add_constraints!(
+        container, HVDCCableOhmsLawConstraint, devices, device_model, network_model,
+    )
+    add_constraints!(
+        container, HVDCVSCConverterPowerConstraint, devices, device_model, network_model,
+    )
+    _maybe_add_reactive_power_constraints!(container, devices, device_model, network_model)
+
+    if _use_linear_loss(F, device_model)
+        ll_devices = _devices_with_linear_loss(devices)
+        if !isempty(ll_devices)
+            _add_abs_value_decomposition_constraints!(
+                container, ll_devices, device_model, network_model,
+                DCLineCurrentFlowVariable, _vsc_shared_i_max,
+            )
+        end
+    end
+
+    add_constraint_dual!(container, sys, device_model)
+    add_feedforward_constraints!(container, device_model, devices)
+    return
+end
+
+# AreaBalancePowerModel warning (consistent with other two-terminal formulations).
+function construct_device!(
+    ::OptimizationContainer,
+    ::PSY.System,
+    ::ArgumentConstructStage,
+    ::DeviceModel{PSY.TwoTerminalVSCLine, <:AbstractTwoTerminalVSCFormulation},
+    ::NetworkModel{AreaBalancePowerModel},
+)
+    @warn "AreaBalancePowerModel doesn't model individual line flows for PSY.TwoTerminalVSCLine. Arguments not built"
+    return
+end
+
+function construct_device!(
+    ::OptimizationContainer,
+    ::PSY.System,
+    ::ModelConstructStage,
+    ::DeviceModel{PSY.TwoTerminalVSCLine, <:AbstractTwoTerminalVSCFormulation},
+    ::NetworkModel{AreaBalancePowerModel},
+)
+    @warn "AreaBalancePowerModel doesn't model individual line flows for PSY.TwoTerminalVSCLine. Model not built"
+    return
+end

src/common_models/quadratic_converter_loss.jl

@@ -0,0 +1,124 @@
+# Shared helpers for quadratic / two-term converter losses
+#   loss(I) = a * I^2 + b * |I| + c
+# Used by multi-terminal InterconnectingConverter formulations
+# (Bin2QuadraticLossConverter, QuadraticLossConverter) and two-terminal
+# HVDCTwoTerminalVSC formulations.
+
+#########################################
+######## Loss-curve introspection #######
+#########################################
+
+_get_quadratic_term(loss_fn::PSY.QuadraticCurve) = PSY.get_quadratic_term(loss_fn)
+_get_quadratic_term(loss_fn) = 0.0
+
+# Whether the formulation wants the b*|I| linear-loss term (which requires
+# decomposing the current into positive/negative parts with a direction binary).
+_use_linear_loss(::Type{Bin2QuadraticLossConverter}, _) = true
+_use_linear_loss(::Type{QuadraticLossConverter}, model) =
+    get_attribute(model, "use_linear_loss")
+_use_linear_loss(::Type{HVDCTwoTerminalVSCBin2}, _) = true
+_use_linear_loss(::Type{HVDCTwoTerminalVSC}, model) =
+    get_attribute(model, "use_linear_loss")
+
+# Per-device test: does this device have a nonzero linear loss term anywhere?
+# Dispatched on device type because different PSY devices store loss curves on
+# different fields.
+_has_linear_loss(d::PSY.InterconnectingConverter) =
+    !iszero(PSY.get_proportional_term(PSY.get_loss_function(d)))
+_has_linear_loss(d::PSY.TwoTerminalVSCLine) =
+    !iszero(PSY.get_proportional_term(PSY.get_converter_loss_from(d))) ||
+    !iszero(PSY.get_proportional_term(PSY.get_converter_loss_to(d)))
+
+function _devices_with_linear_loss(devices)
+    return [d for d in devices if _has_linear_loss(d)]
+end
+
+#########################################
+######## Loss expression builder ########
+#########################################
+
+# Returns the JuMP expression  a*i_sq + b*(i_pos + i_neg) + c
+# for a single (device, time). The b*(i_pos+i_neg) term is included only when
+# the formulation has opted into the linear-loss path AND b is nonzero for
+# this specific device.
+function _quadratic_converter_loss_expr(
+    a, b, c, i_sq_t, i_pos_t, i_neg_t; use_linear_loss::Bool,
+)
+    loss = a * i_sq_t + c
+    if use_linear_loss && !iszero(b)
+        loss += b * (i_pos_t + i_neg_t)
+    end
+    return loss
+end
+
+#########################################
+####### Abs-value decomposition #########
+#########################################
+
+# Adds the three variables (PositiveCurrent, NegativeCurrent, CurrentDirection)
+# that decompose a signed current variable into  i = i^+ - i^-  with a binary
+# direction indicator. Called from the ArgumentConstructStage.
+function _add_abs_value_decomposition_variables!(
+    container::OptimizationContainer,
+    devices,
+    ::DeviceModel{D, F},
+) where {D <: PSY.Device, F}
+    add_variables!(container, PositiveCurrent, devices, F)
+    add_variables!(container, NegativeCurrent, devices, F)
+    add_variables!(container, CurrentDirection, devices, F)
+    return
+end
+
+# Adds the three constraints implementing  i = i^+ - i^-  with the big-M
+# direction binary bounds. The CurrentAbsoluteValueConstraint container is
+# created internally with three meta-tagged sub-containers ("", "pos_ub",
+# "neg_ub"). Caller passes the parent current variable type and a function
+# `d -> i_max` so device-specific bound lookups stay device-specific.
+function _add_abs_value_decomposition_constraints!(
+    container::OptimizationContainer,
+    devices,
+    ::DeviceModel{D, F},
+    ::NetworkModel{<:AbstractPowerModel},
+    parent_var_type::Type{<:VariableType},
+    i_max_getter::Function,
+) where {D <: PSY.Device, F}
+    time_steps = get_time_steps(container)
+    names = [PSY.get_name(d) for d in devices]
+    jump_model = get_jump_model(container)
+    i_var = get_variable(container, parent_var_type, D)
+    i_pos_var = get_variable(container, PositiveCurrent, D)
+    i_neg_var = get_variable(container, NegativeCurrent, D)
+    i_dir_var = get_variable(container, CurrentDirection, D)
+
+    abs_val_const = add_constraints_container!(
+        container, CurrentAbsoluteValueConstraint, D, names, time_steps,
+    )
+    pos_ub_const = add_constraints_container!(
+        container, CurrentAbsoluteValueConstraint, D, names, time_steps;
+        meta = "pos_ub",
+    )
+    neg_ub_const = add_constraints_container!(
+        container, CurrentAbsoluteValueConstraint, D, names, time_steps;
+        meta = "neg_ub",
+    )
+
+    for d in devices
+        name = PSY.get_name(d)
+        i_max = i_max_getter(d)
+        for t in time_steps
+            abs_val_const[name, t] = JuMP.@constraint(
+                jump_model,
+                i_var[name, t] == i_pos_var[name, t] - i_neg_var[name, t],
+            )
+            pos_ub_const[name, t] = JuMP.@constraint(
+                jump_model,
+                i_pos_var[name, t] <= i_max * i_dir_var[name, t],
+            )
+            neg_ub_const[name, t] = JuMP.@constraint(
+                jump_model,
+                i_neg_var[name, t] <= i_max * (1 - i_dir_var[name, t]),
+            )
+        end
+    end
+    return
+end