Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
using CSV
using Dates
using DataFrames
using Statistics
using LinearAlgebra
using Graphs
# path to data folder
network = "elec_s_47_ec_lcopt_Co2L-2H"
# construct network
loads = CSV.read("$(@__DIR__)/$network/loads.csv", DataFrame)
lines = CSV.read("$(@__DIR__)/$network/lines.csv", DataFrame)
g = Graph()
add_vertices!(g, length(nodes))
srcs = indexin(lines.bus0, loads.bus)
dsts = indexin(lines.bus1, loads.bus)
for (u,v) in zip(srcs, dsts)
add_edge!(g, u,v)
end
## Time series
# time stamps
snapshots = CSV.read("$(@__DIR__)/$network/snapshots.csv", DataFrame)
snapshots = DateTime.(snapshots.snapshot, "yyyy-mm-dd HH:MM:SS")
# consumption
loads_ts = CSV.read("$(@__DIR__)/$network/loads-p_set.csv", DataFrame)
loads_ts[!, 1] = snapshots
rename!(loads_ts, :Column1 => :snapshot)
# production
generators = CSV.read("$(@__DIR__)/$network/generators.csv", DataFrame)
generators_ts = CSV.read("$(@__DIR__)/$network/generators-p_max_pu.csv", DataFrame)
generators_ts[!, 1] = snapshots
rename!(generators_ts, :Column1 => :snapshot)
## Renewables
wind = generators[occursin.("wind", generators.carrier), :]
wind_ts = generators_ts[!, ["snapshot", wind.name...]]
solar = generators[occursin.("solar", generators.carrier), :]
solar_ts = generators_ts[!, ["snapshot", solar.name...]]
# aggregate over nodes
renewables = vcat(wind, solar)
renewables_ts = generators_ts[!, ["snapshot", wind.name..., solar.name...]]
for b in unique(renewables.bus)
idx = renewables[ renewables.bus .== b, : ].name
agg = sum.(eachrow(renewables_ts[!, idx]))
for c in idx
select!(renewables_ts, Not(c))
end
renewables_ts[!, b] = agg
end
## Calculate averages
avg_load = mean( sum.(eachrow(loads_ts[!, 2:end])) )
avg_prod = mean( sum.(eachrow(generators_ts[!, 2:end])) )
avg_wind = mean( sum.(eachrow(wind_ts[!, 2:end])) )
avg_solar = mean( sum.(eachrow(solar_ts[!, 2:end])) )
avg_renewables = mean( sum.(eachrow(renewables_ts[!, 2:end])) )
println("Penetration")
println("\tWind\t", avg_wind / avg_load * 100, "%")
println("\tSolar\t", avg_solar / avg_load * 100, "%")
println("\tRenewables\t", avg_renewables / avg_load * 100, "%")
println("\tAll sources\t", avg_prod / avg_load * 100, "%")
# [Q]: why is this number so low? how much is going into storage?
## Calculate power injection and transmission distribution factor
# mismatch
mismatch = renewables_ts[!,2:end] .- loads_ts[!,2:end]
Δ⁺ = max.(0, mismatch)
Δ⁻ = max.(0, (-1) .* mismatch)
# backup energy
b = Δ⁺ .- (mean.(eachcol(loads_ts[!,2:end])) / avg_load)' .* Δ⁻
# [Q]: what about first term? here sum(Δ⁺)/sum(Δ⁺) was left out since this would be 0/0
# power injection
p = hcat(DataFrame(:snapshot => snapshots), mismatch .- b)
# power transmission distribution factor
Ω = diagm(1 ./ lines.x)
K = incidence_matrix(g, oriented=true)
B = diagm( ones(nv(g)) ) #???
H = Ω * K' * pinv(B)
# [Q]: what is the nodal succeptance matrix B?