Skip to content

Commit

Permalink
v0.3.2 single phase shunt_susceptance and Results.shadow_prices
Browse files Browse the repository at this point in the history
  • Loading branch information
@NLaws
NLaws committed May 14, 2023
1 parent 9e4ad22 commit 46dee4e
Show file tree
Hide file tree
Showing 5 changed files with 281 additions and 6 deletions.
4 changes: 4 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,9 @@
# BranchFlowModel Changelog

## v0.3.2
- add `shunt_susceptance` times `vsqrd` to flow balance equations in single phase model
- add `Results.shadow_prices`

## v0.3.1
- change JuMP compat to "1" (was "1.9") for better cross-compatibility

Expand Down
4 changes: 2 additions & 2 deletions Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "BranchFlowModel"
uuid = "73c867df-75f8-459f-abd8-059b58de1e18"
authors = ["Nick Laws"]
version = "0.3.1"
version = "0.3.2"

[deps]
CommonOPF = "dad7ea18-2b21-482e-81c1-e84414cc4b03"
Expand All @@ -13,7 +13,7 @@ MetaGraphs = "626554b9-1ddb-594c-aa3c-2596fe9399a5"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[compat]
CommonOPF = "0.3.1"
CommonOPF = "0.3.2"
Graphs = "1.8"
JSON = "0.21"
JuMP = "1"
Expand Down
7 changes: 4 additions & 3 deletions src/model_single_phase.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,7 @@ function constrain_power_balance(m, p::Inputs)
qcon = @constraint(m, [t in 1:p.Ntimesteps],
sum( Qij[string(i*"-"*j), t] for i in i_to_j(j, p) ) -
sum( lij[string(i*"-"*j), t] * xij(i,j,p) for i in i_to_j(j, p) ) +
Qj[j, t] == 0
Qj[j, t] - p.shunt_susceptance[j] * m[:vsqrd][j,t] == 0
)
else
pcon = @constraint(m, [t in 1:p.Ntimesteps],
Expand All @@ -113,8 +113,9 @@ function constrain_power_balance(m, p::Inputs)
)
qcon = @constraint(m, [t in 1:p.Ntimesteps],
sum( Qij[string(i*"-"*j), t] for i in i_to_j(j, p) )
- sum( lij[string(i*"-"*j), t] * xij(i,j,p) for i in i_to_j(j, p) ) +
Qj[j,t] - sum( Qij[string(j*"-"*k), t] for k in j_to_k(j, p) ) == 0
- sum( lij[string(i*"-"*j), t] * xij(i,j,p) for i in i_to_j(j, p) )
+ Qj[j,t] - sum( Qij[string(j*"-"*k), t] for k in j_to_k(j, p) )
- p.shunt_susceptance[j] * m[:vsqrd][j,t]== 0
)
end
m[:loadbalcons][j] = Dict("p" => pcon, "q" => qcon)
Expand Down
7 changes: 6 additions & 1 deletion src/results.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ struct Results <: AbstractResults
current_magnitudes
real_sending_end_powers
reactive_sending_end_powers
shadow_prices
end


Expand All @@ -34,8 +35,12 @@ function Results(m::AbstractModel, p::Inputs{SinglePhase})
lij = get_variable_values(:lij, m, p)
Pij = get_variable_values(:Pij, m, p)
Qij = get_variable_values(:Qij, m, p)
prices = Dict(
j => JuMP.dual.(m[:loadbalcons][j]["p"])
for j in p.busses
)

Results(vs, Pj, Qj, lij, Pij, Qij)
Results(vs, Pj, Qj, lij, Pij, Qij, prices)
end


Expand Down
265 changes: 265 additions & 0 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,271 @@ end

@testset "BranchFlowModel.jl" begin

@testset "Papavasiliou 2018 with shunts" begin
#=
Copied network from paper "Analysis of DLMPs"
and testing some DLMP values here as well as the addition of shunt susceptance values
=#

T = 1
# edges_as_drawn = [
# ("0", "1"), ("1", "2"), ("2", "3"), ("3", "4"), ("4", "5"), ("5", "6"),
# ("8", "7"), ("3", "8"),
# ("8", "9"), ("9", "10"), ("10", "11"),
# ("0", "12"), ("12", "13"), ("13", "14"),
# ]
# edges sorted according to ending bus
# and 7 switched with 8, gives closer voltages
edges = [
("0", "1"), ("1", "2"), ("2", "3"), ("3", "4"), ("4", "5"), ("5", "6"),
("3", "7"), ("7", "8"),
("7", "9"), ("9", "10"), ("10", "11"),
("0", "12"), ("12", "13"), ("13", "14"),
]
linecodes = ["l"*string(i) for i = 1:14]
linelengths = repeat([1.0], length(edges))
phases = repeat([[1]], length(edges))
substation_bus = "0"

