Add auto_open parameter for plotly plotting functions
1 
# * coding: utf8 *


2  
3 
# Copyright (c) 20162021 by University of Kassel and Fraunhofer Institute for Energy Economics


4 
# and Energy System Technology (IEE), Kassel. All rights reserved.


5  
6  1 
import numpy as np 
7  1 
from scipy.stats import chi2 
8  
9  1 
from pandapower.estimation.algorithm.base import (WLSAlgorithm, 
10 
WLSZeroInjectionConstraintsAlgorithm, 

11 
IRWLSAlgorithm) 

12  1 
from pandapower.estimation.algorithm.lp import LPAlgorithm 
13  1 
from pandapower.estimation.algorithm.optimization import OptAlgorithm 
14  1 
from pandapower.estimation.ppc_conversion import pp2eppci, _initialize_voltage 
15  1 
from pandapower.estimation.results import eppci2pp 
16  1 
from pandapower.estimation.util import set_bb_switch_impedance, reset_bb_switch_impedance 
17  
18  1 
try: 
19  1 
import pplog as logging 
20  1 
except ImportError: 
21  1 
import logging 
22  1 
std_logger = logging.getLogger(__name__) 
23  
24  1 
ALGORITHM_MAPPING = {'wls': WLSAlgorithm, 
25 
'wls_with_zero_constraint': WLSZeroInjectionConstraintsAlgorithm, 

26 
'opt': OptAlgorithm, 

27 
'irwls': IRWLSAlgorithm, 

28 
'lp': LPAlgorithm} 

29  1 
ALLOWED_OPT_VAR = {"a", "opt_method", "estimator"} 
30  
31  
32  1 
def estimate(net, algorithm='wls', 
33 
init='flat', tolerance=1e6, maximum_iterations=10, 

34 
calculate_voltage_angles=True, 

35 
zero_injection='aux_bus', fuse_buses_with_bb_switch='all', 

36 
**opt_vars): 

37 
"""


38 
Wrapper function for WLS state estimation.


39 
INPUT:


40 
**net**  The net within this line should be created


41  
42 
**init**  (string) Initial voltage for the estimation. 'flat' sets 1.0 p.u. / 0° for all


43 
buses, 'results' uses the values from *res_bus* if available and 'slack' considers the


44 
slack bus voltage (and optionally, angle) as the initial values. Default is 'flat'


45 
OPTIONAL:


46 
**tolerance**  (float)  When the maximum state change between iterations is less than


47 
tolerance, the process stops. Default is 1e6


48  
49 
**maximum_iterations**  (integer)  Maximum number of iterations. Default is 10


50  
51 
**calculate_voltage_angles**  (boolean)  Take into account absolute voltage angles and phase


52 
shifts in transformers, if init is 'slack'. Default is True


53  
54 
**zero_injection**  (str, iterable, None)  Defines which buses are zero injection bus or the method


55 
to identify zero injection bus, with 'wls_estimator' virtual measurements will be added, with


56 
'wls_estimator with zero constraints' the buses will be handled as constraints


57 
"auto": all bus without p,q measurement, without p, q value (load, sgen...) and aux buses will be


58 
identified as zero injection bus


59 
"aux_bus": only aux bus will be identified as zero injection bus


60 
None: no bus will be identified as zero injection bus


61 
iterable: the iterable should contain index of the zero injection bus and also aux bus will be identified


62 
as zeroinjection bus


63  
64 
**fuse_buses_with_bb_switch**  (str, iterable, None)  Defines how buses with closed bb switches should


65 
be handled, if fuse buses will only fused to one for calculation, if not fuse, an auxiliary bus and


66 
auxiliary line will be automatically added to the network to make the buses with different p,q injection


67 
measurements identifieble


68 
"all": all buses with bbswitches will be fused, the same as the default behaviour in load flow


69 
None: buses with bbswitches and individual p,q measurements will be reconfigurated


70 
by auxiliary elements


71 
iterable: the iterable should contain the index of buses to be fused, the behaviour is contigous e.g.


72 
if one of the bus among the buses connected through bb switch is given, then all of them will still


73 
be fused


74 
OUTPUT:


75 
**successful** (boolean)  Was the state estimation successful?


76 
"""


