diff --git a/Project.toml b/Project.toml index e749a78..0c2127c 100644 --- a/Project.toml +++ b/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" diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl index 27411c1..d0bdfaf 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. @@ -505,24 +507,14 @@ 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 CurrentAbsoluteValueVariable export HVDCFlowDirectionVariable export HVDCLosses -export ConverterDCPower -export ConverterCurrentDirection +export HVDCFromDCVoltage +export HVDCToDCVoltage +export HVDCReactivePowerFromVariable +export HVDCReactivePowerToVariable # Load Variables export ShiftUpActivePowerVariable @@ -770,11 +762,15 @@ export HVDCTwoTerminalLossless export HVDCTwoTerminalDispatch export HVDCTwoTerminalPiecewiseLoss export HVDCTwoTerminalLCC +export HVDCTwoTerminalVSCNLP +export HVDCTwoTerminalVSCLP # Converter Formulations export LosslessConverter export LinearLossConverter -export QuadraticLossConverter +export AbstractQuadraticLossConverter +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 568b072..0e93c03 100644 --- a/src/ac_transmission_models/branch_constructor.jl +++ b/src/ac_transmission_models/branch_constructor.jl @@ -1670,3 +1670,149 @@ 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), + ) + + add_variables!(container, CurrentAbsoluteValueVariable, devices, F) + + 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 = -_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) + + 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, + ) + + _add_abs_value_constraints!( + container, devices, device_model, network_model, + DCLineCurrentFlowVariable, + ) + + 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..06ffd6c --- /dev/null +++ b/src/common_models/quadratic_converter_loss.jl @@ -0,0 +1,93 @@ +# Shared helpers for quadratic / two-term converter losses +# loss(I) = a * I^2 + b * |I| + c +# Used by multi-terminal InterconnectingConverter formulations +# (QuadraticLossConverterMILP, QuadraticLossConverterNLP) and two-terminal +# TwoTerminalVSCLine formulations (HVDCTwoTerminalVSCLP, HVDCTwoTerminalVSCNLP). +# +# `|I|` is represented by an LP surrogate: a single non-negative variable +# `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 ####### +######################################### + +_get_quadratic_term(loss_fn::PSY.QuadraticCurve) = PSY.get_quadratic_term(loss_fn) +_get_quadratic_term(loss_fn) = 0.0 + +######################################### +######## Loss expression builder ######## +######################################### + +""" + _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`. + +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, + i_sq_t::JuMP.AffExpr, + abs_i_t::JuMP.VariableRef, +) + expr = JuMP.AffExpr(c) + 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 + +function _quadratic_converter_loss_expr( + a::Float64, b::Float64, c::Float64, + i_sq_t::JuMP.QuadExpr, + abs_i_t::JuMP.VariableRef, +) + expr = JuMP.QuadExpr(JuMP.AffExpr(c)) + 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 + +######################################### +######## Absolute-value surrogate ####### +######################################### + +function _add_abs_value_constraints!( + container::OptimizationContainer, + devices, + ::DeviceModel{D, F}, + ::NetworkModel{<:AbstractPowerModel}, + parent_var_type::Type{<:VariableType}, +) 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) + abs_i_var = get_variable(container, CurrentAbsoluteValueVariable, D) + + lower_const = add_constraints_container!( + container, CurrentAbsoluteValueConstraint, D, names, time_steps; + meta = "ge_pos", + ) + upper_const = add_constraints_container!( + container, CurrentAbsoluteValueConstraint, D, names, time_steps; + meta = "ge_neg", + ) + + for d in devices + name = PSY.get_name(d) + for t in time_steps + lower_const[name, t] = JuMP.@constraint( + jump_model, abs_i_var[name, t] >= i_var[name, t], + ) + upper_const[name, t] = JuMP.@constraint( + jump_model, abs_i_var[name, t] >= -i_var[name, t], + ) + end + end + return +end diff --git a/src/core/constraints.jl b/src/core/constraints.jl index 2dc7a17..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. @@ -577,19 +613,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 +626,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..97bda90 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 @@ -201,9 +226,21 @@ Linear Loss InterconnectingConverter Model struct LinearLossConverter <: AbstractConverterFormulation end """ -Quadratic Loss InterconnectingConverter Model +Abstract supertype for InterconnectingConverter formulations with quadratic losses. +""" +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 QuadraticLossConverterMILP <: AbstractQuadraticLossConverter end + +""" +Quadratic Loss InterconnectingConverter using exact bilinear (v·i) and quadratic (i²) +products. Requires an NLP-capable solver (e.g., Ipopt). """ -struct QuadraticLossConverter <: AbstractConverterFormulation 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 660ef5e..1667837 100644 --- a/src/core/variables.jl +++ b/src/core/variables.jl @@ -150,12 +150,6 @@ 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}`` @@ -163,95 +157,13 @@ 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}`` -""" -struct ConverterPositiveCurrent <: VariableType end - -""" -Struct to dispatch the creation of DC Converter Negative Term Current Variables for DC formulations -Docs abbreviation: ``i_c^{-,dc}`` -""" -struct ConverterNegativeCurrent <: VariableType end - -""" -Struct to dispatch the creation of DC Converter Binary for Absolute Value Current Variables for DC formulations -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 - +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 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 +struct CurrentAbsoluteValueVariable <: VariableType end ######################################################### ######################################################### @@ -463,6 +375,34 @@ Docs abbreviation: ``z`` """ struct HVDCPiecewiseBinaryLossVariable <: SparseVariableType end +""" +DC-side voltage at the from-terminal of a two-terminal HVDC link. +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 both `HVDCTwoTerminalVSCNLP` and `HVDCTwoTerminalVSCLP` 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 cfa40e1..8cf9cc9 100644 --- a/src/mt_hvdc_models/HVDCsystems.jl +++ b/src/mt_hvdc_models/HVDCsystems.jl @@ -119,65 +119,19 @@ 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 - +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{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 +get_variable_lower_bound(::Type{CurrentAbsoluteValueVariable}, 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{QuadraticLossConverter}, -) - return Dict{String, Any}( - "voltage_segments" => 3, - "current_segments" => 6, - "bilinear_segments" => 10, - "use_linear_loss" => true, - ) -end +get_variable_upper_bound(::Type{CurrentAbsoluteValueVariable}, d::PSY.InterconnectingConverter, ::Type{<:AbstractConverterFormulation}) = PSY.get_max_dc_current(d) #! format: on @@ -422,7 +376,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 +406,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,381 +518,49 @@ 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)) - 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] - ) - end - end - return -end - -function add_constraints!( - container::OptimizationContainer, - ::Type{ConverterMcCormickEnvelopes}, - 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_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 - return -end function add_constraints!( container::OptimizationContainer, ::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) - 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") - if use_linear_loss - pos_current = get_variable(container, ConverterPositiveCurrent, U) - neg_current = get_variable(container, ConverterNegativeCurrent, U) - end + 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, CurrentAbsoluteValueVariable, U) + + 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) - 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 - for t in time_steps - 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 - end - constraint[name, t] = JuMP.@constraint( - get_jump_model(container), - var_ac_power[name, t] == var_dc_power[name, t] - loss - ) - end - end - return -end - -function add_constraints!( - container::OptimizationContainer, - ::Type{T}, - devices::IS.FlattenIteratorWrapper{U}, - ::DeviceModel{U, V}, - ::NetworkModel{<:AbstractPowerModel}, -) where { - T <: CurrentAbsoluteValueConstraint, - U <: PSY.InterconnectingConverter, - V <: QuadraticLossConverter, -} - 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", - ) - - for d in devices - name = PSY.get_name(d) - I_max = PSY.get_max_dc_current(d) + 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 - constraint[name, t] = JuMP.@constraint( - JuMPmodel, - current_var[name, t] == current_var_pos[name, t] - current_var_neg[name, t] - ) - constraint_pos_ub[name, t] = JuMP.@constraint( - JuMPmodel, - current_var_pos[name, t] <= I_max * current_dir[name, t] + loss = _quadratic_converter_loss_expr( + a, b, c, i_sq_expr[name, t], abs_i_var[name, t], ) - constraint_neg_ub[name, t] = JuMP.@constraint( - JuMPmodel, - current_var_neg[name, t] <= I_max * (1 - current_dir[name, t]) + 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::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..743be6a 100644 --- a/src/mt_hvdc_models/hvdcsystems_constructor.jl +++ b/src/mt_hvdc_models/hvdcsystems_constructor.jl @@ -48,183 +48,112 @@ 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, - ) # ν - 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 #### - ##################### +) where {T <: AbstractQuadraticLossConverter} + devices = get_available_components(model, sys) + add_variables!(container, ActivePowerVariable, devices, T) + add_variables!(container, ConverterCurrent, devices, T) + add_variables!(container, CurrentAbsoluteValueVariable, devices, T) 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 end -function construct_device!( +function _voltage_expr_per_converter( container::OptimizationContainer, - sys::PSY.System, - ::ModelConstructStage, - model::DeviceModel{PSY.InterconnectingConverter, QuadraticLossConverter}, - network_model::NetworkModel{<:AbstractActivePowerModel}, + devices, + ipc_names::Vector{String}, + time_steps, ) - devices = get_available_components( - model, - sys, + 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 - add_constraints!( - container, - ConverterPowerCalculationConstraint, - devices, - model, - network_model, - ) - add_constraints!( - container, - ConverterMcCormickEnvelopes, - devices, - model, - network_model, - ) - add_constraints!( - container, - ConverterLossConstraint, - devices, - model, - network_model, - ) - use_linear_loss = get_attribute(model, "use_linear_loss") - if use_linear_loss - add_constraints!( - container, - CurrentAbsoluteValueConstraint, - devices, - model, - network_model, - ) +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 - add_constraints!( - container, - InterpolationVoltageConstraints, - devices, - model, - network_model, - ) - add_constraints!( - container, - InterpolationCurrentConstraints, - devices, - model, - network_model, - ) - add_constraints!( - container, - InterpolationBilinearConstraints, - devices, - model, - network_model, - ) + return v_bounds, i_bounds +end + +_quad_config(::Type{QuadraticLossConverterMILP}) = + IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH) +_quad_config(::Type{QuadraticLossConverterNLP}) = IOM.NoQuadApproxConfig() +_bilinear_config(::Type{QuadraticLossConverterMILP}) = + IOM.Bin2Config(IOM.SolverSOS2QuadConfig(DEFAULT_INTERPOLATION_LENGTH)) +_bilinear_config(::Type{QuadraticLossConverterNLP}) = IOM.NoBilinearApproxConfig() + +function construct_device!( + container::OptimizationContainer, + sys::PSY.System, + ::ModelConstructStage, + 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] + 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) + + quad_cfg, bilin_cfg = _quad_config(T), _bilinear_config(T) + v_sq_expr = IOM._add_quadratic_approx!( + quad_cfg, + container, PSY.InterconnectingConverter, + ipc_names, time_steps, + v_expr, v_bounds, + "v_sq", + ) + i_sq_expr = IOM._add_quadratic_approx!( + quad_cfg, + container, PSY.InterconnectingConverter, + ipc_names, time_steps, + i_var, i_bounds, + "i_sq", + ) + IOM._add_bilinear_approx!( + bilin_cfg, + container, PSY.InterconnectingConverter, + ipc_names, time_steps, + v_sq_expr, i_sq_expr, + v_expr, i_var, + v_bounds, i_bounds, + "vi", + ) + + # 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!( + container, devices, model, network_model, ConverterCurrent, + ) + + 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), + container, devices, model, get_network_formulation(network_model), ) - #add_constraint_dual!(container, sys, model) return end 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/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/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl index 5422b6f..b71ed14 100644 --- a/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl +++ b/src/twoterminal_hvdc_models/TwoTerminalDC_branches.jl @@ -1433,3 +1433,377 @@ 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{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 + +# Warm starts +get_variable_warm_start_value(::Type{DCLineCurrentFlowVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_dc_current(d) +get_variable_warm_start_value(::Type{HVDCReactivePowerFromVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_reactive_power_from(d) +get_variable_warm_start_value(::Type{HVDCReactivePowerToVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_reactive_power_to(d) +get_variable_warm_start_value(::Type{FlowActivePowerFromToVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_active_power_flow(d) +get_variable_warm_start_value(::Type{FlowActivePowerToFromVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = -PSY.get_active_power_flow(d) + +# Active power flow bounds (per-terminal) +get_variable_lower_bound(::Type{FlowActivePowerFromToVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_active_power_limits_from(d).min +get_variable_upper_bound(::Type{FlowActivePowerFromToVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_active_power_limits_from(d).max +get_variable_lower_bound(::Type{FlowActivePowerToFromVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_active_power_limits_to(d).min +get_variable_upper_bound(::Type{FlowActivePowerToFromVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_active_power_limits_to(d).max + +# Reactive power bounds (per-terminal) +get_variable_lower_bound(::Type{HVDCReactivePowerFromVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_reactive_power_limits_from(d).min +get_variable_upper_bound(::Type{HVDCReactivePowerFromVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_reactive_power_limits_from(d).max +get_variable_lower_bound(::Type{HVDCReactivePowerToVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_reactive_power_limits_to(d).min +get_variable_upper_bound(::Type{HVDCReactivePowerToVariable}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_reactive_power_limits_to(d).max + +# DC voltage bounds (per-terminal) +get_variable_lower_bound(::Type{HVDCFromDCVoltage}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_voltage_limits_from(d).min +get_variable_upper_bound(::Type{HVDCFromDCVoltage}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_voltage_limits_from(d).max +get_variable_lower_bound(::Type{HVDCToDCVoltage}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_voltage_limits_to(d).min +get_variable_upper_bound(::Type{HVDCToDCVoltage}, d::PSY.TwoTerminalVSCLine, ::Type{<:AbstractTwoTerminalVSCFormulation}) = PSY.get_voltage_limits_to(d).max + +# Shared cable current bounds — must respect BOTH terminals' I_max ratings. +_vsc_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) + +# 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 + +####################### 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") + + abs_i_var = get_variable(container, CurrentAbsoluteValueVariable, U) + + 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 + 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], abs_i_t, + ) + loss_tf = _quadratic_converter_loss_expr( + a_t, b_t, c_t, i_sq_expr[name, t], abs_i_t, + ) + 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 + +# `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_octagon" => true) +end diff --git a/test/Project.toml b/test/Project.toml index 9445668..43494b2 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"} +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 2b158b3..30acb05 100644 --- a/test/test_device_hvdc.jl +++ b/test/test_device_hvdc.jl @@ -97,33 +97,215 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED end -@testset "HVDC System with Losses Network" begin +@testset "HVDC MILP vs NLP QuadraticLossConverter agreement" begin + # 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) + set_device_model!(template, DeviceModel(Line, StaticBranch)) + set_device_model!(template, TModelHVDCLine, DCLossyLine) + set_device_model!( + template, + DeviceModel(InterconnectingConverter, formulation), + ) + 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 + + 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 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() 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; - attributes = Dict( - "voltage_segments" => 3, - "current_segments" => 3, - "bilinear_segments" => 3, - "use_linear_loss" => true, - ), + set_device_model!( + template, + DeviceModel(InterconnectingConverter, QuadraticLossConverterMILP), ) - set_device_model!(template, ipc_model) set_hvdc_network_model!(template, VoltageDispatchHVDCNetworkModel) - model = - DecisionModel( - template, - sys; - store_variable_names = true, - optimizer = HiGHS_optimizer, - ) + 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, + CurrentAbsoluteValueVariable, + InterconnectingConverter, + ).data, + ) + @test isapprox(abs_i_vals, abs.(i_vals); atol = 1e-6) +end + +############################################################################## +################ Two-Terminal VSC HVDC tests ################################# +############################################################################## + +# Build a small AC test system and replace an AC line with a TwoTerminalVSCLine +# so we have a concrete VSC device to exercise the formulation against. +function _generate_test_vsc_sys(; + g = 50.0, + rating_from = 2.0, + rating_to = 2.0, + loss_a = 0.01, + loss_b = 0.0, + loss_c = 0.0, +) + sys = build_system(PSITestSystems, "c_sys5_uc"; force_build = true) + line = get_component(Line, sys, "1") + remove_component!(sys, line) + + vsc = TwoTerminalVSCLine(; + name = get_name(line), + available = true, + arc = get_arc(line), + active_power_flow = 0.0, + rating = max(rating_from, rating_to), + active_power_limits_from = (min = -rating_from, max = rating_from), + active_power_limits_to = (min = -rating_to, max = rating_to), + g = g, + dc_current = 0.0, + reactive_power_from = 0.0, + dc_voltage_control_from = true, + ac_voltage_control_from = true, + dc_setpoint_from = 0.0, + ac_setpoint_from = 1.0, + converter_loss_from = QuadraticCurve(loss_a, loss_b, loss_c), + max_dc_current_from = 5.0, + rating_from = rating_from, + reactive_power_limits_from = (min = -rating_from, max = rating_from), + power_factor_weighting_fraction_from = 1.0, + voltage_limits_from = (min = 0.95, max = 1.05), + reactive_power_to = 0.0, + dc_voltage_control_to = true, + ac_voltage_control_to = true, + dc_setpoint_to = 0.0, + ac_setpoint_to = 1.0, + converter_loss_to = QuadraticCurve(loss_a, loss_b, loss_c), + max_dc_current_to = 5.0, + rating_to = rating_to, + reactive_power_limits_to = (min = -rating_to, max = rating_to), + power_factor_weighting_fraction_to = 1.0, + voltage_limits_to = (min = 0.95, max = 1.05), + ) + add_component!(sys, vsc) + return sys +end + +function _build_vsc_model(formulation, network, optimizer; sys = _generate_test_vsc_sys()) + template = OperationsProblemTemplate(NetworkModel(network)) + set_device_model!(template, ThermalStandard, ThermalDispatchNoMin) + set_device_model!(template, RenewableDispatch, RenewableFullDispatch) + set_device_model!(template, PowerLoad, StaticPowerLoad) + set_device_model!(template, DeviceModel(Line, StaticBranch)) + set_device_model!(template, DeviceModel(TwoTerminalVSCLine, formulation)) + return DecisionModel( + template, sys; store_variable_names = true, optimizer = optimizer, + ) +end + +# 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 + # 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: 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