The new HELM Powerflow in GridCal

Hello, I am new at the community. I am catching up with the discussion, I look forward to contribute.

I am the authored of the publication that @robbie.morrison just mentioned, it describes how the distributed slack bus model is successfully implemented in HELM. The slack voltage bus is embedded and there are two models for the PV buses, I believe these are what @Zalnd referred to where Q was eliminated or Q is taken as a variable. The distributed slack bus is implemented by adding an extra variable that represents the total active power losses and an extra equation that constrains the active power generated at the slack bus, this equation happens to be part of the same equation for a PV bus.

All this is implemented and validated in the python repository: https://github.com/HELMpy/HELMpy which was publish on august of the last year. This repository contains fully functional implementations of HELM. The repository needs work but I will start working on it.

Do you mean that the B shunt on Bus 11 is 4.4 Mvar?

If so, we are looking to the same network now. Have you tried to reach the 0.99183 load factor after correcting the network data?

I’m almost finishing to implement my formulation in GridCal. I haven’t finished it yet because I started to investigate the admittance matrices. Somehow, the calculated Ybus, Yseries and Yshunt doesn’t obey ‘Ybus = Yseries + Yshunt’. I guess it’s a numerical issue since Iwamoto’s network is ill-conditioned. Could you confirm that the equation hold in your working station?

The following file contains the variable created by the functions ‘compute’ and ‘calc_connectivity’ in ‘snapshot_static_inputs.py’. It was saved using pickle.

store.zip (5.7 KB)

Hello @Zalnd,
From what I understand the relation ‘Ybus = Yseries + Yshunt’ must be satisfied even if the network is ill-conditioned. If there exist phase shifters, the Yseries matrix must be separated in its symmetric and unsymmetric parts, and the relation ’ Y = Yseries_symmetric + Yseries_unsymmetric + Yshunt’ must be satisfied.

I have made a script that calculates this three matrixes, it is part of a module of my repository. It is ready to be applied to the Iwamoto 11 bus.
Hope this may help.

Y matrix.zip (11.8 KB)

Hi,

I’m not questioning from a mathematical point of view. I’m questioning specifically GridCal implementation. The code adapts an “incidence matrix formulation” and all the three terms, Ybus, Yseries and Yshunt are calculated using independent equations. I mean that the formulation doesn’t calculate two terms independently and enforce the equality by calculating the last one using the equality itself. The formulation obtains the three terms using completely independent equations.

Somehow I’m getting 00041723356008991144 as the maximum absolute error from “Ybus - Yseries - Yshunt”. I didn’t have enough time to take a good look at it but the formulation seems to be correct. If the formulation is correct but it leads to a wrong answer, the issue may be numerical.

You may check out the “incidence matrix formulation” I’m talking about in calc_connectivity function: https://github.com/SanPen/GridCal/blob/e789689a6b7674616b0ed2d795d54ba40d9255bd/src/GridCal/Engine/Core/series_static_inputs.py#L36

Anyway, I appreciate your prompt reply!

I found the source of the problem. It’s the GridCal formulation. Considering the Universal Branch Model, I think the “from-from primitive admittance vector” should be calculated by:

(
)
Yff = ( Ys / ( tap * np.conj(tap) ) + GBc / 2.0) / (tap_f * tap_f) LINE 93 of series_static_inputs.py
(
)

The shunt element on HV side shouldn’t be reflected through the “Real Transformer”. Just through the “HV Virtual Transformer”.

It also must be pointed out that virtual taps are completely ignored in Yseries and Yshunt calculation. This doesn’t make any difference to the Iwamoto’s 11 Bus System, though.

I need to carefully look at this from all angles. I am refactoring the engine to remove the generic branch to specialise into lines, transformers, and other branch models like HVDC or VSC. The virtual taps are a thing that have their use for transformers and should keep working after changes.

The way to proceed is to split each device into symmetric, unsymmetric and shunt as Jose Juan points out. A reverse enforcing of Yshunt by a subtraction is not a good practice. I’ll check the change that you propose and report back.

–Edit–
I added a test, to check that Ybus = Yseries + Yshunt, and yes it is equal.

I have noticed that the Yshunt from the output of “numerical_circuit.compute()” is “.Ysh_helm” and updated the main function of helm_power_flow.py a couple days ago. So, this isn’t the source of the issue.

Have you updated series_static_inputs.py since 3.7.0 version? I restate that this series_static_inputs.py fails to calculate coherent admittance matrices.

When looking deeper into this issue, it seems that the Ybus formulation considers the matpower branch model (page 26) and the Yshunt formulation considers this universal branch model. The position of the shunt element on HV side is different. The solution I pointed out above corrects the Ybus formulation in order to represent the universal branch model.

EDIT: I was treating the “series_static_inputs.py” as the “snapshot_static_inputs.py” all along. You have updated the “series_static_inputs.py” since the 3.7.0. But you haven’t updated the “series_static_inputs.py”. Perhaps only the old formulation causes the issue.

EDIT2: The formulation haven’t changed, so the problem still persists. I think you should change your test code to compare the absolute value of the difference. It seems that numpy less function completely ignores the imaginary part.

Yes, the admittance matrix looks the same and everything seems to match.

However, I still can’t figure out the right s_0 vector to work with. According to the Sigma plot the system has a solution. The conformation of the s_0 vector I implemented is quite rudimentary, but still, it worked nice with a total loading factor of about 0.9. You mentioned that your initial load factor was 0.1, so what where the other ones?

In theory, the P-W method should be capable to compute the solution with an error similar to the desired precision. Therefore, in Python I estimate that the maximum mismatch could be somewhat close to 10^{-12} or 10^{-13} due to the Padé approximants reaching a plateau around that point. If you use more than 3 steps, are you able to improve the solution?

Since this is a scientific forum, I posted a test proving that the admitances match. However it seems that you still disagree with such test, so I would like to see code from you testing the hypothesis that you propose.

Furthermore, it has been implied that somehow I changed the code to correct a bug that only you seem to find, but cannot prove. You may check the commits. I promise I did not hack GitHub.

In the master branch, there is a matpower-like branch model enhanced with the virtual taps (model that has been checked against ETAP) However thais is going to change soon, since in a separated branch I have reworked all the internals to get rid of the branch model and replace it with particular models (line, transformer, etc
)

Currently, the admittance splitting is between Series and Shunt, whereas in the future there shall additionally be a split between symmetrical-series, unsymmetrical-series and shunt. That is the only change that is concerning me at the moment since it enables P-W.

On the personal side, I dialogue with people that show value, humbleness and act with respect.

Best regards,
Santiago

@SanPen
I used pip to install GridCal and the snapshot_static_inputs.py in github is different from the snapshot_static_inputs.py in 3.7.0 version installed using pip (8.5 KB). The old version I was using could have been the source of the problem. This was strictly what I meant.

About your test function, I tried to point out that, yes, the Boolean Expression

((island.Ybus - (island.Yseries + diags(island.Ysh_helm))).data < 1e-9).all()

gives the answer “True”, but

(abs((island.Ybus - (island.Yseries + diags(island.Ysh_helm))).data) < 1e-9).all()

doesn’t.

The test function I’m proposing (1.2 KB)

@jfb
I don’t understand why you have to define all the transformations a priori. I think you should find the next transformation based on the maximum mismatch obtained along s real axis.

I agree, you are correct on that. My implementation initializes all the s values at once, despite that I manually select them while computing each step. I am still puzzled to know the value you have chosen in each step. Can you please share them so I can test the system and find out if the mismatch improves?