77  1 
if algorithm not in ALGORITHM_MAPPING: 
78  1 
raise UserWarning("Algorithm {} is not a valid estimator".format(algorithm)) 
79  
80  1 
se = StateEstimation(net, tolerance, maximum_iterations, algorithm=algorithm) 
81  1 
v_start, delta_start = _initialize_voltage(net, init, calculate_voltage_angles) 
82  1 
return se.estimate(v_start=v_start, delta_start=delta_start, 
83 
calculate_voltage_angles=calculate_voltage_angles, 

84 
zero_injection=zero_injection, 

85 
fuse_buses_with_bb_switch=fuse_buses_with_bb_switch, **opt_vars) 

86  
87  
88  1 
def remove_bad_data(net, init='flat', tolerance=1e6, maximum_iterations=10, 
89 
calculate_voltage_angles=True, rn_max_threshold=3.0): 

90 
"""


91 
Wrapper function for bad data removal.


92  
93 
INPUT:


94 
**net**  The net within this line should be created


95  
96 
**init**  (string) Initial voltage for the estimation. 'flat' sets 1.0 p.u. / 0° for all


97 
buses, 'results' uses the values from *res_bus_est* if available and 'slack' considers the


98 
slack bus voltage (and optionally, angle) as the initial values. Default is 'flat'


99  
100 
OPTIONAL:


101 
**tolerance**  (float)  When the maximum state change between iterations is less than


102 
tolerance, the process stops. Default is 1e6


103  
104 
**maximum_iterations**  (integer)  Maximum number of iterations. Default is 10


105  
106 
**calculate_voltage_angles**  (boolean)  Take into account absolute voltage angles and phase


107 
shifts in transformers, if init is 'slack'. Default is True


108  
109 
**rn_max_threshold** (float)  Identification threshold to determine


110 
if the largest normalized residual reflects a bad measurement


111 
(default value of 3.0)


112  
113 
OUTPUT:


114 
**successful** (boolean)  Was the state estimation successful?


115 
"""


116  1 
wls_se = StateEstimation(net, tolerance, maximum_iterations, algorithm="wls") 
117  1 
v_start, delta_start = _initialize_voltage(net, init, calculate_voltage_angles) 
118  1 
return wls_se.perform_rn_max_test(v_start, delta_start, calculate_voltage_angles, 
119 
rn_max_threshold) 

120  
121  
122  1 
def chi2_analysis(net, init='flat', tolerance=1e6, maximum_iterations=10, 
123 
calculate_voltage_angles=True, chi2_prob_false=0.05): 

124 
"""


125 
Wrapper function for the chisquared test.


126  
127 
INPUT:


128 
**net**  The net within this line should be created.


129  
130 
**init**  (string) Initial voltage for the estimation. 'flat' sets 1.0 p.u. / 0° for all


131 
buses, 'results' uses the values from *res_bus_est* if available and 'slack' considers the


132 
slack bus voltage (and optionally, angle) as the initial values. Default is 'flat'


133  
134 
OPTIONAL:


135 
**tolerance**  (float)  When the maximum state change between iterations is less than


136 
tolerance, the process stops. Default is 1e6


137  
138 
**maximum_iterations**  (integer)  Maximum number of iterations. Default is 10


139  
140 
**calculate_voltage_angles**  (boolean)  Take into account absolute voltage angles and phase


141 
shifts in transformers, if init is 'slack'. Default is True


142  
143 
**chi2_prob_false** (float)  probability of error / false alarms


144 
(default value: 0.05)


145  
146 
OUTPUT:


147 
**bad_data_detected** (boolean)  Returns true if bad data has been detected


148 
"""


