diff --git a/.github/workflows/cross-package-test.yml b/.github/workflows/cross-package-test.yml new file mode 100644 index 0000000000..cda8a3ede3 --- /dev/null +++ b/.github/workflows/cross-package-test.yml @@ -0,0 +1,47 @@ +name: CrossPackageTest + +on: + push: + branches: [main] + tags: [v*] + pull_request: + +jobs: + test: + name: Julia v${{ matrix.julia-version }} - ${{ matrix.package_name }} + runs-on: ${{ matrix.os }} + strategy: + matrix: + julia-version: [1] + os: [ubuntu-latest] + package_name: [HydroPowerSimulations, StorageSystemsSimulations] + continue-on-error: true + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.julia-version }} + arch: x64 + - uses: julia-actions/julia-buildpkg@latest + - name: Clone ${{matrix.package_name}} + uses: actions/checkout@v2 + with: + repository: NREL-Sienna/${{matrix.package_name}}.jl + path: downstream + - name: Run the tests + shell: julia --project=downstream {0} + run: | + using Pkg + try + # Force it to use this PR's version of the package + Pkg.develop(PackageSpec(path=".")) # resolver may fail with main deps + Pkg.update() + Pkg.test() # resolver may fail with test time deps + catch err + err isa Pkg.Resolve.ResolverError || rethrow() + # If we can't resolve that means this is incompatible by SemVer, and this is fine. + # It means we marked this as a breaking change, so we don't need to worry about + # mistakenly introducing a breaking change as we have intentionally made one. + @info "Not compatible with this release. No problem." exception=err + exit(0) # Exit immediately, as a success + end diff --git a/Project.toml b/Project.toml index 53a901c19e..79f3213b88 100644 --- a/Project.toml +++ b/Project.toml @@ -38,15 +38,14 @@ DocStringExtensions = "~v0.9" HDF5 = "~0.17" InfrastructureSystems = "2" InteractiveUtils = "1" -JSON = "0.21" JSON3 = "1" JuMP = "1" LinearAlgebra = "1" Logging = "1" MathOptInterface = "1" PowerModels = "^0.21" -PowerNetworkMatrices = "^0.11.1" -PowerSystems = "4" +PowerNetworkMatrices = "^0.11" +PowerSystems = "^4.4" PrettyTables = "2" ProgressMeter = "^1.5" Serialization = "1" diff --git a/docs/src/api/PowerSimulations.md b/docs/src/api/PowerSimulations.md index be88337f83..a00e3bce49 100644 --- a/docs/src/api/PowerSimulations.md +++ b/docs/src/api/PowerSimulations.md @@ -269,7 +269,6 @@ FlowLimitConstraint FlowRateConstraint FlowRateConstraintFromTo FlowRateConstraintToFrom -HVDCLossesAbsoluteValue HVDCPowerBalance NetworkFlowConstraint RateLimitConstraint diff --git a/src/PowerSimulations.jl b/src/PowerSimulations.jl index f4d882d98d..ddca095e63 100644 --- a/src/PowerSimulations.jl +++ b/src/PowerSimulations.jl @@ -412,7 +412,6 @@ import TimeSeries # I/O Imports import DataFrames -import JSON import CSV import HDF5 import PrettyTables diff --git a/src/core/constraints.jl b/src/core/constraints.jl index a36e8257cd..817130ba45 100644 --- a/src/core/constraints.jl +++ b/src/core/constraints.jl @@ -332,23 +332,8 @@ The specified constraint is formulated as: ``` """ struct PhaseAngleControlLimit <: ConstraintType end -""" -Struct to create the constraints that set the losses through a lossy HVDC two-terminal line. - -For more information check [Branch Formulations](@ref PowerSystems.Branch-Formulations). - -The specified constraints are formulated as: - -```math -\\begin{align*} -& f_t^\\text{to-from} - f_t^\\text{from-to} \\le \\ell_t,\\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ -& f_t^\\text{from-to} - f_t^\\text{to-from} \\le \\ell_t,\\quad \\forall t \\in \\{1,\\dots, T\\} -\\end{align*} -``` -""" -struct HVDCLossesAbsoluteValue <: ConstraintType end -struct HVDCDirection <: ConstraintType end struct InterfaceFlowLimit <: ConstraintType end +struct HVDCFlowCalculationConstraint <: ConstraintType end abstract type PowerVariableLimitsConstraint <: ConstraintType end """ diff --git a/src/core/formulations.jl b/src/core/formulations.jl index 1b6ddb109a..f7ced88017 100644 --- a/src/core/formulations.jl +++ b/src/core/formulations.jl @@ -134,6 +134,11 @@ struct HVDCTwoTerminalLossless <: AbstractTwoTerminalDCLineFormulation end Branch type to represent lossy power flow on DC lines """ struct HVDCTwoTerminalDispatch <: AbstractTwoTerminalDCLineFormulation end +""" +Branch type to represent piecewise lossy power flow on two terminal DC lines +""" +struct HVDCTwoTerminalPiecewiseLoss <: AbstractTwoTerminalDCLineFormulation end + # Not Implemented # struct VoltageSourceDC <: AbstractTwoTerminalDCLineFormulation end diff --git a/src/core/initial_conditions.jl b/src/core/initial_conditions.jl index 329cb47c94..8cb845430d 100644 --- a/src/core/initial_conditions.jl +++ b/src/core/initial_conditions.jl @@ -3,7 +3,7 @@ Container for the initial condition data """ mutable struct InitialCondition{ T <: InitialConditionType, - U <: Union{JuMP.VariableRef, Float64}, + U <: Union{JuMP.VariableRef, Float64, Nothing}, } component::PSY.Component value::U @@ -39,16 +39,22 @@ function get_condition( return jump_value(p.value) end +function get_condition( + ::InitialCondition{T, Nothing}, +) where {T <: InitialConditionType} + return nothing +end + get_component(ic::InitialCondition) = ic.component get_value(ic::InitialCondition) = ic.value get_component_name(ic::InitialCondition) = PSY.get_name(ic.component) get_component_type(ic::InitialCondition) = typeof(ic.component) get_ic_type( ::Type{InitialCondition{T, U}}, -) where {T <: InitialConditionType, U <: Union{JuMP.VariableRef, Float64}} = T +) where {T <: InitialConditionType, U <: Union{JuMP.VariableRef, Float64, Nothing}} = T get_ic_type( ::InitialCondition{T, U}, -) where {T <: InitialConditionType, U <: Union{JuMP.VariableRef, Float64}} = T +) where {T <: InitialConditionType, U <: Union{JuMP.VariableRef, Float64, Nothing}} = T """ Stores data to populate initial conditions before the build call diff --git a/src/core/optimization_container.jl b/src/core/optimization_container.jl index e3e7f2ffb7..b94438f9bd 100644 --- a/src/core/optimization_container.jl +++ b/src/core/optimization_container.jl @@ -33,7 +33,12 @@ function get_objective_expression(v::ObjectiveFunction) if iszero(v.variant_terms) return v.invariant_terms else - return JuMP.add_to_expression!(v.variant_terms, v.invariant_terms) + # JuMP doesn't support expression conversion from Affn to QuadExpressions + if isa(v.invariant_terms, JuMP.GenericQuadExpr) + return JuMP.add_to_expression!(v.invariant_terms, v.variant_terms) + else + return JuMP.add_to_expression!(v.variant_terms, v.invariant_terms) + end end end get_sense(v::ObjectiveFunction) = v.sense @@ -1503,7 +1508,8 @@ function _add_initial_condition_container!( else param_type = Float64 end - ini_conds = Vector{InitialCondition{T, param_type}}(undef, length_devices) + ini_type = Union{InitialCondition{T, param_type}, InitialCondition{T, Nothing}} + ini_conds = Vector{ini_type}(undef, length_devices) _assign_container!(container.initial_conditions, ic_key, ini_conds) return ini_conds end @@ -1716,6 +1722,9 @@ function _process_duals(container::OptimizationContainer, lp_optimizer) if JuMP.has_upper_bound(first(variable)) cache[key][:ub] = JuMP.upper_bound.(variable) end + if JuMP.is_fixed(first(variable)) && is_integer_flag + cache[key][:fixed_int_value] = jump_value.(v) + end cache[key][:integer] = is_integer_flag JuMP.fix.(variable, var_cache[key]; force = true) end @@ -1756,6 +1765,9 @@ function _process_duals(container::OptimizationContainer, lp_optimizer) else JuMP.unfix.(variable) JuMP.set_binary.(variable) + if haskey(cache[key], :fixed_int_value) + JuMP.fix.(variable, cache[key][:fixed_int_value]) + end #= Needed if a model has integer variables if haskey(cache[key], :lb) && JuMP.has_lower_bound(first(variable)) JuMP.set_lower_bound.(variable, cache[key][:lb]) diff --git a/src/core/variables.jl b/src/core/variables.jl index 42d5bc0560..3050c31f51 100644 --- a/src/core/variables.jl +++ b/src/core/variables.jl @@ -228,8 +228,36 @@ Docs abbreviation: ``u^\\text{dir}`` """ struct HVDCFlowDirectionVariable <: VariableType end +""" +Struct to dispatch the creation of HVDC Received Flow at From Bus Variables for PWL formulations + +Docs abbreviation: ``x`` +""" +struct HVDCActivePowerReceivedFromVariable <: VariableType end + +""" +Struct to dispatch the creation of HVDC Received Flow at To Bus Variables for PWL formulations + +Docs abbreviation: ``y`` +""" +struct HVDCActivePowerReceivedToVariable <: VariableType end + abstract type SparseVariableType <: VariableType end +""" +Struct to dispatch the creation of HVDC Piecewise Loss Variables + +Docs abbreviation: ``h`` or ``w`` +""" +struct HVDCPiecewiseLossVariable <: SparseVariableType end + +""" +Struct to dispatch the creation of HVDC Piecewise Binary Loss Variables + +Docs abbreviation: ``z`` +""" +struct HVDCPiecewiseBinaryLossVariable <: SparseVariableType end + """ Struct to dispatch the creation of piecewise linear cost variables for objective function @@ -274,6 +302,8 @@ const START_VARIABLES = (HotStartVariable, WarmStartVariable, ColdStartVariable) should_write_resulting_value(::Type{PieceWiseLinearCostVariable}) = false should_write_resulting_value(::Type{PieceWiseLinearBlockOffer}) = false +should_write_resulting_value(::Type{HVDCPiecewiseLossVariable}) = false +should_write_resulting_value(::Type{HVDCPiecewiseBinaryLossVariable}) = false convert_result_to_natural_units(::Type{ActivePowerVariable}) = true convert_result_to_natural_units(::Type{PowerAboveMinimumVariable}) = true convert_result_to_natural_units(::Type{ActivePowerInVariable}) = true diff --git a/src/devices_models/device_constructors/branch_constructor.jl b/src/devices_models/device_constructors/branch_constructor.jl index 9aff973682..8230cef82d 100644 --- a/src/devices_models/device_constructors/branch_constructor.jl +++ b/src/devices_models/device_constructors/branch_constructor.jl @@ -167,7 +167,7 @@ function construct_device!( end function construct_device!( - ::OptimizationContainer, + container::OptimizationContainer, sys::PSY.System, ::ModelConstructStage, model::DeviceModel{<:PSY.ACBranch, StaticBranchUnbounded}, @@ -694,7 +694,6 @@ function construct_device!( add_constraints!(container, FlowRateConstraintFromTo, devices, model, network_model) add_constraints!(container, FlowRateConstraintToFrom, devices, model, network_model) add_constraints!(container, HVDCPowerBalance, devices, model, network_model) - add_constraints!(container, HVDCLossesAbsoluteValue, devices, model, network_model) return end @@ -719,6 +718,7 @@ function construct_device!( HVDCTwoTerminalDispatch(), ) add_variables!(container, HVDCFlowDirectionVariable, devices, HVDCTwoTerminalDispatch()) + add_variables!(container, HVDCLosses, devices, HVDCTwoTerminalDispatch()) add_to_expression!( container, ActivePowerBalance, @@ -754,6 +754,135 @@ function construct_device!( return end +function construct_device!( + container::OptimizationContainer, + sys::PSY.System, + ::ArgumentConstructStage, + model::DeviceModel{T, U}, + network_model::NetworkModel{<:AbstractPTDFModel}, +) where { + T <: TwoTerminalHVDCTypes, + U <: HVDCTwoTerminalPiecewiseLoss, +} + devices = + get_available_components(model, sys) + add_variables!( + container, + HVDCActivePowerReceivedFromVariable, + devices, + HVDCTwoTerminalPiecewiseLoss(), + ) + add_variables!( + container, + HVDCActivePowerReceivedToVariable, + devices, + HVDCTwoTerminalPiecewiseLoss(), + ) + _add_sparse_pwl_loss_variables!(container, devices, model) + add_to_expression!( + container, + ActivePowerBalance, + HVDCActivePowerReceivedFromVariable, + devices, + model, + network_model, + ) + add_to_expression!( + container, + ActivePowerBalance, + HVDCActivePowerReceivedToVariable, + devices, + model, + network_model, + ) + return +end + +function construct_device!( + container::OptimizationContainer, + sys::PSY.System, + ::ModelConstructStage, + model::DeviceModel{T, U}, + network_model::NetworkModel{<:AbstractPTDFModel}, +) where { + T <: TwoTerminalHVDCTypes, + U <: HVDCTwoTerminalPiecewiseLoss, +} + devices = + get_available_components(model, sys) + add_constraints!(container, FlowRateConstraintFromTo, devices, model, network_model) + add_constraints!(container, FlowRateConstraintToFrom, devices, model, network_model) + add_constraints!( + container, + HVDCFlowCalculationConstraint, + devices, + model, + network_model, + ) + return +end + +# TODO: Other models besides PTDF +#= +function construct_device!( + container::OptimizationContainer, + sys::PSY.System, + ::ArgumentConstructStage, + model::DeviceModel{T, HVDCTwoTerminalPiecewiseLoss}, + network_model::NetworkModel{<:PM.AbstractActivePowerModel}, +) where {T <: TwoTerminalHVDCTypes} + devices = + get_available_components(model, sys) + add_variables!( + container, + FlowActivePowerToFromVariable, + devices, + HVDCTwoTerminalDispatch(), + ) + add_variables!( + container, + FlowActivePowerFromToVariable, + devices, + HVDCTwoTerminalDispatch(), + ) + add_variables!(container, HVDCFlowDirectionVariable, devices, HVDCTwoTerminalDispatch()) + add_to_expression!( + container, + ActivePowerBalance, + FlowActivePowerToFromVariable, + devices, + model, + network_model, + ) + add_to_expression!( + container, + ActivePowerBalance, + FlowActivePowerFromToVariable, + devices, + model, + network_model, + ) + add_feedforward_arguments!(container, model, devices) + return +end + +function construct_device!( + container::OptimizationContainer, + sys::PSY.System, + ::ModelConstructStage, + model::DeviceModel{T, HVDCTwoTerminalPiecewiseLoss}, + network_model::NetworkModel{CopperPlatePowerModel}, +) where {T <: TwoTerminalHVDCTypes} + devices = + get_available_components(model, sys) + @warn "CopperPlatePowerModel models with HVDC ignores inter-area losses" + add_constraints!(container, FlowRateConstraintFromTo, devices, model, network_model) + add_constraints!(container, FlowRateConstraintToFrom, devices, model, network_model) + add_constraint_dual!(container, sys, model) + return +end +=# + function construct_device!( container::OptimizationContainer, sys::PSY.System, @@ -765,7 +894,6 @@ function construct_device!( add_constraints!(container, FlowRateConstraintFromTo, devices, model, network_model) add_constraints!(container, FlowRateConstraintToFrom, devices, model, network_model) add_constraints!(container, HVDCPowerBalance, devices, model, network_model) - add_constraints!(container, HVDCDirection, devices, model, network_model) add_constraint_dual!(container, sys, model) add_feedforward_constraints!(container, model, devices) return diff --git a/src/devices_models/devices/AC_branches.jl b/src/devices_models/devices/AC_branches.jl index e6cb1c74dc..48767cd794 100644 --- a/src/devices_models/devices/AC_branches.jl +++ b/src/devices_models/devices/AC_branches.jl @@ -37,6 +37,10 @@ get_variable_binary(::FlowActivePowerSlackLowerBound, ::Type{<:PSY.ACBranch}, :: # These two methods are defined to avoid ambiguities get_variable_binary(::FlowActivePowerSlackUpperBound, ::Type{<:PSY.TwoTerminalHVDCLine}, ::AbstractTwoTerminalDCLineFormulation,) = false get_variable_binary(::FlowActivePowerSlackLowerBound, ::Type{<:PSY.TwoTerminalHVDCLine}, ::AbstractTwoTerminalDCLineFormulation,) = false +get_variable_binary(::HVDCPiecewiseLossVariable, ::Type{<:PSY.TwoTerminalHVDCLine}, ::AbstractTwoTerminalDCLineFormulation,) = false +get_variable_binary(::HVDCActivePowerReceivedFromVariable, ::Type{<:PSY.TwoTerminalHVDCLine}, ::AbstractTwoTerminalDCLineFormulation,) = false +get_variable_binary(::HVDCActivePowerReceivedToVariable, ::Type{<:PSY.TwoTerminalHVDCLine}, ::AbstractTwoTerminalDCLineFormulation,) = false +get_variable_binary(::HVDCPiecewiseBinaryLossVariable, ::Type{<:PSY.TwoTerminalHVDCLine}, ::AbstractTwoTerminalDCLineFormulation,) = true get_variable_upper_bound(::FlowActivePowerSlackUpperBound, ::PSY.ACBranch, ::AbstractBranchFormulation) = nothing get_variable_lower_bound(::FlowActivePowerSlackUpperBound, ::PSY.ACBranch, ::AbstractBranchFormulation) = 0.0 get_variable_upper_bound(::FlowActivePowerSlackLowerBound, ::PSY.ACBranch, ::AbstractBranchFormulation) = nothing @@ -174,6 +178,142 @@ function branch_rate_bounds!( return end +################################## PWL Loss Variables ################################## + +function _check_pwl_loss_model(devices) + first_loss = PSY.get_loss(first(devices)) + first_loss_type = typeof(first_loss) + for d in devices + loss = PSY.get_loss(d) + if !isa(loss, first_loss_type) + error( + "Not all TwoTerminal HVDC lines have the same loss model data. Check that all loss models are LinearCurve or PiecewiseIncrementalCurve", + ) + end + if isa(first_loss, PSY.PiecewiseIncrementalCurve) + len_first_loss = length(PSY.get_slopes(first_loss)) + len_loss = length(PSY.get_slopes(loss)) + if len_first_loss != len_loss + error( + "Different length of PWL segments for TwoTerminal HVDC losses are not supported. Check that all HVDC data have the same amount of PWL segments.", + ) + end + end + end +end + +function _add_dense_pwl_loss_variables!( + container::OptimizationContainer, + devices, + model::DeviceModel{D, HVDCTwoTerminalPiecewiseLoss}, +) where {D <: TwoTerminalHVDCTypes} + # Check if type and length of PWL loss model are the same for all devices + _check_pwl_loss_model(devices) + + # Create Variables + time_steps = get_time_steps(container) + settings = get_settings(container) + formulation = HVDCTwoTerminalPiecewiseLoss() + T = HVDCPiecewiseLossVariable + binary = get_variable_binary(T(), D, formulation) + first_loss = PSY.get_loss(first(devices)) + if isa(first_loss, PSY.LinearCurve) + len_segments = 4 # 2*1 + 2 + elseif isa(first_loss, PSY.PiecewiseIncrementalCurve) + len_segments = 2 * length(PSY.get_slopes(first_loss)) + 2 + else + error("Should not be here") + end + + segments = ["pwl_$i" for i in 1:len_segments] + T = HVDCPiecewiseLossVariable + variable = add_variable_container!( + container, + T(), + D, + [PSY.get_name(d) for d in devices], + segments, + time_steps, + ) + + for t in time_steps, s in segments, d in devices + name = PSY.get_name(d) + variable[name, s, t] = JuMP.@variable( + get_jump_model(container), + base_name = "$(T)_$(D)_{$(name), $(s), $(t)}", + binary = binary + ) + ub = get_variable_upper_bound(T(), d, formulation) + ub !== nothing && JuMP.set_upper_bound(variable[name, s, t], ub) + + lb = get_variable_lower_bound(T(), d, formulation) + lb !== nothing && JuMP.set_lower_bound(variable[name, s, t], lb) + + if get_warm_start(settings) + init = get_variable_warm_start_value(T(), d, formulation) + init !== nothing && JuMP.set_start_value(variable[name, s, t], init) + end + end +end + +# Full Binary +function _add_sparse_pwl_loss_variables!( + container::OptimizationContainer, + devices, + ::DeviceModel{D, HVDCTwoTerminalPiecewiseLoss}, +) where {D <: TwoTerminalHVDCTypes} + # Check if type and length of PWL loss model are the same for all devices + #_check_pwl_loss_model(devices) + + # Create Variables + time_steps = get_time_steps(container) + settings = get_settings(container) + formulation = HVDCTwoTerminalPiecewiseLoss() + T = HVDCPiecewiseLossVariable + binary_T = get_variable_binary(T(), D, formulation) + U = HVDCPiecewiseBinaryLossVariable + binary_U = get_variable_binary(U(), D, formulation) + first_loss = PSY.get_loss(first(devices)) + if isa(first_loss, PSY.LinearCurve) + len_segments = 3 # 2*1 + 1 + elseif isa(first_loss, PSY.PiecewiseIncrementalCurve) + len_segments = 2 * length(PSY.get_slopes(first_loss)) + 1 + else + error("Should not be here") + end + + var_container = lazy_container_addition!(container, T(), D) + var_container_binary = lazy_container_addition!(container, U(), D) + + for d in devices + name = PSY.get_name(d) + for t in time_steps + pwlvars = Array{JuMP.VariableRef}(undef, len_segments) + pwlvars_bin = Array{JuMP.VariableRef}(undef, len_segments) + for i in 1:len_segments + pwlvars[i] = + var_container[(name, i, t)] = JuMP.@variable( + get_jump_model(container), + base_name = "$(T)_$(name)_{pwl_$(i), $(t)}", + binary = binary_T + ) + ub = get_variable_upper_bound(T(), d, formulation) + ub !== nothing && JuMP.set_upper_bound(var_container[name, i, t], ub) + + lb = get_variable_lower_bound(T(), d, formulation) + lb !== nothing && JuMP.set_lower_bound(var_container[name, i, t], lb) + + pwlvars_bin[i] = + var_container_binary[(name, i, t)] = JuMP.@variable( + get_jump_model(container), + base_name = "$(U)_$(name)_{pwl_$(i), $(t)}", + binary = binary_U + ) + end + end + end +end + ################################## Rate Limits constraint_infos ############################ """ diff --git a/src/devices_models/devices/TwoTerminalDC_branches.jl b/src/devices_models/devices/TwoTerminalDC_branches.jl index b5805456d0..29685dbb8a 100644 --- a/src/devices_models/devices/TwoTerminalDC_branches.jl +++ b/src/devices_models/devices/TwoTerminalDC_branches.jl @@ -35,15 +35,21 @@ get_variable_multiplier( ::FlowActivePowerToFromVariable, ::Type{<:PSY.TwoTerminalHVDCLine}, ::AbstractTwoTerminalDCLineFormulation, -) = 1.0 +) = -1.0 function get_variable_multiplier( ::HVDCLosses, d::PSY.TwoTerminalHVDCLine, ::HVDCTwoTerminalDispatch, ) - l1 = PSY.get_loss(d).l1 - l0 = PSY.get_loss(d).l0 + loss = PSY.get_loss(d) + if !isa(loss, PSY.LinearCurve) + error( + "HVDCTwoTerminalDispatch of branch $(PSY.get_name(d)) only accepts LinearCurve for loss models.", + ) + end + l1 = PSY.get_proportional_term(loss) + l0 = PSY.get_constant_term(loss) if l1 == 0.0 && l0 == 0.0 return 0.0 else @@ -105,13 +111,43 @@ get_variable_lower_bound( ::HVDCTwoTerminalDispatch, ) = PSY.get_active_power_limits_to(d).min +get_variable_upper_bound( + ::HVDCActivePowerReceivedFromVariable, + d::PSY.TwoTerminalHVDCLine, + ::AbstractTwoTerminalDCLineFormulation, +) = PSY.get_active_power_limits_from(d).max + +get_variable_lower_bound( + ::HVDCActivePowerReceivedFromVariable, + d::PSY.TwoTerminalHVDCLine, + ::AbstractTwoTerminalDCLineFormulation, +) = PSY.get_active_power_limits_from(d).min + +get_variable_upper_bound( + ::HVDCActivePowerReceivedToVariable, + d::PSY.TwoTerminalHVDCLine, + ::AbstractTwoTerminalDCLineFormulation, +) = PSY.get_active_power_limits_to(d).max + +get_variable_lower_bound( + ::HVDCActivePowerReceivedToVariable, + d::PSY.TwoTerminalHVDCLine, + ::AbstractTwoTerminalDCLineFormulation, +) = PSY.get_active_power_limits_to(d).min + function get_variable_upper_bound( ::HVDCLosses, d::PSY.TwoTerminalHVDCLine, ::HVDCTwoTerminalDispatch, ) - l1 = PSY.get_loss(d).l1 - l0 = PSY.get_loss(d).l0 + loss = PSY.get_loss(d) + if !isa(loss, PSY.LinearCurve) + error( + "HVDCTwoTerminalDispatch of branch $(PSY.get_name(d)) only accepts LinearCurve for loss models.", + ) + end + l1 = PSY.get_proportional_term(loss) + l0 = PSY.get_constant_term(loss) if l1 == 0.0 && l0 == 0.0 return 0.0 else @@ -119,6 +155,18 @@ function get_variable_upper_bound( end end +get_variable_upper_bound( + ::HVDCPiecewiseLossVariable, + d::PSY.TwoTerminalHVDCLine, + ::Union{HVDCTwoTerminalDispatch, HVDCTwoTerminalPiecewiseLoss}, +) = 1.0 + +get_variable_lower_bound( + ::HVDCPiecewiseLossVariable, + d::PSY.TwoTerminalHVDCLine, + ::Union{HVDCTwoTerminalDispatch, HVDCTwoTerminalPiecewiseLoss}, +) = 0.0 + function get_default_time_series_names( ::Type{U}, ::Type{V}, @@ -139,6 +187,159 @@ get_initial_conditions_device_model( ) where {T <: PSY.TwoTerminalHVDCLine, U <: AbstractTwoTerminalDCLineFormulation} = DeviceModel(T, U) +####################################### PWL Constraints ####################################################### + +function _get_range_segments(::PSY.TwoTerminalHVDCLine, loss::PSY.LinearCurve) + return 1:4 +end + +function _get_range_segments(::PSY.TwoTerminalHVDCLine, loss::PSY.PiecewiseIncrementalCurve) + loss_factors = PSY.get_slopes(loss) + return 1:(2 * length(loss_factors) + 2) +end + +function _get_pwl_loss_params(d::PSY.TwoTerminalHVDCLine, loss::PSY.LinearCurve) + from_to_loss_params = Vector{Float64}(undef, 4) + to_from_loss_params = Vector{Float64}(undef, 4) + loss_factor = PSY.get_proportional_term(loss) + P_send0 = PSY.get_constant_term(loss) + P_max_ft = PSY.get_active_power_limits_from(d).max + P_max_tf = PSY.get_active_power_limits_to(d).max + if P_max_ft != P_max_tf + error( + "HVDC Line $(PSY.get_name(d)) has non-symmetrical limits for from and to, that are not supported in the HVDCTwoTerminalPiecewiseLoss formulation", + ) + end + P_sendS = P_max_ft + ### Update Params Vectors ### + from_to_loss_params[1] = -P_sendS - P_send0 + from_to_loss_params[2] = -P_send0 + from_to_loss_params[3] = 0.0 + from_to_loss_params[4] = P_sendS * (1 - loss_factor) + + to_from_loss_params[1] = P_sendS * (1 - loss_factor) + to_from_loss_params[2] = 0.0 + to_from_loss_params[3] = -P_send0 + to_from_loss_params[4] = -P_sendS - P_send0 + + return from_to_loss_params, to_from_loss_params +end + +function _get_pwl_loss_params( + d::PSY.TwoTerminalHVDCLine, + loss::PSY.PiecewiseIncrementalCurve, +) + p_breakpoints = PSY.get_x_coords(loss) + loss_factors = PSY.get_slopes(loss) + len_segments = length(loss_factors) + len_variables = 2 * len_segments + 2 + from_to_loss_params = Vector{Float64}(undef, len_variables) + to_from_loss_params = similar(from_to_loss_params) + P_max_ft = PSY.get_active_power_limits_from(d).max + P_max_tf = PSY.get_active_power_limits_to(d).max + if P_max_ft != P_max_tf + error( + "HVDC Line $(PSY.get_name(d)) has non-symmetrical limits for from and to, that are not supported in the HVDCTwoTerminalPiecewiseLoss formulation", + ) + end + if P_max_ft != last(p_breakpoints) + error( + "Maximum power limit $P_max_ft of HVDC Line $(PSY.get_name(d)) has different value of last breakpoint from Loss data $(last(p_breakpoints)).", + ) + end + ### Update Params Vectors ### + ## Update from 1 to S + for i in 1:len_segments + from_to_loss_params[i] = -p_breakpoints[2 + len_segments - i] - p_breakpoints[1] # for i = 1: P_end, for i = len_segments: P_2 + to_from_loss_params[i] = + p_breakpoints[2 + len_segments - i] * (1 - loss_factors[len_segments + 1 - i]) + end + ## Update from S+1 and S+2 + from_to_loss_params[len_segments + 1] = -p_breakpoints[1] # P_send0 + from_to_loss_params[len_segments + 2] = 0.0 + to_from_loss_params[len_segments + 1] = 0.0 + to_from_loss_params[len_segments + 2] = -p_breakpoints[1] # P_send0 + ## Update from S+3 to 2S+2 + for i in 1:len_segments + from_to_loss_params[2 + len_segments + i] = + p_breakpoints[i + 1] * (1 - loss_factors[i]) + to_from_loss_params[2 + len_segments + i] = -p_breakpoints[i + 1] - p_breakpoints[1] + end + + return from_to_loss_params, to_from_loss_params +end + +function add_constraints!( + container::OptimizationContainer, + ::Type{T}, + devices::Union{Vector{U}, IS.FlattenIteratorWrapper{U}}, + ::DeviceModel{U, HVDCTwoTerminalPiecewiseLoss}, + ::NetworkModel{<:PM.AbstractPowerModel}, +) where {T <: HVDCFlowCalculationConstraint, U <: PSY.TwoTerminalHVDCLine} + var_pwl = get_variable(container, HVDCPiecewiseLossVariable(), U) + var_pwl_bin = get_variable(container, HVDCPiecewiseBinaryLossVariable(), U) + names = PSY.get_name.(devices) + time_steps = get_time_steps(container) + flow_ft = get_variable(container, HVDCActivePowerReceivedFromVariable(), U) + flow_tf = get_variable(container, HVDCActivePowerReceivedToVariable(), U) + + constraint_from_to = + add_constraints_container!(container, T(), U, names, time_steps; meta = "ft") + constraint_to_from = + add_constraints_container!(container, T(), U, names, time_steps; meta = "tf") + constraint_binary = + add_constraints_container!(container, T(), U, names, time_steps; meta = "bin") + for d in devices + name = PSY.get_name(d) + loss = PSY.get_loss(d) + from_to_params, to_from_params = _get_pwl_loss_params(d, loss) + #@show from_to_params + #@show to_from_params + range_segments = 1:(length(from_to_params) - 1) # 1:(2S+1) + for t in time_steps + ## Add Equality Constraints ## + constraint_from_to[name, t] = JuMP.@constraint( + get_jump_model(container), + flow_ft[name, t] == + sum( + var_pwl_bin[name, ix, t] * from_to_params[ix] for + ix in range_segments + ) + sum( + var_pwl[name, ix, t] * (from_to_params[ix + 1] - from_to_params[ix]) for + ix in range_segments + ) + ) + constraint_to_from[name, t] = JuMP.@constraint( + get_jump_model(container), + flow_tf[name, t] == + sum( + var_pwl_bin[name, ix, t] * to_from_params[ix] for + ix in range_segments + ) + sum( + var_pwl[name, ix, t] * (to_from_params[ix + 1] - to_from_params[ix]) for + ix in range_segments + ) + ) + ## Add Binary Bound ### + constraint_binary[name, t] = JuMP.@constraint( + get_jump_model(container), + sum(var_pwl_bin[name, ix, t] for ix in range_segments) == 1.0 + ) + ## Add Bounds for Continuous ## + for ix in range_segments + JuMP.@constraint( + get_jump_model(container), + var_pwl[name, ix, t] <= var_pwl_bin[name, ix, t] + ) + if ix == div(length(range_segments) + 1, 2) + JuMP.fix(var_pwl[name, ix, t], 0.0; force = true) + end + end + end + end + return +end + #################################### Rate Limits Constraints ################################################## function _get_flow_bounds(d::PSY.TwoTerminalHVDCLine) check_hvdc_line_limits_consistency(d) @@ -246,7 +447,12 @@ end function _add_hvdc_flow_constraints!( container::OptimizationContainer, devices::Union{Vector{T}, IS.FlattenIteratorWrapper{T}}, - var::Union{FlowActivePowerFromToVariable, FlowActivePowerToFromVariable}, + var::Union{ + FlowActivePowerFromToVariable, + FlowActivePowerToFromVariable, + HVDCActivePowerReceivedFromVariable, + HVDCActivePowerReceivedToVariable, + }, constraint::Union{FlowRateConstraintFromTo, FlowRateConstraintToFrom}, ) where {T <: PSY.TwoTerminalHVDCLine} time_steps = get_time_steps(container) @@ -305,8 +511,10 @@ function add_constraints!( devices::IS.FlattenIteratorWrapper{U}, ::DeviceModel{U, HVDCTwoTerminalDispatch}, ::NetworkModel{<:PM.AbstractDCPModel}, -) where {T <: Union{FlowRateConstraintToFrom, FlowRateConstraintFromTo}, - U <: PSY.TwoTerminalHVDCLine} +) where { + T <: Union{FlowRateConstraintToFrom, FlowRateConstraintFromTo}, + U <: PSY.TwoTerminalHVDCLine, +} _add_hvdc_flow_constraints!(container, devices, T()) return end @@ -316,49 +524,83 @@ function add_constraints!( ::Type{T}, devices::IS.FlattenIteratorWrapper{U}, ::DeviceModel{U, HVDCTwoTerminalDispatch}, - ::NetworkModel{<:PM.AbstractDCPModel}, -) where {T <: HVDCDirection, U <: PSY.TwoTerminalHVDCLine} - time_steps = get_time_steps(container) - names = [PSY.get_name(d) for d in devices] + ::NetworkModel{<:AbstractPTDFModel}, +) where { + T <: Union{FlowRateConstraintToFrom, FlowRateConstraintFromTo}, + U <: PSY.TwoTerminalHVDCLine, +} + _add_hvdc_flow_constraints!(container, devices, T()) + return +end - tf_var = get_variable(container, FlowActivePowerToFromVariable(), U) - ft_var = get_variable(container, FlowActivePowerFromToVariable(), U) - direction_var = get_variable(container, HVDCFlowDirectionVariable(), U) - - constraint_ft_ub = - add_constraints_container!(container, T(), U, names, time_steps; meta = "ft_ub") - constraint_tf_ub = - add_constraints_container!(container, T(), U, names, time_steps; meta = "tf_ub") - constraint_ft_lb = - add_constraints_container!(container, T(), U, names, time_steps; meta = "ft_lb") - constraint_tf_lb = - add_constraints_container!(container, T(), U, names, time_steps; meta = "tf_lb") +function add_constraints!( + container::OptimizationContainer, + ::Type{T}, + devices::IS.FlattenIteratorWrapper{U}, + model::DeviceModel{U, V}, + network_model::NetworkModel{CopperPlatePowerModel}, +) where { + T <: Union{FlowRateConstraintFromTo, FlowRateConstraintToFrom}, + U <: PSY.TwoTerminalHVDCLine, + V <: HVDCTwoTerminalPiecewiseLoss, +} + inter_network_branches = U[] for d in devices - min_rate_to, max_rate_to = PSY.get_active_power_limits_to(d) - min_rate_from, max_rate_from = PSY.get_active_power_limits_to(d) - name = PSY.get_name(d) - for t in time_steps - constraint_tf_ub[name, t] = JuMP.@constraint( - get_jump_model(container), - tf_var[name, t] <= max_rate_to * (1 - direction_var[name, t]) - ) - constraint_ft_ub[name, t] = JuMP.@constraint( - get_jump_model(container), - ft_var[name, t] <= max_rate_from * (1 - direction_var[name, t]) - ) - constraint_tf_lb[name, t] = JuMP.@constraint( - get_jump_model(container), - direction_var[name, t] * min_rate_to <= tf_var[name, t] + ref_bus_from = get_reference_bus(network_model, PSY.get_arc(d).from) + ref_bus_to = get_reference_bus(network_model, PSY.get_arc(d).to) + if ref_bus_from != ref_bus_to + push!(inter_network_branches, d) + end + end + if !isempty(inter_network_branches) + if T <: FlowRateConstraintFromTo + _add_hvdc_flow_constraints!( + container, + devices, + HVDCActivePowerReceivedFromVariable(), + T(), ) - constraint_ft_lb[name, t] = JuMP.@constraint( - get_jump_model(container), - direction_var[name, t] * min_rate_from <= tf_var[name, t] + else + _add_hvdc_flow_constraints!( + container, + devices, + HVDCActivePowerReceivedToVariable(), + T(), ) end end return end +function add_constraints!( + container::OptimizationContainer, + ::Type{T}, + devices::IS.FlattenIteratorWrapper{U}, + ::DeviceModel{U, V}, + ::NetworkModel{<:AbstractPTDFModel}, +) where { + T <: Union{FlowRateConstraintFromTo, FlowRateConstraintToFrom}, + U <: PSY.TwoTerminalHVDCLine, + V <: HVDCTwoTerminalPiecewiseLoss, +} + if T <: FlowRateConstraintFromTo + _add_hvdc_flow_constraints!( + container, + devices, + HVDCActivePowerReceivedFromVariable(), + T(), + ) + else + _add_hvdc_flow_constraints!( + container, + devices, + HVDCActivePowerReceivedToVariable(), + T(), + ) + end + return +end + function add_constraints!( container::OptimizationContainer, ::Type{HVDCPowerBalance}, @@ -371,6 +613,7 @@ function add_constraints!( tf_var = get_variable(container, FlowActivePowerToFromVariable(), T) ft_var = get_variable(container, FlowActivePowerFromToVariable(), T) direction_var = get_variable(container, HVDCFlowDirectionVariable(), T) + losses = get_variable(container, HVDCLosses(), T) constraint_ft_ub = add_constraints_container!( container, @@ -404,82 +647,72 @@ function add_constraints!( time_steps; meta = "ft_lb", ) - for d in devices - l1 = PSY.get_loss(d).l1 - l0 = PSY.get_loss(d).l0 - name = PSY.get_name(d) - for t in get_time_steps(container) - if l1 == 0.0 && l0 == 0.0 - constraint_tf_ub[PSY.get_name(d), t] = JuMP.@constraint( - get_jump_model(container), - tf_var[name, t] - ft_var[name, t] == 0.0 - ) - constraint_ft_ub[PSY.get_name(d), t] = JuMP.@constraint( - get_jump_model(container), - ft_var[name, t] - tf_var[name, t] == 0.0 - ) - else - constraint_tf_ub[PSY.get_name(d), t] = JuMP.@constraint( - get_jump_model(container), - tf_var[name, t] - ft_var[name, t] <= l1 * tf_var[name, t] - l0 - ) - constraint_ft_ub[PSY.get_name(d), t] = JuMP.@constraint( - get_jump_model(container), - ft_var[name, t] - tf_var[name, t] >= l1 * ft_var[name, t] + l0 - ) - end - constraint_tf_lb[PSY.get_name(d), t] = JuMP.@constraint( - get_jump_model(container), - ft_var[name, t] - tf_var[name, t] >= - -M_VALUE * (1 - direction_var[name, t]) - ) - constraint_ft_lb[PSY.get_name(d), t] = JuMP.@constraint( - get_jump_model(container), - tf_var[name, t] - ft_var[name, t] >= -M_VALUE * (direction_var[name, t]) - ) - end - end - return -end - -function add_constraints!( - container::OptimizationContainer, - ::Type{HVDCLossesAbsoluteValue}, - devices::IS.FlattenIteratorWrapper{T}, - ::DeviceModel{T, <:AbstractTwoTerminalDCLineFormulation}, - ::NetworkModel{<:PM.AbstractDCPModel}, -) where {T <: PSY.TwoTerminalHVDCLine} - time_steps = get_time_steps(container) - names = [PSY.get_name(d) for d in devices] - losses = get_variable(container, HVDCLosses(), T) - tf_var = get_variable(container, FlowActivePowerToFromVariable(), T) - ft_var = get_variable(container, FlowActivePowerFromToVariable(), T) - constraint_tf = add_constraints_container!( + constraint_loss = add_constraints_container!( container, - HVDCLossesAbsoluteValue(), + HVDCPowerBalance(), T, names, time_steps; - meta = "tf", + meta = "loss", ) - constraint_ft = add_constraints_container!( + constraint_loss_aux1 = add_constraints_container!( container, - HVDCLossesAbsoluteValue(), + HVDCPowerBalance(), T, names, time_steps; - meta = "ft", + meta = "loss_aux1", + ) + constraint_loss_aux2 = add_constraints_container!( + container, + HVDCPowerBalance(), + T, + names, + time_steps; + meta = "loss_aux2", ) for d in devices + loss = PSY.get_loss(d) + if !isa(loss, PSY.LinearCurve) + error( + "HVDCTwoTerminalDispatch of branch $(PSY.get_name(d)) only accepts LinearCurve for loss models.", + ) + end + l1 = PSY.get_proportional_term(loss) + l0 = PSY.get_constant_term(loss) name = PSY.get_name(d) + R_min_from, R_max_from = PSY.get_active_power_limits_from(d) + R_min_to, R_max_to = PSY.get_active_power_limits_to(d) for t in get_time_steps(container) - constraint_tf[PSY.get_name(d), t] = JuMP.@constraint( + constraint_tf_ub[PSY.get_name(d), t] = JuMP.@constraint( + get_jump_model(container), + tf_var[name, t] <= R_max_to * direction_var[name, t] + ) + constraint_tf_lb[PSY.get_name(d), t] = JuMP.@constraint( + get_jump_model(container), + tf_var[name, t] >= R_min_to * (1 - direction_var[name, t]) + ) + constraint_ft_ub[PSY.get_name(d), t] = JuMP.@constraint( + get_jump_model(container), + ft_var[name, t] <= R_max_from * (1 - direction_var[name, t]) + ) + constraint_ft_lb[PSY.get_name(d), t] = JuMP.@constraint( + get_jump_model(container), + ft_var[name, t] >= R_min_from * direction_var[name, t] + ) + constraint_loss[PSY.get_name(d), t] = JuMP.@constraint( + get_jump_model(container), + tf_var[name, t] + ft_var[name, t] == losses[name, t] + ) + constraint_loss_aux1[PSY.get_name(d), t] = JuMP.@constraint( get_jump_model(container), - tf_var[name, t] - ft_var[name, t] <= losses[name, t] + losses[name, t] >= + l1 * ft_var[name, t] + l0 - M_VALUE * direction_var[name, t] ) - constraint_ft[PSY.get_name(d), t] = JuMP.@constraint( + constraint_loss_aux2[PSY.get_name(d), t] = JuMP.@constraint( get_jump_model(container), - -tf_var[name, t] + ft_var[name, t] <= losses[name, t] + losses[name, t] >= + l1 * tf_var[name, t] + l0 - M_VALUE * (1 - direction_var[name, t]) ) end end diff --git a/src/devices_models/devices/common/add_to_expression.jl b/src/devices_models/devices/common/add_to_expression.jl index efbda4845c..9c53a40e0e 100644 --- a/src/devices_models/devices/common/add_to_expression.jl +++ b/src/devices_models/devices/common/add_to_expression.jl @@ -260,13 +260,17 @@ function add_to_expression!( for d in devices name = PSY.get_name(d) device_bus_from = PSY.get_arc(d).from + device_bus_to = PSY.get_arc(d).to ref_bus_from = get_reference_bus(network_model, device_bus_from) - for t in get_time_steps(container) - _add_to_jump_expression!( - expression[ref_bus_from, t], - variable[name, t], - get_variable_multiplier(U(), d, W()), - ) + ref_bus_to = get_reference_bus(network_model, device_bus_to) + if ref_bus_from == ref_bus_to + for t in get_time_steps(container) + _add_to_jump_expression!( + expression[ref_bus_from, t], + variable[name, t], + get_variable_multiplier(U(), d, W()), + ) + end end end return @@ -329,9 +333,9 @@ function add_to_expression!( ref_bus_to = get_reference_bus(network_model, PSY.get_arc(d).to) for t in get_time_steps(container) flow_variable = var[PSY.get_name(d), t] - _add_to_jump_expression!(nodal_expr[bus_no_to, t], flow_variable, 1.0) + _add_to_jump_expression!(nodal_expr[bus_no_to, t], flow_variable, -1.0) if ref_bus_from != ref_bus_to - _add_to_jump_expression!(sys_expr[ref_bus_to, t], flow_variable, 1.0) + _add_to_jump_expression!(sys_expr[ref_bus_to, t], flow_variable, -1.0) end end end @@ -375,6 +379,80 @@ function add_to_expression!( return end +""" +PWL implementation to add FromTo branch variables to SystemBalanceExpressions +""" +function add_to_expression!( + container::OptimizationContainer, + ::Type{T}, + ::Type{U}, + devices::IS.FlattenIteratorWrapper{V}, + ::DeviceModel{V, W}, + network_model::NetworkModel{X}, +) where { + T <: ActivePowerBalance, + U <: HVDCActivePowerReceivedFromVariable, + V <: TwoTerminalHVDCTypes, + W <: HVDCTwoTerminalPiecewiseLoss, + X <: AbstractPTDFModel, +} + var = get_variable(container, U(), V) + nodal_expr = get_expression(container, T(), PSY.ACBus) + sys_expr = get_expression(container, T(), _system_expression_type(X)) + radial_network_reduction = get_radial_network_reduction(network_model) + for d in devices + bus_no_from = + PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).from) + ref_bus_to = get_reference_bus(network_model, PSY.get_arc(d).to) + ref_bus_from = get_reference_bus(network_model, PSY.get_arc(d).from) + for t in get_time_steps(container) + flow_variable = var[PSY.get_name(d), t] + _add_to_jump_expression!(nodal_expr[bus_no_from, t], flow_variable, 1.0) + if ref_bus_from != ref_bus_to + _add_to_jump_expression!(sys_expr[ref_bus_from, t], flow_variable, 1.0) + end + end + end + return +end + +""" +PWL implementation to add FromTo branch variables to SystemBalanceExpressions +""" +function add_to_expression!( + container::OptimizationContainer, + ::Type{T}, + ::Type{U}, + devices::IS.FlattenIteratorWrapper{V}, + ::DeviceModel{V, W}, + network_model::NetworkModel{X}, +) where { + T <: ActivePowerBalance, + U <: HVDCActivePowerReceivedToVariable, + V <: TwoTerminalHVDCTypes, + W <: HVDCTwoTerminalPiecewiseLoss, + X <: AbstractPTDFModel, +} + var = get_variable(container, U(), V) + nodal_expr = get_expression(container, T(), PSY.ACBus) + sys_expr = get_expression(container, T(), _system_expression_type(X)) + radial_network_reduction = get_radial_network_reduction(network_model) + for d in devices + bus_no_to = + PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_arc(d).to) + ref_bus_to = get_reference_bus(network_model, PSY.get_arc(d).to) + ref_bus_from = get_reference_bus(network_model, PSY.get_arc(d).from) + for t in get_time_steps(container) + flow_variable = var[PSY.get_name(d), t] + _add_to_jump_expression!(nodal_expr[bus_no_to, t], flow_variable, 1.0) + if ref_bus_from != ref_bus_to + _add_to_jump_expression!(sys_expr[ref_bus_to, t], flow_variable, 1.0) + end + end + end + return +end + """ Default implementation to add branch variables to SystemBalanceExpressions """ @@ -530,15 +608,24 @@ function add_to_expression!( variable = get_variable(container, U(), V) expression = get_expression(container, T(), PSY.ACBus) radial_network_reduction = get_radial_network_reduction(network_model) - for d in devices, t in get_time_steps(container) + for d in devices name = PSY.get_name(d) bus_no_ = PSY.get_number(PSY.get_bus(d)) bus_no = PNM.get_mapped_bus_number(radial_network_reduction, bus_no_) - _add_to_jump_expression!( - expression[bus_no, t], - variable[name, t], - get_variable_multiplier(U(), d, W()), - ) + for t in get_time_steps(container) + if PSY.get_must_run(d) + _add_to_jump_expression!( + expression[bus_no, t], + get_variable_multiplier(U(), d, W()), + ) + else + _add_to_jump_expression!( + expression[bus_no, t], + variable[name, t], + get_variable_multiplier(U(), d, W()), + ) + end + end end return end @@ -558,15 +645,24 @@ function add_to_expression!( } variable = get_variable(container, U(), V) expression = get_expression(container, T(), PSY.ACBus) - for d in devices, t in get_time_steps(container) + for d in devices + name = PSY.get_name(d) bus = PSY.get_bus(d) area_name = PSY.get_name(PSY.get_area(bus)) - name = PSY.get_name(d) - _add_to_jump_expression!( - expression[area_name, t], - variable[name, t], - get_variable_multiplier(U(), d, W()), - ) + for t in get_time_steps(container) + if PSY.get_must_run(d) + _add_to_jump_expression!( + expression[area_name, t], + get_variable_multiplier(U(), d, W()), + ) + else + _add_to_jump_expression!( + expression[area_name, t], + variable[name, t], + get_variable_multiplier(U(), d, W()), + ) + end + end end return end @@ -1251,6 +1347,9 @@ function add_to_expression!( end expression = get_expression(container, T(), V) for d in devices + if PSY.get_must_run(d) + continue + end mult = get_expression_multiplier(U(), T(), d, W()) for t in get_time_steps(container) name = PSY.get_name(d) @@ -1401,7 +1500,7 @@ end function add_to_expression!( container::OptimizationContainer, ::Type{S}, - cost_expression::JuMP.AbstractJuMPScalar, + cost_expression::Union{JuMP.AbstractJuMPScalar, Float64}, component::T, time_period::Int, ) where {S <: CostExpressions, T <: PSY.Component} diff --git a/src/devices_models/devices/common/duration_constraints.jl b/src/devices_models/devices/common/duration_constraints.jl index 954a7e2aa5..12a1f8d7b8 100644 --- a/src/devices_models/devices/common/duration_constraints.jl +++ b/src/devices_models/devices/common/duration_constraints.jl @@ -52,7 +52,10 @@ function device_duration_retrospective!( varstart = get_variable(container, var_types[2], T) varstop = get_variable(container, var_types[3], T) - set_names = [get_component_name(ic) for ic in initial_duration[:, 1]] + set_names = [ + get_component_name(ic) for + ic in initial_duration[:, 1] if !isnothing(get_value(ic)) + ] con_up = add_constraints_container!( container, cons_type, @@ -72,6 +75,7 @@ function device_duration_retrospective!( for t in time_steps for (ix, ic) in enumerate(initial_duration[:, 1]) + isnothing(get_value(ic)) && continue name = get_component_name(ic) # Minimum Up-time Constraint lhs_on = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}(0) @@ -84,10 +88,11 @@ function device_duration_retrospective!( JuMP.add_to_expression!(lhs_on, 1) end con_up[name, t] = - JuMP.@constraint(container.JuMPmodel, lhs_on - varon[name, t] <= 0.0) + JuMP.@constraint(get_jump_model(container), lhs_on - varon[name, t] <= 0.0) end for (ix, ic) in enumerate(initial_duration[:, 2]) + isnothing(get_value(ic)) && continue name = get_component_name(ic) # Minimum Down-time Constraint lhs_off = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}(0) @@ -100,11 +105,12 @@ function device_duration_retrospective!( JuMP.add_to_expression!(lhs_off, 1) end con_down[name, t] = - JuMP.@constraint(container.JuMPmodel, lhs_off + varon[name, t] <= 1.0) + JuMP.@constraint(get_jump_model(container), lhs_off + varon[name, t] <= 1.0) end end return end + @doc raw""" This formulation of the duration constraints looks ahead in the time frame of the model. @@ -165,6 +171,7 @@ function device_duration_look_ahead!( for t in time_steps for (ix, ic) in enumerate(initial_duration[:, 1]) + isnothing(get_value(ic)) && continue name = get_component_name(ic) # Minimum Up-time Constraint lhs_on = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}(0) @@ -177,12 +184,13 @@ function device_duration_look_ahead!( lhs_on += get_value(ic) end con_up[name, t] = JuMP.@constraint( - container.JuMPmodel, + get_jump_model(container), varstop[name, t] * duration_data[ix].up - lhs_on <= 0.0 ) end for (ix, ic) in enumerate(initial_duration[:, 2]) + isnothing(get_value(ic)) && continue name = get_component_name(ic) # Minimum Down-time Constraint lhs_off = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}(0) @@ -195,7 +203,7 @@ function device_duration_look_ahead!( lhs_off += get_value(ic) end con_down[name, t] = JuMP.@constraint( - container.JuMPmodel, + get_jump_model(container), varstart[name, t] * duration_data[ix].down - lhs_off <= 0.0 ) end @@ -278,8 +286,8 @@ function device_duration_parameters!( for t in time_steps for (ix, ic) in enumerate(initial_duration[:, 1]) + isnothing(get_value(ic)) && continue name = get_component_name(ic) - # Minimum Up-time Constraint lhs_on = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}(0) for i in UnitRange{Int}(Int(t - duration_data[ix].up + 1), t) @@ -294,16 +302,20 @@ function device_duration_parameters!( if t <= duration_data[ix].up lhs_on += get_value(ic) con_up[name, t] = JuMP.@constraint( - container.JuMPmodel, + get_jump_model(container), varstop[name, t] * duration_data[ix].up - lhs_on <= 0.0 ) else con_up[name, t] = - JuMP.@constraint(container.JuMPmodel, lhs_on - varon[name, t] <= 0.0) + JuMP.@constraint( + get_jump_model(container), + lhs_on - varon[name, t] <= 0.0 + ) end end for (ix, ic) in enumerate(initial_duration[:, 2]) + isnothing(get_value(ic)) && continue name = get_component_name(ic) # Minimum Down-time Constraint lhs_off = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}(0) @@ -319,12 +331,15 @@ function device_duration_parameters!( if t <= duration_data[ix].down lhs_off += get_value(ic) con_down[name, t] = JuMP.@constraint( - container.JuMPmodel, + get_jump_model(container), varstart[name, t] * duration_data[ix].down - lhs_off <= 0.0 ) else con_down[name, t] = - JuMP.@constraint(container.JuMPmodel, lhs_off + varon[name, t] <= 1.0) + JuMP.@constraint( + get_jump_model(container), + lhs_off + varon[name, t] <= 1.0 + ) end end end @@ -395,6 +410,7 @@ function device_duration_compact_retrospective!( total_time_steps = length(time_steps) for t in time_steps for (ix, ic) in enumerate(initial_duration[:, 1]) + isnothing(get_value(ic)) && continue name = get_component_name(ic) # Minimum Up-time Constraint lhs_on = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}(0) @@ -413,10 +429,11 @@ function device_duration_compact_retrospective!( continue end con_up[name, t] = - JuMP.@constraint(container.JuMPmodel, lhs_on - varon[name, t] <= 0.0) + JuMP.@constraint(get_jump_model(container), lhs_on - varon[name, t] <= 0.0) end for (ix, ic) in enumerate(initial_duration[:, 2]) + isnothing(get_value(ic)) && continue name = get_component_name(ic) # Minimum Down-time Constraint lhs_off = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}(0) @@ -435,7 +452,7 @@ function device_duration_compact_retrospective!( continue end con_down[name, t] = - JuMP.@constraint(container.JuMPmodel, lhs_off + varon[name, t] <= 1.0) + JuMP.@constraint(get_jump_model(container), lhs_off + varon[name, t] <= 1.0) end end for c in [con_up, con_down] diff --git a/src/devices_models/devices/common/objective_function/common.jl b/src/devices_models/devices/common/objective_function/common.jl index 082cb79932..c02b55a51f 100644 --- a/src/devices_models/devices/common/objective_function/common.jl +++ b/src/devices_models/devices/common/objective_function/common.jl @@ -11,6 +11,7 @@ function add_variable_cost!( for d in devices op_cost_data = PSY.get_operation_cost(d) _add_variable_cost_to_objective!(container, U(), d, op_cost_data, V()) + _add_vom_cost_to_objective!(container, U(), d, op_cost_data, V()) end return end @@ -27,6 +28,7 @@ function add_shut_down_cost!( ) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} multiplier = objective_function_multiplier(U(), V()) for d in devices + PSY.get_must_run(d) && continue op_cost_data = PSY.get_operation_cost(d) cost_term = shut_down_cost(op_cost_data, d, V()) iszero(cost_term) && continue @@ -40,7 +42,6 @@ end ################################## ####### Proportional Cost ######## ################################## - function add_proportional_cost!( container::OptimizationContainer, ::U, @@ -53,12 +54,46 @@ function add_proportional_cost!( cost_term = proportional_cost(op_cost_data, U(), d, V()) iszero(cost_term) && continue for t in get_time_steps(container) - _add_proportional_term!(container, U(), d, cost_term * multiplier, t) + exp = _add_proportional_term!(container, U(), d, cost_term * multiplier, t) + add_to_expression!(container, ProductionCostExpression, exp, d, t) end end return end +################################## +########## VOM Cost ############## +################################## + +function _add_vom_cost_to_objective!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + op_cost::PSY.OperationalCost, + ::U, +) where {T <: VariableType, U <: AbstractDeviceFormulation} + variable_cost_data = variable_cost(op_cost, T(), component, U()) + power_units = PSY.get_power_units(variable_cost_data) + vom_cost = PSY.get_vom_cost(variable_cost_data) + multiplier = 1.0 # VOM Cost is always positive + cost_term = PSY.get_proportional_term(vom_cost) + iszero(cost_term) && return + base_power = get_base_power(container) + device_base_power = PSY.get_base_power(component) + cost_term_normalized = get_proportional_cost_per_system_unit( + cost_term, + power_units, + base_power, + device_base_power, + ) + for t in get_time_steps(container) + exp = + _add_proportional_term!(container, T(), d, cost_term_normalized * multiplier, t) + add_to_expression!(container, ProductionCostExpression, exp, d, t) + end + return +end + ################################## ######## OnVariable Cost ######### ################################## @@ -68,15 +103,24 @@ function add_proportional_cost!( ::U, devices::IS.FlattenIteratorWrapper{T}, ::V, -) where {T <: PSY.ThermalGen, U <: OnVariable, V <: AbstractCompactUnitCommitment} +) where {T <: PSY.ThermalGen, U <: OnVariable, V <: AbstractThermalUnitCommitment} multiplier = objective_function_multiplier(U(), V()) for d in devices op_cost_data = PSY.get_operation_cost(d) cost_term = proportional_cost(op_cost_data, U(), d, V()) iszero(cost_term) && continue for t in get_time_steps(container) - exp = _add_proportional_term!(container, U(), d, cost_term * multiplier, t) - add_to_expression!(container, ProductionCostExpression, exp, d, t) + if !PSY.get_must_run(d) + exp = _add_proportional_term!(container, U(), d, cost_term * multiplier, t) + add_to_expression!(container, ProductionCostExpression, exp, d, t) + end + add_to_expression!( + container, + ProductionCostExpression, + cost_term * multiplier, + d, + t, + ) end end return @@ -114,6 +158,7 @@ function _add_start_up_cost_to_objective!( op_cost::Union{PSY.ThermalGenerationCost, PSY.MarketBidCost}, ::U, ) where {T <: VariableType, U <: AbstractDeviceFormulation} + PSY.get_must_run(component) && return cost_term = start_up_cost(op_cost, component, U()) iszero(cost_term) && return multiplier = objective_function_multiplier(T(), U()) @@ -136,6 +181,7 @@ function _add_start_up_cost_to_objective!( op_cost::PSY.ThermalGenerationCost, ::U, ) where {T <: VariableType, U <: ThermalMultiStartUnitCommitment} + PSY.get_must_run(component) && return cost_terms = start_up_cost(op_cost, component, U()) cost_term = cost_terms[MULTI_START_COST_MAP[T]] iszero(cost_term) && return diff --git a/src/devices_models/devices/common/objective_function/market_bid.jl b/src/devices_models/devices/common/objective_function/market_bid.jl index d46025e62f..635f59e07d 100644 --- a/src/devices_models/devices/common/objective_function/market_bid.jl +++ b/src/devices_models/devices/common/objective_function/market_bid.jl @@ -302,21 +302,6 @@ function _add_pwl_term!( device_base_power, ) - compact_status = validate_compact_pwl_data(component, data, base_power) - if !uses_compact_power(component, V()) && compact_status == COMPACT_PWL_STATUS.VALID - error( - "The data provided is not compatible with formulation $V. Use a formulation compatible with Compact Cost Functions", - ) - # data = _convert_to_full_variable_cost(data, component) - elseif uses_compact_power(component, V()) && compact_status != COMPACT_PWL_STATUS.VALID - @warn( - "The cost data provided is not in compact form. Will attempt to convert. Errors may occur." - ) - data = convert_to_compact_variable_cost(data) - else - @debug uses_compact_power(component, V()) compact_status name T V - end - cost_is_convex = PSY.is_convex(data) if !cost_is_convex error("MarketBidCost for component $(name) is non-convex") diff --git a/src/devices_models/devices/common/objective_function/piecewise_linear.jl b/src/devices_models/devices/common/objective_function/piecewise_linear.jl index bc184fb972..d173cd85d8 100644 --- a/src/devices_models/devices/common/objective_function/piecewise_linear.jl +++ b/src/devices_models/devices/common/objective_function/piecewise_linear.jl @@ -54,6 +54,52 @@ end ################# PWL Constraints ################ ################################################## +function _determine_bin_lhs( + container::OptimizationContainer, + sos_status::SOSStatusVariable, + component::T, + period::Int) where {T <: PSY.Component} + name = PSY.get_name(component) + if sos_status == SOSStatusVariable.NO_VARIABLE + return 1.0 + @debug "Using Piecewise Linear cost function but no variable/parameter ref for ON status is passed. Default status will be set to online (1.0)" _group = + LOG_GROUP_COST_FUNCTIONS + + elseif sos_status == SOSStatusVariable.PARAMETER + param = get_default_on_parameter(component) + return get_parameter(container, param, T).parameter_array[name, period] + @debug "Using Piecewise Linear cost function with parameter OnStatusParameter, $T" _group = + LOG_GROUP_COST_FUNCTIONS + elseif sos_status == SOSStatusVariable.VARIABLE + var = get_default_on_variable(component) + return get_variable(container, var, T)[name, period] + @debug "Using Piecewise Linear cost function with variable OnVariable $T" _group = + LOG_GROUP_COST_FUNCTIONS + else + @assert false + end +end + +function _get_bin_lhs( + container::OptimizationContainer, + sos_status::SOSStatusVariable, + component::T, + period::Int) where {T <: PSY.Component} + return _determine_bin_lhs(container, sos_status, component, period) +end + +function _get_bin_lhs( + container::OptimizationContainer, + sos_status::SOSStatusVariable, + component::PSY.ThermalGen, + period::Int) + if PSY.get_must_run(component) + return 1.0 + else + return _determine_bin_lhs(container, sos_status, component, period) + end +end + """ Implement the constraints for PWL variables. That is: @@ -86,26 +132,7 @@ function _add_pwl_constraint!( variables[name, period] == sum(pwl_vars[name, ix, period] * break_points[ix] for ix in 1:len_cost_data) ) - - if sos_status == SOSStatusVariable.NO_VARIABLE - bin = 1.0 - @debug "Using Piecewise Linear cost function but no variable/parameter ref for ON status is passed. Default status will be set to online (1.0)" _group = - LOG_GROUP_COST_FUNCTIONS - - elseif sos_status == SOSStatusVariable.PARAMETER - param = get_default_on_parameter(component) - bin = get_parameter(container, param, T).parameter_array[name, period] - @debug "Using Piecewise Linear cost function with parameter OnStatusParameter, $T" _group = - LOG_GROUP_COST_FUNCTIONS - elseif sos_status == SOSStatusVariable.VARIABLE - var = get_default_on_variable(component) - bin = get_variable(container, var, T)[name, period] - @debug "Using Piecewise Linear cost function with variable OnVariable $T" _group = - LOG_GROUP_COST_FUNCTIONS - else - @assert false - end - + bin = _get_bin_lhs(container, sos_status, component, period) const_normalization_container = lazy_container_addition!( container, PieceWiseLinearCostConstraint(), diff --git a/src/devices_models/devices/common/objective_function/quadratic_curve.jl b/src/devices_models/devices/common/objective_function/quadratic_curve.jl index 4f5a8aea09..686dc0d236 100644 --- a/src/devices_models/devices/common/objective_function/quadratic_curve.jl +++ b/src/devices_models/devices/common/objective_function/quadratic_curve.jl @@ -45,7 +45,15 @@ function _add_quadraticcurve_variable_cost!( proportional_term_per_unit::Vector{Float64}, quadratic_term_per_unit::Vector{Float64}, ) where {T <: VariableType} + lb, ub = PSY.get_active_power_limits(component) for t in get_time_steps(container) + _check_quadratic_monotonicity( + PSY.get_name(component), + quadratic_term_per_unit[t], + proportional_term_per_unit[t], + lb, + ub, + ) _add_quadraticcurve_variable_term_to_model!( container, T(), @@ -66,6 +74,13 @@ function _add_quadraticcurve_variable_cost!( proportional_term_per_unit::Float64, quadratic_term_per_unit::Float64, ) where {T <: VariableType} + lb, ub = PSY.get_active_power_limits(component) + _check_quadratic_monotonicity(PSY.get_name(component), + quadratic_term_per_unit, + proportional_term_per_unit, + lb, + ub, + ) for t in get_time_steps(container) _add_quadraticcurve_variable_term_to_model!( container, @@ -79,6 +94,23 @@ function _add_quadraticcurve_variable_cost!( return end +function _check_quadratic_monotonicity( + name::String, + quad_term::Float64, + linear_term::Float64, + lb::Float64, + ub::Float64, +) + fp_lb = 2 * quad_term * lb + linear_term + fp_ub = 2 * quad_term * ub + linear_term + + if fp_lb < 0 || fp_ub < 0 + @warn "Cost function for component $name is not monotonically increasing in the range [$lb, $ub]. \ + This can lead to unexpected results" + end + return +end + @doc raw""" Adds to the cost function cost terms for sum of variables with common factor to be used for cost expression for optimization_container model. @@ -149,7 +181,7 @@ function _add_variable_cost_to_objective!( } throw( IS.ConflictingInputsError( - "Quadratic Cost Curves are not allowed for Compact formulations", + "Quadratic Cost Curves are not compatible with Compact formulations", ), ) return diff --git a/src/devices_models/devices/common/range_constraint.jl b/src/devices_models/devices/common/range_constraint.jl index a4e0566066..0a7c4b804c 100644 --- a/src/devices_models/devices/common/range_constraint.jl +++ b/src/devices_models/devices/common/range_constraint.jl @@ -104,7 +104,7 @@ function _add_lower_bound_range_constraints_impl!( ci_name = PSY.get_name(device) limits = get_min_max_limits(device, T, W) # depends on constraint type and formulation type con_lb[ci_name, t] = - JuMP.@constraint(container.JuMPmodel, array[ci_name, t] >= limits.min) + JuMP.@constraint(get_jump_model(container), array[ci_name, t] >= limits.min) end return end @@ -126,7 +126,7 @@ function _add_upper_bound_range_constraints_impl!( ci_name = PSY.get_name(device) limits = get_min_max_limits(device, T, W) # depends on constraint type and formulation type con_ub[ci_name, t] = - JuMP.@constraint(container.JuMPmodel, array[ci_name, t] <= limits.max) + JuMP.@constraint(get_jump_model(container), array[ci_name, t] <= limits.max) end return end @@ -246,20 +246,86 @@ function _add_semicontinuous_lower_bound_range_constraints_impl!( ) where {T <: ConstraintType, V <: PSY.Component, W <: AbstractDeviceFormulation} time_steps = get_time_steps(container) names = [PSY.get_name(d) for d in devices] - binary_variables = [OnVariable()] + con_lb = add_constraints_container!(container, T(), V, names, time_steps; meta = "lb") + varbin = get_variable(container, OnVariable(), V) + + for device in devices + ci_name = PSY.get_name(device) + limits = get_min_max_limits(device, T, W) # depends on constraint type and formulation type + for t in time_steps + con_lb[ci_name, t] = JuMP.@constraint( + get_jump_model(container), + array[ci_name, t] >= limits.min * varbin[ci_name, t] + ) + end + end + return +end +function _add_semicontinuous_lower_bound_range_constraints_impl!( + container::OptimizationContainer, + ::Type{T}, + array, + devices::IS.FlattenIteratorWrapper{V}, + ::DeviceModel{V, W}, +) where {T <: ConstraintType, V <: PSY.ThermalGen, W <: AbstractDeviceFormulation} + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] con_lb = add_constraints_container!(container, T(), V, names, time_steps; meta = "lb") + varbin = get_variable(container, OnVariable(), V) - @assert length(binary_variables) == 1 "Expected $(binary_variables) for $U $V $T $W to be length 1" - varbin = get_variable(container, only(binary_variables), V) + for device in devices + ci_name = PSY.get_name(device) + limits = get_min_max_limits(device, T, W) # depends on constraint type and formulation type + if PSY.get_must_run(device) + for t in time_steps + con_lb[ci_name, t] = JuMP.@constraint( + get_jump_model(container), + array[ci_name, t] >= limits.min + ) + end + else + for t in time_steps + con_lb[ci_name, t] = JuMP.@constraint( + get_jump_model(container), + array[ci_name, t] >= limits.min * varbin[ci_name, t] + ) + end + end + end + return +end - for device in devices, t in time_steps +function _add_semicontinuous_upper_bound_range_constraints_impl!( + container::OptimizationContainer, + ::Type{T}, + array, + devices::IS.FlattenIteratorWrapper{V}, + model::DeviceModel{V, W}, +) where {T <: ConstraintType, V <: PSY.ThermalGen, W <: AbstractDeviceFormulation} + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + con_ub = add_constraints_container!(container, T(), V, names, time_steps; meta = "ub") + varbin = get_variable(container, OnVariable(), V) + + for device in devices ci_name = PSY.get_name(device) limits = get_min_max_limits(device, T, W) # depends on constraint type and formulation type - con_lb[ci_name, t] = JuMP.@constraint( - container.JuMPmodel, - array[ci_name, t] >= limits.min * varbin[ci_name, t] - ) + if PSY.get_must_run(device) + for t in time_steps + con_ub[ci_name, t] = JuMP.@constraint( + get_jump_model(container), + array[ci_name, t] <= limits.max + ) + end + else + for t in time_steps + con_ub[ci_name, t] = JuMP.@constraint( + get_jump_model(container), + array[ci_name, t] <= limits.max * varbin[ci_name, t] + ) + end + end end return end @@ -273,18 +339,14 @@ function _add_semicontinuous_upper_bound_range_constraints_impl!( ) where {T <: ConstraintType, V <: PSY.Component, W <: AbstractDeviceFormulation} time_steps = get_time_steps(container) names = [PSY.get_name(d) for d in devices] - binary_variables = [OnVariable()] - con_ub = add_constraints_container!(container, T(), V, names, time_steps; meta = "ub") - - @assert length(binary_variables) == 1 "Expected $(binary_variables) for $U $V $T $W to be length 1" - varbin = get_variable(container, only(binary_variables), V) + varbin = get_variable(container, OnVariable(), V) for device in devices, t in time_steps ci_name = PSY.get_name(device) limits = get_min_max_limits(device, T, W) # depends on constraint type and formulation type con_ub[ci_name, t] = JuMP.@constraint( - container.JuMPmodel, + get_jump_model(container), array[ci_name, t] <= limits.max * varbin[ci_name, t] ) end @@ -375,7 +437,7 @@ function _add_reserve_lower_bound_range_constraints_impl!( ci_name = PSY.get_name(device) limits = get_min_max_limits(device, T, W) con_lb[ci_name, t] = JuMP.@constraint( - container.JuMPmodel, + get_jump_model(container), array[ci_name, t] >= limits.min * (1 - varbin[ci_name, t]) ) end @@ -409,7 +471,7 @@ function _add_reserve_upper_bound_range_constraints_impl!( ci_name = PSY.get_name(device) limits = get_min_max_limits(device, T, W) con_ub[ci_name, t] = JuMP.@constraint( - container.JuMPmodel, + get_jump_model(container), array[ci_name, t] <= limits.max * (1 - varbin[ci_name, t]) ) end @@ -513,7 +575,7 @@ function _add_reserve_lower_bound_range_constraints_impl!( ci_name = PSY.get_name(device) limits = get_min_max_limits(device, T, X) # depends on constraint type and formulation type con_lb[ci_name, t] = JuMP.@constraint( - container.JuMPmodel, + get_jump_model(container), array[ci_name, t] >= limits.min * varbin[ci_name, t] ) end @@ -545,7 +607,7 @@ function _add_reserve_upper_bound_range_constraints_impl!( ci_name = PSY.get_name(device) limits = get_min_max_limits(device, T, X) # depends on constraint type and formulation type con_ub[ci_name, t] = JuMP.@constraint( - container.JuMPmodel, + get_jump_model(container), array[ci_name, t] <= limits.max * varbin[ci_name, t] ) end diff --git a/src/devices_models/devices/common/rateofchange_constraints.jl b/src/devices_models/devices/common/rateofchange_constraints.jl index b5558b56b1..bca3b5ed5f 100644 --- a/src/devices_models/devices/common/rateofchange_constraints.jl +++ b/src/devices_models/devices/common/rateofchange_constraints.jl @@ -195,14 +195,13 @@ function add_semicontinuous_ramp_constraints!( T::Type{<:ConstraintType}, U::Type{S}, devices::IS.FlattenIteratorWrapper{V}, - model::DeviceModel{V, W}, - X::Type{<:PM.AbstractPowerModel}, + ::DeviceModel{V, W}, + ::Type{<:PM.AbstractPowerModel}, ) where { S <: Union{PowerAboveMinimumVariable, ActivePowerVariable}, V <: PSY.Component, W <: AbstractDeviceFormulation, } - parameters = built_for_recurrent_solves(container) time_steps = get_time_steps(container) variable = get_variable(container, U(), V) varstart = get_variable(container, StartVariable(), V) @@ -222,6 +221,7 @@ function add_semicontinuous_ramp_constraints!( add_constraints_container!(container, T(), V, set_name, time_steps; meta = "dn") for ic in initial_conditions_power + component = get_component(ic) name = get_component_name(ic) # This is to filter out devices that dont need a ramping constraint name ∉ set_name && continue @@ -231,26 +231,58 @@ function add_semicontinuous_ramp_constraints!( ic_power = get_value(ic) @debug "add rate_of_change_constraint" name ic_power + if hasmethod(PSY.get_must_run, Tuple{V}) + must_run = PSY.get_must_run(component) + else + must_run = false + end + + if must_run + rhs_up = ramp_limits.up * minutes_per_period + rhd_dn = ramp_limits.down * minutes_per_period + else + rhs_up = + ramp_limits.up * minutes_per_period + power_limits.min * varstart[name, 1] + rhd_dn = + ramp_limits.down * minutes_per_period + power_limits.min * varstop[name, 1] + end + con_up[name, 1] = JuMP.@constraint( get_jump_model(container), expr_up[name, 1] - ic_power <= - ramp_limits.up * minutes_per_period + power_limits.min * varstart[name, 1] + if must_run + ramp_limits.up * minutes_per_period + else + ramp_limits.up * minutes_per_period + power_limits.min * varstart[name, 1] + end ) con_down[name, 1] = JuMP.@constraint( get_jump_model(container), ic_power - expr_dn[name, 1] <= - ramp_limits.down * minutes_per_period + power_limits.min * varstop[name, 1] + if must_run + ramp_limits.down * minutes_per_period + else + ramp_limits.down * minutes_per_period + power_limits.min * varstop[name, 1] + end ) for t in time_steps[2:end] con_up[name, t] = JuMP.@constraint( get_jump_model(container), expr_up[name, t] - variable[name, t - 1] <= - ramp_limits.up * minutes_per_period + power_limits.min * varstart[name, t] + if must_run + ramp_limits.up * minutes_per_period + else + ramp_limits.up * minutes_per_period + power_limits.min * varstart[name, t] + end ) con_down[name, t] = JuMP.@constraint( get_jump_model(container), variable[name, t - 1] - expr_dn[name, t] <= - ramp_limits.down * minutes_per_period + power_limits.min * varstop[name, t] + if must_run + ramp_limits.down * minutes_per_period + else + ramp_limits.down * minutes_per_period + power_limits.min * varstop[name, t] + end ) end end diff --git a/src/devices_models/devices/thermal_generation.jl b/src/devices_models/devices/thermal_generation.jl index 28b12488da..1c99149f83 100644 --- a/src/devices_models/devices/thermal_generation.jl +++ b/src/devices_models/devices/thermal_generation.jl @@ -61,7 +61,7 @@ get_expression_multiplier(::OnStatusParameter, ::ActivePowerRangeExpressionLB, d get_expression_multiplier(::OnStatusParameter, ::ActivePowerBalance, d::PSY.ThermalGen, ::AbstractThermalFormulation) = PSY.get_active_power_limits(d).min #################### Initial Conditions for models ############### -initial_condition_default(::DeviceStatus, d::PSY.ThermalGen, ::AbstractThermalFormulation) = PSY.get_status(d) +initial_condition_default(::DeviceStatus, d::PSY.ThermalGen, ::AbstractThermalFormulation) = PSY.get_status(d) ? 1.0 : 0.0 initial_condition_variable(::DeviceStatus, d::PSY.ThermalGen, ::AbstractThermalFormulation) = OnVariable() initial_condition_default(::DevicePower, d::PSY.ThermalGen, ::AbstractThermalFormulation) = PSY.get_active_power(d) initial_condition_variable(::DevicePower, d::PSY.ThermalGen, ::AbstractThermalFormulation) = ActivePowerVariable() @@ -75,8 +75,10 @@ initial_condition_variable(::InitialTimeDurationOff, d::PSY.ThermalGen, ::Abstra ########################Objective Function################################################## # TODO: Decide what is the cost for OnVariable, if fixed or constant term in variable -#proportional_cost(cost::PSY.ThermalGenerationCost, S::OnVariable, T::PSY.ThermalGen, U::AbstractThermalFormulation) = no_load_cost(cost, S, T, U) -proportional_cost(cost::PSY.ThermalGenerationCost, S::OnVariable, T::PSY.ThermalGen, U::AbstractThermalFormulation) = PSY.get_fixed(cost) +function proportional_cost(cost::PSY.ThermalGenerationCost, S::OnVariable, T::PSY.ThermalGen, U::AbstractThermalFormulation) + return onvar_cost(cost, S, T, U) + PSY.get_constant_term(PSY.get_vom_cost(PSY.get_variable(cost))) + PSY.get_fixed(cost) +end + proportional_cost(cost::PSY.MarketBidCost, ::OnVariable, ::PSY.ThermalGen, ::AbstractThermalFormulation) = PSY.get_no_load_cost(cost) has_multistart_variables(::PSY.ThermalGen, ::AbstractThermalFormulation)=false @@ -109,17 +111,16 @@ variable_cost(cost::PSY.OperationalCost, ::PowerAboveMinimumVariable, ::PSY.Ther """ Theoretical Cost at power output zero. Mathematically is the intercept with the y-axis """ -function no_load_cost(cost::PSY.ThermalGenerationCost, S::OnVariable, d::PSY.ThermalGen, U::AbstractThermalFormulation) - return _no_load_cost(PSY.get_variable(cost), d) +function onvar_cost(cost::PSY.ThermalGenerationCost, S::OnVariable, d::PSY.ThermalGen, U::AbstractThermalFormulation) + return _onvar_cost(PSY.get_variable(cost), d) end -function _no_load_cost(cost_function::PSY.CostCurve{PSY.PiecewisePointCurve}, d::PSY.ThermalGen) - value_curve = PSY.get_value_curve(cost_function) - cost = PSY.get_function_data(value_curve) - return last(first(PSY.get_points(cost))) +function _onvar_cost(cost_function::PSY.CostCurve{PSY.PiecewisePointCurve}, d::PSY.ThermalGen) + # OnVariableCost is included in the Point itself for PiecewisePointCurve + return 0.0 end -function _no_load_cost(cost_function::Union{PSY.CostCurve{PSY.LinearCurve}, PSY.CostCurve{PSY.QuadraticCurve}}, d::PSY.ThermalGen) +function _onvar_cost(cost_function::Union{PSY.CostCurve{PSY.LinearCurve}, PSY.CostCurve{PSY.QuadraticCurve}}, d::PSY.ThermalGen) value_curve = PSY.get_value_curve(cost_function) cost_component = PSY.get_function_data(value_curve) # Always in \$/h @@ -127,13 +128,22 @@ function _no_load_cost(cost_function::Union{PSY.CostCurve{PSY.LinearCurve}, PSY. return constant_term end -function _no_load_cost(cost_function::PSY.FuelCurve{PSY.PiecewisePointCurve}, d::PSY.ThermalGen) - # value_curve = PSY.get_value_curve(cost_function) - # cost = PSY.get_function_data(value_curve) +function _onvar_cost(cost_function::PSY.CostCurve{PSY.PiecewiseIncrementalCurve}, d::PSY.ThermalGen) + # Input at min is used to transform to InputOutputCurve + return 0.0 +end + +function _onvar_cost(cost_function::PSY.FuelCurve{PSY.PiecewisePointCurve}, d::PSY.ThermalGen) + # OnVariableCost is included in the Point itself for PiecewisePointCurve + return 0.0 +end + +function _onvar_cost(cost_function::PSY.FuelCurve{PSY.PiecewiseIncrementalCurve}, d::PSY.ThermalGen) + # Input at min is used to transform to InputOutputCurve return 0.0 end -function _no_load_cost(cost_function::Union{PSY.FuelCurve{PSY.LinearCurve}, PSY.FuelCurve{PSY.QuadraticCurve}}, d::PSY.ThermalGen) +function _onvar_cost(cost_function::Union{PSY.FuelCurve{PSY.LinearCurve}, PSY.FuelCurve{PSY.QuadraticCurve}}, d::PSY.ThermalGen) value_curve = PSY.get_value_curve(cost_function) cost_component = PSY.get_function_data(value_curve) # In Unit/h (unit typically in ) @@ -307,7 +317,7 @@ function add_variable!( devices::U, formulation::AbstractThermalFormulation, ) where { - T <: OnVariable, + T <: Union{OnVariable, StartVariable, StopVariable}, U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, } where {D <: PSY.ThermalGen} @assert !isempty(devices) @@ -319,23 +329,25 @@ function add_variable!( container, variable_type, D, - [PSY.get_name(d) for d in devices], + [PSY.get_name(d) for d in devices if !PSY.get_must_run(d)], time_steps, ) - for t in time_steps, d in devices - name = PSY.get_name(d) - variable[name, t] = JuMP.@variable( - get_jump_model(container), - base_name = "$(T)_$(D)_{$(name), $(t)}", - binary = binary - ) - if get_warm_start(settings) - init = get_variable_warm_start_value(variable_type, d, formulation) - init !== nothing && JuMP.set_start_value(variable[name, t], init) - end + for d in devices if PSY.get_must_run(d) - JuMP.fix(variable[name, t], 1.0; force = true) + continue + end + name = PSY.get_name(d) + for t in time_steps + variable[name, t] = JuMP.@variable( + get_jump_model(container), + base_name = "$(T)_$(D)_{$(name), $(t)}", + binary = binary + ) + if get_warm_start(settings) + init = get_variable_warm_start_value(variable_type, d, formulation) + init !== nothing && JuMP.set_start_value(variable[name, t], init) + end end end @@ -717,26 +729,35 @@ function add_constraints!( for ic in initial_conditions name = PSY.get_name(get_component(ic)) - constraint[name, 1] = JuMP.@constraint( - get_jump_model(container), - varon[name, 1] == get_value(ic) + varstart[name, 1] - varstop[name, 1] - ) - aux_constraint[name, 1] = JuMP.@constraint( - get_jump_model(container), - varstart[name, 1] + varstop[name, 1] <= 1.0 - ) + if !PSY.get_must_run(get_component(ic)) + constraint[name, 1] = JuMP.@constraint( + get_jump_model(container), + varon[name, 1] == get_value(ic) + varstart[name, 1] - varstop[name, 1] + ) + aux_constraint[name, 1] = JuMP.@constraint( + get_jump_model(container), + varstart[name, 1] + varstop[name, 1] <= 1.0 + ) + end end - for t in time_steps[2:end], ic in initial_conditions - name = get_component_name(ic) - constraint[name, t] = JuMP.@constraint( - get_jump_model(container), - varon[name, t] == varon[name, t - 1] + varstart[name, t] - varstop[name, t] - ) - aux_constraint[name, t] = JuMP.@constraint( - get_jump_model(container), - varstart[name, t] + varstop[name, t] <= 1.0 - ) + for ic in initial_conditions + if PSY.get_must_run(get_component(ic)) + continue + else + name = get_component_name(ic) + for t in time_steps[2:end] + constraint[name, t] = JuMP.@constraint( + get_jump_model(container), + varon[name, t] == + varon[name, t - 1] + varstart[name, t] - varstop[name, t] + ) + aux_constraint[name, t] = JuMP.@constraint( + get_jump_model(container), + varstart[name, t] + varstop[name, t] <= 1.0 + ) + end + end end return end @@ -807,14 +828,19 @@ function calculate_aux_variable_value!( time_steps = get_time_steps(container) for ix in eachindex(JuMP.axes(aux_variable_container)[1]) - IS.@assert_op JuMP.axes(aux_variable_container)[1][ix] == - JuMP.axes(on_variable_results)[1][ix] - IS.@assert_op JuMP.axes(aux_variable_container)[1][ix] == - get_component_name(ini_cond[ix]) - on_var = jump_value.(on_variable_results.data[ix, :]) - ini_cond_value = get_condition(ini_cond[ix]) - aux_variable_container.data[ix, :] .= ini_cond_value - sum_on_var = sum(on_var) + # if its nothing it means the thermal unit was on must run + # so there is nothing to do but to add the total number of time steps + # to the count + if isnothing(get_value(ini_cond[ix])) + sum_on_var = time_steps[end] + else + on_var_name = get_component_name(ini_cond[ix]) + ini_cond_value = get_condition(ini_cond[ix]) + # On Var doesn't exist for a unit that has must_run = true + on_var = jump_value.(on_variable_results[on_var_name, :]) + aux_variable_container.data[ix, :] .= ini_cond_value + sum_on_var = sum(on_var) + end if sum_on_var == time_steps[end] # Unit was always on aux_variable_container.data[ix, :] += time_steps elseif sum_on_var == 0.0 # Unit was always off @@ -848,14 +874,18 @@ function calculate_aux_variable_value!( time_steps = get_time_steps(container) for ix in eachindex(JuMP.axes(aux_variable_container)[1]) - IS.@assert_op JuMP.axes(aux_variable_container)[1][ix] == - JuMP.axes(on_variable_results)[1][ix] - IS.@assert_op JuMP.axes(aux_variable_container)[1][ix] == - get_component_name(ini_cond[ix]) - on_var = jump_value.(on_variable_results.data[ix, :]) - ini_cond_value = get_condition(ini_cond[ix]) - aux_variable_container.data[ix, :] .= ini_cond_value - sum_on_var = sum(on_var) + # if its nothing it means the thermal unit was on must_run = true + # so there is nothing to do but continue + if isnothing(get_value(ini_cond[ix])) + sum_on_var = 0.0 + else + on_var_name = get_component_name(ini_cond[ix]) + # On Var doesn't exist for a unit that has must run + on_var = jump_value.(on_variable_results[on_var_name, :]) + ini_cond_value = get_condition(ini_cond[ix]) + aux_variable_container.data[ix, :] .= ini_cond_value + sum_on_var = sum(on_var) + end if sum_on_var == time_steps[end] # Unit was always on aux_variable_container.data[ix, :] .= 0.0 elseif sum_on_var == 0.0 # Unit was always off diff --git a/src/feedforward/feedforward_constraints.jl b/src/feedforward/feedforward_constraints.jl index a85333b85f..656c2deb89 100644 --- a/src/feedforward/feedforward_constraints.jl +++ b/src/feedforward/feedforward_constraints.jl @@ -178,12 +178,17 @@ function _lower_bound_range_with_parameter!( devices::IS.FlattenIteratorWrapper{V}, ) where {V <: PSY.Component} time_steps = axes(constraint_container)[2] - for device in devices, t in time_steps + for device in devices + if hasmethod(PSY.get_must_run, Tuple{V}) + PSY.get_must_run(device) && continue + end name = PSY.get_name(device) - constraint_container[name, t] = JuMP.@constraint( - jump_model, - lhs_array[name, t] >= param_multiplier[name, t] * param_array[name, t] - ) + for t in time_steps + constraint_container[name, t] = JuMP.@constraint( + jump_model, + lhs_array[name, t] >= param_multiplier[name, t] * param_array[name, t] + ) + end end return end @@ -197,12 +202,17 @@ function _upper_bound_range_with_parameter!( devices::IS.FlattenIteratorWrapper{V}, ) where {V <: PSY.Component} time_steps = axes(constraint_container)[2] - for device in devices, t in time_steps + for device in devices + if hasmethod(PSY.get_must_run, Tuple{V}) + PSY.get_must_run(device) && continue + end name = PSY.get_name(device) - constraint_container[name, t] = JuMP.@constraint( - jump_model, - lhs_array[name, t] <= param_multiplier[name, t] * param_array[name, t] - ) + for t in time_steps + constraint_container[name, t] = JuMP.@constraint( + jump_model, + lhs_array[name, t] <= param_multiplier[name, t] * param_array[name, t] + ) + end end return end diff --git a/src/initial_conditions/add_initial_condition.jl b/src/initial_conditions/add_initial_condition.jl index 5dd0dba455..a5d69f0483 100644 --- a/src/initial_conditions/add_initial_condition.jl +++ b/src/initial_conditions/add_initial_condition.jl @@ -5,7 +5,21 @@ function _get_initial_conditions_value( ::V, container::OptimizationContainer, ) where { - T <: InitialCondition{U, Float64}, + T <: InitialCondition{U, Nothing}, + V <: Union{AbstractDeviceFormulation, AbstractServiceFormulation}, + W <: PSY.Component, +} where {U <: InitialConditionType} + return InitialCondition{U, Nothing}(component, nothing) +end + +function _get_initial_conditions_value( + ::Vector{T}, + component::W, + ::U, + ::V, + container::OptimizationContainer, +) where { + T <: Union{InitialCondition{U, Float64}, InitialCondition{U, Nothing}}, V <: Union{AbstractDeviceFormulation, AbstractServiceFormulation}, W <: PSY.Component, } where {U <: InitialConditionType} @@ -16,9 +30,9 @@ function _get_initial_conditions_value( else val = get_initial_condition_value(ic_data, var_type, W)[1, PSY.get_name(component)] end - @debug "Device $(PSY.get_name(component)) initialized DeviceStatus as $var_type" _group = + @debug "Device $(PSY.get_name(component)) initialized $U as $var_type" _group = LOG_GROUP_BUILD_INITIAL_CONDITIONS - return T(component, val) + return InitialCondition{U, Float64}(component, val) end function _get_initial_conditions_value( @@ -28,8 +42,8 @@ function _get_initial_conditions_value( ::V, container::OptimizationContainer, ) where { - T <: InitialCondition{U, JuMP.VariableRef}, - V <: Union{AbstractDeviceFormulation, AbstractServiceFormulation}, + T <: Union{InitialCondition{U, JuMP.VariableRef}, InitialCondition{U, Nothing}}, + V <: AbstractDeviceFormulation, W <: PSY.Component, } where {U <: InitialConditionType} ic_data = get_initial_conditions_data(container) @@ -39,20 +53,22 @@ function _get_initial_conditions_value( else val = get_initial_condition_value(ic_data, var_type, W)[1, PSY.get_name(component)] end - @debug "Device $(PSY.get_name(component)) initialized DeviceStatus as $var_type" _group = + @debug "Device $(PSY.get_name(component)) initialized $U as $var_type" _group = LOG_GROUP_BUILD_INITIAL_CONDITIONS - return T(component, add_jump_parameter(get_jump_model(container), val)) + return InitialCondition{U, JuMP.VariableRef}( + component, + add_jump_parameter(get_jump_model(container), val), + ) end function _get_initial_conditions_value( - ::Vector{T}, + ::Vector{Union{InitialCondition{U, Float64}, InitialCondition{U, Nothing}}}, component::W, ::U, ::V, container::OptimizationContainer, ) where { - T <: InitialCondition{U, Float64}, - V <: Union{AbstractDeviceFormulation, AbstractServiceFormulation}, + V <: AbstractThermalFormulation, W <: PSY.Component, } where {U <: InitialTimeDurationOff} ic_data = get_initial_conditions_data(container) @@ -66,21 +82,20 @@ function _get_initial_conditions_value( val = PSY.get_time_at_status(component) end end - @debug "Device $(PSY.get_name(component)) initialized DeviceStatus as $var_type" _group = + @debug "Device $(PSY.get_name(component)) initialized $U as $var_type" _group = LOG_GROUP_BUILD_INITIAL_CONDITIONS - return T(component, val) + return InitialCondition{U, Float64}(component, val) end function _get_initial_conditions_value( - ::Vector{T}, + ::Vector{Union{InitialCondition{U, JuMP.VariableRef}, InitialCondition{U, Nothing}}}, component::W, ::U, ::V, container::OptimizationContainer, ) where { - T <: InitialCondition{U, JuMP.VariableRef}, - V <: Union{AbstractDeviceFormulation, AbstractServiceFormulation}, - W <: PSY.Component, + V <: AbstractThermalFormulation, + W <: PSY.ThermalGen, } where {U <: InitialTimeDurationOff} ic_data = get_initial_conditions_data(container) var_type = initial_condition_variable(U(), component, V()) @@ -93,21 +108,23 @@ function _get_initial_conditions_value( val = PSY.get_time_at_status(component) end end - @debug "Device $(PSY.get_name(component)) initialized DeviceStatus as $var_type" _group = + @debug "Device $(PSY.get_name(component)) initialized $U as $var_type" _group = LOG_GROUP_BUILD_INITIAL_CONDITIONS - return T(component, add_jump_parameter(get_jump_model(container), val)) + return InitialCondition{U, JuMP.VariableRef}( + component, + add_jump_parameter(get_jump_model(container), val), + ) end function _get_initial_conditions_value( - ::Vector{T}, + ::Vector{Union{InitialCondition{U, Float64}, InitialCondition{U, Nothing}}}, component::W, ::U, ::V, container::OptimizationContainer, ) where { - T <: InitialCondition{U, Float64}, - V <: Union{AbstractDeviceFormulation, AbstractServiceFormulation}, - W <: PSY.Component, + V <: AbstractThermalFormulation, + W <: PSY.ThermalGen, } where {U <: InitialTimeDurationOn} ic_data = get_initial_conditions_data(container) var_type = initial_condition_variable(U(), component, V()) @@ -120,21 +137,20 @@ function _get_initial_conditions_value( val = PSY.get_time_at_status(component) end end - @debug "Device $(PSY.get_name(component)) initialized DeviceStatus as $var_type" _group = + @debug "Device $(PSY.get_name(component)) initialized $U as $var_type" _group = LOG_GROUP_BUILD_INITIAL_CONDITIONS - return T(component, val) + return InitialCondition{U, Float64}(component, val) end function _get_initial_conditions_value( - ::Vector{T}, + ::Vector{Union{InitialCondition{U, JuMP.VariableRef}, InitialCondition{U, Nothing}}}, component::W, ::U, ::V, container::OptimizationContainer, ) where { - T <: InitialCondition{U, JuMP.VariableRef}, - V <: Union{AbstractDeviceFormulation, AbstractServiceFormulation}, - W <: PSY.Component, + V <: AbstractThermalFormulation, + W <: PSY.ThermalGen, } where {U <: InitialTimeDurationOn} ic_data = get_initial_conditions_data(container) var_type = initial_condition_variable(U(), component, V()) @@ -147,9 +163,12 @@ function _get_initial_conditions_value( val = PSY.get_time_at_status(component) end end - @debug "Device $(PSY.get_name(component)) initialized DeviceStatus as $var_type" _group = + @debug "Device $(PSY.get_name(component)) initialized $U as $var_type" _group = LOG_GROUP_BUILD_INITIAL_CONDITIONS - return T(component, add_jump_parameter(get_jump_model(container), val)) + return InitialCondition{U, JuMP.VariableRef}( + component, + add_jump_parameter(get_jump_model(container), val), + ) end function _get_initial_conditions_value( @@ -160,12 +179,12 @@ function _get_initial_conditions_value( container::OptimizationContainer, ) where { T <: InitialCondition{U, JuMP.VariableRef}, - V <: Union{AbstractDeviceFormulation, AbstractServiceFormulation}, + V <: AbstractDeviceFormulation, W <: PSY.Component, } where {U <: InitialEnergyLevel} var_type = initial_condition_variable(U(), component, V()) val = initial_condition_default(U(), component, V()) - @debug "Device $(PSY.get_name(component)) initialized DeviceStatus as $var_type" _group = + @debug "Device $(PSY.get_name(component)) initialized $U as $var_type" _group = LOG_GROUP_BUILD_INITIAL_CONDITIONS return T(component, add_jump_parameter(get_jump_model(container), val)) end @@ -178,12 +197,12 @@ function _get_initial_conditions_value( container::OptimizationContainer, ) where { T <: InitialCondition{U, Float64}, - V <: Union{AbstractDeviceFormulation, AbstractServiceFormulation}, + V <: AbstractDeviceFormulation, W <: PSY.Component, } where {U <: InitialEnergyLevel} var_type = initial_condition_variable(U(), component, V()) val = initial_condition_default(U(), component, V()) - @debug "Device $(PSY.get_name(component)) initialized DeviceStatus as $var_type" _group = + @debug "Device $(PSY.get_name(component)) initialized $U as $var_type" _group = LOG_GROUP_BUILD_INITIAL_CONDITIONS return T(component, val) end @@ -209,3 +228,35 @@ function add_initial_condition!( end return end + +function add_initial_condition!( + container::OptimizationContainer, + components::Union{Vector{T}, IS.FlattenIteratorWrapper{T}}, + ::U, + ::D, +) where { + T <: PSY.ThermalGen, + U <: AbstractThermalFormulation, + D <: Union{InitialTimeDurationOff, InitialTimeDurationOn, DeviceStatus}, +} + if get_rebuild_model(get_settings(container)) && has_container_key(container, D, T) + return + end + + ini_cond_vector = add_initial_condition_container!(container, D(), T, components) + for (ix, component) in enumerate(components) + if PSY.get_must_run(component) + ini_cond_vector[ix] = InitialCondition{D, Nothing}(component, nothing) + else + ini_cond_vector[ix] = + _get_initial_conditions_value( + ini_cond_vector, + component, + D(), + U(), + container, + ) + end + end + return +end diff --git a/src/initial_conditions/calculate_initial_condition.jl b/src/initial_conditions/calculate_initial_condition.jl index 62af678391..be5add50ff 100644 --- a/src/initial_conditions/calculate_initial_condition.jl +++ b/src/initial_conditions/calculate_initial_condition.jl @@ -23,3 +23,10 @@ function set_ic_quantity!( ic.value = var_value return end + +function set_ic_quantity!( + ::InitialCondition{T, Nothing}, + ::Float64, +) where {T <: InitialConditionType} + return +end diff --git a/src/network_models/pm_translator.jl b/src/network_models/pm_translator.jl index 58381f732e..171297715e 100644 --- a/src/network_models/pm_translator.jl +++ b/src/network_models/pm_translator.jl @@ -271,14 +271,14 @@ function get_branch_to_pm( ::Type{<:PM.AbstractDCPModel}, ) PM_branch = Dict{String, Any}( - "loss1" => PSY.get_loss(branch).l1, + "loss1" => PSY.get_proportional_term(PSY.get_loss(branch)), "mp_pmax" => PSY.get_reactive_power_limits_from(branch).max, "model" => 2, "shutdown" => 0.0, "pmaxt" => PSY.get_active_power_limits_to(branch).max, "pmaxf" => PSY.get_active_power_limits_from(branch).max, "startup" => 0.0, - "loss0" => PSY.get_loss(branch).l0, + "loss0" => PSY.get_constant_term(PSY.get_loss(branch)), "pt" => 0.0, "vt" => PSY.get_magnitude(PSY.get_arc(branch).to), "qmaxf" => PSY.get_reactive_power_limits_from(branch).max, @@ -310,14 +310,14 @@ function get_branch_to_pm( ) check_hvdc_line_limits_unidirectional(branch) PM_branch = Dict{String, Any}( - "loss1" => PSY.get_loss(branch).l1, + "loss1" => PSY.get_proportional_term(PSY.get_loss(branch)), "mp_pmax" => PSY.get_reactive_power_limits_from(branch).max, "model" => 2, "shutdown" => 0.0, "pmaxt" => PSY.get_active_power_limits_to(branch).max, "pmaxf" => PSY.get_active_power_limits_from(branch).max, "startup" => 0.0, - "loss0" => PSY.get_loss(branch).l0, + "loss0" => PSY.get_constant_term(PSY.get_loss(branch)), "pt" => 0.0, "vt" => PSY.get_magnitude(PSY.get_arc(branch).to), "qmaxf" => PSY.get_reactive_power_limits_from(branch).max, @@ -348,14 +348,14 @@ function get_branch_to_pm( ::Type{<:PM.AbstractPowerModel}, ) PM_branch = Dict{String, Any}( - "loss1" => PSY.get_loss(branch).l1, + "loss1" => PSY.get_proportional_term(PSY.get_loss(branch)), "mp_pmax" => PSY.get_reactive_power_limits_from(branch).max, "model" => 2, "shutdown" => 0.0, "pmaxt" => PSY.get_active_power_limits_to(branch).max, "pmaxf" => PSY.get_active_power_limits_from(branch).max, "startup" => 0.0, - "loss0" => PSY.get_loss(branch).l0, + "loss0" => PSY.get_constant_term(PSY.get_loss(branch)), "pt" => 0.0, "vt" => PSY.get_magnitude(PSY.get_arc(branch).to), "qmaxf" => PSY.get_reactive_power_limits_from(branch).max, diff --git a/src/operation/initial_conditions_update_in_memory_store.jl b/src/operation/initial_conditions_update_in_memory_store.jl index 7c7793050b..36707f3b41 100644 --- a/src/operation/initial_conditions_update_in_memory_store.jl +++ b/src/operation/initial_conditions_update_in_memory_store.jl @@ -2,12 +2,25 @@ ################## ic updates from store for emulation problems simulation ################# function update_initial_conditions!( - ics::Vector{T}, + ics::T, store::EmulationModelStore, ::Dates.Millisecond, ) where { - T <: InitialCondition{InitialTimeDurationOn, S}, -} where {S <: Union{Float64, JuMP.VariableRef}} + T <: Union{ + Vector{ + Union{ + InitialCondition{InitialTimeDurationOn, Nothing}, + InitialCondition{InitialTimeDurationOn, Float64}, + }, + }, + Vector{ + Union{ + InitialCondition{InitialTimeDurationOn, Nothing}, + InitialCondition{InitialTimeDurationOn, JuMP.VariableRef}, + }, + }, + }, +} for ic in ics var_val = get_value(store, TimeDurationOn(), get_component_type(ic)) set_ic_quantity!(ic, get_last_recorded_value(var_val)[get_component_name(ic)]) @@ -16,12 +29,25 @@ function update_initial_conditions!( end function update_initial_conditions!( - ics::Vector{T}, + ics::T, store::EmulationModelStore, ::Dates.Millisecond, ) where { - T <: InitialCondition{InitialTimeDurationOff, S}, -} where {S <: Union{Float64, JuMP.VariableRef}} + T <: Union{ + Vector{ + Union{ + InitialCondition{InitialTimeDurationOff, Nothing}, + InitialCondition{InitialTimeDurationOff, Float64}, + }, + }, + Vector{ + Union{ + InitialCondition{InitialTimeDurationOff, Nothing}, + InitialCondition{InitialTimeDurationOff, JuMP.VariableRef}, + }, + }, + }, +} for ic in ics var_val = get_value(store, TimeDurationOff(), get_component_type(ic)) set_ic_quantity!(ic, get_last_recorded_value(var_val)[get_component_name(ic)]) @@ -30,12 +56,25 @@ function update_initial_conditions!( end function update_initial_conditions!( - ics::Vector{T}, + ics::T, store::EmulationModelStore, ::Dates.Millisecond, ) where { - T <: InitialCondition{DevicePower, S}, -} where {S <: Union{Float64, JuMP.VariableRef}} + T <: Union{ + Vector{ + Union{ + InitialCondition{DevicePower, Nothing}, + InitialCondition{DevicePower, Float64}, + }, + }, + Vector{ + Union{ + InitialCondition{DevicePower, Nothing}, + InitialCondition{DevicePower, JuMP.VariableRef}, + }, + }, + }, +} for ic in ics var_val = get_value(store, ActivePowerVariable(), get_component_type(ic)) set_ic_quantity!(ic, get_last_recorded_value(var_val)[get_component_name(ic)]) @@ -44,12 +83,25 @@ function update_initial_conditions!( end function update_initial_conditions!( - ics::Vector{T}, + ics::T, store::EmulationModelStore, ::Dates.Millisecond, ) where { - T <: InitialCondition{DeviceStatus, S}, -} where {S <: Union{Float64, JuMP.VariableRef}} + T <: Union{ + Vector{ + Union{ + InitialCondition{DeviceStatus, Nothing}, + InitialCondition{DeviceStatus, Float64}, + }, + }, + Vector{ + Union{ + InitialCondition{DeviceStatus, Nothing}, + InitialCondition{DeviceStatus, JuMP.VariableRef}, + }, + }, + }, +} for ic in ics var_val = get_value(store, OnVariable(), get_component_type(ic)) set_ic_quantity!(ic, get_last_recorded_value(var_val)[get_component_name(ic)]) @@ -58,12 +110,25 @@ function update_initial_conditions!( end function update_initial_conditions!( - ics::Vector{T}, + ics::T, store::EmulationModelStore, ::Dates.Millisecond, ) where { - T <: InitialCondition{DeviceAboveMinPower, S}, -} where {S <: Union{Float64, JuMP.VariableRef}} + T <: Union{ + Vector{ + Union{ + InitialCondition{DeviceAboveMinPower, Nothing}, + InitialCondition{DeviceAboveMinPower, Float64}, + }, + }, + Vector{ + Union{ + InitialCondition{DeviceAboveMinPower, Nothing}, + InitialCondition{DeviceAboveMinPower, JuMP.VariableRef}, + }, + }, + }, +} for ic in ics var_val = get_value(store, PowerAboveMinimumVariable(), get_component_type(ic)) @@ -72,6 +137,7 @@ function update_initial_conditions!( return end +#= Unused without the AGC model enabled function update_initial_conditions!( ics::Vector{T}, store::EmulationModelStore, @@ -85,14 +151,28 @@ function update_initial_conditions!( end return end +=# function update_initial_conditions!( - ics::Vector{T}, + ics::T, store::EmulationModelStore, ::Dates.Millisecond, ) where { - T <: InitialCondition{InitialEnergyLevel, S}, -} where {S <: Union{Float64, JuMP.VariableRef}} + T <: Union{ + Vector{ + Union{ + InitialCondition{InitialEnergyLevel, Nothing}, + InitialCondition{InitialEnergyLevel, Float64}, + }, + }, + Vector{ + Union{ + InitialCondition{InitialEnergyLevel, Nothing}, + InitialCondition{InitialEnergyLevel, JuMP.VariableRef}, + }, + }, + }, +} for ic in ics var_val = get_value(store, EnergyVariable(), get_component_type(ic)) set_ic_quantity!(ic, get_last_recorded_value(var_val)[get_component_name(ic)]) diff --git a/src/operation/operation_model_interface.jl b/src/operation/operation_model_interface.jl index d2157ee967..bc89287be1 100644 --- a/src/operation/operation_model_interface.jl +++ b/src/operation/operation_model_interface.jl @@ -114,6 +114,10 @@ function solve_impl!(model::OperationModel) tss = replace("$(ts)", ":" => "_") model_export_path = joinpath(model_output_dir, "exported_$(model_name)_$(tss).json") serialize_optimization_model(container, model_export_path) + write_lp_file( + get_jump_model(container), + replace(model_export_path, ".json" => ".lp"), + ) end status = solve_impl!(container, get_system(model)) diff --git a/src/parameters/add_parameters.jl b/src/parameters/add_parameters.jl index 6261ded49c..76f1ef2095 100644 --- a/src/parameters/add_parameters.jl +++ b/src/parameters/add_parameters.jl @@ -338,6 +338,52 @@ function _add_parameters!( return end +function _add_parameters!( + container::OptimizationContainer, + ::T, + key::VariableKey{U, D}, + model::DeviceModel{D, W}, + devices::V, +) where { + T <: OnStatusParameter, + U <: OnVariable, + V <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, + W <: AbstractThermalFormulation, +} where {D <: PSY.ThermalGen} + @debug "adding" T D U _group = LOG_GROUP_OPTIMIZATION_CONTAINER + names = [PSY.get_name(device) for device in devices if !PSY.get_must_run(device)] + time_steps = get_time_steps(container) + parameter_container = add_param_container!(container, T(), D, key, names, time_steps) + jump_model = get_jump_model(container) + for d in devices + if PSY.get_must_run(d) + continue + end + name = PSY.get_name(d) + if get_variable_warm_start_value(U(), d, W()) === nothing + inital_parameter_value = 0.0 + else + inital_parameter_value = get_variable_warm_start_value(U(), d, W()) + end + for t in time_steps + set_multiplier!( + parameter_container, + get_parameter_multiplier(T(), d, W()), + name, + t, + ) + set_parameter!( + parameter_container, + jump_model, + inital_parameter_value, + name, + t, + ) + end + end + return +end + function _add_parameters!( container::OptimizationContainer, ::T, diff --git a/src/services_models/service_slacks.jl b/src/services_models/service_slacks.jl index 0ec2464fda..9d530e3e40 100644 --- a/src/services_models/service_slacks.jl +++ b/src/services_models/service_slacks.jl @@ -1,7 +1,7 @@ function reserve_slacks!( container::OptimizationContainer, service::T, -) where {T <: PSY.Reserve} +) where {T <: Union{PSY.Reserve, PSY.ReserveNonSpinning}} time_steps = get_time_steps(container) variable = add_variable_container!( container, diff --git a/src/simulation/initial_condition_update_simulation.jl b/src/simulation/initial_condition_update_simulation.jl index 4f4a46559e..d360753019 100644 --- a/src/simulation/initial_condition_update_simulation.jl +++ b/src/simulation/initial_condition_update_simulation.jl @@ -22,12 +22,25 @@ function update_initial_conditions!( end function update_initial_conditions!( - ics::Vector{T}, + ics::T, state::SimulationState, model_resolution::Dates.Millisecond, ) where { - T <: InitialCondition{InitialTimeDurationOn, S}, -} where {S <: Union{Float64, JuMP.VariableRef}} + T <: Union{ + Vector{ + Union{ + InitialCondition{InitialTimeDurationOn, Nothing}, + InitialCondition{InitialTimeDurationOn, Float64}, + }, + }, + Vector{ + Union{ + InitialCondition{InitialTimeDurationOn, Nothing}, + InitialCondition{InitialTimeDurationOn, JuMP.VariableRef}, + }, + }, + }, +} for ic in ics var_val = get_system_state_value(state, TimeDurationOn(), get_component_type(ic)) state_resolution = get_data_resolution( @@ -42,13 +55,27 @@ function update_initial_conditions!( end function update_initial_conditions!( - ics::Vector{T}, + ics::T, state::SimulationState, model_resolution::Dates.Millisecond, ) where { - T <: InitialCondition{InitialTimeDurationOff, S}, -} where {S <: Union{Float64, JuMP.VariableRef}} + T <: Union{ + Vector{ + Union{ + InitialCondition{InitialTimeDurationOff, Nothing}, + InitialCondition{InitialTimeDurationOff, Float64}, + }, + }, + Vector{ + Union{ + InitialCondition{InitialTimeDurationOff, Nothing}, + InitialCondition{InitialTimeDurationOff, JuMP.VariableRef}, + }, + }, + }, +} for ic in ics + isnothing(get_value(ic)) && continue var_val = get_system_state_value(state, TimeDurationOff(), get_component_type(ic)) state_resolution = get_data_resolution( get_system_state_data(state, TimeDurationOff(), get_component_type(ic)), @@ -62,19 +89,36 @@ function update_initial_conditions!( end function update_initial_conditions!( - ics::Vector{T}, + ics::T, state::SimulationState, ::Dates.Millisecond, ) where { - T <: InitialCondition{DevicePower, S}, -} where {S <: Union{Float64, JuMP.VariableRef}} + T <: Union{ + Vector{ + Union{ + InitialCondition{DevicePower, Nothing}, + InitialCondition{DevicePower, Float64}, + }, + }, + Vector{ + Union{ + InitialCondition{DevicePower, Nothing}, + InitialCondition{DevicePower, JuMP.VariableRef}, + }, + }, + }, +} for ic in ics comp_name = get_component_name(ic) comp_type = get_component_type(ic) - status_val = get_system_state_value(state, OnVariable(), comp_type)[comp_name] + comp = get_component(ic) + if hasmethod(PSY.get_must_run, Tuple{comp_type}) && PSY.get_must_run(comp) + status_val = 1.0 + else + status_val = get_system_state_value(state, OnVariable(), comp_type)[comp_name] + end var_val = get_system_state_value(state, ActivePowerVariable(), comp_type)[comp_name] if !isapprox(status_val, 0.0; atol = ABSOLUTE_TOLERANCE) - comp = get_component(ic) min = PSY.get_active_power_limits(comp).min max = PSY.get_active_power_limits(comp).max if var_val <= max && var_val >= min @@ -101,13 +145,27 @@ function update_initial_conditions!( end function update_initial_conditions!( - ics::Vector{T}, + ics::T, state::SimulationState, ::Dates.Millisecond, ) where { - T <: InitialCondition{DeviceStatus, S}, -} where {S <: Union{Float64, JuMP.VariableRef}} + T <: Union{ + Vector{ + Union{ + InitialCondition{DeviceStatus, Nothing}, + InitialCondition{DeviceStatus, Float64}, + }, + }, + Vector{ + Union{ + InitialCondition{DeviceStatus, Nothing}, + InitialCondition{DeviceStatus, JuMP.VariableRef}, + }, + }, + }, +} for ic in ics + isnothing(get_value(ic)) && continue var_val = get_system_state_value(state, OnVariable(), get_component_type(ic)) set_ic_quantity!(ic, var_val[get_component_name(ic)]) end @@ -115,12 +173,25 @@ function update_initial_conditions!( end function update_initial_conditions!( - ics::Vector{T}, + ics::T, state::SimulationState, ::Dates.Millisecond, ) where { - T <: InitialCondition{DeviceAboveMinPower, S}, -} where {S <: Union{Float64, JuMP.VariableRef}} + T <: Union{ + Vector{ + Union{ + InitialCondition{DeviceAboveMinPower, Nothing}, + InitialCondition{DeviceAboveMinPower, Float64}, + }, + }, + Vector{ + Union{ + InitialCondition{DeviceAboveMinPower, Nothing}, + InitialCondition{DeviceAboveMinPower, JuMP.VariableRef}, + }, + }, + }, +} for ic in ics var_val = get_system_state_value( state, @@ -133,12 +204,25 @@ function update_initial_conditions!( end function update_initial_conditions!( - ics::Vector{T}, + ics::T, state::SimulationState, ::Dates.Millisecond, ) where { - T <: InitialCondition{InitialEnergyLevel, S}, -} where {S <: Union{Float64, JuMP.VariableRef}} + T <: Union{ + Vector{ + Union{ + InitialCondition{InitialEnergyLevel, Nothing}, + InitialCondition{InitialEnergyLevel, Float64}, + }, + }, + Vector{ + Union{ + InitialCondition{InitialEnergyLevel, Nothing}, + InitialCondition{InitialEnergyLevel, JuMP.VariableRef}, + }, + }, + }, +} for ic in ics var_val = get_system_state_value(state, EnergyVariable(), get_component_type(ic)) set_ic_quantity!(ic, var_val[get_component_name(ic)]) diff --git a/src/utils/file_utils.jl b/src/utils/file_utils.jl index 9ea6116bba..ce219a4886 100644 --- a/src/utils/file_utils.jl +++ b/src/utils/file_utils.jl @@ -3,7 +3,7 @@ Return a decoded JSON file. """ function read_json(filename::AbstractString) open(filename, "r") do io - JSON.parse(io) + JSON3.read(io) end end @@ -28,7 +28,7 @@ end function read_file_hashes(path) data = open(joinpath(path, IS.HASH_FILENAME), "r") do io - JSON.parse(io) + JSON3.read(io) end return data["files"] diff --git a/src/utils/jump_utils.jl b/src/utils/jump_utils.jl index e610a91ae5..129414c0d9 100644 --- a/src/utils/jump_utils.jl +++ b/src/utils/jump_utils.jl @@ -6,7 +6,7 @@ function add_jump_parameter(jump_model::JuMP.Model, val::Number) end function write_data(base_power::Float64, save_path::String) - JSON.write(joinpath(save_path, "base_power.json"), JSON.json(base_power)) + JSON3.write(joinpath(save_path, "base_power.json"), JSON3.json(base_power)) return end @@ -315,6 +315,13 @@ function serialize_jump_optimization_model(jump_model::JuMP.Model, save_path::St return end +function write_lp_file(jump_model::JuMP.Model, save_path::String) + MOF_model = MOPFM(; format = MOI.FileFormats.FORMAT_LP) + MOI.copy_to(MOF_model, JuMP.backend(jump_model)) + MOI.write_to_file(MOF_model, save_path) + return +end + # check_conflict_status functions can't be tested on CI because free solvers don't support IIS function check_conflict_status( jump_model::JuMP.Model, diff --git a/src/utils/recorder_events.jl b/src/utils/recorder_events.jl index 755818c522..57d3704141 100644 --- a/src/utils/recorder_events.jl +++ b/src/utils/recorder_events.jl @@ -60,8 +60,8 @@ end function InitialConditionUpdateEvent( simulation_time, ic::InitialCondition, - previous_value::Float64, - model_name, + previous_value::Union{Nothing, Float64}, + model_name::Symbol, ) return InitialConditionUpdateEvent( IS.RecorderEventCommon("InitialConditionUpdateEvent"), @@ -69,8 +69,8 @@ function InitialConditionUpdateEvent( string(get_ic_type(ic)), string(get_component_type(ic)), get_component_name(ic), - get_condition(ic), - previous_value, + isnothing(get_condition(ic)) ? 1e8 : get_condition(ic), + isnothing(previous_value) ? 1e8 : previous_value, string(model_name), ) end diff --git a/test/performance/performance_test.jl b/test/performance/performance_test.jl index d387796ccf..9c97cce33d 100644 --- a/test/performance/performance_test.jl +++ b/test/performance/performance_test.jl @@ -22,6 +22,11 @@ try sys_rts_rt = build_system(PSISystems, "modified_RTS_GMLC_RT_sys") sys_rts_realization = build_system(PSISystems, "modified_RTS_GMLC_realization_sys") + for sys in [sys_rts_da, sys_rts_rt, sys_rts_realization] + g = get_component(ThermalStandard, sys, "121_NUCLEAR_1") + set_must_run!(g, true) + end + for i in 1:2 template_uc = ProblemTemplate( NetworkModel( @@ -32,8 +37,8 @@ try ), ) - set_device_model!(template_uc, ThermalMultiStart, ThermalCompactUnitCommitment) - set_device_model!(template_uc, ThermalStandard, ThermalCompactUnitCommitment) + set_device_model!(template_uc, ThermalMultiStart, ThermalStandardUnitCommitment) + set_device_model!(template_uc, ThermalStandard, ThermalStandardUnitCommitment) set_device_model!(template_uc, RenewableDispatch, RenewableFullDispatch) set_device_model!(template_uc, PowerLoad, StaticPowerLoad) set_device_model!(template_uc, DeviceModel(Line, StaticBranch)) diff --git a/test/runtests.jl b/test/runtests.jl index f6eaa26fe9..9c0568af49 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,9 +2,11 @@ include("includes.jl") # Code Quality Tests import Aqua -Aqua.test_unbound_args(PowerSimulations) Aqua.test_undefined_exports(PowerSimulations) Aqua.test_ambiguities(PowerSimulations) +Aqua.test_stale_deps(PowerSimulations) +Aqua.test_persistent_tasks(PowerSimulations) +Aqua.test_unbound_args(PowerSimulations) const LOG_FILE = "power-simulations-test.log" diff --git a/test/test_device_branch_constructors.jl b/test/test_device_branch_constructors.jl index 118fda9c60..289b2915a9 100644 --- a/test/test_device_branch_constructors.jl +++ b/test/test_device_branch_constructors.jl @@ -312,18 +312,18 @@ end available = true, active_power_flow = 0.0, # Force the flow in the opposite direction for testing purposes - active_power_limits_from = (min = -2.0, max = -2.0), - active_power_limits_to = (min = -3.0, max = 2.0), + active_power_limits_from = (min = -2.0, max = 2.0), + active_power_limits_to = (min = -2.0, max = 2.0), reactive_power_limits_from = (min = -1.0, max = 1.0), reactive_power_limits_to = (min = -1.0, max = 1.0), arc = get_arc(line), - loss = (l0 = 0.00, l1 = 0.00), + loss = LinearCurve(0.0), ) add_component!(sys_5, hvdc) for net_model in [DCPPowerModel, PTDFPowerModel] @testset "$net_model" begin - PSY.set_loss!(hvdc, (l0 = 0.0, l1 = 0.0)) + PSY.set_loss!(hvdc, PSY.LinearCurve(0.0)) template_uc = ProblemTemplate( NetworkModel(net_model; use_slacks = true), ) @@ -331,7 +331,7 @@ end set_device_model!(template_uc, ThermalStandard, ThermalBasicUnitCommitment) set_device_model!(template_uc, RenewableDispatch, FixedOutput) set_device_model!(template_uc, PowerLoad, StaticPowerLoad) - set_device_model!(template_uc, DeviceModel(Line, StaticBranchUnbounded)) + set_device_model!(template_uc, DeviceModel(Line, StaticBranchBounds)) set_device_model!( template_uc, DeviceModel(TwoTerminalHVDCLine, HVDCTwoTerminalLossless), @@ -424,14 +424,14 @@ end @test all( isapprox.( hvdc_ft_no_loss_values[!, "1"], - hvdc_tf_no_loss_values[!, "1"]; + -hvdc_tf_no_loss_values[!, "1"]; atol = 1e-3, ), ) @test isapprox(no_loss_total_gen, ref_total_gen; atol = 0.1) - PSY.set_loss!(hvdc, (l0 = 0.1, l1 = 0.005)) + PSY.set_loss!(hvdc, PSY.LinearCurve(0.005, 0.1)) model_wl = DecisionModel( template_uc, diff --git a/test/test_device_thermal_generation_constructors.jl b/test/test_device_thermal_generation_constructors.jl index c7387981c3..cdb3a38813 100644 --- a/test/test_device_thermal_generation_constructors.jl +++ b/test/test_device_thermal_generation_constructors.jl @@ -885,26 +885,39 @@ end @testset "Test Must Run ThermalGen" begin sys_5 = build_system(PSITestSystems, "c_sys5_uc") template_uc = - ProblemTemplate(NetworkModel(PTDFPowerModel; PTDF_matrix = PTDF(sys_5))) - set_device_model!(template_uc, ThermalStandard, ThermalBasicUnitCommitment) - set_device_model!(template_uc, RenewableDispatch, FixedOutput) + ProblemTemplate(NetworkModel(CopperPlatePowerModel)) + set_device_model!(template_uc, ThermalStandard, ThermalStandardUnitCommitment) + #set_device_model!(template_uc, RenewableDispatch, FixedOutput) set_device_model!(template_uc, PowerLoad, StaticPowerLoad) - set_device_model!(template_uc, DeviceModel(Line, StaticBranch)) + set_device_model!(template_uc, DeviceModel(Line, StaticBranchUnbounded)) # Set Must Run the most expensive one: Sundance sundance = get_component(ThermalStandard, sys_5, "Sundance") set_must_run!(sundance, true) - model = DecisionModel( - template_uc, - sys_5; - name = "UC", - optimizer = HiGHS_optimizer, - system_to_file = false, - ) + for rebuild in [true, false] + model = DecisionModel( + template_uc, + sys_5; + name = "UC", + optimizer = HiGHS_optimizer, + system_to_file = false, + store_variable_names = true, + rebuild_model = rebuild, + ) - solve!(model; output_dir = mktempdir()) - ptdf_vars = get_variable_values(OptimizationProblemResults(model)) - on = ptdf_vars[PowerSimulations.VariableKey{OnVariable, ThermalStandard}("")] - on_sundance = on[!, "Sundance"] - @test all(isapprox.(on_sundance, 1.0)) + solve!(model; output_dir = mktempdir()) + ptdf_vars = get_variable_values(OptimizationProblemResults(model)) + power = + ptdf_vars[PowerSimulations.VariableKey{ActivePowerVariable, ThermalStandard}( + "", + )] + on = ptdf_vars[PowerSimulations.VariableKey{OnVariable, ThermalStandard}("")] + start = ptdf_vars[PowerSimulations.VariableKey{StartVariable, ThermalStandard}("")] + stop = ptdf_vars[PowerSimulations.VariableKey{StopVariable, ThermalStandard}("")] + power_sundance = power[!, "Sundance"] + @test all(power_sundance .>= 1.0) + for v in [on, start, stop] + @test "Sundance" ∉ names(v) + end + end end diff --git a/test/test_network_constructors.jl b/test/test_network_constructors.jl index 28e40011c8..2b90bb6480 100644 --- a/test/test_network_constructors.jl +++ b/test/test_network_constructors.jl @@ -204,7 +204,7 @@ end test_obj_values = IdDict{System, Float64}( c_sys5 => 342000.0, c_sys14 => 142000.0, - c_sys14_dc => 143000.0, + c_sys14_dc => 135000.0, ) for (ix, sys) in enumerate(systems) template = get_thermal_dispatch_template_network(DCPPowerModel) diff --git a/test/test_simulation_execute.jl b/test/test_simulation_execute.jl index 22b8cbb251..d22a2efaba 100644 --- a/test/test_simulation_execute.jl +++ b/test/test_simulation_execute.jl @@ -42,7 +42,7 @@ end test_path = joinpath(folder, "consecutive", "problems", "ED", "optimization_model_exports") @test ispath(test_path) - @test length(readdir(test_path)) == 2 + @test length(readdir(test_path)) == 4 end function test_2_stage_decision_models_with_feedforwards(in_memory)