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.sparse import csr_matrix, vstack, hstack ``` 8 1 ```from scipy.sparse.linalg import spsolve ``` 9 10 1 ```from pandapower.estimation.algorithm.estimator import BaseEstimatorIRWLS, get_estimator ``` 11 1 ```from pandapower.estimation.algorithm.matrix_base import BaseAlgebra, \ ``` 12 ``` BaseAlgebraZeroInjConstraints ``` 13 1 ```from pandapower.estimation.idx_bus import ZERO_INJ_FLAG, P, P_STD, Q, Q_STD ``` 14 1 ```from pandapower.estimation.ppc_conversion import ExtendedPPCI ``` 15 1 ```from pandapower.pypower.idx_bus import bus_cols ``` 16 17 1 ```try: ``` 18 1 ``` import pplog as logging ``` 19 1 ```except ImportError: ``` 20 1 ``` import logging ``` 21 1 ```std_logger = logging.getLogger(__name__) ``` 22 23 1 ```__all__ = ["WLSAlgorithm", "WLSZeroInjectionConstraintsAlgorithm", "IRWLSAlgorithm"] ``` 24 25 26 1 ```class BaseAlgorithm: ``` 27 1 ``` def __init__(self, tolerance, maximum_iterations, logger=std_logger): ``` 28 1 ``` self.tolerance = tolerance ``` 29 1 ``` self.max_iterations = maximum_iterations ``` 30 1 ``` self.logger = logger ``` 31 1 ``` self.successful = False ``` 32 1 ``` self.iterations = None ``` 33 34 ``` # Parameters for estimate ``` 35 1 ``` self.eppci = None ``` 36 1 ``` self.pp_meas_indices = None ``` 37 38 1 ``` def check_observability(self, eppci: ExtendedPPCI, z): ``` 39 ``` # Check if observability criterion is fulfilled and the state estimation is possible ``` 40 1 ``` if len(z) < 2 * eppci["bus"].shape[0] - 1: ``` 41 1 ``` self.logger.error("System is not observable (cancelling)") ``` 42 1 ``` self.logger.error("Measurements available: %d. Measurements required: %d" % ``` 43 ``` (len(z), 2 * eppci["bus"].shape[0] - 1)) ``` 44 1 ``` raise UserWarning("Measurements available: %d. Measurements required: %d" % ``` 45 ``` (len(z), 2 * eppci["bus"].shape[0] - 1)) ``` 46 47 1 ``` def check_result(self, current_error, cur_it): ``` 48 ``` # print output for results ``` 49 1 ``` if current_error <= self.tolerance: ``` 50 1 ``` self.successful = True ``` 51 1 ``` self.logger.debug("State Estimation successful ({:d} iterations)".format(cur_it)) ``` 52 ``` else: ``` 53 1 ``` self.successful = False ``` 54 1 ``` self.logger.debug("State Estimation not successful ({:d}/{:d} iterations)".format(cur_it, ``` 55 ``` self.max_iterations)) ``` 56 57 1 ``` def initialize(self, eppci: ExtendedPPCI): ``` 58 ``` # Check observability ``` 59 1 ``` self.eppci = eppci ``` 60 1 ``` self.pp_meas_indices = eppci.pp_meas_indices ``` 61 1 ``` self.check_observability(eppci, eppci.z) ``` 62 63 1 ``` def estimate(self, eppci: ExtendedPPCI, **kwargs): ``` 64 ``` # Must be implemented individually!! ``` 65 0 ``` pass ``` 66 67 68 1 ```class WLSAlgorithm(BaseAlgorithm): ``` 69 1 ``` def __init__(self, tolerance, maximum_iterations, logger=std_logger): ``` 70 1 ``` super(WLSAlgorithm, self).__init__(tolerance, maximum_iterations, logger) ``` 71 72 ``` # Parameters for Bad data detection ``` 73 1 ``` self.R_inv = None ``` 74 1 ``` self.Gm = None ``` 75 1 ``` self.r = None ``` 76 1 ``` self.H = None ``` 77 1 ``` self.hx = None ``` 78 79 1 ``` def estimate(self, eppci: ExtendedPPCI, **kwargs): ``` 80 1 ``` self.initialize(eppci) ``` 81 ``` # matrix calculation object ``` 82 1 ``` sem = BaseAlgebra(eppci) ``` 83 84 1 ``` current_error, cur_it = 100., 0 ``` 85 ``` # invert covariance matrix ``` 86 1 ``` r_inv = csr_matrix(np.diagflat(1 / eppci.r_cov ** 2)) ``` 87 1 ``` E = eppci.E ``` 88 1 ``` while current_error > self.tolerance and cur_it < self.max_iterations: ``` 89 1 ``` self.logger.debug("Starting iteration {:d}".format(1 + cur_it)) ``` 90 1 ``` try: ``` 91 ``` # residual r ``` 92 1 ``` r = csr_matrix(sem.create_rx(E)).T ``` 93 94 ``` # jacobian matrix H ``` 95 1 ``` H = csr_matrix(sem.create_hx_jacobian(E)) ``` 96 97 ``` # gain matrix G_m ``` 98 ``` # G_m = H^t * R^-1 * H ``` 99 1 ``` G_m = H.T * (r_inv * H) ``` 100 101 ``` # state vector difference d_E ``` 102 ``` # d_E = G_m^-1 * (H' * R^-1 * r) ``` 103 1 ``` d_E = spsolve(G_m, H.T * (r_inv * r)) ``` 104 105 ``` # Update E with d_E ``` 106 1 ``` E += d_E.ravel() ``` 107 1 ``` eppci.update_E(E) ``` 108 109 ``` # prepare next iteration ``` 110 1 ``` cur_it += 1 ``` 111 1 ``` current_error = np.max(np.abs(d_E)) ``` 112 1 ``` self.logger.debug("Current error: {:.7f}".format(current_error)) ``` 113 0 ``` except np.linalg.linalg.LinAlgError: ``` 114 0 ``` self.logger.error("A problem appeared while using the linear algebra methods." ``` 115 ``` "Check and change the measurement set.") ``` 116 0 ``` return False ``` 117 118 ``` # check if the estimation is successfull ``` 119 1 ``` self.check_result(current_error, cur_it) ``` 120 1 ``` if self.successful: ``` 121 ``` # store variables required for chi^2 and r_N_max test: ``` 122 1 ``` self.R_inv = r_inv.toarray() ``` 123 1 ``` self.Gm = G_m.toarray() ``` 124 1 ``` self.r = r.toarray() ``` 125 1 ``` self.H = H.toarray() ``` 126 ``` # create h(x) for the current iteration ``` 127 1 ``` self.hx = sem.create_hx(eppci.E) ``` 128 1 ``` return eppci ``` 129 130 131 1 ```class WLSZeroInjectionConstraintsAlgorithm(BaseAlgorithm): ``` 132 1 ``` def estimate(self, eppci: ExtendedPPCI, **kwargs): ``` 133 ``` # state vector built from delta, |V| and zero injections ``` 134 ``` # Find pq bus with zero p,q and shunt admittance ``` 135 1 ``` if not np.any(eppci["bus"][:, bus_cols + ZERO_INJ_FLAG]): ``` 136 0 ``` raise UserWarning("Network has no bus with zero injections! Please use WLS instead!") ``` 137 1 ``` zero_injection_bus = np.argwhere(eppci["bus"][:, bus_cols + ZERO_INJ_FLAG] == True).ravel() ``` 138 1 ``` eppci["bus"][np.ix_(zero_injection_bus, [bus_cols + P, bus_cols + P_STD, bus_cols + Q, bus_cols + Q_STD])] = np.NaN ``` 139 ``` # Withn pq buses with zero injection identify those who have also no p or q measurement ``` 140 1 ``` p_zero_injections = zero_injection_bus ``` 141 1 ``` q_zero_injections = zero_injection_bus ``` 142 1 ``` new_states = np.zeros(len(p_zero_injections) + len(q_zero_injections)) ``` 143 144 1 ``` num_bus = eppci["bus"].shape[0] ``` 145 146 ``` # matrix calculation object ``` 147 1 ``` sem = BaseAlgebraZeroInjConstraints(eppci) ``` 148 149 1 ``` current_error, cur_it = 100., 0 ``` 150 1 ``` r_inv = csr_matrix((np.diagflat(1 / eppci.r_cov) ** 2)) ``` 151 1 ``` E = eppci.E ``` 152 ``` # update the E matrix ``` 153 1 ``` E_ext = np.r_[eppci.E, new_states] ``` 154 155 1 ``` while current_error > self.tolerance and cur_it < self.max_iterations: ``` 156 1 ``` self.logger.debug("Starting iteration {:d}".format(1 + cur_it)) ``` 157 1 ``` try: ``` 158 1 ``` c_x = sem.create_cx(E, p_zero_injections, q_zero_injections) ``` 159 160 ``` # residual r ``` 161 1 ``` r = csr_matrix(sem.create_rx(E)).T ``` 162 1 ``` c_rxh = csr_matrix(c_x).T ``` 163 164 ``` # jacobian matrix H ``` 165 1 ``` H_temp = sem.create_hx_jacobian(E) ``` 166 1 ``` C_temp = sem.create_cx_jacobian(E, p_zero_injections, q_zero_injections) ``` 167 1 ``` H, C = csr_matrix(H_temp), csr_matrix(C_temp) ``` 168 169 ``` # gain matrix G_m ``` 170 ``` # G_m = H^t * R^-1 * H ``` 171 1 ``` G_m = H.T * (r_inv * H) ``` 172 173 ``` # building a new gain matrix for new constraints. ``` 174 1 ``` A_1 = vstack([G_m, C]) ``` 175 1 ``` c_ax = hstack([C, np.zeros((C.shape[0], C.shape[0]))]) ``` 176 1 ``` c_xT = c_ax.T ``` 177 1 ``` M_tx = csr_matrix(hstack((A_1, c_xT))) # again adding to the new gain matrix ``` 178 1 ``` rhs = H.T * (r_inv * r) # original right hand side ``` 179 1 ``` C_rhs = vstack((rhs, -c_rxh)) # creating the righ hand side with new constraints ``` 180 181 ``` # state vector difference d_E and update E ``` 182 1 ``` d_E_ext = spsolve(M_tx, C_rhs) ``` 183 1 ``` E_ext += d_E_ext.ravel() ``` 184 1 ``` E = E_ext[:E.shape[0]] ``` 185 1 ``` eppci.update_E(E) ``` 186 187 ``` # prepare next iteration ``` 188 1 ``` cur_it += 1 ``` 189 1 ``` current_error = np.max(np.abs(d_E_ext[:len(eppci.non_slack_buses) + num_bus])) ``` 190 1 ``` self.logger.debug("Current error: {:.7f}".format(current_error)) ``` 191 0 ``` except np.linalg.linalg.LinAlgError: ``` 192 0 ``` self.logger.error("A problem appeared while using the linear algebra methods." ``` 193 ``` "Check and change the measurement set.") ``` 194 0 ``` return False ``` 195 196 ``` # check if the estimation is successfull ``` 197 1 ``` self.check_result(current_error, cur_it) ``` 198 1 ``` return eppci ``` 199 200 201 1 ```class IRWLSAlgorithm(BaseAlgorithm): ``` 202 1 ``` def estimate(self, eppci: ExtendedPPCI, estimator="wls", **kwargs): ``` 203 1 ``` self.initialize(eppci) ``` 204 205 ``` # matrix calculation object ``` 206 1 ``` sem = get_estimator(BaseEstimatorIRWLS, estimator)(eppci, **kwargs) ``` 207 208 1 ``` current_error, cur_it = 100., 0 ``` 209 1 ``` E = eppci.E ``` 210 1 ``` while current_error > self.tolerance and cur_it < self.max_iterations: ``` 211 1 ``` self.logger.debug("Starting iteration {:d}".format(1 + cur_it)) ``` 212 1 ``` try: ``` 213 ``` # residual r ``` 214 1 ``` r = csr_matrix(sem.create_rx(E)).T ``` 215 216 ``` # jacobian matrix H ``` 217 1 ``` H = csr_matrix(sem.create_hx_jacobian(E)) ``` 218 219 ``` # gain matrix G_m ``` 220 ``` # G_m = H^t * Phi * H ``` 221 1 ``` phi = csr_matrix(sem.create_phi(E)) ``` 222 1 ``` G_m = H.T * (phi * H) ``` 223 224 ``` # state vector difference d_E and update E ``` 225 1 ``` d_E = spsolve(G_m, H.T * (phi * r)) ``` 226 1 ``` E += d_E.ravel() ``` 227 1 ``` eppci.update_E(E) ``` 228 229 ``` # prepare next iteration ``` 230 1 ``` cur_it += 1 ``` 231 1 ``` current_error = np.max(np.abs(d_E)) ``` 232 1 ``` self.logger.debug("Current error: {:.7f}".format(current_error)) ``` 233 0 ``` except np.linalg.linalg.LinAlgError: ``` 234 0 ``` self.logger.error("A problem appeared while using the linear algebra methods." ``` 235 ``` "Check and change the measurement set.") ``` 236 0 ``` return False ``` 237 238 ``` # check if the estimation is successfull ``` 239 1 ``` self.check_result(current_error, cur_it) ``` 240 ``` # update V/delta ``` 241 1 ``` return eppci ```

Read our documentation on viewing source code .