149  1 
wls_se = StateEstimation(net, tolerance, maximum_iterations, algorithm="wls") 
150  1 
v_start, delta_start = _initialize_voltage(net, init, calculate_voltage_angles) 
151  1 
return wls_se.perform_chi2_test(v_start, delta_start, calculate_voltage_angles, 
152 
chi2_prob_false) 

153  
154  
155  1 
class StateEstimation: 
156 
"""


157 
Any user of the estimation module only needs to use the class state_estimation. It contains all


158 
relevant functions to control and operator the module. Two functions are used to configure the


159 
system according to the users needs while one function is used for the actual estimation


160 
process.


161 
"""


162  
163  1 
def __init__(self, net, tolerance=1e6, maximum_iterations=10, algorithm='wls', logger=None, recycle=False): 
164  1 
self.logger = logger 
165  1 
if self.logger is None: 
166  1 
self.logger = std_logger 
167 
# self.logger.setLevel(logging.DEBUG)


168  1 
self.net = net 
169  1 
self.solver = ALGORITHM_MAPPING[algorithm](tolerance, 
170 
maximum_iterations, self.logger) 

171  1 
self.ppc = None 
172  1 
self.eppci = None 
173  1 
self.recycle = recycle 
174  
175 
# variables for chi^2 / rn_max tests


176  1 
self.delta = None 
177  1 
self.bad_data_present = None 
178  
179  1 
def estimate(self, v_start='flat', delta_start='flat', calculate_voltage_angles=True, 
180 
zero_injection=None, fuse_buses_with_bb_switch='all', **opt_vars): 

181 
"""


182 
The function estimate is the main function of the module. It takes up to three input


183 
arguments: v_start, delta_start and calculate_voltage_angles. The first two are the initial


184 
state variables for the estimation process. Usually they can be initialized in a


185 
"flatstart" condition: All voltages being 1.0 pu and all voltage angles being 0 degrees.


186 
In this case, the parameters can be left at their default values (None). If the estimation


187 
is applied continuously, using the results from the last estimation as the starting


188 
condition for the current estimation can decrease the amount of iterations needed to


189 
estimate the current state. The third parameter defines whether all voltage angles are


190 
calculated absolutely, including phase shifts from transformers. If only the relative


191 
differences between buses are required, this parameter can be set to False. Returned is a


192 
boolean value, which is true after a successful estimation and false otherwise.


193 
The resulting complex voltage will be written into the pandapower network. The result


194 
fields are found res_bus_est of the pandapower network.


195 
INPUT:


196 
**net**  The net within this line should be created


197 
**v_start** (np.array, shape=(1,), optional)  Vector with initial values for all


198 
voltage magnitudes in p.u. (sorted by bus index)


199 
**delta_start** (np.array, shape=(1,), optional)  Vector with initial values for all


200 
voltage angles in degrees (sorted by bus index)


201 
OPTIONAL:


202 
**tolerance**  (float)  When the maximum state change between iterations is less than


203 
tolerance, the process stops. Default is 1e6


204  
205 
**maximum_iterations**  (integer)  Maximum number of iterations. Default is 10


206  
207 
**calculate_voltage_angles**  (boolean)  Take into account absolute voltage angles and phase


208 
shifts in transformers, if init is 'slack'. Default is True


209 


210 
**zero_injection**  (str, iterable, None)  Defines which buses are zero injection bus or the method


211 
to identify zero injection bus, with 'wls_estimator' virtual measurements will be added, with


212 
'wls_estimator with zero constraints' the buses will be handled as constraints


213 
"auto": all bus without p,q measurement, without p, q value (load, sgen...) and aux buses will be


214 
identified as zero injection bus


215 
"aux_bus": only aux bus will be identified as zero injection bus


216 
None: no bus will be identified as zero injection bus


217 
iterable: the iterable should contain index of the zero injection bus and also aux bus will be identified


218 
as zeroinjection bus


219  
220 
**fuse_buses_with_bb_switch**  (str, iterable, None)  Defines how buses with closed bb switches should


221 
be handled, if fuse buses will only fused to one for calculation, if not fuse, an auxiliary bus and


222 
auxiliary line will be automatically added to the network to make the buses with different p,q injection


223 
measurements identifieble


224 
"all": all buses with bbswitches will be fused, the same as the default behaviour in load flow


225 
None: buses with bbswitches and individual p,q measurements will be reconfigurated


226 
by auxiliary elements


227 
iterable: the iterable should contain the index of buses to be fused, the behaviour is contigous e.g.


228 
if one of the bus among the buses connected through bb switch is given, then all of them will still


229 
be fused


230 
OUTPUT:


231 
**successful** (boolean)  True if the estimation process was successful


232 
Optional estimation variables:


233 
The bus power injections can be accessed with *se.s_node_powers* and the estimated


234 
values corresponding to the (noisy) measurement values with *se.hx*. (*hx* denotes h(x))


235 
EXAMPLE:


236 
success = estimate(np.array([1.0, 1.0, 1.0]), np.array([0.0, 0.0, 0.0]))


237 
"""


