Skip to content

Commit 2a80737

Browse files
prezaei85prezaei85
committed
update to check the kernel restart.
1 parent b5a2380 commit 2a80737

12 files changed

+674
-242
lines changed

‎ACPowerFlow/src/FindSubGraphs.jl

+14-6
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,8 @@
1-
function FindSubsGraphs(nodes,links)
1+
function FindSubGraphs(nodes,links)
22
# find connected components in a graph
33
nodes = vec(nodes) # make sure nodes are in vector format and not matrix
44
N = length(nodes)
5+
nbr = size(links,1)
56
bus_i = zeros(Int64,maximum(nodes))
67
bus_i[nodes] = 1:N
78
F = bus_i[links[:,1]]
@@ -10,8 +11,9 @@ function FindSubsGraphs(nodes,links)
1011

1112
grNo = 1
1213
graphNos = zeros(Int64,N)
14+
linkNos = zeros(Int64,nbr)
1315
next = 1
14-
while !isempty(next)
16+
while next != 0
1517
included = falses(N)
1618
included[next] = true
1719
oldLen = 0
@@ -22,8 +24,14 @@ function FindSubsGraphs(nodes,links)
2224
end
2325
graphNos[included] = grNo
2426
grNo = grNo+1
25-
next = find(graphNos.==0)
27+
next = findfirst(graphNos.==0,1)
2628
end
27-
nSubGraphs = grNo-1
28-
return graphNos
29-
end
29+
# find where each link is
30+
for i = 1:nbr
31+
this_busi = F[i]
32+
linkNos[i] = graphNos[this_busi]
33+
end
34+
35+
return graphNos, linkNos
36+
end
37+

‎ACPowerFlow/src/acpf.jl

