From f838cb634d1cbc1f9e3bdaf3631c7eb68c4514d5 Mon Sep 17 00:00:00 2001 From: jjospina Date: Tue, 3 Sep 2024 10:51:01 -0600 Subject: [PATCH] ADD: linear storage model to be used with all formulations based on a predefined problem specification. --- CHANGELOG.md | 2 + src/PowerModelsITD.jl | 2 + src/core/constraint_storage_linear.jl | 38 + src/core/objective_storage.jl | 4 +- src/prob/opfitd_storage_linear.jl | 1178 +++++++++++++++++++++++++ test/opfitd_storage.jl | 2 +- test/opfitd_storage_linear.jl | 137 +++ test/runtests.jl | 1 + 8 files changed, 1361 insertions(+), 3 deletions(-) create mode 100644 src/core/constraint_storage_linear.jl create mode 100755 src/prob/opfitd_storage_linear.jl create mode 100644 test/opfitd_storage_linear.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index f3b29f3..44404c1 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,8 @@ ## staged - Updated to use the new `NonlinearExpr` syntax introduced in JuMP v1.15. In most cases, there should be no user-visible changes. `Breaking change`: PowerModels v0.21.2+, PowerModelsDistribution v0.16+, and JuMP v1.23.2+ should be used to support new compatibility (enforced in `Project.toml`). This upgrade signifcantly decreases the time to build large-scale models. +- Added a new problem formulation `opfitd_storage_linear.jl` that removes the complementary_nl constraint by making rs=0 and xs=0, allows the use of a linear storage model for T&D co-optimization using nonlinear formulations. +- Fixed issue with functions in `objective_storage.jl` not transforming correctly the cost. ## v0.9.2 diff --git a/src/PowerModelsITD.jl b/src/PowerModelsITD.jl index d0d941f..7f7dfb1 100755 --- a/src/PowerModelsITD.jl +++ b/src/PowerModelsITD.jl @@ -48,6 +48,7 @@ module PowerModelsITD include("core/objective_dmld_simple.jl") include("core/objective_storage.jl") include("core/solution.jl") + include("core/constraint_storage_linear.jl") include("data_model/transformations.jl") @@ -67,6 +68,7 @@ module PowerModelsITD include("prob/opfitd_oltc.jl") include("prob/opfitd_dmld.jl") include("prob/opfitd_storage.jl") + include("prob/opfitd_storage_linear.jl") # This must come last to support automated export. include("core/export.jl") diff --git a/src/core/constraint_storage_linear.jl b/src/core/constraint_storage_linear.jl new file mode 100644 index 0000000..2cde6ed --- /dev/null +++ b/src/core/constraint_storage_linear.jl @@ -0,0 +1,38 @@ +""" + constraint_storage_losses_linear(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default) + +Template function for storage loss constraints for linear model - transmission. +""" +function constraint_storage_losses_linear(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default) + storage = _PM.ref(pm, nw, :storage, i) + + # force storage model to be linear + storage["r"] = 0.0 + storage["x"] = 0.0 + + _PM.constraint_storage_losses(pm, nw, i, storage["storage_bus"], storage["r"], storage["x"], storage["p_loss"], storage["q_loss"]) +end + + +""" + constraint_mc_storage_losses_linear(pm::_PMD.AbstractUnbalancedPowerModel, i::Int; nw::Int=_PMD.nw_id_default)::Nothing + +Template function for storage loss constraints for linear model - distribution. +""" +function constraint_mc_storage_losses_linear(pmd::_PMD.AbstractUnbalancedPowerModel, i::Int; nw::Int=_PMD.nw_id_default)::Nothing + storage = _PMD.ref(pmd, nw, :storage, i) + + # force storage model to be linear in nonlinear formulations + storage["r"] = 0.0 + storage["x"] = 0.0 + + _PMD.constraint_mc_storage_losses(pmd, nw, i, storage["storage_bus"], storage["connections"], storage["r"], storage["x"], storage["p_loss"], storage["q_loss"]) + nothing +end + + +"" +function constraint_mc_storage_losses_linear(pmd::_PMD.AbstractUBFModels, i::Int; nw::Int=_PMD.nw_id_default)::Nothing + _PMD.constraint_mc_storage_losses(pmd, i; nw=nw) + nothing +end diff --git a/src/core/objective_storage.jl b/src/core/objective_storage.jl index 383cf65..c0ac49c 100644 --- a/src/core/objective_storage.jl +++ b/src/core/objective_storage.jl @@ -125,7 +125,7 @@ function _objective_itd_min_fuel_cost_polynomial_linquad_storage(pmitd::Abstract for (n, nw_ref) in _PMD.nws(pmd) for (i, strg) in nw_ref[:storage] dsch = _PMD.var(pmd, n, :sd, i) # get discharge power value - strg_cost_dollar_per_pu = strg["cost"][1]#*nw_ref[:settings]["sbase_default"] # convert from $/kWh -> $/pu + strg_cost_dollar_per_pu = strg["cost"][1]*nw_ref[:settings]["sbase_default"] # convert from $/kWh -> $/pu strg_cost_dollar_per_pu = round(strg_cost_dollar_per_pu, digits=4) pmd_strg_cost[(n,i)] = strg_cost_dollar_per_pu*dsch # compute discharge cost end @@ -260,7 +260,7 @@ function _objective_itd_min_fuel_cost_polynomial_nl_storage(pmitd::AbstractPower for (n, nw_ref) in _PMD.nws(pmd) for (i, strg) in nw_ref[:storage] dsch = _PMD.var(pmd, n, :sd, i) # get discharge power value - strg_cost_dollar_per_pu = strg["cost"][1]#*nw_ref[:settings]["sbase_default"] # convert from $/kWh -> $/pu + strg_cost_dollar_per_pu = strg["cost"][1]*nw_ref[:settings]["sbase_default"] # convert from $/kWh -> $/pu strg_cost_dollar_per_pu = round(strg_cost_dollar_per_pu, digits=4) pmd_strg_cost[(n,i)] = strg_cost_dollar_per_pu*dsch # compute discharge cost end diff --git a/src/prob/opfitd_storage_linear.jl b/src/prob/opfitd_storage_linear.jl new file mode 100755 index 0000000..85df75d --- /dev/null +++ b/src/prob/opfitd_storage_linear.jl @@ -0,0 +1,1178 @@ +# Definitions for solving the integrated T&D opf problem with storage opf dispatch with linear storage model. + +""" + function solve_opfitd_storage_linear( + pm_file, + pmd_file, + pmitd_file, + pmitd_type, + optimizer; + solution_processors::Vector{<:Function}=Function[], + pmitd_ref_extensions::Vector{<:Function}=Vector{Function}([]), + eng2math_passthrough::Dict{String,Vector{String}}=Dict{String,Vector{String}}(), + make_si::Bool=true, + auto_rename::Bool=false, + solution_model::String="eng", + kwargs... + ) + +Solve Integrated T&D Optimal Power Flow with Storage OPF Dispatch with Linear Storage Model. +""" +function solve_opfitd_storage_linear(pm_file, pmd_file, pmitd_file, pmitd_type, optimizer; solution_processors::Vector{<:Function}=Function[], pmitd_ref_extensions::Vector{<:Function}=Vector{Function}([]), eng2math_passthrough::Dict{String,Vector{String}}=Dict{String,Vector{String}}(), make_si::Bool=true, auto_rename::Bool=false, solution_model::String="eng", kwargs...) + + if isempty(eng2math_passthrough) + eng2math_passthrough = Dict("storage"=>["cost"]) # by default, pass the eng2math passthrough + else + eng2math_pass_strg = "storage"=>["cost"] + push!(eng2math_passthrough, eng2math_pass_strg) + end + + return solve_model(pm_file, pmd_file, pmitd_file, pmitd_type, optimizer, build_opfitd_storage_linear; solution_processors=solution_processors, pmitd_ref_extensions=pmitd_ref_extensions, eng2math_passthrough=eng2math_passthrough, make_si=make_si, auto_rename=auto_rename, solution_model=solution_model, kwargs...) +end + + +""" + function solve_opfitd_storage_linear( + pmitd_data::Dict{String,<:Any}, + pmitd_type, + optimizer; + solution_processors::Vector{<:Function}=Function[], + pmitd_ref_extensions::Vector{<:Function}=Vector{Function}([]), + eng2math_passthrough::Dict{String,Vector{String}}=Dict{String,Vector{String}}(), + make_si::Bool=true, + solution_model::String="eng", + kwargs... + ) + +Solve Integrated T&D Optimal Power Flow with Storage OPF Dispatch with Linear Storage Model. +""" +function solve_opfitd_storage_linear(pmitd_data::Dict{String,<:Any}, pmitd_type, optimizer; solution_processors::Vector{<:Function}=Function[], pmitd_ref_extensions::Vector{<:Function}=Vector{Function}([]), eng2math_passthrough::Dict{String,Vector{String}}=Dict{String,Vector{String}}(), make_si::Bool=true, solution_model::String="eng", kwargs...) + + if isempty(eng2math_passthrough) + eng2math_passthrough = Dict("storage"=>["cost"]) # by default, pass the eng2math passthrough + else + eng2math_pass_strg = "storage"=>["cost"] + push!(eng2math_passthrough, eng2math_pass_strg) + end + + return solve_model(pmitd_data, pmitd_type, optimizer, build_opfitd_storage_linear; solution_processors=solution_processors, pmitd_ref_extensions=pmitd_ref_extensions, eng2math_passthrough=eng2math_passthrough, make_si=make_si, solution_model=solution_model, kwargs...) +end + + +""" + function solve_mn_opfitd_storage_linear( + pm_file, + pmd_file, + pmitd_file, + pmitd_type, + optimizer; + solution_processors::Vector{<:Function}=Function[], + pmitd_ref_extensions::Vector{<:Function}=Vector{Function}([]), + eng2math_passthrough::Dict{String,Vector{String}}=Dict{String,Vector{String}}(), + make_si::Bool=true, + auto_rename::Bool=false, + solution_model::String="eng", + kwargs... + ) + +Solve Multinetwork Integrated T&D Optimal Power Flow with Storage OPF Dispatch with Linear Storage Model. +""" +function solve_mn_opfitd_storage_linear(pm_file, pmd_file, pmitd_file, pmitd_type, optimizer; solution_processors::Vector{<:Function}=Function[], pmitd_ref_extensions::Vector{<:Function}=Vector{Function}([]), eng2math_passthrough::Dict{String,Vector{String}}=Dict{String,Vector{String}}(), make_si::Bool=true, auto_rename::Bool=false, solution_model::String="eng", kwargs...) + + if isempty(eng2math_passthrough) + eng2math_passthrough = Dict("storage"=>["cost"]) # by default, pass the eng2math passthrough + else + eng2math_pass_strg = "storage"=>["cost"] + push!(eng2math_passthrough, eng2math_pass_strg) + end + + return solve_model(pm_file, pmd_file, pmitd_file, pmitd_type, optimizer, build_mn_opfitd_storage_linear; multinetwork=true, solution_processors=solution_processors, pmitd_ref_extensions=pmitd_ref_extensions, eng2math_passthrough=eng2math_passthrough, make_si=make_si, auto_rename=auto_rename, solution_model=solution_model, kwargs...) +end + + +""" + function solve_mn_opfitd_storage_linear( + pmitd_data::Dict{String,<:Any} + pmitd_type, + optimizer; + solution_processors::Vector{<:Function}=Function[], + pmitd_ref_extensions::Vector{<:Function}=Vector{Function}([]), + eng2math_passthrough::Dict{String,Vector{String}}=Dict{String,Vector{String}}(), + make_si::Bool=true, + solution_model::String="eng", + kwargs... + ) + +Solve Multinetwork Integrated T&D Optimal Power Flow with Storage OPF Dispatch with Linear Storage Model. +""" +function solve_mn_opfitd_storage_linear(pmitd_data::Dict{String,<:Any}, pmitd_type, optimizer; solution_processors::Vector{<:Function}=Function[], pmitd_ref_extensions::Vector{<:Function}=Vector{Function}([]), eng2math_passthrough::Dict{String,Vector{String}}=Dict{String,Vector{String}}(), make_si::Bool=true, solution_model::String="eng", kwargs...) + + if isempty(eng2math_passthrough) + eng2math_passthrough = Dict("storage"=>["cost"]) # by default, pass the eng2math passthrough + else + eng2math_pass_strg = "storage"=>["cost"] + push!(eng2math_passthrough, eng2math_pass_strg) + end + + return solve_model(pmitd_data, pmitd_type, optimizer, build_mn_opfitd_storage_linear; multinetwork=true, solution_processors=solution_processors, pmitd_ref_extensions=pmitd_ref_extensions, eng2math_passthrough=eng2math_passthrough, make_si=make_si, solution_model=solution_model, kwargs...) +end + + +""" + function build_opfitd_storage_linear( + pmitd::AbstractPowerModelITD + ) +Constructor for Integrated T&D Optimal Power Flow with Storage OPF Dispatch with Linear Storage Model. +""" +function build_opfitd_storage_linear(pmitd::AbstractPowerModelITD) + + # Get Models + pm_model = _get_powermodel_from_powermodelitd(pmitd) + pmd_model = _get_powermodeldistribution_from_powermodelitd(pmitd) + + # PM(Transmission) Variables + _PM.variable_bus_voltage(pm_model) + _PM.variable_gen_power(pm_model) + _PM.variable_branch_power(pm_model) + _PM.variable_dcline_power(pm_model) + _PM.variable_storage_power(pm_model) + + # PMD(Distribution) Variables + _PMD.variable_mc_bus_voltage(pmd_model) + _PMD.variable_mc_branch_power(pmd_model) + _PMD.variable_mc_transformer_power(pmd_model) + _PMD.variable_mc_switch_power(pmd_model) + _PMD.variable_mc_generator_power(pmd_model) + _PMD.variable_mc_load_power(pmd_model) + _PMD.variable_mc_storage_power(pmd_model) + + # PMITD (Boundary) Variables + variable_boundary_power(pmitd) + + # --- PM(Transmission) Constraints --- + _PM.constraint_model_voltage(pm_model) + + # reference buses (this only needs to happen for pm(transmission)) + for i in _PM.ids(pm_model, :ref_buses) + _PM.constraint_theta_ref(pm_model, i) + end + + for i in _PM.ids(pm_model, :storage) + _PM.constraint_storage_state(pm_model, i) + constraint_storage_losses_linear(pm_model, i) + _PM.constraint_storage_thermal_limit(pm_model, i) + end + + # PM branches + for i in _PM.ids(pm_model, :branch) + _PM.constraint_ohms_yt_from(pm_model, i) + _PM.constraint_ohms_yt_to(pm_model, i) + + _PM.constraint_voltage_angle_difference(pm_model, i) + + _PM.constraint_thermal_limit_from(pm_model, i) + _PM.constraint_thermal_limit_to(pm_model, i) + end + + # PM DC lines + for i in _PM.ids(pm_model, :dcline) + _PM.constraint_dcline_power_losses(pm_model, i) + end + + # ------------------------------------------------- + # --- PMD(Distribution) Constraints --- + _PMD.constraint_mc_model_voltage(pmd_model) + + # generators should be constrained before KCL, or Pd/Qd undefined + for id in _PMD.ids(pmd_model, :gen) + _PMD.constraint_mc_generator_power(pmd_model, id) + end + + # loads should be constrained before KCL, or Pd/Qd undefined + for id in _PMD.ids(pmd_model, :load) + _PMD.constraint_mc_load_power(pmd_model, id) + end + + for i in _PMD.ids(pmd_model, :storage) + _PMD.constraint_storage_state(pmd_model, i) + constraint_mc_storage_losses_linear(pmd_model, i) + _PMD.constraint_mc_storage_thermal_limit(pmd_model, i) + end + + for i in _PMD.ids(pmd_model, :branch) + _PMD.constraint_mc_ohms_yt_from(pmd_model, i) + _PMD.constraint_mc_ohms_yt_to(pmd_model, i) + + _PMD.constraint_mc_voltage_angle_difference(pmd_model, i) + + _PMD.constraint_mc_thermal_limit_from(pmd_model, i) + _PMD.constraint_mc_thermal_limit_to(pmd_model, i) + _PMD.constraint_mc_ampacity_from(pmd_model, i) + _PMD.constraint_mc_ampacity_to(pmd_model, i) + end + + for i in _PMD.ids(pmd_model, :switch) + _PMD.constraint_mc_switch_state(pmd_model, i) + _PMD.constraint_mc_switch_thermal_limit(pmd_model, i) + _PMD.constraint_mc_switch_ampacity(pmd_model, i) + end + + for i in _PMD.ids(pmd_model, :transformer) + _PMD.constraint_mc_transformer_power(pmd_model, i) + end + + # ------------------------------------------------- + # --- PMITD(T&D) INDEPENDENT Constraints ---------- + + for i in ids(pmitd, :boundary) + constraint_boundary_power(pmitd, i) + constraint_boundary_voltage_magnitude(pmitd, i) + constraint_boundary_voltage_angle(pmitd, i) + end + + # ------------------------------------------------- + # --- PMITD(T&D) KCL Constraints ---------- + # Note: Both of these need to consider flow on boundaries if bus is connected to boundary + boundary_buses_transmission = Vector{Int}() # vector to store the boundary buses transmission + boundary_buses_distribution = Vector{Int}() # vector to store the boundary buses distribution + for j in ids(pmitd, :boundary) + boundary_pmitd = ref(pmitd, nw_id_default, :boundary, j) + bus_pm = boundary_pmitd["f_bus"] + bus_pmd = boundary_pmitd["t_bus"] + push!(boundary_buses_transmission, bus_pm) + push!(boundary_buses_distribution, bus_pmd) + end + # Convert to Julia Set - Note: membership checks are faster in sets (vs. vectors) in Julia + boundary_buses_transmission_set = Set(boundary_buses_transmission) + boundary_buses_distribution_set = Set(boundary_buses_distribution) + + # # ---- Transmission Power Balance --- + for i in _PM.ids(pm_model, :bus) + if i in boundary_buses_transmission_set + constraint_transmission_power_balance_boundary(pmitd, i) + else + _PM.constraint_power_balance(pm_model, i) + end + end + + # ---- Distribution Power Balance --- + for i in _PMD.ids(pmd_model, :bus) + if i in boundary_buses_distribution_set + constraint_distribution_power_balance_boundary(pmitd, i) + else + _PMD.constraint_mc_power_balance(pmd_model, i) + end + end + + # ------------------------------------------------- + # --- PMITD(T&D) Cost Functions ------------------- + objective_itd_min_fuel_cost_storage(pmitd) + +end + + +""" + function build_opfitd_storage_linear( + pmitd::AbstractIVRPowerModelITD + ) +Constructor for Integrated T&D Optimal Power Flow in current-voltage (IV) variable space with Storage OPF Dispatch with Linear Storage Model. +""" +function build_opfitd_storage_linear(pmitd::AbstractIVRPowerModelITD) + + @error "IVR-IVRU formulation not yet supported for storage problems." + throw(error()) + # TODO + +end + + +""" + function build_opfitd_storage_linear( + pmitd::AbstractBFPowerModelITD + ) +Constructor for Integrated T&D Optimal Power Flow for BF Models with Storage OPF Dispatch with Linear Storage Model. +""" +function build_opfitd_storage_linear(pmitd::AbstractBFPowerModelITD) + + # Get Models + pm_model = _get_powermodel_from_powermodelitd(pmitd) + pmd_model = _get_powermodeldistribution_from_powermodelitd(pmitd) + + # PM(Transmission) Variables + _PM.variable_bus_voltage(pm_model) + _PM.variable_gen_power(pm_model) + _PM.variable_branch_power(pm_model) + _PM.variable_branch_current(pm_model) + _PM.variable_dcline_power(pm_model) + _PM.variable_storage_power(pm_model) + + + # PMD(Distribution) Variables + _PMD.variable_mc_bus_voltage(pmd_model) + _PMD.variable_mc_branch_current(pmd_model) + _PMD.variable_mc_branch_power(pmd_model) + _PMD.variable_mc_transformer_power(pmd_model) + _PMD.variable_mc_switch_power(pmd_model) + _PMD.variable_mc_generator_power(pmd_model) + _PMD.variable_mc_load_power(pmd_model) + _PMD.variable_mc_storage_power(pmd_model) + + # PMITD (Boundary) Variables + variable_boundary_power(pmitd) + + # --- PM(Transmission) Constraints --- + _PM.constraint_model_current(pm_model) + + # reference buses (this only needs to happen for pm(transmission)) + for i in _PM.ids(pm_model, :ref_buses) + _PM.constraint_theta_ref(pm_model, i) + end + + for i in _PM.ids(pm_model, :storage) + _PM.constraint_storage_state(pm_model, i) + constraint_storage_losses_linear(pm_model, i) + _PM.constraint_storage_thermal_limit(pm_model, i) + end + + # PM branches + for i in _PM.ids(pm_model, :branch) + _PM.constraint_power_losses(pm_model, i) + _PM.constraint_voltage_magnitude_difference(pm_model, i) + + _PM.constraint_voltage_angle_difference(pm_model, i) + + _PM.constraint_thermal_limit_from(pm_model, i) + _PM.constraint_thermal_limit_to(pm_model, i) + end + + # PM DC lines + for i in _PM.ids(pm_model, :dcline) + _PM.constraint_dcline_power_losses(pm_model, i) + end + + # ------------------------------------------------- + # --- PMD(Distribution) Constraints --- + _PMD.constraint_mc_model_current(pmd_model) + + # generators should be constrained before KCL, or Pd/Qd undefined + for id in _PMD.ids(pmd_model, :gen) + _PMD.constraint_mc_generator_power(pmd_model, id) + end + + # loads should be constrained before KCL, or Pd/Qd undefined + for id in _PMD.ids(pmd_model, :load) + _PMD.constraint_mc_load_power(pmd_model, id) + end + + for i in _PMD.ids(pmd_model, :storage) + _PMD.constraint_storage_state(pmd_model, i) + constraint_mc_storage_losses_linear(pmd_model, i) + _PMD.constraint_mc_storage_thermal_limit(pmd_model, i) + end + + for i in _PMD.ids(pmd_model, :branch) + _PMD.constraint_mc_power_losses(pmd_model, i) + _PMD.constraint_mc_model_voltage_magnitude_difference(pmd_model, i) + + _PMD.constraint_mc_voltage_angle_difference(pmd_model, i) + + _PMD.constraint_mc_thermal_limit_from(pmd_model, i) + _PMD.constraint_mc_thermal_limit_to(pmd_model, i) + _PMD.constraint_mc_ampacity_from(pmd_model, i) + _PMD.constraint_mc_ampacity_to(pmd_model, i) + end + + for i in _PMD.ids(pmd_model, :switch) + _PMD.constraint_mc_switch_state(pmd_model, i) + _PMD.constraint_mc_switch_thermal_limit(pmd_model, i) + _PMD.constraint_mc_switch_ampacity(pmd_model, i) + end + + + for i in _PMD.ids(pmd_model, :transformer) + _PMD.constraint_mc_transformer_power(pmd_model, i) + end + + + # ------------------------------------------------- + # --- PMITD(T&D) INDEPENDENT Constraints ---------- + + for i in ids(pmitd, :boundary) + constraint_boundary_power(pmitd, i) + constraint_boundary_voltage_magnitude(pmitd, i) + constraint_boundary_voltage_angle(pmitd, i) + end + + # ------------------------------------------------- + # --- PMITD(T&D) KCL Constraints ---------- + # Note: Both of these need to consider flow on boundaries if bus is connected to boundary + boundary_buses_transmission = Vector{Int}() # vector to store the boundary buses transmission + boundary_buses_distribution = Vector{Int}() # vector to store the boundary buses distribution + for j in ids(pmitd, :boundary) + boundary_pmitd = ref(pmitd, nw_id_default, :boundary, j) + bus_pm = boundary_pmitd["f_bus"] + bus_pmd = boundary_pmitd["t_bus"] + push!(boundary_buses_transmission, bus_pm) + push!(boundary_buses_distribution, bus_pmd) + end + # Convert to Julia Set - Note: membership checks are faster in sets (vs. vectors) in Julia + boundary_buses_transmission_set = Set(boundary_buses_transmission) + boundary_buses_distribution_set = Set(boundary_buses_distribution) + + # # ---- Transmission Power Balance --- + for i in _PM.ids(pm_model, :bus) + if i in boundary_buses_transmission_set + constraint_transmission_power_balance_boundary(pmitd, i) + else + _PM.constraint_power_balance(pm_model, i) + end + end + + # ---- Distribution Power Balance --- + for i in _PMD.ids(pmd_model, :bus) + if i in boundary_buses_distribution_set + constraint_distribution_power_balance_boundary(pmitd, i) + else + _PMD.constraint_mc_power_balance(pmd_model, i) + end + end + + # ------------------------------------------------- + # --- PMITD(T&D) Cost Functions ------------------- + objective_itd_min_fuel_cost_storage(pmitd) + +end + + +# -- Combined (Hybrid) Formulations +""" + function build_opfitd_storage_linear( + pmitd::AbstractLNLBFPowerModelITD + ) +Constructor for Integrated T&D Optimal Power Flow for L/NL to BF with Storage OPF Dispatch with Linear Storage Model. +""" +function build_opfitd_storage_linear(pmitd::AbstractLNLBFPowerModelITD) + + # Get Models + pm_model = _get_powermodel_from_powermodelitd(pmitd) + pmd_model = _get_powermodeldistribution_from_powermodelitd(pmitd) + + # PM(Transmission) Variables + _PM.variable_bus_voltage(pm_model) + _PM.variable_gen_power(pm_model) + _PM.variable_branch_power(pm_model) + _PM.variable_dcline_power(pm_model) + _PM.variable_storage_power(pm_model) + + + # PMD(Distribution) Variables + _PMD.variable_mc_bus_voltage(pmd_model) + _PMD.variable_mc_branch_current(pmd_model) + _PMD.variable_mc_branch_power(pmd_model) + _PMD.variable_mc_transformer_power(pmd_model) + _PMD.variable_mc_switch_power(pmd_model) + _PMD.variable_mc_generator_power(pmd_model) + _PMD.variable_mc_load_power(pmd_model) + _PMD.variable_mc_storage_power(pmd_model) + + # PMITD (Boundary) Variables + variable_boundary_power(pmitd) + + # --- PM(Transmission) Constraints --- + _PM.constraint_model_voltage(pm_model) + + # reference buses (this only needs to happen for pm(transmission)) + for i in _PM.ids(pm_model, :ref_buses) + _PM.constraint_theta_ref(pm_model, i) + end + + for i in _PM.ids(pm_model, :storage) + _PM.constraint_storage_state(pm_model, i) + constraint_storage_losses_linear(pm_model, i) + _PM.constraint_storage_thermal_limit(pm_model, i) + end + + # PM branches + for i in _PM.ids(pm_model, :branch) + _PM.constraint_ohms_yt_from(pm_model, i) + _PM.constraint_ohms_yt_to(pm_model, i) + + _PM.constraint_voltage_angle_difference(pm_model, i) + + _PM.constraint_thermal_limit_from(pm_model, i) + _PM.constraint_thermal_limit_to(pm_model, i) + end + + # PM DC lines + for i in _PM.ids(pm_model, :dcline) + _PM.constraint_dcline_power_losses(pm_model, i) + end + + # ------------------------------------------------- + # --- PMD(Distribution) Constraints --- + _PMD.constraint_mc_model_current(pmd_model) + + # generators should be constrained before KCL, or Pd/Qd undefined + for id in _PMD.ids(pmd_model, :gen) + _PMD.constraint_mc_generator_power(pmd_model, id) + end + + # loads should be constrained before KCL, or Pd/Qd undefined + for id in _PMD.ids(pmd_model, :load) + _PMD.constraint_mc_load_power(pmd_model, id) + end + + for i in _PMD.ids(pmd_model, :storage) + _PMD.constraint_storage_state(pmd_model, i) + constraint_mc_storage_losses_linear(pmd_model, i) + _PMD.constraint_mc_storage_thermal_limit(pmd_model, i) + end + + for i in _PMD.ids(pmd_model, :branch) + _PMD.constraint_mc_power_losses(pmd_model, i) + _PMD.constraint_mc_model_voltage_magnitude_difference(pmd_model, i) + + _PMD.constraint_mc_voltage_angle_difference(pmd_model, i) + + _PMD.constraint_mc_thermal_limit_from(pmd_model, i) + _PMD.constraint_mc_thermal_limit_to(pmd_model, i) + _PMD.constraint_mc_ampacity_from(pmd_model, i) + _PMD.constraint_mc_ampacity_to(pmd_model, i) + end + + for i in _PMD.ids(pmd_model, :switch) + _PMD.constraint_mc_switch_state(pmd_model, i) + _PMD.constraint_mc_switch_thermal_limit(pmd_model, i) + _PMD.constraint_mc_switch_ampacity(pmd_model, i) + end + + for i in _PMD.ids(pmd_model, :transformer) + _PMD.constraint_mc_transformer_power(pmd_model, i) + end + + # ------------------------------------------------- + # --- PMITD(T&D) INDEPENDENT Constraints ---------- + + for i in ids(pmitd, :boundary) + constraint_boundary_power(pmitd, i) + constraint_boundary_voltage_magnitude(pmitd, i) + constraint_boundary_voltage_angle(pmitd, i) + end + + # ------------------------------------------------- + # --- PMITD(T&D) KCL Constraints ---------- + # Note: Both of these need to consider flow on boundaries if bus is connected to boundary + boundary_buses_transmission = Vector{Int}() # vector to store the boundary buses transmission + boundary_buses_distribution = Vector{Int}() # vector to store the boundary buses distribution + for j in ids(pmitd, :boundary) + boundary_pmitd = ref(pmitd, nw_id_default, :boundary, j) + bus_pm = boundary_pmitd["f_bus"] + bus_pmd = boundary_pmitd["t_bus"] + push!(boundary_buses_transmission, bus_pm) + push!(boundary_buses_distribution, bus_pmd) + end + # Convert to Julia Set - Note: membership checks are faster in sets (vs. vectors) in Julia + boundary_buses_transmission_set = Set(boundary_buses_transmission) + boundary_buses_distribution_set = Set(boundary_buses_distribution) + + # # ---- Transmission Power Balance --- + for i in _PM.ids(pm_model, :bus) + if i in boundary_buses_transmission_set + constraint_transmission_power_balance_boundary(pmitd, i) + else + _PM.constraint_power_balance(pm_model, i) + end + end + + # ---- Distribution Power Balance --- + for i in _PMD.ids(pmd_model, :bus) + if i in boundary_buses_distribution_set + constraint_distribution_power_balance_boundary(pmitd, i) + else + _PMD.constraint_mc_power_balance(pmd_model, i) + end + end + + # ------------------------------------------------- + # --- PMITD(T&D) Cost Functions ------------------- + objective_itd_min_fuel_cost_storage(pmitd) + +end + + +# ---------------------------------------------------------------------------------------- +# --- Multinetwork OPFITD Problem Specifications +# ---------------------------------------------------------------------------------------- + +""" + function build_mn_opfitd_storage_linear( + pmitd::AbstractPowerModelITD + ) +Constructor for Multinetwork Integrated T&D Optimal Power Flow with Storage OPF Dispatch with Linear Storage Model. +""" +function build_mn_opfitd_storage_linear(pmitd::AbstractPowerModelITD) + + # Get Models + pm_model = _get_powermodel_from_powermodelitd(pmitd) + pmd_model = _get_powermodeldistribution_from_powermodelitd(pmitd) + + for (n, network) in nws(pmitd) + # PM(Transmission) Variables + _PM.variable_bus_voltage(pm_model, nw=n) + _PM.variable_gen_power(pm_model, nw=n) + _PM.variable_branch_power(pm_model, nw=n) + _PM.variable_dcline_power(pm_model, nw=n) + _PM.variable_storage_power(pm_model, nw=n) + + # PMD(Distribution) Variables + _PMD.variable_mc_bus_voltage(pmd_model; nw=n) + _PMD.variable_mc_branch_power(pmd_model; nw=n) + _PMD.variable_mc_transformer_power(pmd_model; nw=n) + _PMD.variable_mc_switch_power(pmd_model; nw=n) + _PMD.variable_mc_generator_power(pmd_model; nw=n) + _PMD.variable_mc_load_power(pmd_model; nw=n) + _PMD.variable_mc_storage_power(pmd_model; nw=n) + + # PMITD (Boundary) Variables + variable_boundary_power(pmitd; nw=n) + + # --- PM(Transmission) Constraints --- + _PM.constraint_model_voltage(pm_model, nw=n) + + # reference buses (this only needs to happen for pm(transmission)) + for i in _PM.ids(pm_model, :ref_buses, nw=n) + _PM.constraint_theta_ref(pm_model, i, nw=n) + end + + for i in _PM.ids(pm_model, :storage, nw=n) + constraint_storage_losses_linear(pm_model, i, nw=n) + _PM.constraint_storage_thermal_limit(pm_model, i, nw=n) + end + + # PM branches + for i in _PM.ids(pm_model, :branch, nw=n) + _PM.constraint_ohms_yt_from(pm_model, i, nw=n) + _PM.constraint_ohms_yt_to(pm_model, i, nw=n) + + _PM.constraint_voltage_angle_difference(pm_model, i, nw=n) + + _PM.constraint_thermal_limit_from(pm_model, i, nw=n) + _PM.constraint_thermal_limit_to(pm_model, i, nw=n) + end + + # PM DC lines + for i in _PM.ids(pm_model, :dcline, nw=n) + _PM.constraint_dcline_power_losses(pm_model, i, nw=n) + end + + # ------------------------------------------------- + # --- PMD(Distribution) Constraints --- + _PMD.constraint_mc_model_voltage(pmd_model; nw=n) + + # generators should be constrained before KCL, or Pd/Qd undefined + for i in _PMD.ids(pmd_model, n, :gen) + _PMD.constraint_mc_generator_power(pmd_model, i; nw=n) + end + + # loads should be constrained before KCL, or Pd/Qd undefined + for i in _PMD.ids(pmd_model, n, :load) + _PMD.constraint_mc_load_power(pmd_model, i; nw=n) + end + + for i in _PMD.ids(pmd_model, n, :storage) + constraint_mc_storage_losses_linear(pmd_model, i; nw=n) + _PMD.constraint_mc_storage_thermal_limit(pmd_model, i; nw=n) + end + + for i in _PMD.ids(pmd_model, n, :branch) + _PMD.constraint_mc_ohms_yt_from(pmd_model, i; nw=n) + _PMD.constraint_mc_ohms_yt_to(pmd_model, i; nw=n) + + _PMD.constraint_mc_voltage_angle_difference(pmd_model, i; nw=n) + + _PMD.constraint_mc_thermal_limit_from(pmd_model, i; nw=n) + _PMD.constraint_mc_thermal_limit_to(pmd_model, i; nw=n) + _PMD.constraint_mc_ampacity_from(pmd_model, i; nw=n) + _PMD.constraint_mc_ampacity_to(pmd_model, i; nw=n) + end + + for i in _PMD.ids(pmd_model, n, :switch) + _PMD.constraint_mc_switch_state(pmd_model, i; nw=n) + _PMD.constraint_mc_switch_thermal_limit(pmd_model, i; nw=n) + _PMD.constraint_mc_switch_ampacity(pmd_model, i; nw=n) + end + + for i in _PMD.ids(pmd_model, n, :transformer) + _PMD.constraint_mc_transformer_power(pmd_model, i; nw=n) + end + + # ------------------------------------------------- + # --- PMITD(T&D) INDEPENDENT Constraints ---------- + + for i in ids(pmitd, :boundary; nw=n) + constraint_boundary_power(pmitd, i; nw=n) + constraint_boundary_voltage_magnitude(pmitd, i; nw=n) + constraint_boundary_voltage_angle(pmitd, i; nw=n) + end + + # ------------------------------------------------- + # --- PMITD(T&D) KCL Constraints ---------- + # Note: Both of these need to consider flow on boundaries if bus is connected to boundary + boundary_buses_transmission = Vector{Int}() # vector to store the boundary buses transmission + boundary_buses_distribution = Vector{Int}() # vector to store the boundary buses distribution + for j in ids(pmitd, :boundary; nw=n) + boundary_pmitd = ref(pmitd, n, :boundary, j) + bus_pm = boundary_pmitd["f_bus"] + bus_pmd = boundary_pmitd["t_bus"] + push!(boundary_buses_transmission, bus_pm) + push!(boundary_buses_distribution, bus_pmd) + end + # Convert to Julia Set - Note: membership checks are faster in sets (vs. vectors) in Julia + boundary_buses_transmission_set = Set(boundary_buses_transmission) + boundary_buses_distribution_set = Set(boundary_buses_distribution) + + # # ---- Transmission Power Balance --- + for i in _PM.ids(pm_model, :bus, nw=n) + if i in boundary_buses_transmission_set + constraint_transmission_power_balance_boundary(pmitd, i; nw_pmitd=n) + else + _PM.constraint_power_balance(pm_model, i, nw=n) + end + end + + # ---- Distribution Power Balance --- + for i in _PMD.ids(pmd_model, n, :bus) + if i in boundary_buses_distribution_set + constraint_distribution_power_balance_boundary(pmitd, i; nw_pmitd=n) + else + _PMD.constraint_mc_power_balance(pmd_model, i; nw=n) + end + end + + end + + # --- PM energy storage state constraint --- + network_ids_pm = sort(collect(_PM.nw_ids(pm_model))) + + n_1_pm = network_ids_pm[1] + for i in _PM.ids(pm_model, :storage, nw=n_1_pm) + _PM.constraint_storage_state(pm_model, i, nw=n_1_pm) + end + + for n_2_pm in network_ids_pm[2:end] + for i in _PM.ids(pm_model, :storage, nw=n_2_pm) + _PM.constraint_storage_state(pm_model, i, n_1_pm, n_2_pm) + end + n_1_pm = n_2_pm + end + + # --- PMD energy storage state constraint --- + network_ids_pmd = sort(collect(_PMD.nw_ids(pmd_model))) + + n_1_pmd = network_ids_pmd[1] + + for i in _PMD.ids(pmd_model, :storage; nw=n_1_pmd) + _PMD.constraint_storage_state(pmd_model, i; nw=n_1_pmd) + end + + for n_2_pmd in network_ids_pmd[2:end] + for i in _PMD.ids(pmd_model, :storage; nw=n_2_pmd) + _PMD.constraint_storage_state(pmd_model, i, n_1_pmd, n_2_pmd) + end + + n_1_pmd = n_2_pmd + end + + # ------------------------------------------------- + # --- PMITD(T&D) Cost Functions ------------------- + objective_itd_min_fuel_cost_storage(pmitd) + +end + + +""" + function build_mn_opfitd_storage_linear( + pmitd::AbstractIVRPowerModelITD + ) +Constructor for Multinetwork Integrated T&D Optimal Power Flow in current-voltage (IV) variable space with Storage OPF Dispatch with Linear Storage Model. +""" +function build_mn_opfitd_storage_linear(pmitd::AbstractIVRPowerModelITD) + + @error "IVR-IVRU formulation not yet supported for multinetwork storage problems." + throw(error()) + # TODO +end + + +""" + function build_mn_opfitd_storage_linear( + pmitd::AbstractBFPowerModelITD + ) +Constructor for Multinetwork Integrated T&D Optimal Power Flow for BF Models with Storage OPF Dispatch with Linear Storage Model. +""" +function build_mn_opfitd_storage_linear(pmitd::AbstractBFPowerModelITD) + + # Get Models + pm_model = _get_powermodel_from_powermodelitd(pmitd) + pmd_model = _get_powermodeldistribution_from_powermodelitd(pmitd) + + for (n, network) in nws(pmitd) + # PM(Transmission) Variables + _PM.variable_bus_voltage(pm_model, nw=n) + _PM.variable_gen_power(pm_model, nw=n) + _PM.variable_branch_power(pm_model, nw=n) + _PM.variable_branch_current(pm_model, nw=n) + _PM.variable_dcline_power(pm_model, nw=n) + _PM.variable_storage_power(pm_model, nw=n) + + # PMD(Distribution) Variables + _PMD.variable_mc_bus_voltage(pmd_model; nw=n) + _PMD.variable_mc_branch_current(pmd_model; nw=n) + _PMD.variable_mc_branch_power(pmd_model; nw=n) + _PMD.variable_mc_switch_power(pmd_model; nw=n) + _PMD.variable_mc_transformer_power(pmd_model; nw=n) + _PMD.variable_mc_generator_power(pmd_model; nw=n) + _PMD.variable_mc_load_power(pmd_model; nw=n) + _PMD.variable_mc_storage_power(pmd_model; nw=n) + + # PMITD (Boundary) Current Variables + variable_boundary_power(pmitd; nw=n) + + # --- PM(Transmission) Constraints --- + _PM.constraint_model_current(pm_model; nw=n) + + # reference buses (this only needs to happen for pm(transmission)) + for i in _PM.ids(pm_model, :ref_buses, nw=n) + _PM.constraint_theta_ref(pm_model, i, nw=n) + end + + for i in _PM.ids(pm_model, :storage, nw=n) + constraint_storage_losses_linear(pm_model, i, nw=n) + _PM.constraint_storage_thermal_limit(pm_model, i, nw=n) + end + + # PM branches + for i in _PM.ids(pm_model, :branch, nw=n) + _PM.constraint_power_losses(pm_model, i, nw=n) + _PM.constraint_voltage_magnitude_difference(pm_model, i, nw=n) + + _PM.constraint_voltage_angle_difference(pm_model, i, nw=n) + + _PM.constraint_thermal_limit_from(pm_model, i, nw=n) + _PM.constraint_thermal_limit_to(pm_model, i, nw=n) + end + + # PM DC lines + for i in _PM.ids(pm_model, :dcline, nw=n) + _PM.constraint_dcline_power_losses(pm_model, i, nw=n) + end + + # ------------------------------------------------- + # --- PMD(Distribution) Constraints --- + _PMD.constraint_mc_model_current(pmd_model; nw=n) + + # generators should be constrained before KCL, or Pd/Qd undefined + for i in _PMD.ids(pmd_model, n, :gen) + _PMD.constraint_mc_generator_power(pmd_model, i; nw=n) + end + + # loads should be constrained before KCL, or Pd/Qd undefined + for i in _PMD.ids(pmd_model, n, :load) + _PMD.constraint_mc_load_power(pmd_model, i; nw=n) + end + + for i in _PMD.ids(pmd_model, n, :storage) + constraint_mc_storage_losses_linear(pmd_model, i; nw=n) + _PMD.constraint_mc_storage_thermal_limit(pmd_model, i; nw=n) + end + + for i in _PMD.ids(pmd_model, n, :branch) + _PMD.constraint_mc_power_losses(pmd_model, i; nw=n) + _PMD.constraint_mc_model_voltage_magnitude_difference(pmd_model, i; nw=n) + + _PMD.constraint_mc_voltage_angle_difference(pmd_model, i; nw=n) + + _PMD.constraint_mc_thermal_limit_from(pmd_model, i; nw=n) + _PMD.constraint_mc_thermal_limit_to(pmd_model, i; nw=n) + _PMD.constraint_mc_ampacity_from(pmd_model, i; nw=n) + _PMD.constraint_mc_ampacity_to(pmd_model, i; nw=n) + end + + for i in _PMD.ids(pmd_model, n, :switch) + _PMD.constraint_mc_switch_state(pmd_model, i; nw=n) + _PMD.constraint_mc_switch_thermal_limit(pmd_model, i; nw=n) + _PMD.constraint_mc_switch_ampacity(pmd_model, i; nw=n) + end + + for i in _PMD.ids(pmd_model, n, :transformer) + _PMD.constraint_mc_transformer_power(pmd_model, i; nw=n) + end + + # ------------------------------------------------- + # --- PMITD(T&D) INDEPENDENT Constraints ---------- + + for i in ids(pmitd, :boundary; nw=n) + constraint_boundary_power(pmitd, i; nw=n) + constraint_boundary_voltage_magnitude(pmitd, i; nw=n) + constraint_boundary_voltage_angle(pmitd, i; nw=n) + end + + # ------------------------------------------------- + # --- PMITD(T&D) KCL Constraints ---------- + # Note: Both of these need to consider flow on boundaries if bus is connected to boundary + boundary_buses_transmission = Vector{Int}() # vector to store the boundary buses transmission + boundary_buses_distribution = Vector{Int}() # vector to store the boundary buses distribution + for j in ids(pmitd, :boundary; nw=n) + boundary_pmitd = ref(pmitd, n, :boundary, j) + bus_pm = boundary_pmitd["f_bus"] + bus_pmd = boundary_pmitd["t_bus"] + push!(boundary_buses_transmission, bus_pm) + push!(boundary_buses_distribution, bus_pmd) + end + # Convert to Julia Set - Note: membership checks are faster in sets (vs. vectors) in Julia + boundary_buses_transmission_set = Set(boundary_buses_transmission) + boundary_buses_distribution_set = Set(boundary_buses_distribution) + + # # ---- Transmission Power Balance --- + for i in _PM.ids(pm_model, :bus, nw=n) + if i in boundary_buses_transmission_set + constraint_transmission_power_balance_boundary(pmitd, i; nw_pmitd=n) + else + _PM.constraint_power_balance(pm_model, i, nw=n) + end + end + + # ---- Distribution Power Balance --- + for i in _PMD.ids(pmd_model, n, :bus) + if i in boundary_buses_distribution_set + constraint_distribution_power_balance_boundary(pmitd, i; nw_pmitd=n) + else + _PMD.constraint_mc_power_balance(pmd_model, i; nw=n) + end + end + + end + + # --- PM energy storage state constraint --- + network_ids_pm = sort(collect(_PM.nw_ids(pm_model))) + + n_1_pm = network_ids_pm[1] + for i in _PM.ids(pm_model, :storage, nw=n_1_pm) + _PM.constraint_storage_state(pm_model, i, nw=n_1_pm) + end + + for n_2_pm in network_ids_pm[2:end] + for i in _PM.ids(pm_model, :storage, nw=n_2_pm) + _PM.constraint_storage_state(pm_model, i, n_1_pm, n_2_pm) + end + n_1_pm = n_2_pm + end + + # --- PMD energy storage state constraint --- + network_ids_pmd = sort(collect(_PMD.nw_ids(pmd_model))) + + n_1_pmd = network_ids_pmd[1] + + for i in _PMD.ids(pmd_model, :storage; nw=n_1_pmd) + _PMD.constraint_storage_state(pmd_model, i; nw=n_1_pmd) + end + + for n_2_pmd in network_ids_pmd[2:end] + for i in _PMD.ids(pmd_model, :storage; nw=n_2_pmd) + _PMD.constraint_storage_state(pmd_model, i, n_1_pmd, n_2_pmd) + end + + n_1_pmd = n_2_pmd + end + + # ------------------------------------------------- + # --- PMITD(T&D) Cost Functions ------------------- + objective_itd_min_fuel_cost_storage(pmitd) + +end + + +""" + function build_mn_opfitd_storage_linear( + pmitd::AbstractLNLBFPowerModelITD + ) +Constructor for Multinetwork Integrated T&D Optimal Power Flow for L/NL to BF with Storage OPF Dispatch with Linear Storage Model. +""" +function build_mn_opfitd_storage_linear(pmitd::AbstractLNLBFPowerModelITD) + + # Get Models + pm_model = _get_powermodel_from_powermodelitd(pmitd) + pmd_model = _get_powermodeldistribution_from_powermodelitd(pmitd) + + for (n, network) in nws(pmitd) + # PM(Transmission) Variables + _PM.variable_bus_voltage(pm_model, nw=n) + _PM.variable_gen_power(pm_model, nw=n) + _PM.variable_branch_power(pm_model, nw=n) + _PM.variable_dcline_power(pm_model, nw=n) + _PM.variable_storage_power(pm_model, nw=n) + + # PMD(Distribution) Variables + _PMD.variable_mc_bus_voltage(pmd_model; nw=n) + _PMD.variable_mc_branch_current(pmd_model; nw=n) + _PMD.variable_mc_branch_power(pmd_model; nw=n) + _PMD.variable_mc_transformer_power(pmd_model; nw=n) + _PMD.variable_mc_switch_power(pmd_model; nw=n) + _PMD.variable_mc_generator_power(pmd_model; nw=n) + _PMD.variable_mc_load_power(pmd_model; nw=n) + _PMD.variable_mc_storage_power(pmd_model; nw=n) + + # PMITD (Boundary) Variables + variable_boundary_power(pmitd; nw=n) + + # --- PM(Transmission) Constraints --- + _PM.constraint_model_voltage(pm_model, nw=n) + + # reference buses (this only needs to happen for pm(transmission)) + for i in _PM.ids(pm_model, :ref_buses, nw=n) + _PM.constraint_theta_ref(pm_model, i, nw=n) + end + + for i in _PM.ids(pm_model, :storage, nw=n) + constraint_storage_losses_linear(pm_model, i, nw=n) + _PM.constraint_storage_thermal_limit(pm_model, i, nw=n) + end + + # PM branches + for i in _PM.ids(pm_model, :branch, nw=n) + _PM.constraint_ohms_yt_from(pm_model, i, nw=n) + _PM.constraint_ohms_yt_to(pm_model, i, nw=n) + + _PM.constraint_voltage_angle_difference(pm_model, i, nw=n) + + _PM.constraint_thermal_limit_from(pm_model, i, nw=n) + _PM.constraint_thermal_limit_to(pm_model, i, nw=n) + end + + # PM DC lines + for i in _PM.ids(pm_model, :dcline, nw=n) + _PM.constraint_dcline_power_losses(pm_model, i, nw=n) + end + + # ------------------------------------------------- + # --- PMD(Distribution) Constraints --- + _PMD.constraint_mc_model_current(pmd_model; nw=n) + + # generators should be constrained before KCL, or Pd/Qd undefined + for i in _PMD.ids(pmd_model, n, :gen) + _PMD.constraint_mc_generator_power(pmd_model, i; nw=n) + end + + # loads should be constrained before KCL, or Pd/Qd undefined + for i in _PMD.ids(pmd_model, n, :load) + _PMD.constraint_mc_load_power(pmd_model, i; nw=n) + end + + for i in _PMD.ids(pmd_model, n, :storage) + constraint_mc_storage_losses_linear(pmd_model, i; nw=n) + _PMD.constraint_mc_storage_thermal_limit(pmd_model, i; nw=n) + end + + for i in _PMD.ids(pmd_model, n, :branch) + _PMD.constraint_mc_power_losses(pmd_model, i; nw=n) + _PMD.constraint_mc_model_voltage_magnitude_difference(pmd_model, i; nw=n) + + _PMD.constraint_mc_voltage_angle_difference(pmd_model, i; nw=n) + + _PMD.constraint_mc_thermal_limit_from(pmd_model, i; nw=n) + _PMD.constraint_mc_thermal_limit_to(pmd_model, i; nw=n) + _PMD.constraint_mc_ampacity_from(pmd_model, i; nw=n) + _PMD.constraint_mc_ampacity_to(pmd_model, i; nw=n) + end + + for i in _PMD.ids(pmd_model, n, :switch) + _PMD.constraint_mc_switch_state(pmd_model, i; nw=n) + _PMD.constraint_mc_switch_thermal_limit(pmd_model, i; nw=n) + _PMD.constraint_mc_switch_ampacity(pmd_model, i; nw=n) + end + + for i in _PMD.ids(pmd_model, n, :transformer) + _PMD.constraint_mc_transformer_power(pmd_model, i; nw=n) + end + + # ------------------------------------------------- + # --- PMITD(T&D) INDEPENDENT Constraints ---------- + + for i in ids(pmitd, :boundary; nw=n) + constraint_boundary_power(pmitd, i; nw=n) + constraint_boundary_voltage_magnitude(pmitd, i; nw=n) + constraint_boundary_voltage_angle(pmitd, i; nw=n) + end + + # ------------------------------------------------- + # --- PMITD(T&D) KCL Constraints ---------- + # Note: Both of these need to consider flow on boundaries if bus is connected to boundary + boundary_buses_transmission = Vector{Int}() # vector to store the boundary buses transmission + boundary_buses_distribution = Vector{Int}() # vector to store the boundary buses distribution + for j in ids(pmitd, :boundary; nw=n) + boundary_pmitd = ref(pmitd, n, :boundary, j) + bus_pm = boundary_pmitd["f_bus"] + bus_pmd = boundary_pmitd["t_bus"] + push!(boundary_buses_transmission, bus_pm) + push!(boundary_buses_distribution, bus_pmd) + end + # Convert to Julia Set - Note: membership checks are faster in sets (vs. vectors) in Julia + boundary_buses_transmission_set = Set(boundary_buses_transmission) + boundary_buses_distribution_set = Set(boundary_buses_distribution) + + # # ---- Transmission Power Balance --- + for i in _PM.ids(pm_model, :bus, nw=n) + if i in boundary_buses_transmission_set + constraint_transmission_power_balance_boundary(pmitd, i; nw_pmitd=n) + else + _PM.constraint_power_balance(pm_model, i, nw=n) + end + end + + # ---- Distribution Power Balance --- + for i in _PMD.ids(pmd_model, n, :bus) + if i in boundary_buses_distribution_set + constraint_distribution_power_balance_boundary(pmitd, i; nw_pmitd=n) + else + _PMD.constraint_mc_power_balance(pmd_model, i; nw=n) + end + end + + end + + # --- PM energy storage state constraint --- + network_ids_pm = sort(collect(_PM.nw_ids(pm_model))) + + n_1_pm = network_ids_pm[1] + for i in _PM.ids(pm_model, :storage, nw=n_1_pm) + _PM.constraint_storage_state(pm_model, i, nw=n_1_pm) + end + + for n_2_pm in network_ids_pm[2:end] + for i in _PM.ids(pm_model, :storage, nw=n_2_pm) + _PM.constraint_storage_state(pm_model, i, n_1_pm, n_2_pm) + end + n_1_pm = n_2_pm + end + + # --- PMD energy storage state constraint --- + network_ids_pmd = sort(collect(_PMD.nw_ids(pmd_model))) + + n_1_pmd = network_ids_pmd[1] + + for i in _PMD.ids(pmd_model, :storage; nw=n_1_pmd) + _PMD.constraint_storage_state(pmd_model, i; nw=n_1_pmd) + end + + for n_2_pmd in network_ids_pmd[2:end] + for i in _PMD.ids(pmd_model, :storage; nw=n_2_pmd) + _PMD.constraint_storage_state(pmd_model, i, n_1_pmd, n_2_pmd) + end + + n_1_pmd = n_2_pmd + end + + # ------------------------------------------------- + # --- PMITD(T&D) Cost Functions ------------------- + objective_itd_min_fuel_cost_storage(pmitd) +end diff --git a/test/opfitd_storage.jl b/test/opfitd_storage.jl index ed32540..4242261 100644 --- a/test/opfitd_storage.jl +++ b/test/opfitd_storage.jl @@ -160,7 +160,7 @@ # with storage cost problem pmitd_result_strg = solve_opfitd_storage(pmitd_data, pmitd_type, ipopt) - @test isapprox(pmitd_result_strg["objective"], 17976.959438016635+(pmitd_result_strg["solution"]["it"]["pmd"]["storage"]["3bus_bal_battery.s1"]["sd"]*strg_cost/100000); atol = 1e-3) + @test isapprox(pmitd_result_strg["objective"], 17978.20946; atol = 1e-3) end diff --git a/test/opfitd_storage_linear.jl b/test/opfitd_storage_linear.jl new file mode 100644 index 0000000..998f63c --- /dev/null +++ b/test/opfitd_storage_linear.jl @@ -0,0 +1,137 @@ +@info "running ITD storage tests for linear model" + +@testset "test/opfitd_storage_linear.jl" begin + + @testset "solve_opfitd_storage: Balanced case5-case3 With Battery ACP-ACPU " begin + pm_file = joinpath(dirname(trans_path), "case5_withload.m") + pmd_file = joinpath(dirname(dist_path), "case3_balanced_withBattery.dss") + pmitd_file = joinpath(dirname(bound_path), "case5_case3_bal_battery.json") + pmitd_type = NLPowerModelITD{ACPPowerModel, ACPUPowerModel} + pmitd_data = parse_files(pm_file, pmd_file, pmitd_file) + + # cost to assign to energy storage + # Units $/kWh + strg_cost = 0.025 + + # add cost to storages in PMD + for (st_name, st_data) in pmitd_data["it"]["pmd"]["storage"] + st_data["cost"] = strg_cost + end + + # with storage cost problem + pmitd_result_strg = solve_opfitd_storage_linear(pmitd_data, pmitd_type, ipopt) + + end + + @testset "solve_opfitd_storage with Transmission storage: Balanced case5-case3 With Battery ACP-ACPU " begin + pm_file = joinpath(dirname(trans_path), "case5_withload_strg.m") + pmd_file = joinpath(dirname(dist_path), "case3_balanced_withBattery.dss") + pmitd_file = joinpath(dirname(bound_path), "case5_case3_bal_battery.json") + pmitd_type = NLPowerModelITD{ACPPowerModel, ACPUPowerModel} + pmitd_result_strg = solve_opfitd_storage_linear(pm_file, pmd_file, pmitd_file, pmitd_type, ipopt) + @test pmitd_result_strg["termination_status"] == LOCALLY_SOLVED + end + + @testset "solve_mn_opfitd_storage with Transmission storage: Multinetwork Balanced case5-case3 With Battery ACP-ACPU " begin + pm_file = joinpath(dirname(trans_path), "case5_withload_strg.m") + pmd_file = joinpath(dirname(dist_path), "case3_balanced_withBattery_mn.dss") + pmitd_file = joinpath(dirname(bound_path), "case5_case3_bal_battery_mn.json") + pmitd_type = NLPowerModelITD{ACPPowerModel, ACPUPowerModel} + pmitd_result_strg = solve_mn_opfitd_storage_linear(pm_file, pmd_file, pmitd_file, pmitd_type, ipopt) + @test pmitd_result_strg["termination_status"] == LOCALLY_SOLVED + end + + @testset "solve_mn_opfitd_storage with Transmission storage: Multinetwork Balanced case5-case3 With Battery ACR-ACRU " begin + pm_file = joinpath(dirname(trans_path), "case5_withload_strg.m") + pmd_file = joinpath(dirname(dist_path), "case3_balanced_withBattery_mn.dss") + pmitd_file = joinpath(dirname(bound_path), "case5_case3_bal_battery_mn.json") + pmitd_type = NLPowerModelITD{ACRPowerModel, ACRUPowerModel} + pmitd_result_strg = solve_mn_opfitd_storage_linear(pm_file, pmd_file, pmitd_file, pmitd_type, ipopt) + @test pmitd_result_strg["termination_status"] == LOCALLY_SOLVED + end + + @testset "solve_mn_opfitd_storage with Transmission storage: Multinetwork Balanced case5-case3 With Battery ACR-FBSUBF " begin + pm_file = joinpath(dirname(trans_path), "case5_withload_strg.m") + pmd_file = joinpath(dirname(dist_path), "case3_balanced_withBattery_mn.dss") + pmitd_file = joinpath(dirname(bound_path), "case5_case3_bal_battery_mn.json") + pmitd_type = NLBFPowerModelITD{ACRPowerModel, FBSUBFPowerModel} + pmitd_result_strg = solve_mn_opfitd_storage_linear(pm_file, pmd_file, pmitd_file, pmitd_type, ipopt) + @test pmitd_result_strg["termination_status"] == LOCALLY_SOLVED + end + + @testset "solve_mn_opfitd_storage with Transmission storage: Multinetwork Balanced case5-case3 With Battery ACR-FOTR " begin + pm_file = joinpath(dirname(trans_path), "case5_withload_strg.m") + pmd_file = joinpath(dirname(dist_path), "case3_balanced_withBattery_mn.dss") + pmitd_file = joinpath(dirname(bound_path), "case5_case3_bal_battery_mn.json") + pmitd_type = NLFOTPowerModelITD{ACRPowerModel, FOTRUPowerModel} + pmitd_result_strg = solve_mn_opfitd_storage_linear(pm_file, pmd_file, pmitd_file, pmitd_type, ipopt) + @test pmitd_result_strg["termination_status"] == LOCALLY_SOLVED + end + + @testset "solve_mn_opfitd_storage with Transmission storage: Multinetwork Balanced case5-case3 With Battery ACP-FOTP " begin + pm_file = joinpath(dirname(trans_path), "case5_withload_strg.m") + pmd_file = joinpath(dirname(dist_path), "case3_balanced_withBattery_mn.dss") + pmitd_file = joinpath(dirname(bound_path), "case5_case3_bal_battery_mn.json") + pmitd_type = NLFOTPowerModelITD{ACPPowerModel, FOTPUPowerModel} + pmitd_result_strg = solve_mn_opfitd_storage_linear(pm_file, pmd_file, pmitd_file, pmitd_type, ipopt) + @test pmitd_result_strg["termination_status"] == LOCALLY_SOLVED + end + + @testset "solve_mn_opfitd_storage with Transmission storage: Multinetwork Balanced case5-case3 With Battery BFA-LinDist3Flow " begin + pm_file = joinpath(dirname(trans_path), "case5_withload_strg.m") + pmd_file = joinpath(dirname(dist_path), "case3_balanced_withBattery_mn.dss") + pmitd_file = joinpath(dirname(bound_path), "case5_case3_bal_battery_mn.json") + pmitd_type = BFPowerModelITD{BFAPowerModel, LinDist3FlowPowerModel} + pmitd_result_strg = solve_mn_opfitd_storage_linear(pm_file, pmd_file, pmitd_file, pmitd_type, ipopt) + @test pmitd_result_strg["termination_status"] == LOCALLY_SOLVED + end + + @testset "solve_mn_opfitd_storage: Multinetwork Balanced case5-case3 With Battery in Transmission ACP-ACPU " begin + pm_file = joinpath(dirname(trans_path), "case5_withload_strg.m") + pmd_file = joinpath(dirname(dist_path), "case3_unbalanced_withoutgen_mn_diff3.dss") + pmitd_file = joinpath(dirname(bound_path), "case5_case3_unbal_nogen_mn_diff.json") + pmitd_type = NLPowerModelITD{ACPPowerModel, ACPUPowerModel} + pmitd_data = parse_files(pm_file, pmd_file, pmitd_file; multinetwork=true) + + # cost to assign to energy storage + # Units $/kWh + strg_cost = 0.0025 + + # add cost to storages in PM + for (nw_id, nw_data) in pmitd_data["it"]["pm"]["nw"] + for (st_name, st_data) in nw_data["storage"] + st_data["cost"] = strg_cost + end + end + + # varying load + pm_residential_prof_p = [0.293391082 0.139045672 0.469951253 0.432787275 0.392755291 0.548138118 0.466415402 0.862109886 0.794416422 0.643809168 0.709608241 0.704126411 0.963765953 0.575280879 0.611590304 0.488662676 0.430361821 0.87041275 1.306424968 1.307958207 0.928544462 0.572729047 0.552861937 0.346550725] + pm_residential_prof_q = [0.247949186 0.171635443 0.183848151 0.227990471 0.192748178 0.319480381 0.180830891 0.415389708 0.224002485 0.322713139 0.241254927 0.300147051 0.298347195 0.186199635 0.22251068 0.303541769 0.185154905 0.432356924 0.425774291 0.50190239 0.294529233 0.320796731 0.14920804 0.136169507] + + pm_community_prof_p = [0.284287776 0.268623055 0.227327982 0.397346384 0.235818756 0.20823994 0.435025687 0.569052345 0.837116674 0.895062747 0.966074521 1.01396914 0.87299946 0.900887326 0.92727811 0.839239282 1.068426425 1.066425748 1.392200763 1.268324779 1.400386576 1.218758216 1.001009647 0.538582071] + pm_community_prof_q = [0.124401652 0.149265062 0.112556201 0.226246144 0.112703236 0.145539091 0.24876825 0.304559383 0.412502283 0.360553174 0.361034872 0.481935002 0.515125962 0.507287921 0.398147224 0.463007061 0.507275554 0.534884845 0.508290756 0.531864342 0.549012302 0.519189403 0.452179964 0.218622959] + + pm_commercial_prof_p = [0.375187827 0.286661472 0.289096564 0.22256313 0.271830914 0.292079786 0.236396231 0.93119696 1.056699043 1.049660067 1.009438919 1.082846718 1.059419722 1.039691666 1.025920791 1.032047161 1.085377745 1.04104745 1.022642982 0.924980776 0.319601874 0.347188983 0.289086358 0.276695314] + pm_commercial_prof_q = [0.167332967 0.299127427 0.279260644 0.207282394 0.159990914 0.219440511 0.178642786 0.559090945 0.514151587 0.540382353 0.531009283 0.518096858 0.514654825 0.591549766 0.525647471 0.549013224 0.520947523 0.557888204 0.553033554 0.511838587 0.161558414 0.257372041 0.20559509 0.295150131] + + for (nw, nw_data) in pmitd_data["it"]["pm"]["nw"] + for (load, load_data) in nw_data["load"] + if (load == "1") + load_data["pd"] = load_data["pd"]*pm_commercial_prof_p[parse(Int64,nw)] + load_data["qd"] = load_data["qd"]*pm_commercial_prof_q[parse(Int64,nw)] + elseif (load == "2") + load_data["pd"] = load_data["pd"]*pm_community_prof_p[parse(Int64,nw)] + load_data["qd"] = load_data["qd"]*pm_community_prof_q[parse(Int64,nw)] + elseif (load == "3") + load_data["pd"] = load_data["pd"]*pm_residential_prof_p[parse(Int64,nw)] + load_data["qd"] = load_data["qd"]*pm_residential_prof_q[parse(Int64,nw)] + end + end + end + + pmitd_result_strg = solve_mn_opfitd_storage(pmitd_data, pmitd_type, ipopt) + @test isapprox(pmitd_result_strg["solution"]["it"]["pm"]["nw"]["10"]["storage"]["1"]["sd"]*100, 37.8596; atol = 1e-1) + + end + +end diff --git a/test/runtests.jl b/test/runtests.jl index 09d5f0f..14ad6d6 100755 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -42,6 +42,7 @@ PowerModelsITD.silence!() include("opfitd_dmld.jl") include("opfitd_solution.jl") include("opfitd_storage.jl") + include("opfitd_storage_linear.jl") include("solve_x.jl") include("pfitd_mn.jl") end