Merge pull request #9 from marcopau/AF-WLS

Implementation of AF-WLS

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

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)

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,13 +555,18 @@ 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])
         self.E_init = np.r_[self.delta_init[self.non_slack_bus_mask], self.v_init]
         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