e2nIEE / pandapower
1
# -*- coding: utf-8 -*-
2

3
# Copyright (c) 2016-2021 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=1e-6, 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 1e-6
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 zero-injection 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 bb-switches will be fused, the same as the default behaviour in load flow
69
            None: buses with bb-switches 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=1e-6, 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 1e-6
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=1e-6, maximum_iterations=10,
123
                  calculate_voltage_angles=True, chi2_prob_false=0.05):
124
    """
125
    Wrapper function for the chi-squared 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 1e-6
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=1e-6, 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
        "flat-start" 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 1e-6
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 zero-injection 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 bb-switches will be fused, the same as the default behaviour in load flow
225
            None: buses with bb-switches 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 0
                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 < 1e-10 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 0
            except np.linalg.linalg.LinAlgError:
434 0
                self.logger.error("A problem appeared while using the linear algebra methods."
435
                                  "Check and change the measurement set.")
436 0
                return False
437

438 1
            self.logger.debug("rN_max identification threshold: %.2f" % rn_max_threshold)
439 1
            num_iterations += 1
440

441 0
        return False

Read our documentation on viewing source code .

Loading