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 .

Loading