Pload = Dict(
"1" => [0.7936],
"2" => [0.0],
"3" => [0.0201],
"4" => [0.0173],
"5" => [0.0291],
"6" => [0.0219],
"7" => [-0.1969], # uncontrolled generator
"8" => [0.0235],
"9" => [0.0229],
"10" => [0.0217],
"11" => [0.0132],
"12" => [0.6219],
"13" => [0.0014],
"14" => [0.0224],
)

Qload = Dict(
"1" => [0.01855],
"2" => [0.0],
"3" => [0.0084],
"4" => [0.0043],
"5" => [0.0073],
"6" => [0.0055],
"7" => [0.0019],
"8" => [0.0059],
"9" => [0.0142],
"10" => [0.0065],
"11" => [0.0033],
"12" => [0.1291],
"13" => [0.0008],
"14" => [0.0083],
)

Zdict = Dict(
"l1" => Dict("rmatrix"=> [0.001], "xmatrix"=> [0.12], "nphases"=> 1),
"l2" => Dict("rmatrix"=> [0.0883], "xmatrix"=> [0.1262], "nphases"=> 1),
"l3" => Dict("rmatrix"=> [0.1384], "xmatrix"=> [0.1978], "nphases"=> 1),
"l4" => Dict("rmatrix"=> [0.0191], "xmatrix"=> [0.0273], "nphases"=> 1),
"l5" => Dict("rmatrix"=> [0.0175], "xmatrix"=> [0.0251], "nphases"=> 1),
"l6" => Dict("rmatrix"=> [0.0482], "xmatrix"=> [0.0689], "nphases"=> 1),
"l7" => Dict("rmatrix"=> [0.0523], "xmatrix"=> [0.0747], "nphases"=> 1),
"l8" => Dict("rmatrix"=> [0.0407], "xmatrix"=> [0.0582], "nphases"=> 1),
"l9" => Dict("rmatrix"=> [0.01], "xmatrix"=> [0.0143], "nphases"=> 1),
"l10" => Dict("rmatrix"=> [0.0241], "xmatrix"=> [0.0345], "nphases"=> 1),
"l11" => Dict("rmatrix"=> [0.0103], "xmatrix"=> [0.0148], "nphases"=> 1),
"l12" => Dict("rmatrix"=> [0.001], "xmatrix"=> [0.12], "nphases"=> 1),
"l13" => Dict("rmatrix"=> [0.1559], "xmatrix"=> [0.1119], "nphases"=> 1),
"l14" => Dict("rmatrix"=> [0.0953], "xmatrix"=> [0.0684], "nphases"=> 1),
)
v0 = 1.0
shunts = Dict(
"0" => 0.0,
"1" => 1.1,
"2" => 2.8,
"3" => 2.4,
"4" => 0.4,
"5" => 0.8,
"6" => 0.6,
"7" => 0.6,
"8" => 1.2,
"9" => 0.4,
"10" => 0.4,
"11" => 0.1,
"12" => 0.1,
"13" => 0.2,
"14" => 0.1,
)
shunts = Dict(k => v*1e-3 for (k,v) in shunts)

p = Inputs(
edges,
linecodes,
linelengths,
phases,
substation_bus;
Pload=Pload,
Qload=Qload,
Sbase=1,
Vbase=1,
Zdict=Zdict,
v0=v0,
v_lolim=0.9,
v_uplim=1.1,
Ntimesteps=T,
P_up_bound=1e4,
Q_up_bound=1e4,
P_lo_bound=-1e4,
Q_lo_bound=-1e4,
Isquared_up_bounds=Dict{String, Float64}(),
relaxed=true, # TODO does the unrelaxed model match unrelaxed
shunt_susceptance=shunts,
);
# TODO check_soc_inequalities
# TODO line limits
m = Model(ECOS.Optimizer)
set_optimizer_attribute(m, "maxit", 10000)
set_optimizer_attribute(m, "verbose", 0)
set_optimizer_attribute(m, "reltol", 1e-7)
build_model!(m,p)

# put in the generator at bus 11
b = "11"
@variable(m, 0.4 >= pgen11 >= 0)
@variable(m, 0.4 >= qgen11 >= 0)
JuMP.delete.(m, m[:injectioncons][b]["p"])
m[:injectioncons][b]["p"] = @constraint(m,
m[:Pj][b, 1] == pgen11 - p.Pload[b][1]
)
JuMP.delete.(m, m[:injectioncons][b]["q"])
m[:injectioncons][b]["q"] = @constraint(m,
m[:Qj][b, 1] == qgen11 - p.Qload[b][1]
)

@objective(m, Min,
m[:Pj][p.substation_bus, 1] * 50 + pgen11 * 10
)
optimize!(m)
r = Results(m,p)