+83-58
Original file line numberDiff line numberDiff line change
@@ -9,30 +9,40 @@ function acpf!(ps::Caseps;
99
nodes = ps.bus[:,:ID]
1010
br_status = ps.branch[:,:status]
1111
links = ps.branch[br_status, [:from, :to]]
12-
graphNos = FindSubsGraphs(nodes,links)
12+
graphNos, linkNos = FindSubGraphs(nodes,links)
1313
Vmag = ps.bus[:Vmag] # note that Vmag and theta are passed by reference, so they are in fact the same as ps.bus[:Vmag] and ps.bus[:Vang]
14-
ps.bus[:Vang] = 0.0
15-
theta = convert(Array,ps.bus[:Vang])
14+
ps.bus[:,:Vang] = 0.0
15+
theta = ps.bus[:Vang]
1616
n_subgraph = maximum(graphNos)
17-
one_converge_flag = false
17+
Ybus, Yf, Yt = getYbus(ps)
18+
Sbus, Sd, Sg = getSbus(ps)
19+
Sd = sum(Sd,2)
20+
D = ps.bus_i[ps.shunt[:,:bus]]
21+
Sd_bus = zeros(Complex{Float64},nBus); Sd_bus[D] = Sd
22+
F = ps.bus_i[ps.branch[:,:from]]
23+
T = ps.bus_i[ps.branch[:,:to]]
24+
not_converged = Int64[] # list of the islands that did not converge
1825
for i = 1:n_subgraph
26+
# find buses in this island
1927
busi_sub = find(graphNos .== i)
28+
# find branches in this island
29+
bri_sub = find(linkNos .== i)
2030
mismatch, x0, ix_Vmag_sub, ix_theta_sub, pq_sub, ref_sub =
2131
init_mismatch(ps, busi_sub, verbose, PartFactFlag, PartFact, PolarFlag)
2232
## Solve power flow
2333
# newton-raphson
24-
x, converged = nrsolve(mismatch, x0, verbose=verbose, linesearch="exact");
34+
x, converged = nrsolve(mismatch, x0, verbose=verbose, linesearch="exact")
2535

2636
if converged
27-
one_converge_flag = true # just check if at least one converged
2837
# save the results
2938
if PolarFlag
30-
idx = [busi_sub][pq_sub]
31-
Vmag[idx] = x[ix_Vmag_sub]
32-
idx = [busi_sub][!ref_sub]
33-
theta[idx] = x[ix_theta_sub]
34-
idx = [busi_sub][ref_sub]
35-
theta[idx] = 0
39+
Vmag_sub = Vmag[busi_sub]
40+
Vmag_sub[pq_sub] = x[ix_Vmag_sub]
41+
Vmag[busi_sub] = Vmag_sub # this updates ps.bus[:,:Vmag]
42+
theta_sub = theta[busi_sub]
43+
theta_sub[!ref_sub] = x[ix_theta_sub]
44+
theta_sub[ref_sub] = 0.
45+
theta[busi_sub] = theta_sub # this updates ps.bus[:,:Vang]
3646
else
3747
# ignore this for now
3848
#=
@@ -44,57 +54,72 @@ function acpf!(ps::Caseps;
4454
theta = atan(f./e)
4555
=#
4656
end
47-
elseif verbose
48-
println("Power flow did not converge on island $i of $n_subgraph.")
49-
end
50-
end
51-
if one_converge_flag
52-
# make sure that only the converged islands have updated data!
53-
54-
# voltage results
55-
V = Vmag.*exp(im*theta)
56-
D = ps.bus_i[ps.shunt[:,:bus]]
57-
Ybus, Yf, Yt = getYbus(ps)
58-
Sbus, Sd, Sg = getSbus(ps)
59-
Sd = sum(Sd,2)
60-
Sd_bus = zeros(Complex{Float64},nBus); Sd_bus[D] = Sd
61-
# branch results
62-
If = Yf*V # branch status is accounted for in Yf
63-
It = Yt*V # branch status is accounted for in Yt
64-
F = ps.bus_i[ps.branch[:,:from]]
65-
T = ps.bus_i[ps.branch[:,:to]]
66-
Sf = V[F] .* conj(If)
67-
St = V[T] .* conj(It)
68-
ps.branch[:,:ImagF] = abs(If)
69-
ps.branch[:,:ImagT] = abs(It)
70-
ps.branch[:,:Pf] = real(Sf) * ps.baseMVA
71-
ps.branch[:,:Qf] = imag(Sf) * ps.baseMVA
72-
ps.branch[:,:Pt] = real(St) * ps.baseMVA
73-
ps.branch[:,:Qt] = imag(St) * ps.baseMVA
74-
# update ps.gen
75-
Sg_bus = V.*conj(Ybus*V) + Sd_bus
76-
Pg_bus = real(Sg_bus)
77-
Qg_bus = imag(Sg_bus)
78-
for gi = 1:size(ps.gen,1)
79-
gen_bus = ps.gen[gi,1]
80-
gen_bus_i = ps.bus_i[gen_bus]
81-
if Vmag[gen_bus_i] == 0
82-
ps.gen[gi,:P] = 0
83-
ps.gen[gi,:Q] = 0
84-
else
85-
gens_at_bus = ps.gen[:,1] .== gen_bus
86-
if sum(gens_at_bus) == 1
87-
Pg_ratio = 1
57+
# update results in ps only for this island
58+
V_sub = Vmag_sub.*exp(im*theta_sub)
59+
V = zeros(Complex{Float64},nBus)
60+
V[busi_sub] = V_sub
61+
V_sub = Vmag_sub.*exp(im*theta_sub)
62+
Yf_sub = Yf[bri_sub,busi_sub]
63+
Yt_sub = Yt[bri_sub,busi_sub]
64+
If_sub = Yf_sub*V_sub
65+
It_sub = Yt_sub*V_sub
66+
F_sub = F[bri_sub]
67+
Sf_sub = V[F_sub] .* conj(If_sub)
68+
T_sub = T[bri_sub]
69+
St_sub = V[T_sub] .* conj(It_sub)
70+
ps.branch[bri_sub,:ImagF] = abs(If_sub)
71+
ps.branch[bri_sub,:ImagT] = abs(It_sub)
72+
ps.branch[bri_sub,:Pf] = real(Sf_sub) * ps.baseMVA
73+
ps.branch[bri_sub,:Qf] = imag(Sf_sub) * ps.baseMVA
74+
ps.branch[bri_sub,:Pt] = real(St_sub) * ps.baseMVA
75+
ps.branch[bri_sub,:Qt] = imag(St_sub) * ps.baseMVA
76+
# update ps.gen
77+
Ybus_sub = Ybus[busi_sub,busi_sub]
78+
Sg_bus_sub = V_sub.*conj(Ybus_sub*V_sub) + Sd_bus[busi_sub]
79+
Pg_bus_sub = real(Sg_bus_sub)
80+
Qg_bus_sub = imag(Sg_bus_sub)
81+
gen_busID = ps.gen[:,:bus]
82+
gen_busi = ps.bus_i[gen_busID]
83+
for (j,this_bus_i) = enumerate(busi_sub)
84+
gi = (gen_busi .== this_bus_i)
85+
if V[this_bus_i] == 0.
86+
ps.gen[gi,:P] = 0.
87+
ps.gen[gi,:Q] = 0.
8888
else
89-
Pg_ratio = ps.gen[gi,:Pmax]/sum(ps.gen[gens_at_bus,:Pmax])
89+
if sum(gi) == 1
90+
Pg_ratio = 1
91+
else
92+
Pg_ratio = ps.gen[gi,:Pmax]/sum(ps.gen[gi,:Pmax])
93+
end
94+
ps.gen[gi,:P] = Pg_bus_sub[j] * Pg_ratio * ps.baseMVA
95+
ps.gen[gi,:Q] = Qg_bus_sub[j] * Pg_ratio * ps.baseMVA
9096
end
91-
ps.gen[gi,:P] = Pg_bus[gen_bus_i] * Pg_ratio * ps.baseMVA
92-
ps.gen[gi,:Q] = Qg_bus[gen_bus_i] * Pg_ratio * ps.baseMVA
9397
end
94-
end
98+
else
99+
if verbose
100+
println("Power flow did not converge on island $i of $n_subgraph.")
101+
end
102+
not_converged = [not_converged, i]
103+
# update ps with 0 values for non-converged islands
104+
Vmag[busi_sub] = 0.
105+
theta[busi_sub] = 0.
106+
ps.branch[bri_sub,:ImagF] = 0.
107+
ps.branch[bri_sub,:ImagT] = 0.
108+
ps.branch[bri_sub,:Pf] = 0.
109+
ps.branch[bri_sub,:Qf] = 0.
110+
ps.branch[bri_sub,:Pt] = 0.
111+
ps.branch[bri_sub,:Qt] = 0.
112+
gen_busID = ps.gen[:,:bus]
113+
gen_busi = ps.bus_i[gen_busID]
114+
for j = busi_sub
115+
gi = (gen_busi .== j)
116+
ps.gen[gi,:P] = 0.
117+
ps.gen[gi,:Q] = 0.
118+
end
119+
end
95120
end
96121

97-
return
122+
return not_converged
98123

99124
#=
100125
# solve the unconstrained optimization problem 1/2 * g(x)' * g(x) using levenberg-marquardt

‎ACPowerFlow/src/mismatch_polar.jl

+4-3
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,4 @@
1-
using Debug
2-
@debug function mismatch_polar(x::Vector{Float64},
1+
function mismatch_polar(x::Vector{Float64},
32
Ybus::SparseMatrixCSC{Complex{Float64},Int64},
43
Vmag::Vector{Float64},
54
Sg_bus::Vector{Complex{Float64}},
@@ -156,4 +155,6 @@ using Debug
156155
end
157156

158157
return g, dg_dx, S
159-
end
158+
end
159+
160+

‎ACPowerFlow/src/nrsolve.jl

+1
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,7 @@ function nrsolve(mismatch::Function,
3131
Jg = Jac.'*g;
3232
max_Jg = maximum(Jg)
3333
g2 = sqrt(sum(g.^2))
34+
# @bp
3435

3536
# choose the search direction
3637
p = -(Jac\g)

‎ACSIMSEP/src/ACSIMSEP.jl

+9-53
Original file line numberDiff line numberDiff line change
@@ -3,8 +3,12 @@ module ACSIMSEP
33
using Makeps, ACPowerFlow
44

55
export get_left_eigv,
6-
find_Sm,
7-
find_sens_fact
6+
find_sens_fact,
7+
find_beta
8+
9+
# find the sensitivity factors to laod shedding
10+
include("find_sens_fact.jl")
11+
# include("find_sens_fact_old.jl")
812

913
function get_left_eigv(J::SparseMatrixCSC)
1014
# gives the left eigen vector of the zero eigen value for the jacobian and checks for some errors
@@ -35,13 +39,13 @@ function get_left_eigv(J::SparseMatrixCSC)
3539
return w
3640
end
3741

38-
function find_Sm(mismatch::Function, x0::Vector{Float64}; verbose::Bool = true)
39-
# This function finds the minimum distance in parameter space from
42+
function find_beta(mismatch::Function, x0::Vector{Float64}; verbose::Bool = true)
43+
# This function finds the minimum distance in parameter space (beta) from
4044
# an unsolvable point to the solvable region (based on Overbye's paper)
4145
x, converged = nrsolve(mismatch, x0, verbose=verbose, linesearch="exact")
4246
if converged == true
4347
dist = 0
44-
println("The initial power flow converged.")
48+
println("Warning: The initial power flow converged.")
4549
end
4650
g(x) = mismatch(x)[1]
4751
Jac(x) = mismatch(x, NeedJac=true)[2]
@@ -73,54 +77,6 @@ function find_Sm(mismatch::Function, x0::Vector{Float64}; verbose::Bool = true)
7377
return x, Jac(x), beta
7478
end
7579

76-
function find_sens_fact(ps::Caseps, PartFactFlag; verbose = true)
77-
mismatch, x0, Vmag, pq, ref, Yf, Yt, ix_Vmag, ix_theta = init_mismatch(ps, true, PartFactFlag, true);
78-
# find the minimum distance to solvability
79-
xm, Jm, beta = find_Sm(mismatch,x0,verbose = verbose);
80-
# get the left eigen vector of Jacobian associated with zero eigen value
81-
w_m = get_left_eigv(Jm)
82-
# resolve the direction of w_m (it needs to be outward the solvable region)
83-
g(x) = mismatch(x)[1]
84-
if dot(g(xm),w_m) < 0
85-
w_m = -w_m
86-
end
87-
nbus = size(ps.bus,1)
88-
Sbus, = getSbus(ps)
89-
if PartFactFlag
90-
# beta_P: the sensitivity factors when reducing real power at buses
91-
beta_P = w_m[1:nbus]
92-
93-
# beta_Q: the sensitivity factors when reducing reactive power at buses
94-
beta_Q = zeros(nbus)
95-
beta_Q[pq] = w_m[nbus+1:end]
96-
else
97-
# beta_P: the sensitivity factors when reducing real power at buses
98-
beta_P = zeros(nbus)
99-
beta_P[!ref] = w_m[1:nbus-1]
100-
101-
# beta_Q: the sensitivity factors when reducing reactive power at buses
102-
beta_Q = zeros(nbus)
103-
beta_Q[pq] = w_m[nbus:end]
104-
end
105-
# beta_S: the sensitivity factors when reducing complex power at buses while maintaining power factor
106-
# only considers buses with positive load
107-
Pbus = real(Sbus)
108-
Qbus = imag(Sbus)
109-
LoadBusID = (Pbus.<=0) & (Qbus.<=0) & !((Pbus.==0) & (Qbus.==0)) # finds buses with at least one non-zero P/Q
110-
LoadBusID = find(LoadBusID)
111-
power_factor = zeros(nbus)
112-
power_factor[LoadBusID] = -Pbus[LoadBusID]./sqrt(Pbus[LoadBusID].^2 + Qbus[LoadBusID].^2)
113-
beta_S = zeros(nbus)
114-
beta_S[LoadBusID] = beta_P[LoadBusID].*power_factor[LoadBusID] + beta_Q[LoadBusID].*sqrt(1-power_factor[LoadBusID].^2)
115-
116-
# beta_C: the sensitivity factors when switching a capacitor at buses assuming there was no capacitor before
117-
beta_C = zeros(nbus)
118-
Vmag[pq] = xm[ix_Vmag]
119-
beta_C[pq] = Vmag[pq].^2
120-
beta_C = beta_C.*beta_Q
121-
122-
return beta, beta_P, beta_Q, beta_S, beta_C
123-
end
12480

12581
end
12682

‎ACSIMSEP/src/find_sens_fact.jl

+75
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,75 @@
1+
function find_sens_fact(ps::Caseps, PartFactFlag, not_converged; verbose = true)
2+
nBus = size(ps.bus,1)
3+
nodes = ps.bus[:,:ID]
4+
br_status = ps.branch[:,:status]
5+
links = ps.branch[br_status, [:from, :to]]
6+
graphNos, linkNos = FindSubGraphs(nodes,links)
7+
n_subgraph = maximum(graphNos)
8+
beta_P = zeros(nBus); beta_Q = zeros(nBus)
9+
beta_S = zeros(nBus); beta_C = zeros(nBus)
10+
Sbus, = getSbus(ps)
11+
Vmag = ps.bus[:,:Vmag]
12+
C = Constps()
13+
for i = not_converged
14+
# find buses in this island
15+
busi_sub = find(graphNos .== i)
16+
nBus_sub = length(busi_sub)
17+
beta_P_sub = zeros(nBus_sub); beta_Q_sub = zeros(nBus_sub)
18+
beta_S_sub = zeros(nBus_sub); beta_C_sub = zeros(nBus_sub)
19+
# find branches in this island
20+
bri_sub = find(linkNos .== i)
21+
mismatch, x0, ix_Vmag_sub, ix_theta_sub, pq_sub, ref_sub =
22+
init_mismatch(ps, busi_sub, verbose, PartFactFlag, Float64[], true)
23+
# find the minimum distance to solvability
24+
xm, Jm, beta = find_beta(mismatch,x0,verbose = verbose);
25+
26+
# get the left eigen vector of Jacobian associated with zero eigen value
27+
w_m = get_left_eigv(Jm)
28+
# resolve the direction of w_m (it needs to be outward the solvable region)
29+
g(x) = mismatch(x)[1]
30+
if dot(g(xm),w_m) < 0 # resolve the direction of the eigen vector
31+
w_m = -w_m
32+
end
33+
34+
if PartFactFlag
35+
# beta_P: the sensitivity factors when reducing real power at buses
36+
beta_P_sub = w_m[1:nBus_sub]
37+
38+
# beta_Q: the sensitivity factors when reducing reactive power at buses
39+
beta_Q_sub[pq_sub] = w_m[nBus_sub+1:end]
40+
else
41+
# beta_P: the sensitivity factors when reducing real power at buses
42+
beta_P_sub[!ref_sub] = w_m[1:nBus_sub-1]
43+
44+
# beta_Q: the sensitivity factors when reducing reactive power at buses
45+
beta_Q_sub[pq_sub] = w_m[nBus_sub:end]
46+
end
47+
# beta_S: the sensitivity factors when reducing complex power at buses while maintaining power factor
48+
# only considers buses with positive load
49+
Sbus_sub = Sbus[busi_sub]
50+
Pbus_sub = real(Sbus_sub)
51+
Qbus_sub = imag(Sbus_sub)
52+
load_busi_sub = (Pbus_sub.<=0) & (Qbus_sub.<=0) & !((Pbus_sub.==0) & (Qbus_sub.==0)) # finds buses with at least one non-zero P/Q
53+
load_busi_sub = find(load_busi_sub)
54+
power_factor = zeros(nBus_sub)
55+
power_factor[load_busi_sub] = -Pbus_sub[load_busi_sub]./sqrt(Pbus_sub[load_busi_sub].^2 + Qbus_sub[load_busi_sub].^2)
56+
beta_S_sub[load_busi_sub] = beta_P_sub[load_busi_sub].*power_factor[load_busi_sub] +
57+
beta_Q_sub[load_busi_sub].*sqrt(1-power_factor[load_busi_sub].^2)
58+
59+
# beta_C: the sensitivity factors when switching a capacitor at buses assuming there was no capacitor before
60+
Vmag_sub = Vmag[busi_sub]
61+
Vmag_sub[pq_sub] = xm[ix_Vmag_sub]
62+
beta_C_sub[pq_sub] = Vmag_sub[pq_sub].^2
63+
beta_C_sub = beta_C_sub.*beta_Q_sub
64+
65+
# update the original list
66+
beta_P[busi_sub] = beta_P_sub
67+
beta_Q[busi_sub] = beta_Q_sub
68+
beta_S[busi_sub] = beta_S_sub
69+
beta_C[busi_sub] = beta_C_sub
70+
end
71+
72+
return beta, beta_P, beta_Q, beta_S, beta_C
73+
end
74+
75+

0 commit comments

Comments
 (0)