238 
# check if all parameter are allowed


239  1 
for var_name in opt_vars.keys(): 
240  1 
if var_name not in ALLOWED_OPT_VAR: 
241 
self.logger.warning("Caution! %s is not allowed as parameter" % var_name \ 

242 
+ " for estimate and will be ignored!") 

243  
244  1 
if self.net is None: 
245  1 
raise UserWarning("SE Component was not initialized with a network.") 
246  
247 
# change the configuration of the pp net to avoid auto fusing of buses connected


248 
# through bb switch with elements on each bus if this feature enabled


249  1 
bus_to_be_fused = None 
250  1 
if fuse_buses_with_bb_switch != 'all' and not self.net.switch.empty: 
251  1 
if isinstance(fuse_buses_with_bb_switch, str): 
252  1 
raise UserWarning("fuse_buses_with_bb_switch parameter is not correctly initialized") 
253  1 
elif hasattr(fuse_buses_with_bb_switch, '__iter__'): 
254  1 
bus_to_be_fused = fuse_buses_with_bb_switch 
255  1 
set_bb_switch_impedance(self.net, bus_to_be_fused) 
256  
257  1 
self.net, self.ppc, self.eppci = pp2eppci(self.net, v_start=v_start, delta_start=delta_start, 
258 
calculate_voltage_angles=calculate_voltage_angles, 

259 
zero_injection=zero_injection, ppc=self.ppc, eppci=self.eppci) 

260  
261 
# Estimate voltage magnitude and angle with the given estimator


262  1 
self.eppci = self.solver.estimate(self.eppci, **opt_vars) 
263  
264  1 
if self.solver.successful: 
265  1 
self.net = eppci2pp(self.net, self.ppc, self.eppci) 
266 
else: 

267  1 
self.logger.warning("Estimation failed! Pandapower network failed to update!") 
268  
269 
# clear the aux elements and calculation results created for the substitution of bb switches


270  1 
if fuse_buses_with_bb_switch != 'all' and not self.net.switch.empty: 
271  1 
reset_bb_switch_impedance(self.net) 
272  
273 
# if recycle is not wished, reset ppc, ppci


274  1 
if not self.recycle: 
275  1 
self.ppc, self.eppci = None, None 
276  1 
return self.solver.successful 
277  
278  1 
def perform_chi2_test(self, v_in_out=None, delta_in_out=None, 
279 
calculate_voltage_angles=True, chi2_prob_false=0.05): 