r.real_power_injections["0"]
# 1.064072, looking for 1.063

#=
Increase load on 11 s.t. generator can't export
and test if all prices > 50
=#
m = Model(ECOS.Optimizer)
set_optimizer_attribute(m, "maxit", 100000)
set_optimizer_attribute(m, "verbose", 0)
set_optimizer_attribute(m, "reltol", 1e-6)

p.Pload["11"][1] = 0.4 # this is not enough to get prices > 50 b/c 7 is injecting
p.Pload["7"][1] = 0 # so set load at 7 to zero

build_model!(m,p)

# put in the generator at bus 11
b = "11"
@variable(m, 0.4 >= pgen11 >= 0)
@variable(m, 0.4 >= qgen11 >= 0)
JuMP.delete.(m, m[:injectioncons][b]["p"])
m[:injectioncons][b]["p"] = @constraint(m,
m[:Pj][b, 1] == pgen11 - p.Pload[b][1]
)
JuMP.delete.(m, m[:injectioncons][b]["q"])
m[:injectioncons][b]["q"] = @constraint(m,
m[:Qj][b, 1] == qgen11 - p.Qload[b][1]
)

@objective(m, Min,
m[:Pj][p.substation_bus, 1] * 50 + pgen11 * 10
)
optimize!(m)
r = Results(m,p)
@test all((price[1] >= 49.999 for price in values(r.shadow_prices)))

#=
now with voltage limit preventing full gen the gen should set the price
=#
p.v_uplim = 1.01 # pgen11 not curtailed at 1.03 ?!
p.Pload["11"][1] = 0.0132 # back to original value
p.relaxed = false
m = Model(Ipopt.Optimizer)
set_optimizer_attribute(m, "print_level", 0)

build_model!(m,p)

# put in the generator at bus 11
b = "11"
@variable(m, 0.4 >= pgen11 >= 0)
@variable(m, 0.4 >= qgen11 >= 0)
JuMP.delete.(m, m[:injectioncons][b]["p"])
m[:injectioncons][b]["p"] = @constraint(m,
m[:Pj][b, 1] == pgen11 - p.Pload[b][1]
)
JuMP.delete.(m, m[:injectioncons][b]["q"])
m[:injectioncons][b]["q"] = @constraint(m,
m[:Qj][b, 1] == qgen11 - p.Qload[b][1]
)

@objective(m, Min,
m[:Pj][p.substation_bus, 1] * 50 + pgen11 * 10
)
optimize!(m)
r = Results(m,p)
@test r.shadow_prices["11"][1] ≈ 10.0


#=
is the DLMP lower with DER not net injecting?
i.e. is the total DSO cost lower when paying DLMP > LMP?
=#

p.v_uplim = 1.05 # pgen11 not curtailed at 1.03 ?!
p.Pload["11"][1] = 0.1
p.relaxed = false
m = Model(Ipopt.Optimizer)
set_optimizer_attribute(m, "print_level", 0)

build_model!(m,p)
@objective(m, Min,
m[:Pj][p.substation_bus, 1] * 50
)
optimize!(m)
prices_no_der = Dict(
j => JuMP.dual.(m[:loadbalcons][j]["p"][1])
for j in p.busses
) # TODO put this as option in Results?
cost_no_der = objective_value(m)


m = Model(Ipopt.Optimizer)
set_optimizer_attribute(m, "print_level", 0)
build_model!(m,p)
# put in the generator at bus 11
b = "11"
@variable(m, 0.05 >= pgen11 >= 0)
@variable(m, 0.05 >= qgen11 >= 0)
JuMP.delete.(m, m[:injectioncons][b]["p"])
m[:injectioncons][b]["p"] = @constraint(m,
m[:Pj][b, 1] == pgen11 - p.Pload[b][1]
)
JuMP.delete.(m, m[:injectioncons][b]["q"])
m[:injectioncons][b]["q"] = @constraint(m,
m[:Qj][b, 1] == qgen11 - p.Qload[b][1]
)

@objective(m, Min,
m[:Pj][p.substation_bus, 1] * 50 + pgen11 * 10
)
optimize!(m)
r = Results(m,p)
cost_with_der = value(m[:Pj][p.substation_bus, 1]) * 50 + r.shadow_prices["11"][1] * value(pgen11)
@test cost_with_der < cost_no_der
@test r.shadow_prices["11"][1] < prices_no_der["11"][1]

end


@testset "Results" begin
# p = singlephase38linesInputs(); # TODO use this once CommonOPF updated
Expand Down

2 comments on commit 46dee4e

@NLaws
Copy link
Owner Author

@NLaws NLaws commented on 46dee4e May 14, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/83569

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.3.2 -m "<description of version>" 46dee4e5070162671276891999261854cb580ad9
git push origin v0.3.2

Please sign in to comment.