Skip to content

Commit

Permalink
Merge pull request #8 from marcopau/Trafo-phase-shift-issue
Browse files Browse the repository at this point in the history 
WLS developments with 50Hz
  • Loading branch information
@marcopau
marcopau authored Dec 18, 2024
2 parents 71b7e3a + a56e378 commit 14e7a7c
Show file tree
Hide file tree
Showing 5 changed files with 117 additions and 43 deletions.
2 changes: 1 addition & 1 deletion pandapower/create.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5413,7 +5413,7 @@ def create_dcline(net, from_bus, to_bus, p_mw, loss_percent, loss_mw, vm_from_pu


def create_measurement(net, meas_type, element_type, value, std_dev, element, side=None,
check_existing=True, index=None, name=None, **kwargs):
check_existing=False, index=None, name=None, **kwargs):
"""
Creates a measurement, which is used by the estimation module. Possible types of measurements
are: v, p, q, i, va, ia
Expand Down
24 changes: 17 additions & 7 deletions pandapower/estimation/algorithm/base.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@

import numpy as np
from scipy.sparse import csr_matrix, vstack, hstack
from scipy.sparse.linalg import spsolve
from scipy.sparse.linalg import spsolve, norm, inv

from pandapower.estimation.algorithm.estimator import BaseEstimatorIRWLS, get_estimator
from pandapower.estimation.algorithm.matrix_base import BaseAlgebra, \
Expand Down Expand Up @@ -86,7 +86,7 @@ def estimate(self, eppci: ExtendedPPCI, **kwargs):

current_error, cur_it = 100., 0
# invert covariance matrix
eppci.r_cov[eppci.r_cov<(10**(-6))] = 10**(-6)
eppci.r_cov[eppci.r_cov<(10**(-5))] = 10**(-5)
r_inv = csr_matrix(np.diagflat(1 / eppci.r_cov ** 2))
E = eppci.E
while current_error > self.tolerance and cur_it < self.max_iterations:
Expand All @@ -109,20 +109,30 @@ def estimate(self, eppci: ExtendedPPCI, **kwargs):
# gain matrix G_m
# G_m = H^t * R^-1 * H
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 = G_m^-1 * (H' * R^-1 * r)
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.35:
d_E = d_E*0.35/current_error

# Update E with d_E
E += d_E.ravel()
eppci.update_E(E)

# 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))
# 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 with current measurements
if cur_it == 0 and eppci.any_i_meas:
Expand All @@ -139,7 +149,7 @@ def estimate(self, eppci: ExtendedPPCI, **kwargs):
# check if the estimation is successfull
self.check_result(current_error, cur_it)
self.iterations = cur_it
self.obj_func = obj_func
# 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()
Expand Down
101 changes: 79 additions & 22 deletions pandapower/estimation/ppc_conversion.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ def _init_ppc(net, v_start, delta_start, calculate_voltage_angles):
if np.any(net.trafo.shift_degree):
vm_backup = ppci["bus"][:, 7].copy()
pq_backup = ppci["bus"][:, [2, 3]].copy()
ppci["bus"][:, [2, 3]] = 0.
# ppci["bus"][:, [2, 3]] = 0.
ppci = _run_dc_pf(ppci)
ppci["bus"][:, 7] = vm_backup
ppci["bus"][:, [2, 3]] = pq_backup
Expand Down Expand Up @@ -180,15 +180,29 @@ def _add_measurements_to_ppci(net, ppci, zero_injection, algorithm):
if meas_type in ("p", "q"):
# Convert injection reference to consumption reference (P, Q)
this_meas.value *= -1
unique_bus_positions = np.unique(bus_positions)
if len(unique_bus_positions) < len(bus_positions):
std_logger.debug("P,Q Measurement duplication will be automatically merged!")
for bus in unique_bus_positions:
this_meas_on_bus = this_meas.iloc[np.argwhere(bus_positions == bus).ravel(), :]
bus_append[bus, BUS_MEAS_PPCI_IX[meas_type]["VALUE"]] = this_meas_on_bus.value.sum()
bus_append[bus, BUS_MEAS_PPCI_IX[meas_type]["STD"]] = this_meas_on_bus.std_dev.max()
unique_bus_positions = np.unique(bus_positions)
if len(unique_bus_positions) < len(bus_positions):
std_logger.debug("P,Q Measurement duplication will be automatically merged!")
for bus in unique_bus_positions:
this_meas_on_bus = this_meas.iloc[np.argwhere(bus_positions == bus).ravel(), :]
element_positions = this_meas_on_bus["element"].values.astype(np.int64)
unique_element_positions = np.unique(element_positions)
if meas_type in ("v", "va"):
merged_value, merged_std_dev = merge_measurements(this_meas_on_bus.value, this_meas_on_bus.std_dev)
bus_append[bus, BUS_MEAS_PPCI_IX[meas_type]["VALUE"]] = merged_value
bus_append[bus, BUS_MEAS_PPCI_IX[meas_type]["STD"]] = merged_std_dev
bus_append[bus, BUS_MEAS_PPCI_IX[meas_type]["IDX"]] = this_meas_on_bus.index[0]
continue
continue
for element in unique_element_positions:
this_meas_on_element = this_meas_on_bus.iloc[np.argwhere(element_positions == element).ravel(), :]
merged_value, merged_std_dev = merge_measurements(this_meas_on_element.value, this_meas_on_element.std_dev)
this_meas_on_bus.loc[this_meas_on_element.index[0], ["value", "std_dev"]] = [merged_value, merged_std_dev]
this_meas_on_bus.loc[this_meas_on_element.index[1:], ["value", "std_dev"]] = [0, 0]
sum_value, sum_std_dev = sum_measurements(this_meas_on_bus.value, this_meas_on_bus.std_dev)
bus_append[bus, BUS_MEAS_PPCI_IX[meas_type]["VALUE"]] = sum_value
bus_append[bus, BUS_MEAS_PPCI_IX[meas_type]["STD"]] = sum_std_dev
bus_append[bus, BUS_MEAS_PPCI_IX[meas_type]["IDX"]] = this_meas_on_bus.index[0]
continue

bus_append[bus_positions, BUS_MEAS_PPCI_IX[meas_type]["VALUE"]] = this_meas.value.values
bus_append[bus_positions, BUS_MEAS_PPCI_IX[meas_type]["STD"]] = this_meas.std_dev.values
Expand All @@ -197,12 +211,12 @@ def _add_measurements_to_ppci(net, ppci, zero_injection, algorithm):
# add zero injection measurement and labels defined in parameter zero_injection
bus_append = _add_zero_injection(net, ppci, bus_append, zero_injection)
# add virtual measurements for artificial buses, which were created because
# of an open line switch. p/q are 0. and std dev is 1. (small value)
# of an open line switch. p/q are 0. and std dev is 1e-6. (small value)
new_in_line_buses = np.setdiff1d(np.arange(ppci["bus"].shape[0]), map_bus[map_bus >= 0])
bus_append[new_in_line_buses, 2] = 0.
bus_append[new_in_line_buses, 3] = 1.
bus_append[new_in_line_buses, 3] = 1e-6
bus_append[new_in_line_buses, 4] = 0.
bus_append[new_in_line_buses, 5] = 1.
bus_append[new_in_line_buses, 5] = 1e-6

# add 15 columns to mpc[branch] for Im_from, Im_from std dev, Im_to, Im_to std dev,
# P_from, P_from std dev, P_to, P_to std dev, Q_from, Q_from std dev, Q_to, Q_to std dev,
Expand All @@ -224,15 +238,19 @@ def _add_measurements_to_ppci(net, ppci, zero_injection, algorithm):
net[br_type][BR_SIDE[br_type][br_side]+"_bus"]
[this_meas.element]).values]
ix_side = br_map[meas_this_side.element.values].values
branch_append[ix_side,
BR_MEAS_PPCI_IX[(meas_type, br_side)]["VALUE"]] =\
meas_this_side.value.values
branch_append[ix_side,
BR_MEAS_PPCI_IX[(meas_type, br_side)]["STD"]] =\
meas_this_side.std_dev.values
branch_append[ix_side,
BR_MEAS_PPCI_IX[(meas_type, br_side)]["IDX"]] =\
meas_this_side.index.values
unique_ix_side = np.unique(ix_side)
if len(unique_ix_side) < len(ix_side):
for branch in unique_ix_side:
this_meas_on_branch = meas_this_side.iloc[np.argwhere(ix_side == branch).ravel(), :]
merged_value, merged_std_dev = merge_measurements(this_meas_on_branch.value, this_meas_on_branch.std_dev)
branch_append[branch, BR_MEAS_PPCI_IX[(meas_type, br_side)]["VALUE"]] = merged_value
branch_append[branch, BR_MEAS_PPCI_IX[(meas_type, br_side)]["STD"]] = merged_std_dev
branch_append[branch, BR_MEAS_PPCI_IX[(meas_type, br_side)]["IDX"]] = this_meas_on_branch.index[0]
continue

branch_append[ix_side, BR_MEAS_PPCI_IX[(meas_type, br_side)]["VALUE"]] = meas_this_side.value.values
branch_append[ix_side, BR_MEAS_PPCI_IX[(meas_type, br_side)]["STD"]] = meas_this_side.std_dev.values
branch_append[ix_side, BR_MEAS_PPCI_IX[(meas_type, br_side)]["IDX"]] = meas_this_side.index.values

# Add measurements for trafo3w
if map_trafo3w is not None:
Expand All @@ -250,6 +268,31 @@ def _add_measurements_to_ppci(net, ppci, zero_injection, algorithm):
ix_hv = [map_trafo3w[t]['hv'] for t in meas_hv.element.values]
ix_mv = [map_trafo3w[t]['mv'] for t in meas_mv.element.values]
ix_lv = [map_trafo3w[t]['lv'] for t in meas_lv.element.values]
unique_ix_hv = np.unique(ix_hv)
unique_ix_mv = np.unique(ix_mv)
unique_ix_lv = np.unique(ix_lv)
if len(unique_ix_hv) < len(ix_hv):
for branch in unique_ix_hv:
this_meas_on_branch = meas_hv.iloc[np.argwhere(ix_hv == branch).ravel(), :]
merged_value, merged_std_dev = merge_measurements(this_meas_on_branch.value, this_meas_on_branch.std_dev)
branch_append[branch, BR_MEAS_PPCI_IX[(meas_type, "f")]["VALUE"]] = merged_value
branch_append[branch, BR_MEAS_PPCI_IX[(meas_type, "f")]["STD"]] = merged_std_dev
branch_append[branch, BR_MEAS_PPCI_IX[(meas_type, "f")]["IDX"]] = this_meas_on_branch.index[0]
if len(unique_ix_mv) < len(ix_mv):
for branch in unique_ix_mv:
this_meas_on_branch = meas_mv.iloc[np.argwhere(ix_mv == branch).ravel(), :]
merged_value, merged_std_dev = merge_measurements(this_meas_on_branch.value, this_meas_on_branch.std_dev)
branch_append[branch, BR_MEAS_PPCI_IX[(meas_type, "t")]["VALUE"]] = merged_value
branch_append[branch, BR_MEAS_PPCI_IX[(meas_type, "t")]["STD"]] = merged_std_dev
branch_append[branch, BR_MEAS_PPCI_IX[(meas_type, "t")]["IDX"]] = this_meas_on_branch.index[0]
if len(unique_ix_lv) < len(ix_lv):
for branch in unique_ix_lv:
this_meas_on_branch = meas_lv.iloc[np.argwhere(ix_lv == branch).ravel(), :]
merged_value, merged_std_dev = merge_measurements(this_meas_on_branch.value, this_meas_on_branch.std_dev)
branch_append[branch, BR_MEAS_PPCI_IX[(meas_type, "t")]["VALUE"]] = merged_value
branch_append[branch, BR_MEAS_PPCI_IX[(meas_type, "t")]["STD"]] = merged_std_dev
branch_append[branch, BR_MEAS_PPCI_IX[(meas_type, "t")]["IDX"]] = this_meas_on_branch.index[0]
continue
branch_append[ix_hv, BR_MEAS_PPCI_IX[(meas_type, "f")]["VALUE"]] = meas_hv.value.values
branch_append[ix_hv, BR_MEAS_PPCI_IX[(meas_type, "f")]["STD"]] = meas_hv.std_dev.values
branch_append[ix_hv, BR_MEAS_PPCI_IX[(meas_type, "f")]["IDX"]] = meas_hv.index.values
Expand All @@ -272,6 +315,20 @@ def _add_measurements_to_ppci(net, ppci, zero_injection, algorithm):
ppci["branch"][:, branch_cols: branch_cols + branch_cols_se] = branch_append
return ppci

def merge_measurements(value, std_dev):
weight = np.divide(1, np.square(std_dev))
merged_variance = np.divide(1, weight.sum())
merged_std_dev = np.sqrt(merged_variance)
weighted_value = np.multiply(value, weight)
merged_value = np.multiply(weighted_value.sum(), merged_variance)
return merged_value, merged_std_dev

def sum_measurements(value, std_dev):
sum_values = value.values.sum()
variance = np.square(std_dev.values)
sum_variance = variance.sum()
sum_std_dev = np.sqrt(sum_variance)
return sum_values, sum_std_dev

def _add_zero_injection(net, ppci, bus_append, zero_injection):
"""
Expand All @@ -294,7 +351,7 @@ def _add_zero_injection(net, ppci, bus_append, zero_injection):
if isinstance(zero_injection, str):
if zero_injection == 'auto':
# identify bus without elements and pq measurements as zero injection
zero_inj_bus_mask = (ppci["bus"][:, 1] == 1) & (ppci["bus"][:, 2:6] == 0).all(axis=1) & \
zero_inj_bus_mask = (ppci["bus"][:, 1] == 1) & (ppci["bus"][:, 2:4] == 0).all(axis=1) & \
np.isnan(bus_append[:, P:(Q_STD + 1)]).all(axis=1)
bus_append[zero_inj_bus_mask, ZERO_INJ_FLAG] = True
elif zero_injection != "aux_bus":
Expand Down
31 changes: 19 additions & 12 deletions pandapower/estimation/results.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,18 +51,25 @@ def _extract_result_ppci_to_pp(net, ppc, ppci):
net.res_bus_est.loc[merged_bus_idx, 'p_mw'] = 0
net.res_bus_est.loc[merged_bus_idx, "q_mvar"] = 0
# add shunt power because the injection at the node computed via Ybus is only the extra injection on top of the shunt
if ~net["shunt"].empty:
for i in range(net["shunt"].shape[0]):
bus = net.shunt.bus.iloc[i]
Sn = complex(net.shunt.p_mw.iloc[i],net.shunt.q_mvar.iloc[i])*net.shunt.step.iloc[i]
Ysh = Sn / (net.shunt.vn_kv.iloc[i]**2)
V = net["res_bus_est"].loc[bus,"vm_pu"]*net["bus"].loc[bus,"vn_kv"]
Sinj = Ysh*(V**2)
net["res_bus_est"].loc[bus,"p_mw"] += Sinj.real
net["res_bus_est"].loc[bus,"q_mvar"] += Sinj.imag
net["res_shunt_est"].loc[net["shunt"].loc[:,"bus"]==bus,"p_mw"] = Sinj.real
net["res_shunt_est"].loc[net["shunt"].loc[:,"bus"]==bus,"q_mvar"] = Sinj.imag
net["res_shunt_est"].loc[net["shunt"].loc[:,"bus"]==bus,"vm_pu"] = net["res_bus_est"].loc[bus,"vm_pu"]
for element in ["shunt", "ward", "xward"]:
if ~net[element].empty:
for i in range(net[element].shape[0]):
bus = net[element].bus.iloc[i]
if element == "shunt":
Sn = complex(net[element].p_mw.iloc[i],net[element].q_mvar.iloc[i])*net[element].step.iloc[i]
Ysh = Sn / (net[element].vn_kv.iloc[i]**2)
else:
Sn = complex(net[element].pz_mw.iloc[i],net[element].qz_mvar.iloc[i])
Ysh = Sn / (net.bus.loc[bus,"vn_kv"]**2)
V = net["res_bus_est"].loc[bus,"vm_pu"]*net["bus"].loc[bus,"vn_kv"]
Sinj = Ysh*(V**2)
net["res_bus_est"].loc[bus,"p_mw"] += Sinj.real
net["res_bus_est"].loc[bus,"q_mvar"] += Sinj.imag
if element == "shunt":
element_res_est = "res_" + element + "_est"
net[element_res_est].loc[net[element].loc[:,"bus"]==bus,"p_mw"] = Sinj.real
net[element_res_est].loc[net[element].loc[:,"bus"]==bus,"q_mvar"] = Sinj.imag
net[element_res_est].loc[net[element].loc[:,"bus"]==bus,"vm_pu"] = net["res_bus_est"].loc[bus,"vm_pu"]
return net


Expand Down
2 changes: 1 addition & 1 deletion pandapower/results_branch.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -415,7 +415,7 @@ def _get_trafo_results_3ph(net, ppc0, ppc1, ppc2, I012_f, V012_f, I012_t, V012_t
Iabc_sum += abs(I_branch_abc[:, x])

# Loads
load_index = np.where(net.asymmetric_load['bus'] == lv_bus)[0]
load_index = net.asymmetric_load.index[net.asymmetric_load['bus'] == lv_bus]
if len(load_index > 0):
S_load_abc = abs(np.array([
np.array(net.res_asymmetric_load_3ph['p_a_mw'][load_index]
Expand Down

0 comments on commit 14e7a7c

Please sign in to comment.