Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

State Estimation modifications for WLS #2414

Draft
wants to merge 6 commits into
base: develop
Choose a base branch
from

Conversation

marcopau
Copy link

@marcopau marcopau commented Oct 7, 2024

Fixed some bugs to state estimation code, like:

  • power injection results
  • shunts handling
  • current magnitude measurements handling

@vogt31337 vogt31337 marked this pull request as draft October 7, 2024 13:11
@Heiner92
Copy link

Great initiative! I hope these changes will improve the results I'm getting in distribution system state estimation where current measurements are usually only available as magnitudes.

root1 = ds_find(ar, bus1)
root2 = ds_find(ar, bus2)
if root1 == root2:
return
if bus_is_pv[root2]:
if bus_is_active[root2] and ~bus_is_pv[root1]: # if root2 is an active bus (load or gen) and root1 is not a PV bus, it will merge root1 into root2

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here we should actually also make sure that root2 is not external grid or reference bus as well

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If root2 is an external grid it is still fine to merge root1 into root2. You want to keep the PV bus, and in this case root2 (which is the PV bus) would be kept.

if net._fused_bb_switches.any():
# Note: this might seem a little odd - first constructing a pp to ppc mapping without
# fused busses and then update the entries. The alternative (to construct the final
# mapping at once) would require to determine how to fuse busses and which busses
# are not part of any opened bus-bus switches first. It turns out, that the latter takes
# quite some time in the average usecase, where #busses >> #bus-bus switches.

# create mask arrays to distinguish different grid elements

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Redundant code

# one pv bus. Get bus from pv_buses_in_set
map_to = bus_lookup[pv_buses_in_set.pop()]
ref_bus = pv_buses_in_set.pop()

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

not a very robust way since if pandapower decides to not have the first element as REF, this wont work and will be difficult to debug

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Explain what "if pandapower decides to not have the first element as REF" means. Feel free also to provide an example where this does not work, so I can look at it.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if the pv_buses_in_set first element which you are popping is not REF bus because of breaking change from pandapower

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In pv_buses_in_set you have the list of pv buses for the considered set. What I am popping is a PV bus. That bus is then taken as bus to be kept later when doing the merging among multiple buses.

@@ -222,11 +220,12 @@ def _build_pp_gen(net, ppc, f, t):

# set bus values for generator buses
ppc["bus"][gen_buses[ppc["bus"][gen_buses, BUS_TYPE] != REF], BUS_TYPE] = PV

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why here all other buses are set to PV

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here there is no change. Gen buses are always modelled as PV nodes in pandapower, sgens as PQ.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes i am aware of that, the question was why here are the PQ buses exclusively excluded

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function only deals with gen_buses, and it assigns the label "PV" to them. Other buses are handled somewhere else.


vb = range(len(V))
diagV = sparse((V, (vb, vb)))
diagVnorm = sparse((V / abs(V), (vb, vb)))

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here there should be try except block because division with zero can occur due to abs(V), it will help enormously if throughindices one could also find the pandapower element (in dataframe) which leads to this exception , for example NaNor zero g_us_per_km attribute for lines, etc

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A division by 0 when dividing by the voltage magnitude? The second part of the comment (looking at the pandapower elements that can bring a division by zero) is not something that has to be addressed in a function to calculate the Jacobian. There are other functions and routines in pandapower that do that.

Copy link

@AnkurArohi AnkurArohi Oct 14, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes we have this many times for networks which have problems in attribute and is very central as then the algorithm does not work. I am unaware of other functions that does that

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see, but then the problem in that case is with the network data and with its processing, not with this pull request or with this specific function, so I guess we should shift the discussion somewhere else.

ib = range(len(If))
idxf = abs(If) == 0
idxt = abs(It) == 0
diagIfnorm = sparse((conj(If) / abs(If), (ib, ib)))

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same here with the try and except block

Copy link
Author

@marcopau marcopau Oct 14, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is a routine in the WLS to avoid using the current magnitude measurements at the first iteration, when currents are initialized as zero. In the other iterations it is quite unlikely to have branch currents exactly equal to 0. If that case happens, I believe we should have some more elaborated routine to deal with it, not just a try and except. It is something that requires some thinking and that would come in one of the next round of updates. For the moment, the priority was to provide a Jacobian for current magnitude measurements that works, at least in most of the cases.

idxt = abs(It) == 0
diagIfnorm = sparse((conj(If) / abs(If), (ib, ib)))
diagItnorm = sparse((conj(It) / abs(It), (ib, ib)))
diagIfnorm[idxf,idxf] = 0

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why this is equated to zero?, please comment

dita_dv = (np.angle(1e-5 * dit_dv + It.reshape((-1, 1))) - np.angle(It.reshape((-1, 1))))/1e-5
difm_dth, difm_dv, ditm_dth, ditm_dv = dIbr_dV_new(self.eppci['branch'], self.Yf, self.Yt, V)
difm_dth, difm_dv, ditm_dth, ditm_dv = map(lambda m: m.toarray(), (difm_dth, difm_dv, ditm_dth, ditm_dv))
difa_dth, difa_dv, dita_dth, dita_dv = 0*difm_dth, 0*difm_dv, 0*ditm_dth, 0*ditm_dv

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why are all these values zero?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

At the moment, we are not yet supporting current phase angles for state estimation. I will add a line of comment there.

@@ -83,6 +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)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This value of 10-6 should ideally be parameterizable

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is an option I am considering, but I am not yet sure if to put it as an optional parameter for the user or not (I see some pros and cons). In any case, this will be definitely improved in the next round of updates, for the moment it was urgent to avoid having 0 or close to 0 values in the covariance matrix.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is central to the feasibility of the solution and to evaluate how good the estimate is and where are problems in the network, so one could change it and then simulate multiple usecases

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants