From 17f3abd98f0b779d823dce48b7e1101cd08524ec Mon Sep 17 00:00:00 2001 From: Tucker Oddleifson Date: Mon, 21 Oct 2024 22:32:04 -0600 Subject: [PATCH] Combined the extend.jl and power_flow.jl files --- src/REopt.jl | 1 - .../electric_utility_constraints.jl | 4 +- src/lindistflow/extend.jl | 155 ------------------ src/lindistflow/power_flow.jl | 150 ++++++++++++++++- 4 files changed, 151 insertions(+), 159 deletions(-) delete mode 100644 src/lindistflow/extend.jl diff --git a/src/REopt.jl b/src/REopt.jl index 0f0cc21b8..5cfe121c4 100644 --- a/src/REopt.jl +++ b/src/REopt.jl @@ -197,7 +197,6 @@ include("outagesim/outage_simulator.jl") include("outagesim/backup_reliability.jl") include("lindistflow/power_flow.jl") -include("lindistflow/extend.jl") include("mpc/results.jl") include("mpc/model.jl") diff --git a/src/constraints/electric_utility_constraints.jl b/src/constraints/electric_utility_constraints.jl index c62e92816..31e712663 100644 --- a/src/constraints/electric_utility_constraints.jl +++ b/src/constraints/electric_utility_constraints.jl @@ -3,7 +3,7 @@ function add_export_constraints(m, p; _n="") ##Constraint (8e): Production export and curtailment no greater than production if string(p.s.site.node) != p.s.settings.facilitymeter_node - print("\n Updated***Test 4, Adding constraint 8e to node $(p.s.site.node)") + print("\n Updated***4000 timesteps, Adding constraint 8e to node $(p.s.site.node)") @constraint(m, [t in p.techs.elec, ts in p.time_steps_with_grid], p.production_factor[t,ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t,ts] >= sum(m[Symbol("dvProductionToGrid"*_n)][t, u, ts] for u in p.export_bins_by_tech[t]) + @@ -11,7 +11,7 @@ function add_export_constraints(m, p; _n="") ) else print("\n Not adding constraint 8e to the facility meter node, node $(p.s.site.node)") - TempVector = collect(201:8760) # implemented only for initial testing; TODO: remove this after testing + TempVector = collect(4001:8760) # implemented only for initial testing; TODO: remove this after testing @constraint(m, [t in p.techs.elec, ts in p.time_steps_with_grid], m[Symbol("dvRatedProduction"*_n)][t,ts] == 0) @constraint(m, [t in p.techs.elec, ts in TempVector], sum(m[Symbol("dvProductionToGrid"*_n)][t, u, ts] for u in p.export_bins_by_tech[t]) == 0) @constraint(m, [t in p.techs.elec, ts in p.time_steps_with_grid], m[Symbol("dvCurtail"*_n)][t, ts] == 0) diff --git a/src/lindistflow/extend.jl b/src/lindistflow/extend.jl deleted file mode 100644 index 72e51b189..000000000 --- a/src/lindistflow/extend.jl +++ /dev/null @@ -1,155 +0,0 @@ - -# Additional code for interfacing between the power_flow.jl and other REopt code - -function build_power_flow!(m::JuMP.AbstractModel, p::PowerFlowInputs, ps::Array{REoptInputs{Scenario}, 1}; - make_import_export_complementary::Bool=true - ) - power_flow_add_variables(m, p) - add_expressions(m, ps) - constrain_power_balance(m, p) - constrain_substation_voltage(m, p) - create_line_variables(m, p) - constrain_loads(m, p, ps) - # Note: the constrain_KVL(m, p) function is called in the microgrid.jl file - if make_import_export_complementary - add_complementary_constraints(m, ps) - end -end - -#2. Add power flow constraints to m, -# - set Pⱼ's = -1 * (dvGridPurchase_j - dvProductionToGrid_j) -function add_expressions(m::JuMP.AbstractModel, ps::Array{REoptInputs{Scenario}, 1}) - for p in ps - _n = string("_", p.s.site.node) - if string(p.s.site.node) != p.s.settings.facilitymeter_node - print("\n Applying total export equation for node $(p.s.site.node)") - m[Symbol("TotalExport"*_n)] = @expression(m, [ts in p.time_steps], - sum( - m[Symbol("dvProductionToGrid"*_n)][t,u,ts] - for t in p.techs.elec, u in p.export_bins_by_tech[t] - ) - + sum(m[Symbol("dvStorageToGrid"*_n)][b,ts] for b in p.s.storage.types.all )# added this line to include battery export in the total export - ) - else - # set the total node export to 0 for the facility meter grid, because that node has no techs - # also, for that node, the dvProductionToGrid is used for the grid export benefits and set to the powerflow of the substation line when flow on that line is negative - print("\n Setting the total export for node $(p.s.site.node) to zero. This node is the facility meter node.") - - dv = "TotalExport_"*p.s.settings.facilitymeter_node - m[Symbol(dv)] = @variable(m, [p.time_steps], base_name=dv, lower_bound =0) - - @constraint(m, [t in p.time_steps], m[Symbol("TotalExport_"*p.s.settings.facilitymeter_node)][t] == 0) - - end - end -end - -# TODO add complementary constraint to UL for dvProductionToGrid_ and dvGridPurchase_ (don't want it in LL s.t. it stays linear) -function add_complementary_constraints(m::JuMP.AbstractModel, ps::Array{REoptInputs{Scenario}, 1}) - for p in ps - _n = string("_", p.s.site.node) - #if string(p.s.site.node) != "15" # the complementary constraint has an affect when it is applied to the facility meter node - #print("\n Adding the complementary constraint to node $(p.s.site.node)") - for (i, e) in zip(m[Symbol("dvGridPurchase"*_n)], m[Symbol("TotalExport"*_n)]) - @constraint(m, - [i, e] in MOI.SOS1([1.0, 2.0]) - ) - end - #else - # print("\n Not adding the complementary constraint to the facility meter node, node $(p.s.site.node)") - #end - end -end - - -function constrain_loads(m::JuMP.AbstractModel, p::PowerFlowInputs, ps::Array{REoptInputs{Scenario}, 1}) - - Pⱼ = m[:Pⱼ] - Qⱼ = m[:Qⱼ] - - reopt_nodes = [p.s.site.node for p in ps] - print("\n The REopt nodes being applied in extend.jl are: $(reopt_nodes) \n") - - # Note: positive values are injections - #print("\n p.Pload is: $(p.Pload) ") - - for j in p.busses - if j in keys(p.Pload) - #print("\n the j variable is: $(j) ") - if parse(Int, j) in reopt_nodes - if j != "15" - #print("\n j is not 15") - @constraint(m, [t in 1:p.Ntimesteps], - Pⱼ[j,t] == 1e3/p.Sbase * ( # 1e3 b/c REopt values in kW - m[Symbol("TotalExport_" * j)][t] - - m[Symbol("dvGridPurchase_" * j)][t] - ) - ) - else - print("\n j is 15 and the j variable is: $(j)") - @constraint(m, [t in 1:p.Ntimesteps], Pⱼ["15",t] == 0) - end - else - # Power balance for nodes that have loads defined through ldf, but not in the REopt inputs - @constraint(m, [t in 1:p.Ntimesteps], - Pⱼ[j,t] == -p.Pload[j][t] - ) - end - elseif j != p.substation_bus - @constraint(m, [t in 1:p.Ntimesteps], - Pⱼ[j,t] == 0 - ) - end - - if j in keys(p.Qload) - if parse(Int, j) in reopt_nodes - if j != "15" - #print("\n (reactive power) j is not 15") - @constraint(m, [t in 1:p.Ntimesteps], - Qⱼ[j,t] == 1e3/p.Sbase * p.pf * ( # 1e3 b/c REopt values in kW - m[Symbol("TotalExport_" * j)][t] - - m[Symbol("dvGridPurchase_" * j)][t] - ) - ) - else - print("\n (reactive power) j is 15 and the j variable is: $(j)") - @constraint(m, [t in 1:p.Ntimesteps], Qⱼ["15",t] == 0) - end - else - @constraint(m, [t in 1:p.Ntimesteps], - Qⱼ[j,t] == -p.Qload[j][t] - ) - end - elseif j != p.substation_bus - @constraint(m, [t in 1:p.Ntimesteps], - Qⱼ[j,t] == 0 - ) - end - end - p.Nequality_cons += 2 * (p.Nnodes - 1) * p.Ntimesteps - - # Constrain loads on the transformers - P = m[:Pᵢⱼ] - Q = m[:Qᵢⱼ] - all_transformers = [] - # Define decision variable for the transformer maximum kVa - for i in keys(p.transformers) - if p.transformers[i]["Transformer Side"] == "downstream" - push!(all_transformers, i) - end - end - print("\n The all_transformers variable is: ") - print(all_transformers) - @variable(m, transformer_max_kva[all_transformers] >= 0) - - # Apply the transformer max kva to the constraints - for i in keys(p.transformers) - if p.transformers[i]["Transformer Side"] == "downstream" - DirectlyUpstreamNode = i_to_j(i, p) - transformer_internal_line = string(DirectlyUpstreamNode[1])*"-"*string(i) - @constraint(m, [T in 1:p.Ntimesteps], P[transformer_internal_line, T] + Q[transformer_internal_line, T] <= ((m[:transformer_max_kva][i]*1000)/p.Sbase)) - @constraint(m, [T in 1:p.Ntimesteps], P[transformer_internal_line, T] + Q[transformer_internal_line, T] >= -((m[:transformer_max_kva][i]*1000)/p.Sbase)) - end - end -end - diff --git a/src/lindistflow/power_flow.jl b/src/lindistflow/power_flow.jl index 213bf6afd..3ac145457 100644 --- a/src/lindistflow/power_flow.jl +++ b/src/lindistflow/power_flow.jl @@ -1,4 +1,6 @@ - +#= +Acknowledgement: This code is based on code from the LinDistFlow.jl package +=# """ mutable struct PowerFlowInputs @@ -183,6 +185,152 @@ PowerFlowInputs( ) end +function build_power_flow!(m::JuMP.AbstractModel, p::PowerFlowInputs, ps::Array{REoptInputs{Scenario}, 1}; + make_import_export_complementary::Bool=true) + + power_flow_add_variables(m, p) + add_expressions(m, ps) + constrain_power_balance(m, p) + constrain_substation_voltage(m, p) + create_line_variables(m, p) + constrain_loads(m, p, ps) + # Note: the constrain_KVL(m, p) function is called in the microgrid.jl file + if make_import_export_complementary + add_complementary_constraints(m, ps) + end +end + +#2. Add power flow constraints to m, +# - set Pⱼ's = -1 * (dvGridPurchase_j - dvProductionToGrid_j) +function add_expressions(m::JuMP.AbstractModel, ps::Array{REoptInputs{Scenario}, 1}) +for p in ps + _n = string("_", p.s.site.node) + if string(p.s.site.node) != p.s.settings.facilitymeter_node + print("\n Applying total export equation for node $(p.s.site.node)") + m[Symbol("TotalExport"*_n)] = @expression(m, [ts in p.time_steps], + sum( + m[Symbol("dvProductionToGrid"*_n)][t,u,ts] + for t in p.techs.elec, u in p.export_bins_by_tech[t] + ) + + sum(m[Symbol("dvStorageToGrid"*_n)][b,ts] for b in p.s.storage.types.all )# added this line to include battery export in the total export + ) + else + # set the total node export to 0 for the facility meter grid, because that node has no techs + # also, for that node, the dvProductionToGrid is used for the grid export benefits and set to the powerflow of the substation line when flow on that line is negative + print("\n Setting the total export for node $(p.s.site.node) to zero. This node is the facility meter node.") + + dv = "TotalExport_"*p.s.settings.facilitymeter_node + m[Symbol(dv)] = @variable(m, [p.time_steps], base_name=dv, lower_bound =0) + + @constraint(m, [t in p.time_steps], m[Symbol("TotalExport_"*p.s.settings.facilitymeter_node)][t] == 0) + + end +end +end + +# TODO add complementary constraint to UL for dvProductionToGrid_ and dvGridPurchase_ (don't want it in LL s.t. it stays linear) +function add_complementary_constraints(m::JuMP.AbstractModel, ps::Array{REoptInputs{Scenario}, 1}) +for p in ps + _n = string("_", p.s.site.node) + print("\n Adding the complementary constraint to node $(p.s.site.node)") + for (i, e) in zip(m[Symbol("dvGridPurchase"*_n)], m[Symbol("TotalExport"*_n)]) + @constraint(m, + [i, e] in MOI.SOS1([1.0, 2.0]) + ) + end +end +end + + +function constrain_loads(m::JuMP.AbstractModel, p::PowerFlowInputs, ps::Array{REoptInputs{Scenario}, 1}) + +Pⱼ = m[:Pⱼ] +Qⱼ = m[:Qⱼ] + +reopt_nodes = [p.s.site.node for p in ps] + +# Note: positive values are injections + +for j in p.busses + if j in keys(p.Pload) + #print("\n Debugging: the j variable is: $(j) ") + if parse(Int, j) in reopt_nodes + if j != "15" + #print("\n Debuggin: j is not 15") + @constraint(m, [t in 1:p.Ntimesteps], + Pⱼ[j,t] == 1e3/p.Sbase * ( # 1e3 b/c REopt values in kW + m[Symbol("TotalExport_" * j)][t] + - m[Symbol("dvGridPurchase_" * j)][t] + ) + ) + else + print("\n j is 15 and the j variable is: $(j)") + @constraint(m, [t in 1:p.Ntimesteps], Pⱼ["15",t] == 0) + end + else + # This constraint is for the power balance for nodes that have loads defined through ldf, but not in the REopt inputs + @constraint(m, [t in 1:p.Ntimesteps], + Pⱼ[j,t] == -p.Pload[j][t] + ) + end + elseif j != p.substation_bus + @constraint(m, [t in 1:p.Ntimesteps], + Pⱼ[j,t] == 0 + ) + end + + if j in keys(p.Qload) + if parse(Int, j) in reopt_nodes + if j != "15" + #print("\n Debuggin: (reactive power) j is not 15") + @constraint(m, [t in 1:p.Ntimesteps], + Qⱼ[j,t] == 1e3/p.Sbase * p.pf * ( # 1e3 b/c REopt values in kW + m[Symbol("TotalExport_" * j)][t] + - m[Symbol("dvGridPurchase_" * j)][t] + ) + ) + else + print("\n (reactive power) j is 15 and the j variable is: $(j)") + @constraint(m, [t in 1:p.Ntimesteps], Qⱼ["15",t] == 0) + end + else + @constraint(m, [t in 1:p.Ntimesteps], + Qⱼ[j,t] == -p.Qload[j][t] + ) + end + elseif j != p.substation_bus + @constraint(m, [t in 1:p.Ntimesteps], + Qⱼ[j,t] == 0 + ) + end +end +p.Nequality_cons += 2 * (p.Nnodes - 1) * p.Ntimesteps + +# Constrain loads on the transformers +P = m[:Pᵢⱼ] +Q = m[:Qᵢⱼ] +all_transformers = [] +# Define decision variable for the transformer maximum kVa +for i in keys(p.transformers) + if p.transformers[i]["Transformer Side"] == "downstream" + push!(all_transformers, i) + end +end +print("\n The all_transformers variable is: ") +print(all_transformers) +@variable(m, transformer_max_kva[all_transformers] >= 0) + +# Apply the transformer max kva to the constraints +for i in keys(p.transformers) + if p.transformers[i]["Transformer Side"] == "downstream" + DirectlyUpstreamNode = i_to_j(i, p) + transformer_internal_line = string(DirectlyUpstreamNode[1])*"-"*string(i) + @constraint(m, [T in 1:p.Ntimesteps], P[transformer_internal_line, T] + Q[transformer_internal_line, T] <= ((m[:transformer_max_kva][i]*1000)/p.Sbase)) + @constraint(m, [T in 1:p.Ntimesteps], P[transformer_internal_line, T] + Q[transformer_internal_line, T] >= -((m[:transformer_max_kva][i]*1000)/p.Sbase)) + end +end +end + function power_flow_add_variables(m, p::PowerFlowInputs) T = 1:p.Ntimesteps