280 
"""


281 
The function perform_chi2_test performs a Chi^2 test for bad data and topology error


282 
detection. The function can be called with the optional input arguments v_in_out and


283 
delta_in_out. Then, the Chi^2 test is performed after calling the function estimate using


284 
them as input arguments. It can also be called without these arguments if it is called


285 
from the same object with which estimate had been called beforehand. Then, the Chi^2 test is


286 
performed for the states estimated by the funtion estimate and the result, the existence of bad data,


287 
is given back as a boolean. As a optional argument the probability


288 
of a false measurement can be provided additionally. For bad data detection, the function


289 
perform_rn_max_test is more powerful and should be the function of choice. For topology


290 
error detection, however, perform_chi2_test should be used.


291  
292 
INPUT:


293 
**v_in_out** (np.array, shape=(1,), optional)  Vector with initial values for all


294 
voltage magnitudes in p.u. (sorted by bus index)


295  
296 
**delta_in_out** (np.array, shape=(1,), optional)  Vector with initial values for all


297 
voltage angles in degrees (sorted by bus index)


298  
299 
OPTIONAL:


300 
**calculate_voltage_angles**  (boolean)  Take into account absolute voltage angles and phase


301 
shifts in transformers, if init is 'slack'. Default is True


302  
303 
**chi2_prob_false** (float)  probability of error / false alarms (standard value: 0.05)


304  
305 
OUTPUT:


306 
**successful** (boolean)  True if bad data has been detected


307  
308 
EXAMPLE:


309 
perform_chi2_test(np.array([1.0, 1.0, 1.0]), np.array([0.0, 0.0, 0.0]), 0.97)


310  
311 
"""


312 
# perform SE


313  1 
self.estimate(v_in_out, delta_in_out, calculate_voltage_angles) 
314  
315 
# Performance index J(hx)


316  1 
J = np.dot(self.solver.r.T, np.dot(self.solver.R_inv, self.solver.r)) 
317  
318 
# Number of measurements


319  1 
m = len(self.net.measurement) 
320  
321 
# Number of state variables (the 1 is due to the reference bus)


322  1 
n = len(self.solver.eppci.v) + len(self.solver.eppci.delta)  1 
323  
324 
# Chi^2 test threshold


325  1 
test_thresh = chi2.ppf(1  chi2_prob_false, m  n) 
326  
327 
# Print results


328  1 
self.logger.debug("Result of Chi^2 test:") 
329  1 
self.logger.debug("Number of measurements: %d" % m) 
330  1 
self.logger.debug("Number of state variables: %d" % n) 
331  1 
self.logger.debug("Performance index: %.2f" % J) 
332  1 
self.logger.debug("Chi^2 test threshold: %.2f" % test_thresh) 
333  
334  1 
if J <= test_thresh: 
335  1 
self.bad_data_present = False 
336  1 
self.logger.debug("Chi^2 test passed. No bad data or topology error detected.") 
337 
else: 

338  1 
self.bad_data_present = True 
339  1 
self.logger.debug("Chi^2 test failed. Bad data or topology error detected.") 
340  
341  1 
if self.solver.successful: 
342  1 
return self.bad_data_present 
343  
344  1 
def perform_rn_max_test(self, v_in_out=None, delta_in_out=None, 
345 
calculate_voltage_angles=True, rn_max_threshold=3.0): 

346 
"""


347 
The function perform_rn_max_test performs a largest normalized residual test for bad data


348 
identification and removal. It takes two input arguments: v_in_out and delta_in_out.


349 
These are the initial state variables for the combined estimation and bad data


350 
identification and removal process. They can be initialized as described above, e.g.,


351 
using a "flat" start. In an iterative process, the function performs a state estimation,


352 
identifies a bad data measurement, removes it from the set of measurements


353 
(only if the rn_max threshold is violated by the largest residual of all measurements,


354 
which can be modified), performs the state estimation again,


355 
and so on and so forth until no further bad data measurements are detected.


356  
357 
INPUT:


358 
**v_in_out** (np.array, shape=(1,), optional)  Vector with initial values for all


359 
voltage magnitudes in p.u. (sorted by bus index)


360  
361 
**delta_in_out** (np.array, shape=(1,), optional)  Vector with initial values for all


362 
voltage angles in degrees (sorted by bus index)


363  
364 
OPTIONAL:


365 
**calculate_voltage_angles**  (boolean)  Take into account absolute voltage angles and phase


366 
shifts in transformers, if init is 'slack'. Default is True


367  
368 
**rn_max_threshold** (float)  Identification threshold to determine


369 
if the largest normalized residual reflects a bad measurement


370 
(standard value of 3.0)


371  
372 
OUTPUT:


373 
**successful** (boolean)  True if all bad data could be removed


374  
375 
EXAMPLE:


376 
perform_rn_max_test(np.array([1.0, 1.0, 1.0]), np.array([0.0, 0.0, 0.0]), 5.0, 0.05)


377  
378 
"""


