-
Notifications
You must be signed in to change notification settings - Fork 465
/
run_newton_raphson_pf.py
232 lines (185 loc) · 10.2 KB
/
run_newton_raphson_pf.py
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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
# -*- coding: utf-8 -*-
# Copyright (c) 2016-2022 by University of Kassel and Fraunhofer Institute for Energy Economics
# and Energy System Technology (IEE), Kassel. All rights reserved.
from time import perf_counter
from numpy import flatnonzero as find, r_, zeros, argmax, setdiff1d, union1d, any, int32, sum as np_sum, abs as np_abs
from pandapower.pf.ppci_variables import _get_pf_variables_from_ppci, _store_results_from_pf_in_ppci
from pandapower.pf.run_dc_pf import _run_dc_pf
from pandapower.pypower.bustypes import bustypes
from pandapower.pypower.idx_bus import BUS_I, PD, QD, BUS_TYPE, PQ, GS, BS, SL_FAC as SL_FAC_BUS
from pandapower.pypower.idx_gen import PG, QG, QMAX, QMIN, GEN_BUS, GEN_STATUS, SL_FAC
from pandapower.pypower.makeSbus import makeSbus
from pandapower.pypower.makeYbus import makeYbus as makeYbus_pypower
from pandapower.pypower.newtonpf import newtonpf
from pandapower.pypower.pfsoln import _update_v
from pandapower.pypower.pfsoln import pfsoln as pfsoln_pypower
try:
from pandapower.pf.makeYbus_numba import makeYbus as makeYbus_numba
from pandapower.pf.pfsoln_numba import pfsoln as pfsoln_numba, pf_solution_single_slack
numba_installed = True
except ImportError:
numba_installed = False
try:
from lightsim2grid.newtonpf import newtonpf_new as newton_ls
except ImportError:
newton_ls = None
def _run_newton_raphson_pf(ppci, options):
"""
Runs a Newton-Raphson power flow.
INPUT
ppci (dict) - the "internal" ppc (without out ot service elements and sorted elements)
options(dict) - options for the power flow
"""
t0 = perf_counter()
# we cannot run DC pf before running newton with distributed slack because the slacks come pre-solved after the DC pf
if isinstance(options["init_va_degree"], str) and options["init_va_degree"] == "dc":
if options['distributed_slack']:
pg_copy = ppci['gen'][:, PG].copy()
pd_copy = ppci['bus'][:, PD].copy()
ppci = _run_dc_pf(ppci)
ppci['gen'][:, PG] = pg_copy
ppci['bus'][:, PD] = pd_copy
else:
ppci = _run_dc_pf(ppci)
if options["enforce_q_lims"]:
ppci, success, iterations, bus, gen, branch = _run_ac_pf_with_qlims_enforced(ppci, options)
else:
ppci, success, iterations = _run_ac_pf_without_qlims_enforced(ppci, options)
# update data matrices with solution store in ppci
bus, gen, branch = ppci_to_pfsoln(ppci, options)
et = perf_counter() - t0
ppci = _store_results_from_pf_in_ppci(ppci, bus, gen, branch, success, iterations, et)
return ppci
def ppci_to_pfsoln(ppci, options):
internal = ppci["internal"]
if options["only_v_results"]:
# time series relevant hack which ONLY saves V from ppci
_update_v(internal["bus"], internal["V"])
return internal["bus"], internal["gen"], internal["branch"]
else:
# reads values from internal ppci storage to bus, gen, branch and returns it
if options['distributed_slack']:
# consider buses with non-zero slack weights as if they were slack buses,
# and gens with non-zero slack weights as if they were reference machines
# this way, the function pfsoln will extract results for distributed slack gens, too
# also, the function pfsoln will extract results for the PQ buses for xwards
gens_with_slack_weights = find(internal["gen"][:, SL_FAC] != 0)
# gen_buses_with_slack_weights = internal["gen"][gens_with_slack_weights, GEN_BUS].astype(int32)
buses_with_slack_weights = internal["bus"][find(internal["bus"][:, SL_FAC_BUS] != 0), BUS_I].astype(int32)
# buses_with_slack_weights = union1d(gen_buses_with_slack_weights, buses_with_slack_weights)
ref = union1d(internal["ref"], buses_with_slack_weights)
ref_gens = union1d(internal["ref_gens"], gens_with_slack_weights)
else:
ref = internal["ref"]
ref_gens = internal["ref_gens"]
_, pfsoln = _get_numba_functions(ppci, options)
result_pfsoln = pfsoln(internal["baseMVA"], internal["bus"], internal["gen"], internal["branch"], internal["Ybus"],
internal["Yf"], internal["Yt"], internal["V"], ref, ref_gens)
return result_pfsoln
def _get_Y_bus(ppci, options, makeYbus, baseMVA, bus, branch):
recycle = options["recycle"]
if isinstance(recycle, dict) and not recycle["trafo"] and ppci["internal"]["Ybus"].size:
Ybus, Yf, Yt = ppci["internal"]['Ybus'], ppci["internal"]['Yf'], ppci["internal"]['Yt']
else:
# build admittance matrices
Ybus, Yf, Yt = makeYbus(baseMVA, bus, branch)
return ppci, Ybus, Yf, Yt
def _get_numba_functions(ppci, options):
"""
pfsoln from pypower maybe slow in some cases. This function chooses the fastest for the given pf calculation
"""
if options["numba"] and numba_installed:
makeYbus = makeYbus_numba
shunt_in_net = any(ppci["bus"][:, BS]) or any(ppci["bus"][:, GS])
# faster pfsoln function if only one slack is in the grid and no gens
pfsoln = pf_solution_single_slack if ppci["gen"].shape[0] == 1 \
and not options["voltage_depend_loads"] \
and not options['distributed_slack'] \
and not shunt_in_net \
else pfsoln_numba
else:
makeYbus = makeYbus_pypower
pfsoln = pfsoln_pypower
return makeYbus, pfsoln
def _store_internal(ppci, internal_storage):
# internal storage is a dict with the variables to store in net["_ppc"]["internal"]
for key, val in internal_storage.items():
ppci["internal"][key] = val
return ppci
def _get_Sbus(ppci, recycle=None):
baseMVA, bus, gen = ppci["baseMVA"], ppci["bus"], ppci["gen"]
if not isinstance(recycle, dict) or "Sbus" not in ppci["internal"]:
return makeSbus(baseMVA, bus, gen)
if recycle["bus_pq"] or recycle["gen"]:
return makeSbus(baseMVA, bus, gen)
return ppci["internal"]["Sbus"]
def _run_ac_pf_without_qlims_enforced(ppci, options):
makeYbus, pfsoln = _get_numba_functions(ppci, options)
baseMVA, bus, gen, branch, ref, pv, pq, _, _, V0, ref_gens = _get_pf_variables_from_ppci(ppci)
ppci, Ybus, Yf, Yt = _get_Y_bus(ppci, options, makeYbus, baseMVA, bus, branch)
# compute complex bus power injections [generation - load]
Sbus = _get_Sbus(ppci, options["recycle"])
# run the newton power flow
if options["lightsim2grid"]:
V, success, iterations, J, Vm_it, Va_it = newton_ls(Ybus.tocsc(), Sbus, V0, ref, pv, pq, ppci, options)
else:
V, success, iterations, J, Vm_it, Va_it = newtonpf(Ybus, Sbus, V0, ref, pv, pq, ppci, options)
# keep "internal" variables in memory / net["_ppc"]["internal"] -> needed for recycle.
ppci = _store_internal(ppci, {"J": J, "Vm_it": Vm_it, "Va_it": Va_it, "bus": bus, "gen": gen, "branch": branch,
"baseMVA": baseMVA, "V": V, "pv": pv, "pq": pq, "ref": ref, "Sbus": Sbus,
"ref_gens": ref_gens, "Ybus": Ybus, "Yf": Yf, "Yt": Yt})
return ppci, success, iterations
def _run_ac_pf_with_qlims_enforced(ppci, options):
baseMVA, bus, gen, branch, ref, pv, pq, on, _, V0, ref_gens = _get_pf_variables_from_ppci(ppci)
qlim = options["enforce_q_lims"]
limited = [] # list of indices of gens @ Q lims
fixedQg = zeros(gen.shape[0]) # Qg of gens at Q limits
while True:
ppci, success, iterations = _run_ac_pf_without_qlims_enforced(ppci, options)
bus, gen, branch = ppci_to_pfsoln(ppci, options)
# find gens with violated Q constraints
gen_status = gen[:, GEN_STATUS] > 0
qg_max_lim = gen[:, QG] > gen[:, QMAX]
qg_min_lim = gen[:, QG] < gen[:, QMIN]
mx = setdiff1d(find(gen_status & qg_max_lim), ref_gens)
mn = setdiff1d(find(gen_status & qg_min_lim), ref_gens)
if len(mx) > 0 or len(mn) > 0: # we have some Q limit violations
# one at a time?
if qlim == 2: # fix largest violation, ignore the rest
k = argmax(r_[gen[mx, QG] - gen[mx, QMAX],
gen[mn, QMIN] - gen[mn, QG]])
if k > len(mx):
mn = mn[k - len(mx)]
mx = []
else:
mx = mx[k]
mn = []
# save corresponding limit values
fixedQg[mx] = gen[mx, QMAX]
fixedQg[mn] = gen[mn, QMIN]
mx = r_[mx, mn].astype(int)
# convert to PQ bus
gen[mx, QG] = fixedQg[mx] # set Qg to binding
for i in range(len(mx)): # [one at a time, since they may be at same bus]
gen[mx[i], GEN_STATUS] = 0 # temporarily turn off gen,
bi = gen[mx[i], GEN_BUS].astype(int) # adjust load accordingly,
bus[bi, [PD, QD]] = (bus[bi, [PD, QD]] - gen[mx[i], [PG, QG]])
# if len(ref) > 1 and any(bus[gen[mx, GEN_BUS].astype(int), BUS_TYPE] == REF):
# raise ValueError('Sorry, pandapower cannot enforce Q '
# 'limits for slack buses in systems '
# 'with multiple slacks.')
changed_gens = gen[mx, GEN_BUS].astype(int)
bus[setdiff1d(changed_gens, ref), BUS_TYPE] = PQ # & set bus type to PQ
# update bus index lists of each type of bus
ref, pv, pq = bustypes(bus, gen)
limited = r_[limited, mx].astype(int)
else:
break # no more generator Q limits violated
if len(limited) > 0:
# restore injections from limited gens [those at Q limits]
gen[limited, QG] = fixedQg[limited] # restore Qg value,
for i in range(len(limited)): # [one at a time, since they may be at same bus]
bi = gen[limited[i], GEN_BUS].astype(int) # re-adjust load,
bus[bi, [PD, QD]] = bus[bi, [PD, QD]] + gen[limited[i], [PG, QG]]
gen[limited[i], GEN_STATUS] = 1 # and turn gen back on
return ppci, success, iterations, bus, gen, branch