diff --git a/pandapower/estimation/algorithm/base.py b/pandapower/estimation/algorithm/base.py index f5f252089..dc8d6a2c0 100644 --- a/pandapower/estimation/algorithm/base.py +++ b/pandapower/estimation/algorithm/base.py @@ -37,7 +37,8 @@ def __init__(self, tolerance, maximum_iterations, logger=std_logger): def check_observability(self, eppci: ExtendedPPCI, z): # Check if observability criterion is fulfilled and the state estimation is possible - if len(z) < 2 * eppci["bus"].shape[0] - 1: + num_slacks = sum(~eppci.non_slack_bus_mask) + if len(z) < 2 * eppci["bus"].shape[0] - num_slacks: self.logger.error("System is not observable (cancelling)") self.logger.error("Measurements available: %d. Measurements required: %d" % (len(z), 2 * eppci["bus"].shape[0] - 1)) @@ -272,3 +273,100 @@ def estimate(self, eppci: ExtendedPPCI, estimator="wls", **kwargs): self.check_result(current_error, cur_it) # update V/delta return eppci + + +class AFWLSAlgorithm(BaseAlgorithm): + def __init__(self, tolerance, maximum_iterations, logger=std_logger): + super(AFWLSAlgorithm, self).__init__(tolerance, maximum_iterations, logger) + + # Parameters for Bad data detection + self.R_inv = None + self.Gm = None + self.r = None + self.H = None + self.hx = None + self.iterations = None + self.obj_func = None + logging.basicConfig(level=logging.DEBUG) + + def estimate(self, eppci: ExtendedPPCI, **kwargs): + self.initialize(eppci) + # matrix calculation object + sem = BaseAlgebra(eppci) + + current_error, cur_it = 100., 0 + # invert covariance matrix + eppci.r_cov[eppci.r_cov<(10**(-5))] = 10**(-5) + r_inv = csr_matrix(np.diagflat(1 / eppci.r_cov ** 2)) + E = eppci.E + num_clusters = len(self.eppci["clusters"]) + while current_error > self.tolerance and cur_it < self.max_iterations: + self.logger.debug("Starting iteration {:d}".format(1 + cur_it)) + try: + # residual r + r = csr_matrix(sem.create_rx(E)).T + + # jacobian matrix H + H = csr_matrix(sem.create_hx_jacobian(E)) + + if cur_it == 0 and eppci.any_i_meas: + idx = eppci.idx_non_imeas + r_inv = r_inv[idx,:][:,idx] + r = r[idx,:] + H = H[idx,:] + + # gain matrix G_m + G_m = H.T * (r_inv * H) + norm_G = norm(G_m, np.inf) + norm_invG = norm(inv(G_m), np.inf) + cond = norm_G*norm_invG + if cond > 10**18: + self.logger.warning("WARNING: Gain matrix is ill-conditioned: {:.2E}".format(cond)) + + # state vector difference d_E + d_E = spsolve(G_m, H.T * (r_inv * r)) + + # # Scaling of Delta_X to avoid divergence due o ill-conditioning and + # # operating conditions far from starting state variables + # current_error = np.max(np.abs(d_E)) + # if current_error > 0.25: + # d_E = d_E*0.25/current_error + + # Update E with d_E + E += d_E.ravel() + # eppci.update_E(E1) + + # log data + current_error = np.max(np.abs(d_E)) + obj_func = (r.T*r_inv*r)[0,0] + self.logger.debug("Current delta_x: {:.7f}".format(current_error)) + self.logger.debug("Current objective function value: {:.1f}".format(obj_func)) + + # Restore full weighting matrix + if cur_it == 0 and eppci.any_i_meas: + r_inv = csr_matrix(np.diagflat(1 / eppci.r_cov ** 2)) + + # prepare next iteration + cur_it += 1 + + except np.linalg.linalg.LinAlgError: + self.logger.error("A problem appeared while using the linear algebra methods." + "Check and change the measurement set.") + return False + + # check if the estimation is successfull + self.check_result(current_error, cur_it) + self.iterations = cur_it + self.obj_func = obj_func + if self.successful: + # store variables required for chi^2 and r_N_max test: + self.R_inv = r_inv.toarray() + self.Gm = G_m.toarray() + self.r = r.toarray() + self.H = H.toarray() + # split voltage and allocation factor variables + E1 = E[:-num_clusters] + E2 = E[-num_clusters:] + eppci.update_E(E1) + eppci.clusters = E2 + return eppci \ No newline at end of file diff --git a/pandapower/estimation/algorithm/matrix_base.py b/pandapower/estimation/algorithm/matrix_base.py index ac1bfdd74..de3b778ed 100644 --- a/pandapower/estimation/algorithm/matrix_base.py +++ b/pandapower/estimation/algorithm/matrix_base.py @@ -52,7 +52,13 @@ def create_rx(self, E): def create_hx(self, E): f_bus, t_bus = self.fb, self.tb - V = self.eppci.E2V(E) + if self.eppci.algorithm == "af-wls": + num_clusters = len(self.eppci["clusters"]) + E1 = E[:-num_clusters] + E2 = E[-num_clusters:] + else: + E1 = E + V = self.eppci.E2V(E1) Sfe = V[f_bus] * np.conj(self.Yf * V) Ste = V[t_bus] * np.conj(self.Yt * V) Sbuse = V * np.conj(self.Ybus * V) @@ -64,25 +70,40 @@ def create_hx(self, E): np.imag(Ste), np.abs(V)] - if self.any_i_meas or self.any_degree_meas: - va = np.angle(V) - Ife = self.Yf * V - ifem = np.abs(Ife) - ifea = np.angle(Ife) - Ite = self.Yt * V - item = np.abs(Ite) - itea = np.angle(Ite) + # if self.any_i_meas or self.any_degree_meas: + va = np.angle(V) + Ife = self.Yf * V + ifem = np.abs(Ife) + ifea = np.angle(Ife) + Ite = self.Yt * V + item = np.abs(Ite) + itea = np.angle(Ite) + hx = np.r_[hx, + va, + ifem, + item, + ifea, + itea] + + if self.eppci.algorithm == "af-wls": + Pbuse2 = np.sum(np.multiply(E2,self.eppci["rated_power_clusters"][:,:num_clusters]),axis=1) + Qbuse2 = np.sum(np.multiply(E2,self.eppci["rated_power_clusters"][:,num_clusters:2*num_clusters]),axis=1) hx = np.r_[hx, - va, - ifem, - item, - ifea, - itea] + np.real(Sbuse)-Pbuse2, + np.imag(Sbuse)-Qbuse2, + E2] + return hx[self.non_nan_meas_selector] def create_hx_jacobian(self, E): # Using sparse matrix in creation sub-jacobian matrix - V = self.eppci.E2V(E) + if self.eppci.algorithm == "af-wls": + num_clusters = len(self.eppci["clusters"]) + E1 = E[:-num_clusters] + else: + E1 = E + + V = self.eppci.E2V(E1) dSbus_dth, dSbus_dv = self._dSbus_dv(V) dSf_dth, dSf_dv, dSt_dth, dSt_dv = self._dSbr_dv(V) @@ -106,29 +127,56 @@ def create_hx_jacobian(self, E): jac = np.r_[s_jac, vm_jac] - if self.any_i_meas or self.any_degree_meas: - dva_dth, dva_dv = self._dvabus_dV(V) - va_jac = np.c_[dva_dth, dva_dv] - difm_dth, difm_dv, ditm_dth, ditm_dv,\ - difa_dth, difa_dv, dita_dth, dita_dv = self._dimiabr_dV(V) - im_jac_th = np.r_[difm_dth, - ditm_dth] - im_jac_v = np.r_[difm_dv, - ditm_dv] - ia_jac_th = np.r_[difa_dth, - dita_dth] - ia_jac_v = np.r_[difa_dv, - dita_dv] - - im_jac = np.c_[im_jac_th, im_jac_v] - ia_jac = np.c_[ia_jac_th, ia_jac_v] + # if self.any_i_meas or self.any_degree_meas: + dva_dth, dva_dv = self._dvabus_dV(V) + va_jac = np.c_[dva_dth, dva_dv] + difm_dth, difm_dv, ditm_dth, ditm_dv,\ + difa_dth, difa_dv, dita_dth, dita_dv = self._dimiabr_dV(V) + im_jac_th = np.r_[difm_dth, + ditm_dth] + im_jac_v = np.r_[difm_dv, + ditm_dv] + ia_jac_th = np.r_[difa_dth, + dita_dth] + ia_jac_v = np.r_[difa_dv, + dita_dv] + + im_jac = np.c_[im_jac_th, im_jac_v] + ia_jac = np.c_[ia_jac_th, ia_jac_v] + + jac = np.r_[jac, + va_jac, + im_jac, + ia_jac] + + if self.eppci.algorithm == "af-wls": + p_eq_bal_jac_E1 = hstack((dSbus_dth.real, dSbus_dv.real)).toarray() + q_eq_bal_jac_E1 = hstack((dSbus_dth.imag, dSbus_dv.imag)).toarray() + af_vmeas_E1 = np.zeros((num_clusters,jac.shape[1])) + + jac_E2 = np.zeros((jac.shape[0],num_clusters)) + p_eq_bal_jac_E2 = - self.eppci["rated_power_clusters"][:,:num_clusters] + q_eq_bal_jac_E2 = - self.eppci["rated_power_clusters"][:,num_clusters:2*num_clusters] + af_vmeas_E2 = np.identity(num_clusters) jac = np.r_[jac, - va_jac, - im_jac, - ia_jac] + p_eq_bal_jac_E1, + q_eq_bal_jac_E1, + af_vmeas_E1] + + jac_E2 = np.r_[jac_E2, + p_eq_bal_jac_E2, + q_eq_bal_jac_E2, + af_vmeas_E2] + + jac = jac[self.non_nan_meas_selector, :][:, self.delta_v_bus_selector] + jac_E2 = jac_E2[self.non_nan_meas_selector, :][:] + jac = np.c_[jac, jac_E2] + + else: + jac = jac[self.non_nan_meas_selector, :][:, self.delta_v_bus_selector] - return jac[self.non_nan_meas_selector, :][:, self.delta_v_bus_selector] + return jac def _dSbus_dv(self, V): dSbus_dv, dSbus_dth = dSbus_dV(self.Ybus, V) diff --git a/pandapower/estimation/ppc_conversion.py b/pandapower/estimation/ppc_conversion.py index 79619336f..86867871d 100644 --- a/pandapower/estimation/ppc_conversion.py +++ b/pandapower/estimation/ppc_conversion.py @@ -313,6 +313,42 @@ def _add_measurements_to_ppci(net, ppci, zero_injection, algorithm): ppci["branch"] = np.hstack((ppci["branch"], branch_append)) else: ppci["branch"][:, branch_cols: branch_cols + branch_cols_se] = branch_append + + # Add rated power information for AF-WLS estimator + if algorithm == 'af-wls': + cluster_list_loads = net.load["type"].unique() + cluster_list_gen = net.sgen["type"].unique() + cluster_list_tot = np.concatenate((cluster_list_loads, cluster_list_gen), axis=0) + ppci["clusters"] = cluster_list_tot + num_clusters = len(cluster_list_tot) + num_buses = ppci["bus"].shape[0] + ppci["rated_power_clusters"] = np.zeros([num_buses, 4*num_clusters]) + for var in ["load", "sgen"]: + in_service = net[var]["in_service"] + active_elements = net[var][in_service] + bus = net._pd2ppc_lookups["bus"][active_elements.bus].astype(int) + P = active_elements.p_mw.values/ppci["baseMVA"] + Q = active_elements.q_mvar.values/ppci["baseMVA"] + if var == 'load': + P *= -1 + Q *= -1 + cluster = active_elements.type.values + if (bus >= ppci["bus"].shape[0]).any(): + std_logger.warning("Loads or sgen defined in pp-grid do not exist in ppci, will be deleted!") + P = P[bus < ppci["bus"].shape[0]] + Q = Q[bus < ppci["bus"].shape[0]] + cluster = cluster[bus < ppci["bus"].shape[0]] + bus = bus[bus < ppci["bus"].shape[0]] + for k in range(num_clusters): + cluster[cluster == cluster_list_tot[k]] = k + cluster = cluster.astype(int) + for i in range(len(P)): + bus_i, cluster_i, P_i, Q_i = bus[i], cluster[i], P[i], Q[i] + ppci["rated_power_clusters"][bus_i, cluster_i] += P_i + ppci["rated_power_clusters"][bus_i, cluster_i + num_clusters] += Q_i + ppci["rated_power_clusters"][bus_i, cluster_i + 2*num_clusters] += abs(0.03*P_i) # std dev cluster variability hardcoded, think how to change it + ppci["rated_power_clusters"][bus_i, cluster_i + 3*num_clusters] += abs(0.03*Q_i) # std dev cluster variability hardcoded, think how to change it + return ppci def merge_measurements(value, std_dev): @@ -416,6 +452,10 @@ def _build_measurement_vectors(ppci, update_meas_only=False): np.zeros(sum(i_degree_line_t_not_nan)) )).astype(bool) idx_non_imeas = np.flatnonzero(~imag_meas) + if ppci.algorithm == "af-wls": + balance_eq_meas = np.zeros(ppci["rated_power_clusters"].shape[0]).astype(np.float64) + af_vmeas = 0.4*np.ones(len(ppci["clusters"])) + z = np.concatenate((z, balance_eq_meas[ppci.non_slack_bus_mask], balance_eq_meas[ppci.non_slack_bus_mask], af_vmeas)) if not update_meas_only: # conserve the pandapower indices of measurements in the ppci order @@ -462,6 +502,14 @@ def _build_measurement_vectors(ppci, update_meas_only=False): any_degree_meas = np.any(np.r_[v_degree_bus_not_nan, i_degree_line_f_not_nan, i_degree_line_t_not_nan]) + if ppci.algorithm == "af-wls": + num_clusters = len(ppci["clusters"]) + P_balance_dev_std = np.sqrt(np.sum(np.square(ppci["rated_power_clusters"][:,2*num_clusters:3*num_clusters]),axis=1)) + Q_balance_dev_std = np.sqrt(np.sum(np.square(ppci["rated_power_clusters"][:,3*num_clusters:4*num_clusters]),axis=1)) + af_vmeas_dev_std = 0.15*np.ones(len(ppci["clusters"])) + r_cov = np.concatenate((r_cov, P_balance_dev_std[ppci.non_slack_bus_mask], Q_balance_dev_std[ppci.non_slack_bus_mask], af_vmeas_dev_std)) + meas_mask = np.concatenate((meas_mask, ppci.non_slack_bus_mask, ppci.non_slack_bus_mask, np.ones(len(ppci["clusters"])))) + return z, pp_meas_indices, r_cov, meas_mask, any_i_meas, any_degree_meas, idx_non_imeas else: return z @@ -472,7 +520,7 @@ def pp2eppci(net, v_start=None, delta_start=None, algorithm='wls', ppc=None, eppci=None): if isinstance(eppci, ExtendedPPCI): eppci.data = _add_measurements_to_ppci(net, eppci.data, zero_injection, algorithm) - eppci.update_meas() + eppci.update_meas(algorithm) return net, ppc, eppci else: # initialize ppc @@ -481,13 +529,14 @@ def pp2eppci(net, v_start=None, delta_start=None, # add measurements to ppci structure # Finished converting pandapower network to ppci ppci = _add_measurements_to_ppci(net, ppci, zero_injection, algorithm) - return net, ppc, ExtendedPPCI(ppci) + return net, ppc, ExtendedPPCI(ppci, algorithm) class ExtendedPPCI(UserDict): - def __init__(self, ppci): + def __init__(self, ppci, algorithm): """Initialize ppci object with measurements.""" self.data = ppci + self.algorithm = algorithm # Measurement relevant parameters self.z = None @@ -497,7 +546,6 @@ def __init__(self, ppci): self.non_nan_meas_selector = None self.any_i_meas = False self.any_degree_meas = False - self._initialize_meas() # check slack bus self.non_slack_buses = np.argwhere(ppci["bus"][:, idx_bus.BUS_TYPE] != 3).ravel() @@ -507,6 +555,9 @@ def __init__(self, ppci): np.ones(self.non_slack_bus_mask.shape[0], dtype=bool)].ravel() self.delta_v_bus_selector = np.flatnonzero(self.delta_v_bus_mask) + # Iniialize measurements + self._initialize_meas() + # Initialize state variable self.v_init = ppci["bus"][:, idx_bus.VM] self.delta_init = np.radians(ppci["bus"][:, idx_bus.VA]) @@ -514,6 +565,8 @@ def __init__(self, ppci): self.v = self.v_init.copy() self.delta = self.delta_init.copy() self.E = self.E_init.copy() + if algorithm == "af-wls": + self.E = np.concatenate((self.E, np.full(ppci["clusters"].shape,0.5))) def _initialize_meas(self): # calculate relevant vectors from ppci measurements diff --git a/pandapower/estimation/state_estimation.py b/pandapower/estimation/state_estimation.py index 109672d67..0f33bbe8f 100644 --- a/pandapower/estimation/state_estimation.py +++ b/pandapower/estimation/state_estimation.py @@ -9,7 +9,8 @@ from pandapower.estimation.algorithm.base import (WLSAlgorithm, WLSZeroInjectionConstraintsAlgorithm, - IRWLSAlgorithm) + IRWLSAlgorithm, + AFWLSAlgorithm) from pandapower.estimation.algorithm.lp import LPAlgorithm from pandapower.estimation.algorithm.optimization import OptAlgorithm from pandapower.estimation.ppc_conversion import pp2eppci, _initialize_voltage @@ -26,7 +27,8 @@ 'wls_with_zero_constraint': WLSZeroInjectionConstraintsAlgorithm, 'opt': OptAlgorithm, 'irwls': IRWLSAlgorithm, - 'lp': LPAlgorithm} + 'lp': LPAlgorithm, + 'af-wls': AFWLSAlgorithm} ALLOWED_OPT_VAR = {"a", "opt_method", "estimator"} @@ -272,6 +274,8 @@ def estimate(self, v_start='flat', delta_start='flat', calculate_voltage_angles= if self.solver.successful: self.net = eppci2pp(self.net, self.ppc, self.eppci) + if self.algorithm == "af-wls": + self.net["res_cluster_est"] = self.eppci.clusters else: self.logger.warning("Estimation failed! Pandapower network failed to update!") @@ -283,7 +287,7 @@ def estimate(self, v_start='flat', delta_start='flat', calculate_voltage_angles= if not self.recycle: self.ppc, self.eppci = None, None - if algorithm == "wls": + if algorithm == "wls" or algorithm == "af-wls": now = datetime.now() se_results = { "success": self.solver.successful,