379  1 
num_iterations = 0 
380  
381  1 
while num_iterations <= 10: 
382 
# Estimate the state with bad data identified in previous iteration


383 
# removed from set of measurements:


384  1 
self.estimate(v_in_out, delta_in_out, calculate_voltage_angles) 
385  
386 
# Try to remove the bad data


387  1 
try: 
388 
# Error covariance matrix:


389  1 
R = np.linalg.inv(self.solver.R_inv) 
390  
391 
# for future debugging: this line's results have changed with the ppc


392 
# overhaul in April 2017 after commit 9ae5b8f42f69ae39f8c8cf (which still works)


393 
# there are differences of < 1e10 for the Omega entries which cause


394 
# the function to work far worse. As of now it is unclear if it's just numerical


395 
# accuracy to blame or an error in the code. a sort in the ppc creation function


396 
# was removed which caused this issue


397 
# Covariance matrix of the residuals: \Omega = S*R = R  H*G^(1)*H^T


398 
# (S is the sensitivity matrix: r = S*e):


399  1 
Omega = R  np.dot(self.solver.H, np.dot(np.linalg.inv(self.solver.Gm), self.solver.H.T)) 
400  
401 
# Diagonalize \Omega:


402  1 
Omega = np.diag(np.diag(Omega)) 
403  
404 
# Compute squareroot (. since some 0.0 produced nans):


405  1 
Omega = np.sqrt(np.absolute(Omega)) 
406  
407  1 
OmegaInv = np.linalg.inv(Omega) 
408  
409 
# Compute normalized residuals (r^N_i = r_i/sqrt{Omega_ii}):


410  1 
rN = np.dot(OmegaInv, np.absolute(self.solver.r)) 
411  
412  1 
if max(rN) <= rn_max_threshold: 
413  1 
self.logger.debug("Largest normalized residual test passed. " 
414 
"No bad data detected.") 

415  1 
return True 
416 
else: 

417  1 
self.logger.debug( 
418 
"Largest normalized residual test failed (%.1f > %.1f)." 

419 
% (max(rN), rn_max_threshold)) 

420  
421 
# Identify bad data: Determine index corresponding to max(rN):


422  1 
idx_rN = np.argsort(rN, axis=0)[1] 
423  
424 
# Determine pandapower index of measurement to be removed:


425  1 
meas_idx = self.solver.pp_meas_indices[idx_rN] 
426  
427 
# Remove bad measurement:


428  1 
self.logger.debug("Removing measurement: %s" 
429 
% self.net.measurement.loc[meas_idx].values[0]) 

430  1 
self.net.measurement.drop(meas_idx, inplace=True) 
431  1 
self.logger.debug("Bad data removed from the set of measurements.") 
432  
433 
except np.linalg.linalg.LinAlgError: 

434 
self.logger.error("A problem appeared while using the linear algebra methods." 

435 
"Check and change the measurement set.") 

436 
return False 

437  
438  1 
self.logger.debug("rN_max identification threshold: %.2f" % rn_max_threshold) 
439  1 
num_iterations += 1 
440  
441 
return False 
Read our documentation on viewing source code .